Distinct frequency bands in the local ﬁeld potential are differently tuned to stimulus drift rate

Distinct frequency bands in the local ﬁeld potential are differently tuned to stimulus drift rate. ﬁeld potential (LFP) recorded with a microelectrode reﬂects the activity of several neural processes, including afferent synaptic inputs, microcircuit-level computations, and spiking activity. Objectively probing their contribution requires a design that allows dissociation between these potential contributors. Earlier reports have shown that the primate lateral geniculate nucleus (LGN) has a higher temporal frequency (drift rate) cutoff than the primary visual cortex (V1), such that at higher drift rates inputs into V1 from the LGN continue to persist, whereas output ceases, permit- ting partial dissociation. Using chronic microelectrode arrays, we recorded spikes and LFP from V1 of passively ﬁxating macaques while presenting sinusoidal gratings drifting over a wide range. We further optimized the gratings to produce strong gamma oscillations, since recent studies in rodent V1 have reported LGN-dependent narrow-band gamma oscillations. Consistent with earlier reports, power in higher LFP frequencies (above ~140 Hz) tracked the population ﬁring rate and were tuned to preferred drift rates similar to those for spikes. Signiﬁcantly, power in the lower (up to ~40 Hz) frequencies increased transiently in the early epoch after stimulus onset, even at high drift rates, and had preferred drift rates higher than for spikes/high gamma. Narrow-band gamma (50–80 Hz) power was not strongly correlated with power in high or low frequencies and had much lower preferred temporal frequencies. Our results demonstrate that distinct frequency bands of the V1 LFP show diverse tuning proﬁles, which may potentially convey different attributes of the underlying neural activity. NEW & NOTEWORTHY In recent years the local ﬁeld potential (LFP) has been increasingly studied, but interpreting its rich fre- quency content has been difﬁcult. We use a stimulus manipulation that generates different tuning proﬁles for low, gamma, and high frequen- cies of the LFP, suggesting contributions from potentially different sources. Our results have possible implications for design of better neural prosthesis systems and brain-machine interfacing applications.


INTRODUCTION
Local field potentials (LFPs) are low-frequency electrical potentials obtained from the extracellular space around a pop-ulation of neurons, extracted by low-pass filtering (typically below 500 Hz) the electrical activity obtained from intracortical recordings. The LFP is an inherently ambiguous measure of neuronal activity, with contributions from many sources, such as synaptic activity in neurons due to inputs arriving from upstream areas Mitzdorf 1985Mitzdorf , 1987, action potentials (Manning et al. 2009;Ray et al. 2008a;Ray and Maunsell 2011;Reimann et al. 2013;Schomburg et al. 2012), calcium spikes (Schiller et al. 2000), intrinsic currents (Llinás 1988), spike afterpotentials (Gustafsson and Pinter 1984;Harada and Takahashi 1983;Reimann et al. 2013), gap junctions (Draguhn et al. 1998;Traub and Bibbig 2000), neuron-glia interactions (Kang et al. 1998), intrinsic oscillations , and ephaptic coupling (Anastassiou et al. 2011). The LFP also often shows stimulus-induced narrow-band oscillations in the gamma range (30 -80 Hz), which, in primates, have been shown to have no consistent relationship with spiking activity (Jia et al. 2013) and are thought to reflect complex excitation-inhibition interactions instead (Buzsáki and Wang 2012). However, recent studies in rodent primary visual cortex (V1) have shown a direct relationship between inputs from lateral geniculate nucleus (LGN) and a narrow-band gamma rhythm around 60 Hz (Saleem et al. 2017;Storchi et al. 2017;Veit et al. 2017). Therefore, a stimulus paradigm in which synaptic (input) vs. spiking (output) activity can potentially be modulated differently can shed light on the neural mechanisms behind the generation of gamma rhythm as well.
In primate V1, such differential modulation is possible, at least partially, with drifting gratings. In light of evidence that V1 neurons do not respond to the full range of temporal frequencies (stimulus drift rates) present in the thalamic input, such that both optimal temporal frequency range and highfrequency cutoff are much higher in the LGN (Hawken et al. 1996) compared with V1 (Foster et al. 1985), an earlier study used such a paradigm to show that the blood oxygen leveldependent (BOLD) signal in functional magnetic resonance imaging (fMRI) is more closely associated with synaptic than spiking activity (Viswanathan and Freeman 2007). However, they only used two drift rates (3 and 30 Hz) and did not study the power in different frequency bands of the LFP in detail. Furthermore, their stimuli were not optimized to generate salient gamma rhythm. In addition, although these studies characterized the temporal frequency tuning of spiking activ-ity, no comparable data are available for the V1 LFP, which could provide meaningful insights into its composition.
Onset of a visual stimulus activates the entire visual pathway. Although the first inputs are expected to be from an upstream area (from LGN in the case of V1), later inputs could also arrive from lateral connections and feedback from downstream areas, which may reflect in the LFP signal. To relate the LFP power at different frequencies to various neural processes, it is therefore essential to generate the time-frequency spectrum of the LFP at high resolution. We showed previously that a matching pursuit (MP) algorithm is particularly useful in resolving the fast transients in the signal (Chandran KS et al. 2016), which is likely associated with the first inputs arriving from LGN.
Here we recorded LFP and spiking activity from V1 of passively fixating macaques while they viewed gratings drifting at eight different frequencies, generated high-resolution time-frequency spectra of the LFP with MP, and studied the relationship between the power in different frequency bands and their temporal frequency tuning preferences. Furthermore, unlike the previous study that used a similar paradigm (Viswanathan and Freeman 2007), gratings were optimized to generate strong gamma oscillations in the LFP.

