Abstract
Populations of cortical neurons generate rhythmic fluctuations in their ongoing spontaneous activity. These fluctuations can be seen in the local field potential (LFP), which reflects summed return currents from synaptic activity in the local population near a recording electrode. The LFP is spectrally broad, and many researchers view this breadth as containing many narrowband oscillatory components that may have distinct functional roles. This view is supported by the observation that the phase of narrowband oscillations is often correlated with cortical excitability and can relate to the timing of spiking activity and the fidelity of sensory evoked responses. Accordingly, researchers commonly tune in to these channels by narrowband filtering the LFP. Alternatively, neural activity may be fundamentally broadband and composed of transient, nonstationary rhythms that are difficult to approximate as oscillations. In this view, the instantaneous state of the broad ensemble relates directly to the excitability of the local population with no particular allegiance to any frequency band. To test between these alternatives, we asked whether the spiking activity of neocortical neurons in marmoset of either sex is better aligned with the phase of the LFP within narrow frequency bands or with a broadband measure. We find that the phase of broadband LFP fluctuations provides a better predictor of spike timing than the phase after filtering in narrow bands. These results challenge the view of the neocortex as a system composed of narrowband oscillators and supports a view in which neural activity fluctuations are intrinsically broadband.
SIGNIFICANCE STATEMENT Research into the dynamical state of neural populations often attributes unique significance to the state of narrowband oscillatory components. However, rhythmic fluctuations in cortical activity are nonstationary and broad spectrum. We find that the timing of spontaneous spiking activity is better captured by the state of broadband fluctuations over any latent oscillatory component. These results suggest narrowband interpretations of rhythmic population activity may be limited, and broader representations may provide higher fidelity in describing moment-to-moment fluctuations in cortical activity.
Introduction
Since the first human electroencephalogram (EEG) recordings by Hans Berger (Berger, 1929), neuroscientists have inferred cortical function from the state of rhythmic fluctuations in neural population activity (Buzsáki and Draguhn, 2004; Wang, 2010). These brain rhythms are believed to arise from return currents generated by large-scale spiking activity in cortical neural populations (Logothetis, 2003; Katzner et al., 2009; Buzsáki et al., 2012). When recorded intracranially with penetrating electrodes, rhythmic activity can be measured in the local field potential (LFP), which typically reflects neural signals arising within ∼250 µm of the electrode tip (Katzner et al., 2009; Lindén et al., 2011). LFP fluctuations are spectrally broad but are often thought to be composed of activity in narrow frequency bands correlated with distinct neural functions (Canolty et al., 2010; Einevoll et al., 2013; Friston et al., 2015). For example, in the visual cortex, alpha band rhythms (8–15 Hz) are thought to reflect feedback processes of suppression (Jensen and Mazaheri, 2010; van Kerkoerle et al., 2014) and have been shown to be attenuated with or modulated by attention (Worden et al., 2000; Busch and VanRullen, 2010). Beta band rhythms (15–30 Hz) have been linked to motor planning (Sanes and Donoghue, 1993; Rubino et al., 2006) and feedback regulation of excitability (Bastos et al., 2015; Friston et al., 2015). Theta band (4–8 Hz) activity has been related to attention (Fiebelkorn and Kastner, 2019), working memory load (Jensen and Tesche, 2002), and hippocampal function (Buzsáki, 2002). Delta band (<4 Hz) activity has been related to sleep and states of arousal (Sanes and Donoghue, 1993; Steriade et al., 2001; McGinley et al., 2015). Higher frequency gamma activity (30–90 Hz) has been linked to local coordination in excitation and inhibition (Brunel and Wang, 2003; Bartos et al., 2007; Buzsáki and Wang, 2012), attention (Fries et al., 2001, 2008; Gregoriou et al., 2009), memory (Pesaran et al., 2002; Colgin et al., 2009; van Vugt et al., 2010; Lundqvist et al., 2018), and perception (Singer and Gray, 1995; Panagiotaropoulos et al., 2012; Misselhorn et al., 2019), and has been used as a surrogate for measuring cortical activation (Crone et al., 2006; Ray et al., 2008a; Merker, 2013). Oscillatory activity can be induced under certain conditions, such as the increased low-frequency power that is observed in the EEG when eyes are closed (Berger, 1929; Geller et al., 2014), optogenetically (Lu et al., 2015; Bitzenhofer et al., 2017; Zutshi et al., 2018), electrically (Contreras et al., 1997; Kirov et al., 2009; Escobar Sanabria et al., 2020), or pharmacologically as in the alpha oscillations that occur in medial prefrontal cortex under propofol-induced anesthesia (Purdon et al., 2013; Flores et al., 2017; Bastos et al., 2021).
It has been proposed that certain frequency bands play a privileged role in routing information among brain areas (Akam and Kullmann, 2010; Bonnefond et al., 2017; Khamechian et al., 2019). The idea that communication between brain areas occurs through oscillatory processes within narrow frequency bands bears similarity to a radio, where signals are broadcast within different frequency bands and a receiver can be tuned to receive them (Hoppensteadt and Izhikevich, 1998). For example, the idea of cross-cortical communication through coherence views synchrony in gamma oscillations as periods of coordination between presynaptic and postsynaptic groups to transmit signals about, for example, an attended stimulus while blocking competing inputs (Fries, 2015). These patterns of gamma-band synchronization are proposed to be regulated across cortical areas by top-down signals within a slower (8–20 Hz) frequency band (Bastos et al., 2015). Other theories suggest that the LFP is composed of multiplexed oscillatory neural signals that are separate streams of information processing (Lisman and Idiart, 1995; Panzeri et al., 2010; Akam and Kullmann, 2014; Tingley et al., 2018). If oscillatory activity in separate frequencies encodes distinct information channels, and the spiking activity of neurons are the fundamental units of information transmission in the nervous system, then the spiking activity of individual neurons should show preferential alignment of their spiking activity to oscillatory rhythms to tune in to a channel of information (Canolty et al., 2010; Belluscio et al., 2012). There is evidence to suggest this can occur, as spikes have been found to preferentially align with the phase of theta (Takahashi et al., 2014; Souza and Tort, 2017; Strüber et al., 2022), alpha (Haegens et al., 2011), gamma (Fries et al., 2001; Womelsdorf et al., 2007; Ray et al., 2008b), and beta (Donoghue et al., 1998; Canolty et al., 2010) frequencies.
An alternative view is that neurons spike with no preference for any particular narrowband frequency. Rather, spiking is modulated by the instantaneous state of fluctuations in the local population, which varies from moment to moment across a broad range of frequencies. Supporting this view is the observation that balanced excitation and inhibition creates fluctuating neural activity patterns in the awake state, which often exhibit 1/fɑ power spectra across a broad range of frequencies (Destexhe et al., 2001; Gao et al., 2017). Studies in humans have found that changes in cognitive state are associated with broad spectral changes in the EEG (Voytek et al., 2015). The membrane potential of individual neurons is correlated with the population fluctuations measured in the instantaneous LFP (Haider et al., 2016), as opposed to any narrowband component, which suggests the broadband LFP is therefore informative about the instantaneous excitability of neurons in the population (Davis et al., 2020). Accordingly, previous work has found that spikes are weakly coupled to all frequencies of the broadband LFP (Martin and Schröder, 2016), and specific interactions in narrowband frequencies may at times be because of spurious artifacts from narrowband filtering (Scheffer-Teixeira and Tort, 2016).
Even when approximately oscillatory activity may be transiently apparent in LFP recordings, it is difficult to describe the phase of neural fluctuations within a narrow range of frequencies because of their nonstationarity (Pesaran et al., 2018). LFP phase is a useful measure for tracking the state of neural fluctuations because it is indicative of the relative transition in the balance of excitation and inhibition with, for example, the falling phase reflecting a transition from inhibition to excitation, and the rising phase transitioning from excitation to inhibition (Atallah and Scanziani, 2009; Poo and Isaacson, 2009; Isaacson and Scanziani, 2011; Teleńczuk et al., 2017). This is in contrast to amplitude measures, which can be ambiguous as the same negative voltage value could reflect neurons becoming more depolarized or more hyperpolarized depending on the signal history. Under this view, one can better characterize the state of neural populations from the phase of broadband fluctuations in LFP activity, and neurons will show preferential alignment of their spiking activity to the broadband signal phase, not to any narrowband oscillatory phase.
To ask whether neuronal spiking is better coupled to narrowband oscillations or broadband fluctuations during waking visual function, we compared spike-phase coupling after filtering the LFP in various filter bands. If the spiking probability of a neuron is phase locked with the LFP within some frequency band, this is evidence that the neuron in question participates, to some degree, in oscillatory activity of the larger ensemble of neurons whose transmembrane currents give rise to that rhythm. If narrowband rhythms do reflect distinct information channels, then the phase of these oscillations should be particularly informative about the excitability of neurons participating in that oscillatory rhythm and therefore the timing of their spontaneous spiking activity. Alternatively, if the excitability of the population is reflected in the phase of the broad spectrum fluctuations, then the spiking activity of neurons should be more poorly predicted by any individual oscillatory component and better predicted by the phase of the broadband LFP. Therefore, in this work we take the magnitude of spike-phase coupling as a direct measure of the degree to which oscillatory activity reflects a discrete information channel.
The ability to test between these alternatives has been limited, however, because the calculation of phase using the Hilbert transform breaks down when the frequency content of a signal is too broad (Le Van Quyen et al., 2001). It had been infeasible to directly compare the relative phase coupling of spiking activity to narrowband or broadband LFP signals without consideration of this potential confound. To overcome this technical limitation, we have developed a measure of phase [generalized phase (GP); Davis et al., 2020], a generalization of the Hilbert transform that can be applied to spectrally broad signals, allowing us to directly compare narrowband and broadband phase estimates of cortical excitability. This enabled us to test whether the timing of spontaneous spiking activity in cortical populations is better aligned with the phase of classically defined narrowband oscillations, similar to channels on a radio, or is more tightly coupled to the phase of the broad ensemble of nonstationary components. In recordings made from the marmoset middle temporal (MT) extrastriate visual cortex, we find that spontaneous spiking is more strongly phase coupled to the broadband LFP than to any individual narrow band. Thus, fluctuations in spontaneous neuronal spiking are not coupled preferentially to individual narrowband oscillations but rather track with the instantaneous fluctuations of neural activity as they change from moment to moment.
Materials and Methods
Electrophysiology recordings
One male (monkey W) and one female (monkey T) marmoset monkey (Callithrix jacchus) were surgically implanted with a headpost for head stabilization and eye tracking. The headpost contained a hollow chamber housing an Omnetics connector for a Utah array (Blackrock Microsystems), which was implanted in a 7 × 10 mm craniotomy over area MT (stereotaxic coordinates, 2 mm anterior, 12 mm dorsal). An 8 × 8 (64 channel, monkey W) and 9 × 9 with alternating channels removed (40 channel, monkey T) Utah array was chronically implanted over area MT using a pneumatic inserter wand. The electrode spacing was 400 µm with a pitch depth of 1.5 mm. The craniotomy was closed with DuraSeal (Integra Life Sciences, monkey W) or DuraGen (Integra Life Sciences, monkey T), and covered with a titanium mesh embedded in dental acrylic. All surgical procedures were performed with the monkeys under general anesthesia in an aseptic environment in compliance with National Institutes of Health guidelines. All experimental methods were approved by the Institutional Animal Care and Use Committee of the Salk Institute for Biological Studies and conformed with National Institutes of Health guidelines. Data used in this study were previously used in Davis et al. (2020).
Marmosets were trained to enter a custom-built marmoset chair that was placed inside a Faraday box with an LCD monitor (ASUS VG248QE) at a distance of 40 cm. The monitor was set to a refresh rate of 100 Hz and gamma corrected with a mean gray luminance of 75 candelas/m2. Electrode voltages were recorded from the Utah arrays using two Intan Technologies RHD2132 amplifiers connected to an Intan Technologies RHD2000 USB interface board. Data were sampled at 30 kHz from all channels. The marmosets were head fixed by a headpost for all recordings. Eye position was measured with an ISCAN CCD infrared camera sampling eye position at 500 Hz. Stimulus presentation and behavioral control was managed through MonkeyLogic (Asaad et al., 2013) in MATLAB. Digital and analog signals were coordinated through National Instruments DAQ cards (catalog #PCI6621) and BNC breakout boxes (catalog #BNC2090A, National Instruments). Neural data were broken into two streams for offline processing of spikes (single-unit and multiunit activity) and LFPs. Spike data were high-pass filtered at 500 Hz, and candidate spike waveforms were defined as exceeding 4 SDs of a sliding 1 s window of ongoing voltage fluctuations. Artifacts were rejected if appearing synchronously (within 0.5 ms) on over a quarter of all recorded channels. Segments of data (1.5 ms) around the time of candidate spikes were selected for spike sorting using principal component analysis through the open source spike sorting software MClust in MATLAB (A. David Redish, University of Minnesota). Sorted units were classified as single- or multiunits, and single units were validated by the presence of a clear refractory period in the autocorrelogram. LFP data were low-pass filtered at 300 Hz and downsampled to 1000 Hz.
Fixation behavior
The marmosets were trained to saccade to a marmoset face to initiate each trial. On the gaze arriving at the face, it disappeared and was replaced with a white fixation point [0.15 degrees of visual angle (DVAs)]. The marmosets held fixation on the fixation point (1.5 DVA tolerance) for a minimum duration (400 ms monkey W, 300 ms monkey T) awaiting the appearance of a drifting Gabor target (4 DVA diameter; appearing 6–7 DVAs eccentricity at 1 of 2 equally eccentric locations in the visual field contralateral to the recording array). Spontaneous data were analyzed from the period of fixation preceding the appearance of a target and excluding the initial 100 ms following fixation initiation. Early fixation breaks (defined by the excursion of the eye position from the fixation window) were excluded from analysis.
Free-viewing natural scenes
Marmosets were head fixed and their gaze monitored as in the previous task. Grayscale versions of naturalistic images (spanning 20–30 DVAs) were randomly interleaved and presented to the monkey. The monkey was free to look at the images, and after 10 s was given a juice reward. Visual activity was analyzed as in the spontaneous fixation data excluding a 250 ms window around the times of saccades. Saccades were defined as velocity peaks exceeding 25° per second. The time of saccade was taken from the peak velocity after threshold crossing. Velocity was calculated from the absolute value of the first numerical derivative of the smoothed vertical and horizontal eye traces (5 ms sliding Gaussian). We excluded from our analysis spikes that occurred from 50 ms before to 200 ms after detected saccades. Multiunit spiking activity from two recording sessions in Monkey T and one session in Monkey W (N = 142 units) were combined and analyzed as there was no significant difference in spike-phase index (SPI) effects between the monkeys (p = 0.10; Wilcoxon rank-sum test).
Spike artifact elimination
To eliminate spike artifacts from the LFP, we applied a despiking algorithm first described in Zanos et al. 2011. The goal of the algorithm is to eliminate the contribution of spike waveforms to the signal that after being downsampled and low-pass filtered constitutes the LFP. The algorithm assumes the LFP is based on the measured wideband voltage trace recorded from the electrode (y), which is composed of a low-frequency signal (the LFP, w), high-frequency spike components
Here, m is the number of spikes for kth neuron k. The high-frequency component of k is the convolution of the spike train
Rather than using a spike-triggered average approach to generate a mean template of the spike waveform, which is subtracted at the time of each spike, the algorithm optimally estimates the local field potential w, each spike waveform
The first assumption is that the LFP is smooth with most of its power in the lower frequencies as follows:
An implementation of this algorithm in MATLAB is available from the original authors (http://apps.mni.mcgill.ca/research/cpack/lfpcode.zip).
Generalized phase
We calculated GP as described previously (Davis et al., 2020). The purpose of GP is to mitigate the breakdown of the analytic signal representation for spectrally broad signals. As an initial step in the GP representation, then, we filter the signal within a wide bandpass (i.e., 5–50 Hz; fourth-order zero-phase Butterworth filter), excluding low-frequency content that contributes to origin offsets in the complex plane that distort the estimate of phase angles for higher frequency signals. We then use the single-sided Fourier transform approach (http://www.fuchs-braun.com/media/d9140c7b3d5004fbffff8007fffffff0.pdf; Marple, 1999) on the wideband signal and compute phase derivatives as finite differences, which are calculated by multiplications in the complex plane (Feldman, 2011; Muller et al., 2014, 2016). High-frequency intrusions appear in the analytic signal representation as complex riding cycles (Feldman, 2011), which manifest as periods of negative frequencies in the analytic signal representation. As a secondary step, we then numerically detect these complex riding cycles (Nc points of negative frequency) and use shape-preserving piecewise cubic interpolation on the next 2Nc points following the detected negative frequency epoch. The resulting representation captures the phase of the largest fluctuation on the recording electrode at any moment in time (Fig. 1f), without the distortions because of the large, low-frequency intrusions or the smaller, high-frequency intrusions characteristic of the 1/f-type fluctuations in cortical LFP (Pereda et al., 1998; Linkenkaer-Hansen et al., 2001; Milstein et al., 2009). All phase estimates of filtered LFP segments were calculated using the GP algorithm.
Spike-phase coupling
Three-second LFP epochs centered on the period of fixation were analyzed during the fixational behavioral task. The LFP segments were filtered (fourth-order zero-phase Butterworth filter with varying filter bandwidths depending on the analysis condition), and spike-phase coupling was calculated over epochs of fixation excluding the initial 100 ms following fixation initiation. The degree of spike-phase coupling was measured as the mean resultant vector length for the LFP phase distribution collected at the time of observed spikes. This measure was calculated using the circ_r function in the Circular Statistics Toolbox for MATLAB (Berens, 2009). The mean resultant vector r of the spike-phase distribution is the normalized sum over complex exponentials of the phase angles ϕ as follows:
Filtered-raw LFP signal to noise ratio
We calculated the signal-to-noise ratio (SNR) in decibels by computing the ratio of the summed squared magnitude of the filtered LFP in either theta (4–8 Hz), alpha (8–15 Hz), beta (15–30 Hz), low gamma (30–50 Hz), or the wideband (5–50 Hz) filter to the summed squared magnitude of the broadband 1–100 Hz LFP. The SNR was calculated over a window corresponding to approximately a single cycle of the mean frequency of each filter band (150, 75, 50, 25, and 50 ms respectively). The tested window was slid by 1/fifth the window width over the entire fixation period. Only spike times that occurred in a window that exceeded −5 dB SNR was included in the SPI calculation for that narrowband filter.
Generalized linear model analysis
To compare the relative predictive power of spike timing between multiple narrow and a single wideband measure of LFP phase (GP), we tested general linear models (GLMs) trained to predict the likelihood of spiking activity. In particular, both GLMs were trained using LFP phases recorded at points in time when spikes occurred, and an equal size sample of LFP phases, selected at random, when no spike occurred. The first model used as predictors the phase at the time of each spike or nonspike for theta (4–8 Hz), alpha (8–15 Hz), beta (15–30 Hz), and low gamma (30–50 Hz) narrowband filtered LFP. The second model used a single predictor, the narrowband beta phase (15–30 Hz), and the third model also used a single predictor, the wideband (4–50 Hz) LFP GP computed on the same training set. To linearize the circular phase variables, we used the sine and cosine of each phase value as separate predictors (Cremers and Klugkist, 2018), resulting in eight predictors for the narrowband model and two predictors for the single narrow and wideband models as follows:
Simulated spike and LFP generation
To generate surrogate spiking and LFP data, we first generated a normal distribution of random frequency values with a mean of 10 Hz and a standard distribution of 1 Hz. We then generated a 100 s sinusoidal signal, whose frequency drifted with random draws from the frequency distribution. In the case where spikes were generated from the phase of this narrowband signal, we first filtered this signal between 8 and 15 Hz and used the phase to generate spike times. We also generated a broadband noise signal generated from a Gaussian distribution with mean of 0 and an SD of 1, whose power spectrum followed a 1/f power law (Kasdin, 1995). In the case where spikes were generated from the phase of the broadband signal, the drifting sinusoidal and pink noise signals were summed in the frequency domain and transformed back into the temporal domain and filtered between 1 and 100 Hz. The combination of the sinusoidal signal and the noise signal made up our surrogate LFP signal, which was identical between the alternative spike generating conditions.
Spike times were generated using a phase-dependent Poisson spike generator. The phase-dependent spiking probability was defined with a circular-linear function across 21 phase bins with a 0% spiking probability at 0 rad phases and a 1% spiking probability at ±π rad phases. At each millisecond in time, a random value was drawn from a Poisson distribution, whose lambda corresponded to the probability of a spike occurring at the phase of either the sinusoidal (narrowband hypothesis) or surrogate LFP (broadband hypothesis) signal at that millisecond. Any drawn value that exceeded zero produced a single spike time. The relative phase-dependent spike probabilities produced irregular spike trains with mean firing rates roughly between 5 and 6 Hz in both conditions. The calculation of spike-phase coupling was performed identically as that in the recorded data. The surrogate LFP was filtered in either narrow or wide filters, and the GP was drawn at the time of each spike to generate spike-phase distributions.
Statistical analysis
Statistical tests used in this study include the parametric pairwise Student's t test, the nonparametric Wilcoxon signed-rank test, and the Wilcoxon rank sum test. Two monkeys were used in this work. No power analyses were performed as the number of monkeys used followed standard conventions to reduce the number of primates required for neuroscience research. All results were consistent across both monkeys and were therefore collapsed for analysis. Individual measurements within N = 20 recording session were averaged, and statistical tests were performed on the averages across recording sessions.
Data Availability
The data that support the findings of this study are available from the corresponding authors on reasonable request. An open-source code repository for the generalized phase algorithm is available from http://mullerlab.github.io.
Results
We measured spike-phase coupling for single- and multiunit spiking activity across traditional narrowband and broadband filtered LFP signals. Spiking activity and LFP data were previously recorded from chronically implanted multielectrode arrays (Utah array, Blackrock Microsystems) in area MT of two common marmosets (Callithrix jacchus; 10 recording sessions in each monkey) as they fixated on a point on an otherwise blank screen (gray background, 75 candela/m2; Fig. 1a), awaiting the appearance of a faint visual target during a challenging visual detection task (Davis et al., 2020). Similar experimental paradigms have been used to study the relationship between prestimulus oscillatory phase and sensory processing and behavioral performance (Busch et al., 2009; Balasubramanian et al., 2020; Zareian et al., 2020). The raw LFP (filtered from 1 to 100 Hz) sporadically exhibited rhythmic fluctuations across a range of time scales, but there was not a clear peak in the power spectral density that would be consistent with a clear and consistent oscillatory component (Fig. 1b).
The LFP during periods of fixation was filtered in classically defined frequency bands—theta (4–8 Hz), alpha (8–15 Hz), beta (15–30 Hz), and low gamma (30–50 Hz)—or in a wideband filter that spanned all of these narrow bands from 5 to 50 Hz (Fig. 1c–g). The bounds of the wideband filter were selected to exclude low-frequency fluctuations (<5 Hz) that are associated with slow changes in arousal (Steriade et al., 2001; Petersen et al., 2003) and high-frequency components that may be contaminated by spiking artifacts and could, therefore, induce spurious spike–LFP correlations (Ray et al., 2008a; Zanos et al., 2011). However, it could be possible spiking artifacts exist at sub-50 Hz frequencies, which could, in principle, bias our estimate of the relationship between spiking activity and LFP phase in our broadband representation. To mitigate this potential confound, we performed a despiking procedure on the data as described in Zanos et al. (2011). This removes spike waveforms from the raw (30 kHz) recorded electrode data signals through spike-waveform subtraction and interpolation before downsampling and filtering into the LFP. Any remaining relationship between the phases of sub-50 Hz activity in any frequency band must therefore be because of an indirect relationship between the population currents that give rise to the LFP and individual neuronal spiking, and not the direct contribution of that spike occurring itself.
If spiking activity is either organized into oscillations, giving rise to narrowband fluctuations in the LFP, or if LFP oscillations reflect population-wide subthreshold fluctuations that modulate the probability of spiking within a particular band, then spikes should tend to be aligned in phase with the LFP within these frequency ranges. Alternatively, if no individual rhythmic component of the LFP dictates the excitability of neurons, but rather the precise, moment-by-moment fluctuations of the LFP reflect the state of the population, we would expect spikes to occur more often at phases of the broadband LFP that correspond to states of depolarization across the local population, regardless of frequency.
To test these competing hypotheses, we measured the phase of each filtered LFP signal at the times of multiunit spiking activity. Phase is conventionally measured for oscillatory or spectrally narrow signals by calculating the analytic signal (Marple, 1999; Feldman, 2011), where instantaneous amplitude and phase can be expressed in polar coordinates and whose real and imaginary parts are related to each other by the Hilbert transform. However, for spectrally broad signals, the standard computational implementations break down (Le Van Quyen et al., 2001). Low frequencies can shift the analytic signal representation by a constant in the complex plane, distorting the estimated phase angle. In addition, high-frequency intrusions introduce complex riding cycles that generate phase reversals and appear as negative frequencies that distort the analytic signal. To address these problems, we introduced an updated approach to the analytic signal representation, termed “generalized phase” (Davis et al., 2020). Briefly, in this approach we first impose a high-pass cutoff on the signal (5 Hz). This step aims to eliminate low-frequency intrusions while also preserving a significant portion of the signal spectrum and minimizing waveform distortion. Second, we identify negative frequencies, which can arise from high-frequency intrusions, and remove them, replacing the phase values with shape-preserving interpolation. This approximates the continuation of the trajectory of the dominant fluctuation. The result, after filtering, is an estimate of phase that tracks with the dominant frequency component of the LFP as it shifts over time (Fig. 2a) while minimizing phase distortions that arise because of narrowband filtering a nonstationary broad spectrum signal such as the raw LFP (Yael et al., 2018). All results reported here, for both broadband and narrowband filtered data, were computed using GP. Low- and high-frequency intrusions are rare in narrowband filtered signals, so for narrowband filtered data, computation of GP should yield very similar phase estimates to those estimated using the Hilbert transform. To confirm this, all analyses were repeated for narrowband filtered signals using the Hilbert transform. As expected, the results were virtually identical. All future mentions of phase therefore refer to the GP of the signal.
The phase of the wideband filtered signal is strongly coupled to the timing of measured multiunit spiking activity (Fig. 2b). We measured an index of the coupling of spikes to each filtered LFP by calculating the mean resultant length of the circular spike-phase distribution. This SPI value ranges from 0 (uniform spike-phase distribution) to 1 (spikes perfectly coupled to a single phase), and for the 5–50 Hz wideband filtered signal, the average SPI was 0.15 ± 0.009 SEM (N = 20 sessions across two monkeys). The wideband filtered SPI was significantly stronger than the coupling observed after filtering in theta (SPI = 0.08 ± 0.005; p < 1 × 10−9; two-tailed paired sample t test), alpha (SPI = 0.07 ± 0.005; p < 1 × 10−10), beta (SPI = 0.11 ± 0.007; p < 1 × 10−11), or gamma (SPI = 0.08 ± 0.009; p < 1 × 10−9) frequency bands (Fig. 2c–f). These results suggest that the instantaneous rhythmic state of neuronal excitability is better reflected in the phase of the ensemble LFP activity rather than in the phase of any particular narrowband subcomponent.
If oscillations reflect information streams analogous to channels on a radio, then it could be the case that some neurons are more coupled to one embedded oscillation, and other neurons are more coupled to a different oscillation, and by collapsing across multiunit activity, the phase-dependence of the spiking activity is diluted for each narrowband filter. If true, we might find stronger spike-phase coupling for the wideband filter across the populations, although individual neurons are best coupled to different narrowband oscillations. To test this, we measured the spike-phase coupling across filters for well-isolated single units in our recordings. We did not find any evidence of differential preference across neurons for narrowband signals. Rather, the majority of neurons had a stronger SPI to the state of the wideband signal compared with theta (78.50%, N = 107 single units; Fig. 3a), alpha (78.50%, Fig. 3b), beta (83.18%, Fig. 3c), or gamma (85.98%, Fig. 3c) filtered signals. Additionally, for the minority of neurons that did show a stronger SPI to a narrowband filtered signal, they were more weakly coupled in general (average SPI = 0.09) and did not show specific preference to any one narrowband frequency. Thus, variable population phase coupling to narrowband signals could not explain why the wideband filtered signals exhibit stronger spike-phase coupling.
One possibility is that spike timing is only governed by each narrowband oscillation when that oscillation is strongly present in the data, and as each oscillation is only transiently present, it is unfair to expect, for example, gamma to predict spike timing when gamma is not present in the data. To test this we restricted our analysis to only count spikes for each frequency band at times when there is strong oscillatory power in that band. To do this we calculated the SNR between the filtered and raw (1–100 Hz) LFP and identified epochs where the narrowband signal exceeded a −5 dB threshold for at least one cycle of the center frequency of the filter bandwidth. Only spikes that occurred during these epochs were included for that narrowband SPI measure. Despite restricting each filter band to spikes that occur when those oscillations are transiently apparent in the data, the wideband measure still captures the strongest SPI values (Fig. 4a; wideband mean SPI = 0.16 ± 0.010 SEM; compared with theta, 0.11 ± 0.010; alpha, 0.08 ± 0.005; beta, 0.12 ± 0.010; and gamma, 0.09 ± 0.005; p < 0.001, Wilcoxon signed-rank test) while also describing a majority of the recorded data (approximate fraction above dB threshold; wideband 92% vs theta 16%; alpha 46%; beta 42%; and low gamma 19%).
Thus, spike timing is better predicted by broadband phase than narrowband phase for any of the bands tested. We next asked how well spike timing could be predicted based on the combination of all four narrowband filtered signals. To test this we constructed a GLM that took as its input the phase values measured over the four narrow band frequencies (spanning 4–50 Hz) at times when a spike occurred and an equal number of randomly drawn times when no spike occurred. The GLM was trained to predict whether a spike occurred, based on the four phases. The model was trained on half the data in each recording session, with the remaining data held out as a test set. The ability of the model to predict spiking was measured using ROC analysis.
We reasoned that if oscillatory activity across the multiple narrow bands drives spiking activity, the four-factor GLM, which has simultaneous access to the phases of all four oscillatory signals, should predict spiking better than a GLM trained to predict spiking based on the phase computed in an individual band (four-factor GLM AUC = 0.578 ± 0.004 SEM; single narrowband GLM AUC = 0.545 ± 0.003 SEM; p = 0.00,009, Wilcoxon signed-rank test). This analysis shows that more information about spiking is present across multiple bands. This is consistent with two different hypotheses. The first is that the narrow bands capture the individual contribution of oscillations that fall within each band, and the four-factor GLM reflects the joint contributions of these oscillatory drivers. An alternative hypothesis is that the processes that drive spiking activity fluctuate over time in their power spectrum, and spiking activity follows these fluctuations over time, regardless of where they travel in frequency. If the first hypothesis is true, the four-factor GLM, which has access to the phase within each band, should perform better than a single-factor broadband GLM, which is provided with a single measure of phase that is blind to the interactions across the same frequency space. If the second hypothesis is true, the single-factor broadband GLM, which uses a measure of phase that tracks with the dominant LFP frequency as it changes over time, should do as well as the four-factor GLM.
To test this, a GLM was trained on the same data that was used to train the four-factor GLM, but instead of providing it with four phases computed within the four narrow bands, it was trained using only a unitary measure of phase—GP applied to the wideband (4–50 Hz) signal as its input—and its ability to predict spiking was measured using the same ROC analysis. As shown in Figure 4b, there was no significant difference in the ability of the combined four narrowband or one wideband GLM to predict spike times as defined by the area under the curve for each the ROC for each session (wideband mean AUC = 0.579 ± 0.005 SEM; p = 0.16, Wilcoxon signed-rank test). Thus, even when combining signals across multiple frequency bands, narrowband filtering adds no information beyond what is already present in the phase of the momentarily dominant fluctuation in the LFP preserved in the wideband representation and as measured using generalized phase.
Our results suggest that spontaneous neuronal spiking in the neocortex is not organized by oscillatory activity but rather is modulated by fluctuations in synaptic activity that can be estimated from the instantaneous phase of the broadband LFP. If true, then SPI values should be correlated with how much filtering alters the LFP phase relative to the raw recorded LFP. To test this, we compared the strength of spike-phase coupling to each bandpass filtered signal with the degree of correlation between the LFP signal before and after filtering (Fig. 4c). There was a significant positive correlation between SPI and the raw-filtered LFP correlation across recording sessions (Pearson's r = 0.65 ± 0.11 95% CI, p < 1 × 10−12), suggesting a direct relationship between spike-phase coupling and how well the filtered LFP tracked with the raw LFP.
If spikes are more coupled to the broadband LFP than any embedded narrowband oscillation, then the optimal filter band for maximizing SPI should be one that is as broad as possible. To test for an optimal filter band, we scanned across a large parameter space varying the lower and upper bounds of the bandpass filter. The lower bound ranged from 1 to 50 Hz, and the upper bound ranged from 5 to 125 Hz with a minimum bandwidth of 4 Hz. Consistent with our prior results, the strongest spike-phase coupling was observed for filters that included the largest width of the signal spectrum, with an exception for the lowest frequencies (Fig. 4d). These results indicate that optimal filters for maximizing spike-phase coupling estimates span from 3 Hz in the lower band and as high as we sampled in the upper band (125 Hz), assuming spike artifacts are effectively removed from high-frequency components in the LFP. If not, a cautious step then is maintaining a low-pass filter, which serves to help mitigate spurious coupling values because of residual spike artifacts in higher frequencies.
Spiking activity can bleed into the LFP, artifactually inflating estimates of spike-phase coupling in high-frequency bands. Spike artifacts may be responsible for some gamma phase relationships with spiking activity, as the contribution of spike artifacts in the LFP had been previously observed down to 50 Hz (Ray and Maunsell, 2011; Zanos et al., 2011). To avoid this, we performed a despiking procedure and examined the consequence of that despiking on SPI estimates. A comparison of SPI values on the same data with and without despiking found that the despiking procedure significantly reduced SPI values for frequency bands that included frequencies below 50 Hz (but not below 15 Hz) such as low gamma (30–50 Hz; not despiked SPI = 0.13, p = 1.65 × 10−7, two-tailed Wilcoxon rank sum test), beta (15–30 Hz; SPI = 0.13, p = 0.026), and the wideband (5–50 Hz; SPI = 0.18, p = 0.009). There was no significant reduction in either alpha (8–15 Hz; SPI = 0.08, p = 0.067) or theta (4–8 Hz; SPI = 0.09, p = 0.190) when we despiked the LFP. These observations are consistent with recent reports of spike artifacts having an impact on spike-LFP synchronization at frequencies as low as 20 Hz (Banaie Boroujeni et al., 2020). These results argue either that the artificial coupling of spiking activity to LFP phase may be present at low frequencies or that despiking techniques are overly liberal in the removal of spike waveforms. Regardless, even if we consider the possibility that the despiking procedure is introducing more noise than it is eliminating, the main result—that the broadband LFP phase produces the strongest SPI values—holds when this technique is not applied, and the raw data are left intact.
The results described so far are limited to spontaneous activity recorded during a period in which animals foveated a fixation point at the center of a blank screen while awaiting the appearance of a faint visual target. Do these findings generalize to more naturalistic viewing conditions? To test this, we calculated the SPI for each frequency band in animals as they freely viewed natural scene images. As the focus here is on intrinsic fluctuations, not the transient responses that are evoked at the time of the saccade, neural activity at the time of the saccade (from 50 ms before and ending 200 ms after saccades) was eliminated from analysis. Consistent with the pattern observed during fixation of a blank screen, the wideband filtered signal produced the strongest SPI values (0.16 ± 0.008; N = 142 multiunits across two sessions in Monkey T and one session in Monkey W), which was significantly stronger than the SPI values measured for theta (0.14 ± 0.007; p < 1 × 10−7, Wilcoxon signed-rank test), alpha (0.13 ± 0.007; p < 1 × 10−16), beta (0.10 ± 0.006; p < 1 × 10−14), and low gamma (0.11 ± 0.006; p < 1 × 10−12). Thus, the spontaneous coupling of spiking activity to broadband fluctuations is not limited to fixating on a blank screen but is apparent during more dynamic active vision.
Although our experimental results suggest spiking activity is better correlated with the instantaneous state of the broadband LFP rather than any individual oscillatory component, the ground truth mechanism relating spiking to rhythmic LFP activity is unknown in our recordings. To explore whether our observations can be explained by the hypothesis that spiking activity is coupled to broadband LFP phase as opposed to a narrowband oscillation, we simulated an LFP signal by combining a narrowband oscillatory fluctuation that consisted of spectral power drifting between 8 and 15 Hz with broad spectrum noise. The power spectral density of this simulated LFP fluctuation was designed to be consistent with the typical 1/f power law observed in cortical recordings in vivo (Miller et al., 2009; Fig. 5d). We then generated spike times from a Poisson spike generator, where the probability was dependent on either the phase of the narrowband 8–15 Hz oscillatory signal (Fig. 5a, hypothesis A) or the phase of the combined narrowband and broad spectrum signals (Fig. 5b, hypothesis B). Spike probability was phase dependent, with spikes most likely to occur near ± π radians and spikes least likely to occur near 0 radians. Importantly, the spectral content of the simulated LFP was identical between the two alternative hypotheses and only the timing of spikes differed between the two conditions (Fig. 5c).
To recover the signal correlated with spike generation, we filtered the simulated LFP in either a narrow bandpass filter from 8 to 15 Hz, or a wide bandpass filter from 5 to 100 Hz. In the case where spikes were correlated with the phase of the narrowband oscillatory component (hypothesis A, blue), the narrowband filtered LFP signal was relatively weakly correlated with the raw simulated LFP (Pearson's r = 0.49; Fig. 6a). However, spike timing was strongly coupled to the phase of the narrowband filtered LFP signal (SPI = 0.34 ± 0.002 SEM; N = 20 simulations). This coupling was significantly stronger than when using the wideband filter to recover spike phases (SPI = 0.16 ± 0.002 SEM; p = < 0.0001, two-tailed Wilcoxon signed-rank test; Fig. 6b). In the case where the spikes were generated from the phase of the broad spectral content of the simulated LFP (hypothesis B, red), the wideband filtered LFP was strongly correlated with the raw simulated LFP (Pearson's r = 0.84, Fig. 6c), and the spike-phase relationship was significantly stronger after filtering in the wideband (SPI = 0.29 ± 0.002 SEM) compared with the spike-phase coupling to the narrowband filtered LFP (SPI = 0.14 ± 0.002 SEM; p = < 0.0001 two-tailed Wilcoxon signed-rank test; Fig. 6d). These results indicate, in principle, that if neurons were coupled to an oscillatory component, then narrowband filtering to extract that oscillation would indeed yield stronger spike-phase coupling than the broadband signal.
We next asked whether a narrowband or broadband spike-correlated signal could reproduce the observed relationship of increasing spike-phase coupling with increasing correlation between the filtered and raw LFP signal. We filtered the signal under various filters (theta, 4–8 Hz; alpha, 8–15 Hz; beta, 15–30 Hz; and wideband, 5–50 Hz) as in the cortical recordings, as well as a broad bandpass from 1 to 100 Hz, and measured both the spike-phase coupling for narrowband and broadband correlated spike generation and the correlation between the filtered and raw LFP signal. In the case where spikes were correlated with the narrowband signal, the best filter was the 8–15 Hz filter (matching the source of the spike-generating signal), followed by the wideband and broadband filters, which each included the spike-generating signal band within its bandwidth but also included a smaller and larger part of the noise spectrum, respectively (Fig. 7a). In the case where spikes were correlated with the broadband signal, the best filter was the broadband filter and decreased as the filters became narrower (Fig. 7b). The narrowband spike source had a weak correlation between the SPI (Pearson's r = 0.27), and the degree of filter raw-signal similarity as the optimal filter was one that eliminated the broadband noise from the simulated LFP. In contrast, the broadband spike source reproduced the strong positive correlation between SPI and the filtered raw LFP similarity observed in our recordings (Pearson's r = 0.92; Fig. 7c). Our results indicate a that model where spikes are coupled to the state of fluctuations in the broad spectral content of the LFP is sufficient to account for our observations in vivo and suggest neuronal spiking is not preferentially coupled to narrowband oscillations.
Discussion
A central goal of systems neuroscience is to understand how brain activity underlies information processing and behavior. Ideally, we would like to record every action potential of every neuron and ask how they relate to one another in the service of behavior, but even with the best available neurophysiological tools—sets of electrode arrays with contacts numbering in the thousands—we can only sample a tiny fraction of the neurons in the brain. Therefore, neurophysiologists typically rely on indirect measures of the activity to estimate the spiking statistics of larger cortical populations. These include LFP, EEG, or MEG, which provide indirect measures of the activity of larger populations of neurons. Rhythmic patterns of activity are often observed in these measures, and it is common to treat these rhythmic patterns as meaningful computational units, potentially serving as independent channels of information processing, or if not independent in the context of cross-frequency phase-amplitude coupling (Munia and Aviyente, 2019), at least functionally dissociable from the signal in which they are embedded (Thut et al., 2012; Einevoll et al., 2013), similar to turning the dial on a radio to receive different streams of information.
One way of thinking about rhythmic dynamics is that the spiking probabilities of the neurons in the larger population covary within some frequency band and that this results in an oscillation, for the example studied here, in the LFP. If so, then by filtering the LFP within that oscillatory band and asking how it relates to some measure of either behavior (e.g., performance on a discrimination task), a neural property (such as spike timing or transmission of information across areas), or its covariation with some behavioral manipulation (e.g., directing attention into or away from the retinotopic locus of the electrode), one can identify the contribution of the oscillation to neural computations or behavior. However, there are some problems with treating neural fluctuations as oscillations. First, neural fluctuations are often only transiently rhythmic in the awake state (Jones, 2016), and even then they are not purely sinusoidal (Cole and Voytek, 2017) as they drift in frequency content from moment to moment with changes in arousal (Vinck et al., 2015), attention (Fries et al., 2001), or sensory input (Henrie and Shapley, 2005). Even in the case when neural fluctuations are strongly rhythmic, we find narrowband filtering captures less of the spike-phase relationship than when maintaining a wideband representation. This may be because the application of narrowband filters to signals that are nonstationary in their frequency content can result in a loss of timing precision in phase estimates (Yael et al., 2018).
The results presented here argue that neurons are not specifically coupled to narrowband oscillatory activity, but rather it is the state of the broadband moment-to-moment fluctuations that are informative of the relative excitability of the local population. This is not to say that rhythms are not apparent in fluctuating dynamics or that they are irrelevant for cortical function. Nor are we suggesting that rhythmic power is limited to what one would expect from stochastic synchronizations in a 1/f noise process. For example, it is not the case that oscillatory rhythms are only as informative as their fraction of the spectral content of broadband fluctuations. We observed that low gamma filtered signals had stronger SPI values than one might expect based on their relative power in the PSD and given how poorly correlated the gamma filtered signals were to the raw LFP. Similarly, the alpha band filtered signals had much more power and were relatively well correlated with the raw LFP yet had weaker SPI values than the beta band filtered signals, which were more poorly correlated with the raw LFP. Indeed, there is variation in the degree to which spikes couple to LFP phase across the five frequency bands studied here. However, that does not imply that those frequency bands are independent information channels, distinct from the rest of the LFP. It is evident that they are not, as we see the strongest SPI values for the broadest frequency bands.
To test what one would expect to see if it were the case spikes preferentially coupled to a narrow set of frequencies, we simulated spike trains generated from the phase of oscillatory signals embedded in an otherwise 1/f noise spectrum (hypothesis A). Under these conditions, we found a stronger SPI to the narrowband filter that best matched the signal underlying the spike generation signal. We also found a reduction in SPI values when the broadband filter was used. This matches what one would intuitively expect from a system composed of an oscillatory signal combined additively with a broad noise. This is the intuition that often underlies narrowband filtering approaches in electrophysiological signal analysis. Although there may be alternative explanations for why a broadband signal produces stronger SPI values in our cortical recordings, the second model, where spikes are fluctuation driven (hypothesis B), was sufficient to account for the spike-LFP coupling relationships observed in the data.
Because the phase of narrowband oscillatory activity does not predict spiking activity as well as the phase of wideband activity, it raises the question on whether and when narrowband filtering is appropriate to study rhythmic spiking dynamics. The use of narrowband filters assumes a frequency resolved signal in the brain that is embedded in noise. As shown by hypothesis A and in Figure 6, when neural activity is strongly coupled to latent oscillatory activity, narrowband filtering is effective at recovering the signal. Therefore, in situations with steady, ongoing oscillatory activity that has low variance in frequency, such as sleep spindles, hippocampal theta, or gamma oscillations because of strong feedforward input, narrowband filtering may better capture spiking. However, if the signal is not known, narrowband filtering imposes an assumption of what is signal and noise that may not be warranted and may yield misleading results. Analytic techniques that allow for the contribution of broader frequency ranges, as used here, may reveal the degree to which results are frequency dependent or filter dependent.
It is important to note the limitations of the present findings. First, all analyses here have focused on spontaneous activity. We cannot generalize the present results to neural data collected under other conditions such as data collected during stimulus-evoked responses. Some narrowband frequency ranges, such as the gamma band, do not exhibit much power in the absence of strong sensory input (Henrie and Shapley, 2005; Ray and Maunsell, 2010). Additional experiments will be needed to determine the degree to which gamma band and generalized phase predict spike timing under these conditions. Further, the majority of the data analyzed here was recorded from the visual cortex in monkeys performing a particular task in which they foveated a fixation spot at the center of a blank screen, awaiting the appearance of a faint visual target. In our spontaneous cortical recordings, which are largely representative of the aperiodic 1/f power law observed in primate visual cortex (Fries et al., 2001; Henrie and Shapley, 2005; Yu and Ferster, 2010), even when oscillations are transiently present, narrowband filtering produces a weaker estimate of the spike–LFP relationship than a wider representation.
The generalized phase approach used here provides a meaningful measure of phase for spectrally broad signals (Davis et al., 2020) and reveals a stronger relationship between broadband LFP fluctuations and spiking probability than could be estimated from any individual narrowband filtered signal. The advantage of GP over narrowband signals is that it follows the moment-to-moment fluctuations in the signal and provides a phase value that generalizes across changes in frequency content. This approach can reveal patterns that would not be clear from an analysis of narrowband oscillations. For example, analysis of broadband measures of phase led to the discovery that the alignment of spontaneous traveling waves of cortical activity with the retinotopic locations of faint visual targets was predictive of the magnitude of evoked activity and perceptual sensitivity (Davis et al., 2020). These effects were only apparent in the data when the state of broadband LFP fluctuations was considered. When filtered in narrow bands, the predictive power of wave phase on behavioral performance was abolished. Consistent with those findings, the results presented here show that at least in the spontaneous waking activity of area MT, the instantaneous state of cortical populations is better estimated from the GP of broadband LFP fluctuations than from any narrowband oscillatory component. These results suggest that the phase of broadband neural fluctuations, rather than any specific narrowband frequency content, is the main influence on spontaneous spiking activity in the cortex.
Footnotes
This work was supported by the Gatsby Charitable Foundation; the Fiona and Sanjay Jha Chair in Neuroscience; the Canadian Institute for Health Research; the Swartz Foundation; National Institutes of Health Grants R01-EY028723, T32 EY020503-06, T32 MH020002-16A, and P30 EY019005; NSF (NeuroNex Grant No. 2015276); and BrainsCAN at Western University through the Canada First Research Excellence Fund. We thank Katie Williams, Sean Adams, Mat LeBlanc, and Tom Franken for contributions to this project.
The authors declare no competing financial interests.
- Correspondence should be addressed to John H. Reynolds at reynolds{at}salk.edu or Zachary W. Davis at zdavis{at}salk.edu