MATERIALS AND METHODS
All animal procedures and experiments in this study were reviewed and approved by the Institutional Animal Ethics Committee of the Indian Institute of Science and conducted in accordance with the guidelines approved by the Committee for the Purpose of Control and Supervision of Experiments on Animals.
Experimental setup. Recordings were performed on two female bonnet monkeys (Macaca radiata; weights 3 and 4.6 kg, aged 13 and 17 yr, respectively). Before training, a titanium head post was implanted under general anesthesia. After the monkey learned the behavioral task, a microelectrode array grid (Blackrock Microsystems) was implanted in the V1 cortex,~12 mm from the occipital ridge and~10 -15 mm from midline. In monkey 1 a 9 ϫ 9 microelectrode array that had 81 active electrodes was implanted in the left V1 cortex, whereas in monkey 2 a 10 ϫ 10 microelectrode array with 96 active electrodes was implanted in the right V1 cortex. The reference wires were either positioned securely over the dura near the array implant or wrapped around titanium screws on the skull near the site of craniotomy. Receptive fields of neurons recorded in monkey 1 were in the lower right quadrant, spanning eccentricities of~3.3°to~5.4°, whereas receptive fields of neurons recorded in monkey 2 were more foveal and in the lower left quadrant, spanning eccentricities of~1.4°t o~1.75°. The microelectrodes were 1 mm long; from estimates of the thickness of primate V1 (Hubel and Wiesel 1977), we expect them to be in layer 4 if the base of the array was flush with the cortical surface such that the individual microelectrodes penetrated orthogonally, or layer 2/3 otherwise. Histology was not performed. Care was taken to minimize pain and discomfort in the animals during all surgical and experimental procedures.
During the experiment, the monkey was seated in a primate chair at~50 cm from the monitor. Stimuli were presented on a gammacorrected LCD monitor (BenQ XL2411Z) at a resolution of 1,280 ϫ 720 and a refresh rate of 100 Hz with custom stimulus presentation software (Lablib, Mac OS X). During the task, the monkey's eye position was continuously monitored with an infrared eye tracking system (ISCAN) and a data acquisition interface (InstruTECH ITC-18; HEKA Elektronik).
Neural recordings were done with the Cerebus Neural Signal Processor (Blackrock Microsystems). During acquisition, the raw signal was band-pass filtered between 0.3 Hz (Butterworth filter, 1st order, analog) and 500 Hz (Butterworth filter, 4th order, digital) and digitized (2 kHz, 16 bit) to obtain the LFP. Spiking activity was obtained by band-pass filtering the raw signal between 250 Hz (Butterworth filter, 4th order, digital) and 7,500 Hz (Butterworth filter, 3rd order, analog) followed by an amplitude threshold set at~5 times the standard deviation of the signal. Concurrently, the raw waveform segments from each electrode crossing the amplitude threshold were digitized at 30 kHz and saved as spike waveforms without any online sorting.
Visual stimuli. We have recently shown that large grating stimuli generate a second gamma rhythm at a lower frequency between 25 and 40 Hz in addition to the traditional ("fast") gamma rhythm in V1 (Murty et al. 2018). This "slow" gamma is salient at midlevel contrasts, low drift rates, and large stimulus size exceeding 4°and is tuned to an orientation that is different from fast gamma (for the 2 monkeys tested here, slow gamma was most salient at~20°and~45°, while fast gamma was maximal at 90°). Here we optimized gratings to produce a salient fast gamma between 50 and 80 Hz but negligible slow gamma.
For each monkey, the radius of the grating was chosen to be just enough to cover the receptive field span of the microelectrode array (monkey 1: 1.8°; monkey 2: 0.8°). The gratings were presented at full contrast, at a spatial frequency of 4 cycles per degree (cpd), and at two orientations (0°and 90°), chosen pseudorandomly; all results are shown for the orientation that produced stronger fast gamma (90°in both monkeys). Similar results were obtained with the 0°orientation (data not shown). These gratings were presented at eight drift rates: 0 (static), 1, 2, 4, 8, 16, 32, and 50 cycles per second (cps).
Each experimental trial began with the presentation of a small central fixation spot (0.1°radius), which the monkey was required to fixate within a 2°square window (not visible to the monkey). After a delay (monkey 1: 1,000 ms; monkey 2: 1,500 ms), two or three stimuli were sequentially presented in each trial. For monkey 1, we presented two or three stimuli with a duration of 800 ms each and an interstimulus interval of 700 ms. For monkey 2, we presented two stimuli with a duration of 1,500 ms each and an interstimulus interval of 1,500 ms. Only trials for which the monkey maintained fixation throughout were used for analysis, and the monkey was rewarded with a juice drop at the end of the trial.
Electrode selection. As in our previous reports, we used only those electrodes that yielded reliable LFP responses and stable estimates of receptive field centers across days, resulting in a total of 78 electrodes for monkey 1 and 39 electrodes for monkey 2, which we called LFP electrodes. Spike sorting was performed off-line on all LFP electrodes. For obtaining "spiking" electrodes, we selected only those LFP electrodes that had a signal-to-noise ratio (SNR; as described in Kelly et al. 2007) of at least 2 and a mean firing rate increase of at least 1 spike/s in the time period from 50 ms to 250 ms after stimulus onset for the static grating compared with baseline. This resulted in 29 (mean SNR ϭ 2.88 Ϯ 0.16; mean increase in firing rate ϭ 15.51 Ϯ 2.43 spikes/s) and 19 (mean SNR ϭ 2.66 Ϯ 0.19; mean increase in firing rate ϭ 7.62 Ϯ 3.44 spikes/s) spiking electrodes in the two monkeys. Since we did not get good isolation (as is typical in chronically implanted microelectrode recordings), we refer to them as multiunits. The ratio of the amplitude of the Fourier transform of the peristimulus time histogram (PSTH) of the firing rate at the drift rate frequency vs. the DC component (F1-to-F0 ratio) was less than unity for all the multiunits in both monkeys (i.e., they behaved like complex cells; Skottun et al. 1991). This is not surprising, since there is a preponderance of complex cells over simple cells in the superficial layers of V1 where our microelectrodes were placed (Kelly et al. 2007), and therefore multiunit activity is expected to resemble the complex cells. Relatively smaller drift rate-locked responses in the LFP (smaller F1 than F0) were not because of an insufficient number of stimulus repeats (~17, as described below), since the LFP spectra from the same sites showed a prominent peak at twice the temporal frequency if instead of drifting the gratings were counterphased (data not shown), consistent with previous studies (Katzner et al. 2009;Ray and Maunsell 2011).
Data processing. Data extraction was performed with custom scripts written in MATLAB (The MathWorks). For each successful trial, we extracted a portion of the signal (monkey 1: 2,048-ms length; monkey 2: 4,096 ms length) around each stimulus onset. Segments with excessive artifacts (sharp transients and high-amplitude drifts) were rejected (~12% and~6% of segments), yielding, on average,~18 and~16 stimulus repeats per condition for the two monkeys. Analyses were performed in an "early" time range extending from 50 ms to 250 ms after stimulus onset (to capture transient effects related to the onset of visual stimulus) and compared to a baseline period defined as the 200-ms time range preceding stimulus onset. The analysis duration was limited to 200 ms to ensure that the influence of feedback from higher cortical areas, which could also have power in the alpha/beta range (Michalareas et al. 2016;van Kerkoerle et al. 2014), was unlikely to be prominent during this period. To obtain a comparable estimate of the steady-state neural activity, we also performed additional analyses in a "late" time range extending from 250 ms to 750 ms after stimulus onset, with baseline defined as the 500-ms time range preceding stimulus onset. The power spectral density of the LFP was similar if a smaller subinterval was used, although the tuning curves for firing rates were noisier for shorter intervals because the average firing rates reduced considerably after~200 ms of stimulus onset.
Spiking analysis. To obtain the PSTH for a spiking electrode, we binned spikes across each trial using a bin width of 10 ms and averaged across trials for each condition. For visualization , the PSTH for each electrode was smoothed with a Gaussian kernel (length ϭ 11, ϭ 10 ms), normalized by the maximum firing activity, and averaged across all spiking electrodes for each condition (see Fig. 2A).
Mean spike counts were obtained for each spiking electrode during baseline and stimulus time ranges for each condition and were separately normalized with respect to the maximum spike count during stimulus period across conditions. Finally, these normalized spike counts were averaged over spiking electrodes (see Fig. 2B).
To obtain a measure of the temporal frequency selectivity of spiking responses, we fit the drift rate responses of each spiking electrode, measured as the change in spike counts from baseline, to a difference of exponentials function (Hawken et al. 1996;Levitt et al. 1994) of the form where r(TF) is the change in spike counts from baseline as a function of drift rate, TF is the drift rate used (0, 1, 2, 4, 8, 16, 32, or 50 cps), and coefficients ␣ 1 , ␣ 2 , ␤ 1 , and ␤ 2 are constrained to be nonnegative. We used the MATLAB function lsqcurvefit to minimize the sum of squared differences between the data points and the fitted values and calculated the preferred temporal frequency of each multiunit (an example is shown in Fig. 2C) as the drift rate at which the spiking response peaked (using a resolution of 0.1 cps). To compute the suitability of the fit (that is, whether the fit matched the data), we calculated where r(TF) is the change in spike counts for each drift rate, r is the mean change in spike counts across all drift rates, and r(TF) is the fitted value corresponding to r(TF). Intuitively, qval signifies the fraction of variance in the data explained by the fit; values of qval close to 1 correspond to very good fits whereas values of qval close to 0 or negative values indicate poor fits to the data. Using the above approach, we obtained the preferred temporal frequency of spikes for both the early and late time periods and used them for characterizing the temporal frequency selectivity of the output of V1 (see Fig. 2D). If the preferred temporal frequency of a spiking electrode was Ͻ1 cps, we classified it as low pass (Fig. 2D); otherwise, it was classified as band pass. Note that statistics of preferred temporal frequency, wherever reported, are calculated based on the preferred temporal frequencies obtained from band-pass electrodes only.
Time-frequency analysis. We used the MP algorithm (Mallat and Zhang 1993) for estimating the time-frequency spectrograms. MP is a greedy, iterative algorithm that approximates a signal using a linear combination of waveforms (called atoms) chosen from an overcomplete dictionary. Briefly, every iteration of MP works on a residue (the as-yet unexplained part of the signal) to decompose it into a new residue and the product of a coefficient and an atom. The initial residue is initialized to the signal to be decomposed. The atom chosen in each iteration is the one that best explains the residue (that is, has the largest inner product with it), and the coefficient is the value of the inner product between the residue and the atom. Mathematical details of the procedure are presented elsewhere (Ray et al. 2008b). This method is particularly suitable to resolve fast transients in the signal, such as the activity related to the early stimulus period (see Chandran KS et al. 2016 for comparison with other methods).
We used an overcomplete dyadic dictionary of Gaussian-modulated sinusoidal functions (Gabor functions). For monkey 1, each stimulus trial had a length of 4,096 points (2,048-ms duration, 2-kHz sampling frequency). In case of monkey 2, we downsampled each stimulus trial by a factor of 2 before running MP, resulting in a trial length of 4,096 points and a sampling frequency of 1 kHz. Using MP, we approximated the signal in each stimulus trial using 500 atoms and reconstructed the single-trial energy spectrogram. The high-resolution spectrograms were appropriately downsampled to obtain a frequency resolution of~1 Hz and a time resolution of 4 ms for both monkeys. For each LFP electrode, we averaged the spectrograms across all trials for each drift rate condition. We obtained the change in power spectrograms ( Fig. 3) by subtracting the log of baseline power (mean over the 200-ms period before stimulus onset) from the log of the raw spectrogram and multiplied by 10 to get units in decibels (dB). Results using multitaper method (Babadi and Brown 2014;Jarvis and Mitra 2001) yielded almost similar results for the late stimulus period but failed to capture the early transient effectively (data not shown).
To obtain frequency band power responses for each drift rate condition, we summed the power values in each of alpha (8 -15 Hz), low-frequency (20 -40 Hz), gamma (50 -80 Hz), and high gamma (140 -250 Hz) ranges in the raw power spectrogram, separately for the early and late time ranges for each electrode. These summed power responses were separately normalized with respect to the maximum during stimulus period across conditions for each electrode and finally averaged over LFP electrodes (Figs. 4A and 5A). Note that in our definition of the high gamma range, we have chosen a higher range (140 -250 Hz) than usual to avoid a harmonic of gamma rhythm that is noticeable in both monkeys (Fig. 3).
We correlated the power in each frequency band across all LFP electrodes with the average firing rate of all spiking electrodes (Figs. 4A and 5A), using it as a proxy for the total output in V1 (Fig. 4A, early; Fig. 5A, late). Specifically, we computed Spearman's rank correlation coefficient () between the normalized firing rate for different drift rates (averaged over all spiking electrodes, n ϭ 29 and 19 for the 2 monkeys) and normalized frequency band power for different drift rates (averaged over all LFP electrodes, n ϭ 78 and 39 for the 2 monkeys), separately for the early and late time ranges. Both the correlation and P values were computed with the function corr in MATLAB, which computes the P value using the exact permutation distributions. Importantly, it should be noted that even though the spectrograms shown in Fig. 3 are baseline subtracted for ease of visualization, we compute the correlations in Figs. 4A and 5A from normalized spiking activity and frequency band power obtained directly during the stimulus period, without any baseline correction.
To directly compare the temporal frequency selectivity of power in different LFP frequency bands with the temporal frequency selectivity of V1 spiking activity, difference of exponentials fit (Eq. 1) and associated quality of fit (Eq. 2) were computed for each LFP electrode by taking r(TF) as the change in power (in dB) from baseline for the frequency range of interest. This was done separately for each frequency range during both early and late epochs, obtaining distributions of preferred temporal frequency values (Fig. 4C, early; Fig. 5C, late). Figure 1A shows the raster plots and PSTH traces for a multiunit in our data (electrode 7, monkey 1). Firing rates were transiently elevated within the first 100 ms after stimulus onset before settling down at a steady lower value. This multiunit had a mean firing rate increase of~29 spikes/s in the stimulus period (defined between 50 and 250 ms after stimulus onset) to the static grating and showed a shallow increase in firing with increasing drift rate of the stimulus, before falling off at higher drift rates. Figure 1B shows a couple of raw LFP traces and the average LFP response. Gratings drifting at low/medium drift rates induced large gamma oscillations that were clearly visible in the individual LFP traces. However, because gamma oscillations in individual trials were not phase locked to the stimulus, the average (evoked) response showed a much weaker oscillatory pattern in the gamma range. Figure 1, C and D, show the raw time-frequency spectrograms and change in power spectrograms, respectively, for the same electrode as in Fig. 1A, constructed with the MP algorithm and averaged across trials for each stimulus condition. There was a transient increase in low-frequency power below~50 Hz, corresponding in time to the transient elevation in firing rate within the first 100 ms after stimulus onset, which was preserved across all stimulus drift rates and was evident in both the raw and baseline-subtracted spectrograms. Figure 2A shows the normalized firing rate PSTH across the population of spiking electrodes, while Fig. 2B shows the mean normalized spike counts across all excited spiking electrodes (n ϭ 29 and 19) for the two monkeys during the early and late stimulus period. For both analysis epochs, preferred drift rates were in the low to mid-range values. Similar trends were observed if the F0 or F1 component of the spiking responses (see MATERIALS AND METHODS for details) was used instead of the total firing rate.

RESULTS
To quantify the temporal frequency tuning of the spiking activity, we fitted a difference of exponentials function to the change in spike counts across conditions for each spiking electrode and used the fits to quantify the preferred temporal frequency across the population (see MATERIALS AND METHODS for details). Figure 2C shows the spike counts and the difference of exponentials fit for the same example multiunit reported in Fig. 1A. This multiunit had a mean preferred drift rate of~14 cps during the early period and~12 cps during the late Fig. 1. Spiking activity and time-frequency spectrograms for a representative electrode (electrode 7, monkey 1). A: raster plots showing spiking activity in individual trials for each drift rate. Ticks in trials for each raster plot have been scaled to obtain the same overall vertical scale. Trial-averaged peristimulus time histogram is overlaid on individual rasters. cps, Cycles per second. B: local field potential traces for 2 example trials (light and dark gray) and the evoked response (black) for each condition. C: time-frequency raw power plots constructed with matching pursuit (MP) for each drift rate, averaged across trials for that condition. D: time-frequency change in power plots showing the difference in power from baseline (200-ms period before stimulus onset), constructed with MP for each condition, averaged across trials for that condition.
period. Across the population of spiking electrodes (Fig. 2D), consistent with earlier reports in V1 to drifting gratings (Foster et al. 1985;Hawken et al. 1996), most multiunits showed a shallow band-pass tuning of firing rate to stimulus drift rate (monkey 1: 23 and 26 of 29 multiunits were band pass in early and late epochs, respectively; monkey 2: 16 and 14 of 19 multiunits). For these band-pass multiunits, the mean preferred temporal frequency was significantly higher in monkey 1 than in monkey 2 during the early [monkey 1: 11.9 Ϯ 1.0 cps, 95% confidence interval (CI) (9.8 cps, 14.0 cps); monkey 2: 9.0 Ϯ 1.1 cps, 95% CI (6.7 cps, 11.3 cps); P ϭ 0.033, 2-sample t-test, single tailed] but not the late [monkey 1: 11.2 Ϯ 1.6 cps, 95% CI (8.0 cps, 14.4 cps); monkey 2: 10.0 Ϯ 1.4 cps, 95% CI (7.1 cps, 13.0 cps); P ϭ 0.31, 2-sample t-test, single tailed] period. Similar results were obtained if medians instead of means (and Wilcoxon rank sum test instead of t-test) were used. Figure 3 shows the time-frequency change in power spectrum of the LFP across the population of LFP electrodes. Different frequency bands indeed showed very different trends with varying drift rates. Consistent with previous studies (Manning et al. 2009;Ray et al. 2008a;Ray and Maunsell 2011;Whittingstall and Logothetis 2009), higher frequencies of the LFP (above~140 Hz; chosen to avoid the harmonic of gamma that is present at~100 -120 Hz) showed a trend similar to the population firing rate, peaking at intermediate and rapidly decaying at high drift rates (also readily visible in the repre-sentative site shown in Fig. 1D). Gamma rhythm (~50 -80 Hz in our data) showed an increase in center frequency with increasing drift rates, consistent with previous results in cats Viana Di Prisco 1997), primates (Friedman-Hill et al. 2000), and humans (Orekhova et al. 2015). Importantly, at high drift rates, even though spiking activity, gamma rhythm, and high gamma activity were all substantially reduced, we observed a strong transient activity between~20 and 40 Hz, peaking at~75 ms after stimulus onset. This initial transient was spectrally quite broadband and coincided in time with the transient elevation in firing rate ( Fig. 2A), but the higherfrequency content of this broadband decreased at higher drift rates, unlike the low-frequency power.
To quantify these results, we computed the total power at alpha (8 -15 Hz), low (20 -40 Hz), gamma (50 -80 Hz), and high gamma (140 -250 Hz) bands during the early stimulus period and compared the normalized power values with the normalized firing rates during the same period (Fig. 4A). The high gamma power as a function of drift rate showed an inverted U-shape similar to the firing rate trace (Fig. 2B), and the two curves were significantly correlated (Fig. 4A). Gamma band power peaked at lower temporal frequencies compared with the average firing rates and showed modest correlations (Fig. 4A). Low-frequency power showed inconsistent trends, with the first monkey showing an increase in power at high drift rates while the second monkey showed almost no change in power with increasing drift rates (Fig. 4A). Alpha band power, which captured part of the transient in the early stimulus period, also showed inconsistent trends (Fig. 4A).
Low correlations of power in a particular frequency band with spiking activity could result either because the power itself is untuned to temporal frequency or because the power does have a temporal frequency tuning but one that is distinct from the V1 output tuning. For example, it is clear that the low correlation of gamma band power with spiking activity is due to their differential tuning. To gain quantitative insights into the temporal frequency tuning of LFP frequencies, we fitted difference of exponentials functions to the change in power (in dB) in the different frequency bands for each LFP electrode, in the same way as for spiking responses. From the fits for each LFP electrode, we computed the preferred temporal frequency for each of the frequency bands. Figure 4B summarizes the fit quality for each of the frequency ranges during the early period, and Fig. 4C shows the estimated preferred temporal frequency values. For comparison, we also show the fit quality and preferred temporal frequencies for spikes (same data as Fig. 2D but shown in a different format). For the high gamma range, the fit quality varied considerably across the population, possibly because the power in this band was low and noisy at some sites, although the average fit quality was decent ( Fig. 4B; monkey 1: 64 Ϯ 3%; monkey 2: 70 Ϯ 3%). Most of the sites were classified as band pass (monkey 1: 65 of 78, monkey 2: 36 of 39), and the mean preferred temporal frequency of these band-pass electrodes was 10.2 Ϯ 0.6 cps (n ϭ 65, 95% CI [9.0 cps, 11.4 cps]) in monkey 1 and 6.0 Ϯ 0.5 cps (n ϭ 36, 95% CI [4.9 cps, 7.0 cps]) in monkey 2, similar to the preferred temporal frequency for spikes (Fig. 2D).
Lower frequencies between 20 and 40 Hz showed more complex dynamics (Fig. 4A). In monkey 1, most electrodes showed almost no change in power at low to intermediate temporal frequencies and a noticeable increase at higher temporal frequencies. Although the resulting mean fit quality was good (Fig. 4B), we could not obtain values for the preferred temporal frequency (within 100 cps) for a majority (77 of 78) of the electrodes (the only band-pass electrode for this monkey showed a preferred temporal frequency of 39.5 cps). In monkey 2, the trend in low-frequency power was noisier and the mean fit quality was lower, with about one-third of the electrodes (13 of 39) showing a band-pass behavior. The mean preferred temporal frequency of these band-pass electrodes was 14.1 Ϯ 3.2 cps (n ϭ 13, 95% CI [7.1 cps, 21.1 cps]), substantially higher than the corresponding values for spikes and high gamma power.
Fit quality was the best in the gamma range for both monkeys ( Fig. 4B; monkey 1: 90 Ϯ 1%; monkey 2: 82 Ϯ 2%), with the majority of electrodes showing low-pass behavior (monkey 1: 44 of 78; monkey 2: 31 of 39; no electrodes showed high-pass behavior). The preferred temporal frequency of the remaining band-pass electrodes was 2.1 Ϯ 0.2 cps (n ϭ 34, 95% CI [1.8 cps, 2.4 cps]) in monkey 1 and 2.0 Ϯ 0.4 cps (n ϭ 8, 95% CI [1.1 cps, 2.8 cps]) in monkey 2, much lower than the preferred temporal frequency values for high gamma, spikes, or low frequencies. Finally, the alpha range showed large variability in the fit quality for both monkeys ( Fig. 4B; monkey 1: 63 Ϯ 3%; monkey 2: 46 Ϯ 3%), with a considerable fraction of electrodes showing a band-pass behavior (monkey 1: 54 of 78, monkey 2: 36 of 39). The mean preferred temporal frequency of the alpha band was 7.1 Ϯ 0.5 cps (n ϭ 54, 95% CI [6.0 cps, 8.2 cps]) in monkey 1, which was lower than the mean For clarity, only the time range from Ϫ250 to 500 ms (0 denotes stimulus onset) has been shown. Thick vertical bars on the y-axis of the leftmost plots indicate the alpha (8 -15 Hz), low-frequency (20 -40 Hz), gamma (50 -80 Hz), and high gamma (140 -250 Hz) ranges used for analysis. Note that monkey 1 showed a higher change in power compared with monkey 2, indicated by the different limits for the color axes. cps, Cycles per second. preferred drift rate of low frequencies, whereas in monkey 2, the mean preferred temporal frequency of alpha was 13.2 Ϯ 0.6 cps (n ϭ 36, 95% CI [12.0 cps, 14.4 cps]), comparable to that preferred by low frequencies.
To test whether our results depended on the choice of the analysis period, we repeated the analysis after computing power between 250 and 750 ms (the "late" stimulus period; Fig. 5). High gamma fit quality (Fig. 5B) was comparable to that during the early period, and a majority of the sites were band pass ( Fig. 5C; monkey 1: 71 of 78, monkey 2: 35 of 39). The mean preferred temporal frequency was 9.5 Ϯ 0.6 cps (n ϭ 71, 95% CI [8.4 cps, 10.7 cps]) in monkey 1 and 9.0 Ϯ 1.1 cps (n ϭ 35, 95% CI [6.7 cps, 11.3 cps]) in monkey 2, again in the range of preferred temporal frequency for spiking activity during the same period.
Unlike the high gamma frequencies, which had a clear tuning peak, the lower frequencies neither had a proper low pass or shallow band pass structure nor had a clearly defined peak in most electrodes. This is not surprising, because the LFP power was itself very low (sometimes below baseline levels) during this late stimulus period. There was a pronounced suppression of power at low to intermediate drift rates, which transitioned to an increase in power at still higher drift rates and a tendency to fall off again at the highest drift rate in many electrodes (Fig. 5A). Not surprisingly, the temporal frequency tuning curves in many electrodes did not fit the change in low-frequency power well, and the tuning fit quality was considerably low (Fig. 5B; monkey 1: 58 Ϯ 2%; monkey 2: 24 Ϯ 3%). In the absence of a single discernible peak, at some sites the difference of exponentials had a bias toward more closely fitting the downward-going trend at the lower end and the preferred temporal frequency estimated from these fits was very low; these sites were classified as low pass (monkey 1: 5 of 78, monkey 2: 26 of 39). A higher fraction of sites showed high-pass behavior for monkey 1 (68 of 78) than for monkey 2 (2 of 39). The mean preferred temporal frequency of the remaining band-pass electrodes was 55.5 Ϯ 15.2 cps (n ϭ 5, 95% CI [13.3 cps, 97.8 cps]) in monkey 1 and 15.7 Ϯ 3.0 cps (n ϭ 11, 95% CI [9.0 cps, 22.3 cps]) in monkey 2. Since the power itself was very low during this analysis period, not much can be inferred from these fits.

DISCUSSION
We presented drifting grating stimuli, tuned to elicit strong gamma responses, to awake fixating macaque monkeys over a wide range of drift rates and characterized the temporal frequency tuning of spiking activity as well as LFP power in  Fig. 1. cps, Cycles per second. B: temporal frequency tuning fit quality (see MATERIALS AND METHODS) for the 4 frequency bands in monkey 1 (light gray) and monkey 2 (dark gray). Error bars indicate SE. Fit quality for individual electrodes (small magenta circles) is overlaid on the bar plots. C: preferred drift rates estimated from the temporal frequency tuning for various frequency bands. Error bars indicate SE. Same convention as in B. Electrodes classified as low pass and electrodes for which preferred drift rates could not be estimated within 100 cps are shown separately (small gray circles), below and above each bar plot, respectively. different frequency ranges. Consistent with earlier studies, there was a reduction in V1 spiking activity at higher drift rates, with mean preferred drift rates around~10 cps. However, we found that power in different frequency bands of the LFP was tuned differently to stimulus drift rates. Consistent with earlier studies, power changes in the high gamma range (140 -250 Hz) correlated strongly with spiking activity and were tuned to similar drift rates, reflecting a signature of output. Gamma-band (50 -80 Hz) center frequency varied systematically with drift rate, as in previous studies, and gamma-band power was tuned to very low drift rates, especially during the sustained stimulus epoch. Significantly, however, lower frequencies of the LFP, up to~40 Hz, were relatively invariant in their tuning to stimulus drift rates and appeared to prefer higher drift rates than spikes or high gamma frequencies, especially in the early stimulus epoch. Although we do not have a direct measure of the input LGN spiking activity in our study (ideally, the average LGN firing rate tuning profile in the transient regime should have matched the transient~20 -40 Hz power tuning profile), the early timing and an almost equally strong response at high drift rates suggest that this could be related to the incoming input. Beyond the initial transient, however, the low-frequency LFP power did not appear to be consistent with the known dynamics of firing rates in LGN or retina. Specifically, the low-frequency LFP power reduced to baseline levels or even lower after~125 ms (a fall of~15 dB from the power at~75 ms), even though no such reduction was observed in V1 firing in our data or LGN firing in previous studies (Hawken et al. 1996). As described below, this could be due to lateral/ feedback inputs that could influence LFP power in the low frequency range. Overall, these results add to the growing body of literature on the relationship between spikes, LFP, and gamma oscillations in V1 (Bastos et al. 2014;Belitski et al. 2008;Burns et al. 2010;Frien et al. 2000;Lashgari et al. 2012;Williams et al. 2004).
Low-frequency power and synaptic inputs. LFPs have traditionally been linked to synaptic activity, which are the inputs, mostly based on theoretical arguments and the fact that action potentials are not strongly related to the ongoing LFP activity (Mitzdorf 1985(Mitzdorf , 1987Nunez and Srinivasan 2006), although this link is rather simplistic, since LFPs also represent a host of nonsynaptic events . Earlier work toward understanding the composition of LFP has addressed this partly through experiments using techniques such as information theory for quantification of the LFP frequency content (Belitski et al. 2008) or using computational studies (Einevoll et al. 2013;Reimann et al. 2013).
A recent study used an input-output dissociation paradigm in the middle temporal area (MT) and the medial superior temporal area (MST) by using moving plaid stimuli, which are formed by superposing two gratings drifting at different orientations (Khawaja et al. 2009). While V1 cortex responds primarily to the individual directions of the two gratings (component selectivity), many MT and MST neurons, which receive input from V1, respond to the overall drift direction of the plaid (pattern selectivity; see Khawaja et al. 2009 for details). Khawaja and colleagues showed that the lower frequency bands of the LFP exhibited signatures of component selectivity similar to their afferent input, whereas the high gamma band showed pattern selectivity like the output spiking activity, consistent with the findings of the present study. Our results extend their work in three ways. First, drifting gratings allowed precise manipulation of gamma rhythm to show that it is distinct from input as well as output. The plaid stimuli used by Khawaja and colleagues did not generate salient gamma oscillations in MT or MST, similar to what has been reported in V1 (Lima et al. 2009). Consequently, these stimuli did not manipulate the gamma rhythm differently from the output, such that the results for the gamma band were similar to high gamma, albeit the correlation was weaker. Second, the population response in LGN and V1 to gratings at varying drift rates is likely to be more homogeneous than the response to plaids moving in varying orientations. More specifically, different neurons in MT and MST are likely to be tuned to different orientations, each receiving tuned input from V1. Given that LFP has a spatial spread of at least~250 m (Dubey and Ray 2016;Lindén et al. 2011;Xing et al. 2009), the summed input of neurons in this area may be less tuned for any particular orientation. On the other hand, in response to gratings at varying drift rates, almost all cells in LGN and V1 are low or band pass, such that the dissociation between input and output at high drift rate will be less influenced by the spatial averaging of the LFP. Finally, use of MP allowed us to parse out the signature of input at a finer time resolution than previous studies, although when averaged over suitable time and frequency ranges power computed with multitaper method also yielded comparable results.
Our approach was motivated by the work of Viswanathan and Freeman (2007), who showed that fMRI BOLD responses primarily reflected synaptic rather than spiking activity. Our work extends their results in several ways. First, they did not focus on different frequencies in the LFP, since their main goal was to relate the inputs to the BOLD signal. Second, since they used only two drift rates (3 and 30 Hz), they could not quantitatively characterize the tuning preferences of spiking activity or different frequency bands of the LFP as done here. Most importantly, the grating stimuli used by Viswanathan and Freeman did not generate salient gamma oscillations, presumably because of the use of very low spatial frequency.
Gamma rhythm and LGN inputs. Recent reports (Saleem et al. 2017;Storchi et al. 2017;Veit et al. 2017) have highlighted a narrow-band gamma rhythm in mouse V1 (in addition to the more traditional stimulus-induced gamma rhythm), whose power and frequency are primarily modulated by luminance (Saleem et al. 2017;Storchi et al. 2017) as well as behavioral state, including locomotion (Niell and Stryker 2010). More importantly, this narrow-band gamma rhythm has been shown to reflect excitatory currents of subcortical origin (Saleem et al. 2017) and therefore directly related to inputs coming into mouse V1. Therefore, it becomes important to examine the gamma rhythm that we observe in light of these new reports.
The gamma rhythm in primates differs from this narrowband gamma in rodents in several ways. First, the luminancemodulated narrow-band gamma in mice is present at baseline and suppressed by stimulus contrast (Saleem et al. 2017), whereas the primate gamma rhythm is weak or absent during spontaneous periods and becomes stronger with stimulus contrast (Henrie and Shapley 2005;Ray and Maunsell 2010). Second, the mouse narrow-band gamma was shown to be concentrated in the input layers of V1 (Saleem et al. 2017), whereas the gamma rhythm in primate V1 is more prominent in the corticocortical output layers (Xing et al. 2012). Third, the mouse narrow-band gamma was shown to be mediated by excitatory currents (Saleem et al. 2017) and inherited from subcortical sources (Saleem et al. 2017;Storchi et al. 2017), whereas the stimulus-induced gamma rhythm in sensory cortices is thought to be mediated by rhythmic inhibition from interneurons to principal cells (Bartos et al. 2007; Wang 2012) and therefore of cortical origin (Whittington et al. 2011). Finally, the luminance-modulated gamma in mice has a narrow bandwidth of~2-5 Hz (Saleem et al. 2017), whereas the stimulus-induced gamma is broader, with a bandwidth of 20 -25 Hz in primates (Hadjipapas et al. 2015;Ray and Maunsell 2010).
Our results clearly indicate that the gamma rhythm in primate V1 is not directly related only to the inputs, extending previous work that has shown that gamma rhythm has no consistent relationship with spiking activity either (Jia et al. 2013;Ray and Maunsell 2011). Instead, our results are consistent with the hypothesis that gamma rhythm depends on excitation-inhibition interactions (Atallah and Scanziani 2009;Brunel and Wang 2003;Buzsáki and Wang 2012), and therefore on complex interactions between different types of inputs coming into a cortical area, although it still remains unclear why gamma power reduces at midlevel drift rates when both putative inputs and outputs are still high.
Like temporal frequency, gratings with varying orientations can also be used to dissociate inputs from outputs, since the LGN inputs are thought to be untuned for orientation (Priebe 2016). We recently studied the dependence of different frequency bands of the LFP on orientation for the two monkeys used in this study (Murty et al. 2018) and observed relatively lower orientation selectivity at low frequencies, which is consistent with earlier results in primates (Jia et al. 2011(Jia et al. , 2013 and in line with the hypothesis that low frequencies represent LGN inputs. Furthermore, the gamma rhythm across sites showed a peculiar tuning preference toward a single orientation, even though spiking activity across sites prefers a wide range of orientations, consistent with earlier primate studies (Berens et al. 2008;Jia et al. 2011), again suggesting that the primate gamma rhythm is not directly related only to the inputs or the outputs. Note that while these results appear different from a study in cats (Katzner et al. 2009), for which even the low frequencies appear to be tuned to orientation (see their Fig.  1E and Fig. 4A/C), this could be just because they used orientation noise stimuli in which grating stimuli of varying spatial frequency and orientations were presented only for 32 ms, to highlight the evoked response.
Importance of analysis interval. The analysis interval was chosen to be short (50 -250 ms) to avoid the influence of lateral and top-down feedback, which could be in the alpha/beta range (Michalareas et al. 2016;van Kerkoerle et al. 2014). Note that because of the use of MP the frequency resolution does not suffer due to the short analysis interval in a traditional way, because the actual decomposition is multiscale, such that the available basis functions have a wide span of time and frequency support (see Chandran KS et al. 2016 for a detailed discussion). For example, although the initial transient increase in LFP power at~75 ms is broadband and appears to spill over to gamma range and beyond, it is still easily distinguishable from the actual gamma rhythm that starts after~100 ms, such that the tuning of low frequency and gamma band computed within 50 -250 ms is remarkably different (compare Fig. 4A, low vs. ␥).
Although the LGN cells are expected to continue to fire throughout the stimulus period, the low-frequency power actually reduced, even below spontaneous levels, at longer stimulus intervals, which appears to be inconsistent with known dynamics of firing activity in the retina, LGN, and cortex. One reason for the low-frequency power dynamics is that the LGN neurons may have higher firing at stimulus onset and a subsequent reduction during the steady state, similar to the firing of the V1 multiunits ( Fig. 2A), and therefore the input may have been more dominant during the onset. In addition, there could be stronger top-down and lateral feedback as well as suppression of the alpha rhythm during the later period, because of which the V1 LFP would transition to more intricate dynamics. The late LFP signal would then possibly exhibit features that reflect microcircuit-level computations, making its interpretation more challenging. Our study does not address this directly; indeed, we think that it will be extremely challenging to tease out influences of afferent activity from lateral as well as neuromodulatory and top-down feedback connections, some or all of which may produce spectrally overlapped signatures.
One issue with the analysis interval very close to the stimulus onset is that the sudden onset causes a leakage of the stimulus energy over a broad frequency range, so that the brain may be responding to this onset differently from the sustained period for which the stimulus energy is localized near the drift rate. However, in that case we would not expect the firing rates in V1 at the initial onset to follow the type of tuning that is observed during the steady-state period, but they were almost identical (Fig. 2B), suggesting that the differences in the distribution of stimulus energy during early and late periods may not play a big role. Because the band-pass tuning of firing rates in V1 is observed very early on before lateral/feedback inputs may become prominent, it is likely that this tuning is inherited from LGN, although comparable recordings from LGN neurons would be needed to completely settle this issue.
Limitations of study. As described above, a limitation of our study is that we do not have an objective measure of the input activity, which would have required comparable recordings from LGN as well. It is difficult to infer the average LGN firing profile from earlier studies that have recorded from LGN (Hawken et al. 1996) or to compare our V1 firing rate profiles from previous studies (Foster et al. 1985;Hawken et al. 1996). First, these studies were done in anesthetized animals, for which the temporal frequency tuning could be different (Alitto et al. 2011). Second, we used the same orientation and spatial frequency for all cells, unlike previous work where such stimulus parameters were optimized for each cell (as explained in MATERIALS AND METHODS, this was done to ensure that the stimuli generate salient gamma rhythm). However, it is unlikely that the main results shown here would depend drastically on the choice of stimulus parameters. For example, the low-frequency component of the LFP was largely untuned for orientation as well as spatial frequency (see Fig. 1 and Fig.  3-1B of Murty et al. 2018), and therefore we would expect the LFP responses at low frequencies to be similar irrespective of whether the choice of parameters was optimal for the multiunit recorded from the same site or not. Similarly, for the two orientations tested here (0°and 90°), gamma power was tuned differently from spikes for varying drift rates, so it is likely that the tuning will remain dissimilar even if the orientation that was optimum for the multiunit was used. Only for the high gamma range might we expect to see an improvement in the correlation between the tuning curves for spikes and LFP power. However, given that the correlations were high even with stimuli that were nonoptimal for the multiunits (0.83 and 0.95 in the early and late stimulus period for monkey 1, 0.9 and 0.93 for monkey 2; Figs. 4A and 5A), only a minor improvement, if at all, is possible.
Another issue is related to volume conduction, especially if it depends on the frequency of the LFP Schroeder 2011, 2015). For example, consider a hypothetical case where both lower and higher frequencies are tuned to orientation, with neighboring brain regions preferring different orientations. If lower frequencies spread more than high frequencies because of volume conduction, low frequencies will appear untuned to orientation just because of greater spatial averaging across different orientations, while high frequencies will remain tuned. However, while this type of frequency-dependent volume conduction can partially explain our results related to orientation, for which well-defined maps exist in V1 (Ohki et al. 2006), it is unlikely to explain our results related to temporal frequencies. First, temporal frequency maps in V1 are less distinct, with at best some evidence of clustering of high temporal frequency-preferring regions scattered among low temporal frequency-preferring regions (Shoham et al. 1997;Sun et al. 2007). More importantly, even if different regions of V1 were tuned to different temporal frequencies, spatial averaging due to volume conduction could reduce the overall tuning preference but not yield a preferred frequency that was much higher (as was the case for low LFP frequencies) or much lower (as was the case for gamma range) than the average preferred temporal frequency of individual sites. However, more experiments, especially with laminar probes, are needed to fully account for the effects of volume conduction.
Although the major results shown here were consistent across the two monkeys, there were interesting animal-specific differences. In particular, the low-frequency power increased at high drift rates for monkey 1 but not for monkey 2. It is unclear why the results are different, although it could be because the microelectrodes were in different cortical layers or the receptive fields had different eccentricities or because of possible differences in the position of the reference wire, although how these factors could influence low-frequency power at high drift rates is unclear. We observed that the firing rates remained relatively higher at 32 cps and 50 cps for monkey 1, which could be simply because this animal had a higher mean preferred temporal frequency for spiking activity and, by extension, a higher low-frequency drift rate cutoff than monkey 2 in both LGN and V1. This could explain why the normalized 20 -40 Hz power for the 32 cps and 50 cps drift rates in monkey 1 appeared to deviate from the values at lower drift rates, driving correlations with averaged firing rates to negative values that would otherwise have been close to 0. Another possibility for why the 50 Hz drift rate generated a stronger response in the LFP could be because this rate is actually half of the monitor refresh rate and therefore takes the highest and lowest luminance values in alternate frames, although why this would affect the power in only the first monkey is not clear. Importantly, despite these animal-specific differences, the major results shown here-presence of a low-frequency transient across temporal frequencies, low correlation between lowfrequency power and firing rates, peak frequency shift in gamma band with drift rate, and high correlation between high gamma power and firing rates-are all consistent across the two monkeys.
The hypothesis that distinct zones of measured neural activity could carry signatures of different cortical processes is attractive from a computational standpoint as well as for brain-computer interface (BCI) research. In the last few years, the LFP has become popular as a candidate neural signal for BCI systems because of its rich frequency content and relative stability across time (Andersen et al. 2014;Hwang and Andersen 2013). The present study represents another step toward characterizing the LFP signal in terms of its frequency ranges and, together with earlier work (Khawaja et al. 2009), shows that the behavior and interpretation of the LFP signal may be surprisingly similar across different visual areas. Our results agree with earlier findings on high gamma being a reliable marker of cortical output and further suggest that BCI applications decoding sensory activity should focus on low to midrange frequencies in the early, transient regime of the LFP to extract information that is different from the output.

GRANTS
This work was supported by the Wellcome Trust/Department of Biotechnology (DBT) India Alliance (500145/Z/09/Z; Intermediate fellowship to S. Ray), a Tata Trusts Grant, and the DBT-IISc Partnership Programme.

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