## Abstract

Neuroscience research spans multiple spatiotemporal scales, from subsecond dynamics of individual neurons to the slow coordination of billions of neurons during resting state and sleep. Here it is shown that a single functional principle—temporal fluctuations in oscillation peak frequency (“frequency sliding”)—can be used as a common analysis approach to bridge multiple scales within neuroscience. Frequency sliding is demonstrated in simulated neural networks and in human EEG data during a visual task. Simulations of biophysically detailed neuron models show that frequency sliding modulates spike threshold and timing variability, as well as coincidence detection. Finally, human resting-state EEG data demonstrate that frequency sliding occurs endogenously and can be used to identify large-scale networks. Frequency sliding appears to be a general principle that regulates brain function on multiple spatial and temporal scales, from modulating spike timing in individual neurons to coordinating large-scale brain networks during cognition and resting state.

## Introduction

The brain is complex and operates over myriad spatial and temporal scales. Multiscale interactions are a defining characteristic of the brain, and are thought to underlie aspects of cognition and consciousness (Varela et al., 2001; Breakspear and Stam, 2005; Le Van Quyen, 2011). Integration of findings across multiple scales is limited by the qualitatively distinct data analyses performed at each scale (e.g., spike rate vs event-related EEG potential vs BOLD fMRI general linear model). Here a novel neural principle is considered (termed “frequency sliding,” explained below) that can bridge findings across multiple spatial and temporal scales of brain function. This principle is based on a fundamental functional feature of neurons, and yet can be applied to neural ensembles and large-scale network connectivity.

The principle is based on the phenomenon that the firing rate of a neuron is proportional to the strength of its input (the so-called “*f–I* curve”; Schwindt et al., 1997; Hong et al., 2008). In the living brain, neurons and neural networks are bombarded with synaptic input that varies over time, and thus their firing rates also vary over time.

Curiously, this phenomenon is not taken into consideration in the majority of EEG and local field potential (LFP) data analyses, which instead assume that the frequency of a neural oscillator is approximately stationary over time, and that power within specified frequency bands is the primary outcome measure (Cohen, 2014). Modifying standard analysis methods allows a time series of instantaneous frequencies to be computed via the temporal derivative of the phase angle time series (Boashash, 1992). The method developed here is simple, hypothesis driven, and robust to noise in the phase angle time series (Fig. 1*A*).

Using this method, it is shown that both artificial and real neural networks exhibit exogenously and endogenously driven temporal fluctuations in peak frequency over time, hereafter termed frequency sliding. These subtle changes in frequency have implications for the spike timing dynamics of individual neurons, such that relatively slower oscillation frequencies decrease spike threshold and increase spike timing variability. These oscillation frequency-induced changes in spike threshold have implications for detecting weak synaptic inputs. Furthermore, co-fluctuations in frequency sliding are shown to be a sensitive and physiologically interpretable measure of functional connectivity in artificial networks, in human task-related EEG data, and in human EEG resting-state data. Together, these findings show that frequency sliding is a useful, novel, and neurophysiologically interpretable analysis tool that can help integrate findings across multiple spatial and temporal scales of neuroscience.

## Materials and Methods

##### Identifying time-varying instantaneous oscillation frequencies (frequency sliding).

The instantaneous frequency of a dynamical oscillating system can be defined as the change in the phase per unit time (Boashash, 1992). For EEG time series, this can be understood as the first temporal derivative of the phase angle time series. To convert to units of hertz, the derivative can be multiplied by the data sampling rate in hertz and then divided by 2π, in other words: Hz* _{t}* =

*s*(φ

*− φ*

_{t}

_{t}_{-1})/2π, where

*s*indicates the data sampling rate and φ

*indicates the phase angle at time*

_{t}*t*. Note that phase angles must first be unwrapped. In MATLAB code, this can be implemented as srate*diff(unwrap(phaseangles))/(2*π), where phaseangles is the phase angle time series vector with sampling rate srate. The procedure is illustrated in Figure 1

*A*. Simulated results shown in Figure 1

*B*illustrate that frequency sliding is not biased toward the center of the frequency band. Indeed, when frequency sliding is large relative to the background activity, accurate frequency sliding can be obtained even outside the filter bands (filter transition widths were set to 15%).

In the literature, there are two methods that have been proposed for estimating and analyzing changes in instantaneous frequencies in EEG. One is called empirical model decomposition (Sweeney-Reed and Nasuto, 2007), which involves repeatedly iterating through the data to identify intrinsic modes. Empirical mode decomposition has been successfully applied to EEG data for, among other applications, brain–computer interfaces (Panoulas et al., 2008) and to identify possible frequency-varying contributions to the event-related potential (Burgess, 2012). Though a useful approach, three of its drawbacks are the following: (1) the frequencies returned by empirical mode decomposition are influenced by noise and by temporal filtering that is applied to the data before the analysis; (2) it is less amenable to hypothesis-driven analyses because each frequency recovered is around one-half of the frequency of the previous intrinsic mode (thus, if one mode is around 20 Hz, the next mode is likely to be ∼10 Hz, thus impeding the ability to examine 15 Hz activity); and (3) the iterative procedure can be time-consuming, thus limiting its usefulness for large datasets involving hundreds of trials per condition, electrode, and subject. Another approach is called frequency-flows analysis (Amor et al., 2005; Rudrauf et al., 2006), which has been used to study, for example, hypersynchronization during epileptic seizures. This method involves estimating instantaneous frequencies from the phase angle time series (similar to empirical mode decomposition and the method described below) based on the result of wavelet convolution. Next, the number of electrodes exhibiting the same time course of instantaneous frequency is computed, providing information regarding large-scale multivariate synchronization over all electrodes. Because the phase angle time series is computed from wavelet convolution, and because wavelets have a Gaussian shape in the frequency domain, many frequency-overlapping wavelets are useful to avoid a potential bias of the estimation of frequency toward the peak frequency of the wavelet. Frequency-flows analysis is well suited for large-scale multivariate network analysis, but is less suitable for a priori spatially restricted and frequency band-limited analyses.

The method developed here follows from these previous approaches, but is more suitable for hypothesis-driven task- or resting-state investigations, is robust to low-pass filtering and noise (including phase slips and other phase impurities, as described below), is computationally efficient and thus quite fast, and can be used to link findings across multiple spatial scales in neuroscience such as human EEG and computational network models. The procedure is to filter the data in an a priori specified band (transition zones of 15% of the lower and upper frequency bounds were used here), apply the Hilbert transform to obtain the analytic representation, extract the phase angle time series, and then convert the phase angle time series to frequencies in hertz as described above. Due to bandpass filtering, the frequency sliding time series is bound by the pass-bands, which prevents the results from being dominated by the frequency with the largest power (for real EEG/LFP data, this is typically lower frequencies) and thus facilitates hypothesis-driven analyses. A plateau-shaped filter should be used to prevent a frequency bias toward the peak frequency of the filter (this could be the case, for example, with Morlet wavelet convolution, because Morlet wavelets have a Gaussian frequency shape). To measure functional connectivity, the frequency sliding time series from two electrodes can be correlated time point by time point. Spearman correlations are preferred because the distribution of frequency sliding is not necessarily Gaussian.

Frequency sliding can also be used as a method of measuring stationarity, which is important for many time series analyses including Granger-based connectivity and autoregressive modeling. Furthermore, phase lag-based connectivity analyses (PLI or dwPLI; Stam et al., 2007; Vinck et al., 2011) are valid only if the frequencies of the two electrodes are approximately matched (otherwise, time-varying phase lags will bias connectivity estimates). Frequency stationarity is achieved for the duration that frequency sliding maintains a relatively stable trajectory.

##### Median filter to attenuate nonphysiological noise spikes.

One signal processing issue that arises when computing frequency sliding is that a small amount of noise or a brief slip in the phase angle time series can result in a large “blip” in the frequency sliding time series due to jumps in the unwrapped phase angle time series, and, in some cases, a briefly negative derivative. An example of the effect of such noise is shown in Figure 1*A*, step 6. This can produce nonphysiological frequency jumps, e.g., from 40 to 500 Hz for only a few milliseconds, or brief negative frequency jumps (e.g., to −500 Hz) in the case of a transient negative derivative. The problem of phase noise becomes worse with frequencies that have relatively low power or when there is a power suppression, leading to decreased signal-to-noise ratio when estimating the phase angle time series.

A good solution to this problem is the median filter. The median filter is a nonparametric signal-processing approach for dealing with noise spikes in images or time series. The idea is to reassign each time point to be the median of a distribution made from surrounding time points. Because the median has low sensitivity to outliers compared with the mean, a median filter will strongly attenuate noise spikes, and will thus outperform a mean-based or convolution-based smoothing filter when the noise spikes are extreme. It is beneficial to apply the median filter several times using different window widths, and then to take the median of the median filters. For EEG time series data, an order of 10 and a maximum window size of 400 ms are recommended.

MATLAB code to implement frequency sliding and the median filter, including more comments on its functioning and parameters, is available online (mikexcohen.com/book).

##### Computational models.

Four computational models were used to examine simulated network dynamics at multiple spatial scales and with varying degrees of biological realism, and using two different network-level mechanisms of producing oscillations.

The first model was based on simplified single-compartment neurons that are defined by the differential equations provided by Izhikevich (2003). Parameters were set to standard values for regular spiking (excitatory) and fast-spiking (inhibitory) neurons (Izhikevich, 2003). There were 800 excitatory neurons and 200 inhibitory neurons, and the model is referred to as the “*E–I*” (excitatory–inhibitory) model. Population activity was estimated by the average membrane voltage of all spiking excitatory cells. Results are averaged over 100 simulations. The model was programmed and run in MATLAB (simulation code is available on the modelDB website, Accession number 154770).

The second model was also based on simplified single-compartment neurons but using the adaptive exponential equations as detailed by Brette and Gerstner, (2005). Standard model parameters were used to simulate regular spiking, fast-spiking, and bursting neurons. The architecture of this model was more complex and more anatomically realistic compared with the *E–I* model. The model contained two “columns,” each comprising three layers that simulated cortical layers 2/3, 4, and 5/6. Layer 2/3 contained 500 regular spiking neurons and 125 fast-spiking neurons, and layers 4 and 5/6 each contained 400 regular spiking neurons, 100 bursting neurons, and 125 fast-spiking neurons. All regular spiking and bursting neurons were excitatory, and all fast-spiking neurons were inhibitory. Connectivity patterns among layers followed published patterns of canonical columnar cortical connectivity (Thomson et al., 2003; Bannister, 2005), and are illustrated in Figure 3*D*. Intracolumnar connectivity probability was set to 2.5% for *E*→*E* and *E*→*I*, and 10% for *I*→*E* and *I*→*I* connections. Fifteen percent of neurons within each layer were considered an ensemble and had a connection probability of 30%. All neurons received background noise generated by 200 Poisson-spiking cells ranging in time-varying firing rates from 0 to 50 Hz. External (“thalamic”) input was simulated by stimulating layer 4 cells in the first column. There was connectivity from bursting and from regular spiking neurons in layer 5/6 of the first column to fast-spiking neurons of layer 2/3 in the second column, thus simulating one expression of feedback connectivity (Thomson et al., 2003).

The LFP was estimated as the sum of all excitatory and inhibitory postsynaptic potentials across all excitatory cells in each layer. The model was simulated using the Brian toolbox in python 2.7 (Goodman and Brette, 2008). Data from 100 simulations were exported to and analyzed in MATLAB using custom-written scripts. Because this network simulated laminar organization and interactions, it is termed the “Layer” model. The purpose of including this model in addition to the *E–I* model was to demonstrate that the relationship between input strength and peak oscillation frequency, and the use of frequency sliding to measure connectivity, is robust to the type of network simulation and to the mechanism of generating oscillations. Python code for this model is available on the modelDB website.

Initial testing of the *E–I* and Layer networks revealed that the main findings reported here (relationship between stimulation intensity and peak oscillation frequency, and measuring connectivity via correlated frequency sliding) were robust to parameter settings.

The third and fourth models were biophysically realistic conductance-based models of single neurons. Unlike the *E–I* and Layer models, which demonstrate that frequency sliding can be used as a mechanism of encoding incoming information, the purpose of the biophysically detailed model neurons was to test whether frequency sliding would have consequences for a neuron receiving synaptic inputs from such networks as the *E–I* or Layer models. The two model neurons were taken from a neocortical layer 5 bursting pyramidal cell (Hay et al., 2011) and a neocortical fast-spiking interneuron (Golomb et al., 2007). The morphology and biophysical parameters were kept from the original code downloaded from modelDB (accession numbers 139653 and 97747). Simulations were performed in NEURON 7.3 (Hines and Carnevale, 1997). The original code was not modified except for the dendritic input, which was a sinusoidal current of varying frequencies.

The model neurons were tested on two measures of functioning: spike timing and coincidence detection. For spike timing, the onset of the first action potential, and the variability of action potentials, on each sine wave peak was measured. “Coincidence detection” is a neural computation that can be conceptualized as an action potential in response to two but not one synaptic inputs. Synaptic inputs to the layer 5 pyramidal cell were modeled as an EPSP-like current injected into two apical dendrites, one 5 ms after the first. The EPSP kinetics were the same as those used by Hay et al. (2011). In separate simulation runs, these two inputs were provided at different phases of the simultaneous sine wave input. The amplitude of the input sine wave was decreased to prevent action potentials resulting from the sinusoidal stimulation, and to provide variability in whether the synaptic inputs could elicit an action potential (that is, action potentials resulted from the combination of two synaptic inputs and a near-peak phase of the input sine wave). To quantify the changes in sensitivity to coincidence detection as a function of the sinusoidal input, the average sinusoid input phase angle from the first action potential, or the number of action potentials elicited, following the two synaptic inputs was computed. Statistical evaluations were performed via linear fits in MATLAB (function regstats).

##### EEG task and resting state.

The human EEG task was designed to mimic input strength in a manner similar to the direct-current and alternating-current inputs in the *E–I* and Layer models. In the EEG task, nine subjects (six male) recruited from the University of Amsterdam community viewed a large gray square on a black background. The square was shown in the left or in the right lower visual hemifield, remained on screen for 4 s, and was followed by a 2 s intertrial interval (only the first 2 s of the stimulus was used for Fig. 6*B*, because the frequency structure became less stable toward the end of the trial). There were four levels of luminance ranging from 31 to 100% of monitor maximum brightness in equal steps (pixel values 80, 138, 197, and 255). In a fifth condition, the luminance oscillated at 1 Hz, with a random phase on each trial. The purpose of having a random phase on each trial was to prevent any potential confound between frequency sliding induced by the oscillating stimulus, and possible task-evoked and stimulus phase-locked frequency sliding.

In a random 25% of trials, the fixation turned purple at trial onset. In these trials, subjects were instructed to count to three in their head slowly and then press a button with their left or right hand, depending on whether the stimulus appeared in the left or right hemifield. The purpose of this manipulation was to engage a visual parietal-motor network involved in time-coordinated manual responses. This was used to test for modulations in functional connectivity between visual and motor areas. Sinusoidal luminance trials were not included in connectivity analyses, because the oscillating luminance may have naturally induced counting (with the “beats” of light). There were 100 trials per condition, and subjects received self-paced rest breaks every 3–4 min.

After the task recording, subjects completed an eyes-closed and an eyes-open rest recording (order was counterbalanced over subjects). Each recording was 2 min long, and subjects were instructed simply to relax and to try not to think about anything in particular.

EEG data were collected from 64 scalp channels and additional electrodes to monitor horizontal eye movements, using BioSemi equipment (see www.biosemi.com for hardware details). The sampling rate was 1024 Hz. The study was approved by the local ethics committee at the University of Amsterdam. The entire experiment lasted around 100 min, which included setup, task, and resting state.

##### EEG data processing and analyses.

EEG data were first high-pass filtered at 0.5 Hz. Task data were cut into epochs of −2 to +5 s peristimulus onset, and resting-state data were cut into nonoverlapping epochs of 2 s. Data were visually inspected, and any epochs containing excessive noise or other artifacts were rejected. Trials with horizontal eye movements were also rejected. Further cleaning was done through independent components analysis using the eeglab toolbox (Delorme and Makeig, 2004). Components reflecting eye blinks or muscle artifacts were visually identified and removed (one to five components were removed per subject).

Statistical thresholding was done through standard procedures for nonparametric permutation testing (Maris and Oostenveld, 2007), in which the mappings between condition label and data were randomly permuted 1000 times, and a *p* value was computed relative to this null distribution. Electrodes showing a monotonic effect of stimulus luminance on peak frequency, or connectivity via correlated frequency sliding, were considered significant if the *p* value was smaller than 0.05. In addition to single-electrode statistical thresholding, a cluster-based threshold (*p* < 0.05, distance of 10 cm) was applied to correct for multiple comparisons over the topography.

For the relationship between frequency sliding and stimulus luminance phase, statistics developed for phase-amplitude coupling were applied (Canolty et al., 2006). Specifically, the following equation was evaluated: |*n*^{−1}Σ*ae ^{ip}*

^{.}|, where

*n*is the number of time points,

*a*is the peak frequency at each time point, and

*p*is the phase of the 1 Hz stimulus luminance. The result of this equation reflects the extent to which peak frequency values are modulated by stimulus luminance phase. Next, this result was compared against a distribution of null hypothesis values, generated by shuffling the peak frequency values, thus generating a

*z*value corresponding to the likelihood of observing a modulation value at least that extreme under the null hypothesis. This procedure was done for each subject individually. Finally, for group-level statistical evaluations, a

*t*test was performed on the

*z*values.

Time-frequency power plots were obtained via convolution with a family of complex Morlet wavelets (Cohen, 2014). Power was converted to a decibel scale relative to a −500 to −200 ms prestimulus period.

Resting-state EEG data involved source reconstruction via linearly constrained minimum variance beamforming (Dalal et al., 2008), using analysis procedures and parameters that we have previously found to be appropriate for 64-channel scalp EEG (Cohen and Ridderinkhof, 2013). The idea of beamforming is to create frequency-specific spatial filters that allow estimation of activity in different brain voxels based on weighted combinations of the scalp-recorded electrode activity. The frequency band-specific time series was projected into 2004 brain-space voxels, and frequency sliding was subsequently computed at each voxel as described above. The standard MNI brain and standard electrode locations were used to compute a single leadfield that was applied to all subjects; for low spatial-resolution EEG, the MNI template brain provides reasonable localization accuracy (Fuchs et al., 2002), although the spatial precision is blurry compared with higher density recordings and subject-specific leadfield models. Large-scale networks were identified by first computing a voxel-by-voxel covariance matrix (separately for each 2 s epoch, then averaged together to increase signal-to-noise ratio), and then applying a principal components analysis via eigenvalue decomposition of that covariance matrix. The principal component weights for each voxel were then projected back to brain space.

Thresholding of principal components maps was done via permutation testing. At each of 200 iterations, the time series from each voxel was cut once at a random time point, and the second half was placed before the first half. Thereafter another principal components analysis was performed. From the distribution of principal component loadings at each voxel, the loading of the true principal components analysis (without permuted time series) was set to zero if it was smaller than 99% of the permutation distribution (thus, *p* < 0.01). Although with principal components analysis the loading at each voxel contributes to the total component structure, thresholding is useful to help visualize the spatial distribution of the components. Results were visualized using MRIcro software.

The relationship between frequency sliding and the power time series during resting state was tested by computing the Spearman correlation coefficient between power and the Gaussian transform of the frequency sliding, **e**^{−(f−f̄)2/4} where *f* refers to frequency sliding. The Gaussian transform was selected based on visual inspection of the nonlinear relationship between resting-state band-specific power and frequency sliding (Fig. 8*A*,*B*).

## Results

### Frequency sliding in artificial neural networks

Figures 2 and 3 illustrate the relationship between input strength and peak oscillation frequency, as well as the principle of frequency sliding, in artificial neural networks. Two networks were simulated using two distinct architectures for generating oscillations: one network relies on interactions between excitatory and inhibitory cells (Tiesinga and Sejnowski, 2009) to generate alpha-band and gamma-band oscillations (*E–I* model; Fig. 2) while the other network generates a theta-alpha-beta oscillation complex, typical of EEG data, due to its laminar organization (Thomson et al., 2003; Layer model; Fig. 3). Both networks generate endogenous frequency band-limited oscillations, even without oscillatory input. Despite their divergent neuron types, architectures, sizes, and intrinsic dynamics, both networks responded similarly to monotonic increases in input strength, by monotonically increasing their peak oscillation frequencies (Figs. 2*A*,*B*, 3*A*,*B*). In contrast to the close correspondence between peak oscillation frequency and input strength, oscillation power did not show a consistent monotonic relationship with input strength (Jia et al., 2013; Figs. 2*B*, 3*B*). Both networks quickly modulated their peak frequencies in response to time-varying input strength (Figs. 2*C*, 3*C*; input oscillations were well below the network oscillation frequencies, ruling out the possibility that the network oscillations were simply entrained to the input frequency).

These two networks were expanded to test whether correlated fluctuations in frequency sliding reflected functional connectivity. Three *E–I* networks received common input and independent noise, and network “1” provided input to network “2” (Fig. 2*D*); correlated frequency sliding correctly identified the simulated pattern of connectivity (Fig. 2*E*). In the Layer model, two columns were simulated with synaptic connections from excitatory neurons of layer 5/6 in column 1 to inhibitory neurons of layer 2/3 in column 2 (Fig. 3*D*; this simulates one type of “feedback” connection; Thomson et al., 2003). Again, correlated frequency sliding correctly identified both the intracolumnar connectivity as well as the layer-specific connectivity from column 1 to column 2 (Fig. 3*E*). The correlation coefficient between frequency sliding in layer 5/6 of column 1 (L5_{1}) and frequency sliding in layer 2/3 of column 2 (L3_{2}) was highly significant (average correlation coefficient across 100 simulations: 0.39, *t*_{(99)} = 15.84, *p* < 0.001). Furthermore, correlated frequency sliding between L5_{1} and L3_{2} was significantly stronger than correlated frequency sliding between L5_{1} and L4_{2} or L5_{2} (average correlation coefficients: 0.18 and 0.19, *t*_{(99)} = 8.30 and 7.37, *p*s < 0.001), and was also significantly stronger than correlated frequency sliding between L3_{2} and L3_{1} or L4_{1} (average correlation coefficients: 0.33 and 0.31, *t*_{(99)} = 2.05 and 3.89, *p*s = 0.0428 and <0.001). The simulations and analyses thus far validate that frequency sliding is a sensitive measure of input strength with high temporal precision, and that correlated frequency sliding can be used to measure functional connectivity between networks.

### Consequences of frequency sliding for single neurons

The next set of simulations was designed to understand what the implications of subtle changes in oscillation frequency might be for individual neurons. Two different biophysically detailed conductance-based neurons were simulated, a morphologically accurate neocortical layer 5 pyramidal cell (Hay et al., 2011), and a neocortical fast-spiking interneuron (Golomb et al., 2007). In the experiment, the neurons received oscillatory input on a dendrite while the somatic membrane voltage potential was recorded (Fig. 4*A*,*D*).

Relatively lower frequency input reduced the neurons' spike threshold, allowing the neurons to fire with weaker synaptic input. In both neurons, the decrease in spike threshold with lower input frequency was highly significant (Fig. 4*B*,*E*; see Fig. 4 for statistics from a linear model). Furthermore, relatively lower frequency input oscillations increased the temporal variability of action potentials (Fig. 4*C*,*F*). In addition to changes in the spike threshold, both neurons emitted action potentials in bursts coinciding with the positive phase of the input sine wave. Slower oscillations resulted in more spikes per burst (data not shown), consistent with previous studies (Kepecs et al., 2002). Spike count per burst is a discrete integer-based information coding scheme, whereas spike timing and variability provide a larger analog space for coding information with millisecond precision, particularly when this information is pooled over an ensemble of neurons (Havenith et al., 2011).

Figure 4 also shows that these two simulated cells were sensitive to different input frequency ranges, suggesting that neuronal morphology and biophysical characteristics result in maximal sensitivity (“tuning”) to particular ranges of input frequencies. This tuning appears to be more subtle than previously suggested (e.g., low vs high frequency; Nowak et al., 1997). This phenomenon would allow different neurons to be modulated by frequency sliding on different timescales.

The layer 5 model neuron was next tested in an experiment on coincidence detection, defined here as an action potential in response to two synaptic inputs when neither input is strong enough to elicit an action potential on its own. Synaptic inputs were given to different synapses (one delayed by 5 ms from the other) at various phases of the sinusoidal current (the current intensity was reduced such that the neuron did not produce somatic action potentials in response to the current alone; see Fig. 5*A* for an overview of this experiment). The results in Figure 5*B* confirm that lower frequency oscillations increased sensitivity to weak input. Specifically, action potentials in response to synaptic input were produced during earlier phases of input current oscillations with lower frequencies (see Fig. 5*B* for statistics). This relationship depended to some extent on synaptic input strength, with the relationship between oscillatory current and coincidence detection being stronger for stronger synaptic inputs (Fig. 5*B*).

This result was complemented by a moderate increase in the average number of action potentials elicited by the two synaptic inputs during lower frequency oscillations (Fig. 5*C*). Thus, two weak synaptic inputs will produce earlier and stronger neural responses when those inputs are embedded in a background of a relatively slower oscillatory current, although this effect depends in part on the strength of the synaptic inputs.

### Frequency sliding in human EEG task data

Having demonstrated that frequency sliding can be encoded in artificial neural networks, functionally group multiple networks, and modulate spike threshold and timing variability of biophysically plausible model neurons, the next step was to test whether frequency sliding occurs in real data and at a spatial scale larger than that of neural ensembles. Scalp EEG data were thus acquired from nine human volunteers performing a simple visual perceptual task, and a subsequent resting state. Subjects passively viewed large gray boxes in the left or right hemifield that varied in luminance (Fig. 6*A*). As expected based on the model results, EEG data demonstrated a monotonic relationship between input strength (stimulus luminance) and EEG peak oscillation frequency, a finding that was observed over stimulus-contralateral posterior electrodes (Fig. 6*B*). Although the task induced an alpha power suppression (Fig. 6*C*), there remained a strong peak in the alpha band during stimulus presentation (Fig. 6*D*).

In one condition, luminance varied sinusoidally over time to mimic the sinusoidal input provided to the *E–I* and Layer models. As expected based on the simulated neural networks, peak frequency in the human EEG was also modulated as a function of the phase of the sinusoidal stimulus luminance (Fig. 6*E*). The size of the peak frequency modulations was relatively small in magnitude compared with that during the static conditions, but the modulation of stimulus-contralateral frequency sliding by stimulus phase was significant across subjects for both hemifields (*t*_{(8)} = 5.59 and 5.5 for left and right hemifields; *p*s < 0.001; Fig. 6*E*, blue lines). This modulation was not significant in the static conditions, which for this analysis acted as a control (*t*_{(8)} = 0.79 and 0.49 for left and right hemifields).

In 25% of trials, the central fixation spot turned purple, and subjects counted to three and then pressed a button (“target trials”). Correlated frequency sliding, seeded from the electrodes identified in Figure 6*B*, revealed increased connectivity in a distributed parietal-motor network in target compared with nontarget trials (Fig. 6*F*). This pattern was observed for both left- and right-hemifield trials in the beta band (20–30 Hz). Qualitatively similar patterns of results were observed in the theta and alpha bands but were not statistically significant.

### Frequency sliding in human EEG resting-state data

The final experiment was to investigate whether endogenous frequency sliding in absence of experimentally induced input could be used to identify large-scale networks. Subjects rested quietly after the task, and frequency band-specific time series were estimated for 2004 brain voxels via beamforming (Dalal et al., 2008; Cohen and Ridderinkhof, 2013). Frequency sliding was then computed at each voxel, and principal components analysis was performed to identify “modes” of large-scale correlated frequency sliding (Fig. 7*A*). Axial and sagittal slices from modes 2–6 are shown in Figure 7*B*. This analysis used alpha-band frequency sliding, but the procedure could be applied to any frequency band.

The resting-state data also provided an opportunity to further investigate frequency sliding dynamics in absence of task effects (Fig. 8). First, frequency sliding showed a nonlinear relationship with simultaneous band-limited power, which was due to increased power variability near the center of each frequency band (Fig. 8*A*,*B*). This relationship was quantified by testing a fit between power and Gaussian-transformed frequency sliding for each subject, and then comparing the coefficients against zero at the group level. The Gaussian relationship was significant for theta-band (*t*_{(8)} = 9.77, *p* < 0.001) and alpha-band (*t*_{(8)} = 5.70, *p* < 0.001) frequency sliding, but not for beta-band (*t*_{(8)} = 1.13, *p* = 0.291) frequency sliding. Second, frequency sliding did not result from broadband shifts in the power spectrum, because frequency sliding in different frequency bands was weakly negative though significantly correlated (theta-alpha frequency sliding average correlation coefficient over subjects = −0.04, *t*_{(8)} = −2.89, *p* = 0.01; theta-beta frequency sliding average correlation coefficient = −0.015, *t*_{(8)} = −0.89, *p* = 0.199; alpha-beta frequency sliding average correlation coefficient = −0.09, *t*_{(8)} = −4.04, *p* = 0.001). Small negative correlations can be expected between neighboring frequency bands due to minor spectral leakage across the frequency boundaries. Finally, frequency sliding was compared with intersite phase clustering, a commonly used metric of frequency band-specific functional connectivity (Cohen, 2014). Correlated frequency sliding was significantly but not perfectly correlated with phase clustering (average correlation coefficients in theta/alpha band: 0.431/0.495, *p*s < 10^{−5}; Fig. 8*C*). This is not a surprising result: Perfect phase synchronization can occur only if there is perfect frequency locking (Amor et al., 2005).

There are many other methods of estimating brain functional connectivity from neural time series data, including power envelope correlations, Granger-based directional coupling methods, cross-frequency coupling, and mutual information. It would be interesting but beyond the scope of this paper to compare correlated frequency sliding with other measures of connectivity. However, given that frequency sliding is defined by the phase angle time series, and that its relationship with frequency band-specific power is tenuous, it is likely that correlated frequency sliding will correspond most closely with other phase-based connectivity methods and less strongly with connectivity methods that are independent of phase.

## Discussion

The present findings confirm previous observations that brain oscillation frequencies can change according to external stimuli and internal behavioral states (Lopes da Silva and Kamp, 1969; Arnolds et al., 1980; Ardid et al., 2010; Ray and Maunsell, 2010; Burns et al., 2011; Roberts et al., 2013), and extend these findings by providing putative neuronal consequences of subtle changes in oscillation frequency, and by demonstrating that frequency sliding is a sensitive measure of functional connectivity. In fact, these observations show that neural time series data violate one assumption of the analyses commonly used to analyze those data (such as Fourier- and wavelet-based methods), which is that that the frequency structure of an oscillator does not change over time, or at least changes more slowly than the time window of analyses (typically, several hundreds of milliseconds; Cohen, 2014). Frequency sliding represents a form of signal nonstationarity that is often overlooked, but may prove important for characterizing brain function and communication.

### Implications of frequency sliding for physiology and neural computation

The biophysically detailed models suggest at least two implications of small shifts in peak frequency within a frequency band for neurophysiological mechanisms of computation. First, neural networks may encode input intensity (both tonic as well as time varying) as changes in peak frequency of their oscillation, rather than in the power of those oscillations. Second, the modeling results suggest that neurons will emit somatic action potentials with lower levels of dendritic input when that input arrives with a relatively lower frequency oscillation. This results in part from slower input having a longer voltage integration time, and in part from voltage-gated Na^{+} and K^{+} currents (Azouz and Gray, 2000). Similar findings have been linked to the slope of an input (Azouz and Gray, 2000); the present results suggest that such input slopes may reflect oscillations (indeed, the slope of part of an oscillation simply reflects its frequency). The implications of this change in spike threshold are evident in detecting temporally coincidental synaptic inputs (Fig. 5), which is one of the key computational functions of pyramidal cells (London and Häusser, 2005).

Furthermore, relatively lower frequencies allow a longer time window in which action potentials can occur, thus increasing spike timing variability. The present findings show that the range of frequencies that shape spike timing variability is smaller than previously suspected (Nowak et al., 1997). Together, these findings suggest a novel hypothesis: frequency sliding can be a means of modulating neural excitability—a gain-control mechanism. Specifically, if increasing peak frequency causes neurons to require stronger input before firing, relatively faster oscillations thus facilitate cautious but accurate responding. In contrast, if decreasing peak frequency causes neurons to respond to weaker input, relatively slower oscillations facilitate fast but potentially noise-driven responding. This novel prediction could be empirically tested in humans and nonhuman animals, for example, during perceptual discrimination tasks with luminance-oscillating stimuli, or by varying the frequency of a transcranial alternating current stimulation.

Does the brain use frequency sliding to encode/decode information? This is an attractive possibility because frequency sliding is a property of networks and of individual neurons, changes dynamically over hundreds of milliseconds (corresponding to the time frame of many cognitive and perceptual processes), is maintained within canonical frequency bands (that is, frequency sliding can occur independently in theta, alpha, gamma, etc.), and correlates with external stimulus properties. On the other hand, it is also possible that frequency sliding does not contain any unique information beyond that already available in average spike rate. In this case, frequency sliding is a useful index of the input to a neural network and a novel measure of functional connectivity that is particularly useful when spike rate information is not available, such as in LFP or magnetoencephalography(MEG)/EEG. Additional empirical research is required to determine whether frequency sliding itself directly contributes to neural computation, or whether it is useful mainly as an analysis tool that can bridge research findings across computational models, *in vivo* and *in vitro* animal physiology, and human cognitive and clinical electrophysiology.

### Implications of frequency sliding for EEG/LFP

It is important to develop and use EEG data analysis techniques that are inspired by and can be linked to known neurophysiological properties of neurons and networks of neurons (Cohen and Gulbinaite, 2014). The well known relationship between input intensity and neuron firing rate provides a solid foundation for developing data analysis approaches that can be used and meaningfully interpreted across multiple spatial and temporal scales, ranging from individual neurons to large-scale networks of millions or billions of neurons spanning tens of centimeters. At present, power within a frequency band is the most commonly interpreted feature of EEG time-frequency results (Cohen, 2014). Though clearly a useful and widely used measure of neural activity, it has an unclear relationship with network input strength (Figs. 2, 3, 8; Jia et al., 2013), and may be more related to across-network coherence rather than within-network activity (Musall et al., 2014). Furthermore, the lack of stationarity in peak oscillation frequencies over time suggests that what previously may have been considered noise or uncertainty in the frequency precision of time-frequency results may instead reflect frequency sliding, which is not easily detected by standard time-frequency analysis approaches and which would thus decrease frequency precision in wavelet- or short-time FFT-based analyses.

The EEG task data confirmed predictions of the computational models: alpha-band peak frequency tracked the strength of the visual input, both in the static and in the dynamic conditions. It is unclear why the sinusoidal modulations in the dynamic condition were of smaller magnitude compared with the static conditions. It is possible that the luminance modulation was too weak or too easy to ignore (it was a task-irrelevant manipulation), or that the neural subnetworks most strongly affected by this manipulation were relatively small compared with the overall task-responsive networks. The EEG task data also confirmed that correlated frequency sliding can be used as a measure of task-related functional connectivity. The physiological mechanisms that produce such precise long-range co-modulations of peak frequency may be driven in part by pyramidal cell long-range inputs onto inhibitory interneurons, which in turn can precisely set the timing of network oscillations (Bush and Sejnowski, 1996; Chow et al., 1998; Maex and De Schutter, 2003).

The resting-state data showed that frequency sliding is an endogenous property of neural networks and does not require exogenous experimental manipulation. Resting-state networks are typically defined by correlated fluctuations over a slower timescale (seconds to minutes; Cabral et al., 2014), using the fMRI BOLD response or slow fluctuations in MEG frequency band-specific power fluctuations. As such, it is difficult to directly compare frequency sliding-defined resting-state networks with fMRI- or MEG-defined resting-state networks. Nonetheless, the networks shown in Figure 7*B* show some similarities to fMRI-defined resting-state networks, including positive and negative component loadings in the same network (respectively, red and blue colors), and networks with a spatial distribution at anterior and posterior medial areas. Although the BOLD response exhibits slow oscillations (Bluhm et al., 2007), it is unclear whether frequency sliding in hemodynamics (if it occurs) can be interpreted in a similar fashion as frequency sliding in neurophysiological data.

### Limitations

First, the precise relationship between frequency sliding and spiking activity is not fully understood, and may depend on neuron type (as was shown for two example neurons studied here) and additional factors such as background activity and synaptic strengths. Second, additional biophysical features of neurons may be modulated by frequency sliding but were not investigated here. For example, dendritic calcium waves propagate with a speed proportional to the input (Loewenstein and Sompolinsky, 2003). Voltage-gated synaptic dynamics including NMDA receptors or Na^{+} and K^{+} conductances may also be modulated by frequency sliding (Azouz and Gray, 2000, 2003). Third, source reconstruction of low-dimensional (<100 electrodes) EEG data involves spatial blurring that limits the spatial precision of the results. Furthermore, correlated frequency sliding does not address potential issues of volume conduction (although spatial filters such as beamforming or the scalp Laplacian, and condition comparisons, attenuate volume conduction artifacts). It is therefore likely that smaller subnetworks of correlated frequency sliding would be observed with high spatial resolution recordings, such as 300+ sensor MEG/EEG, or intracranial EEG in epilepsy patients.

### Conclusions

In conclusion, frequency sliding appears to be a fundamental principle of neural function that modulates spike threshold of individual neurons, encodes input strength in neural ensembles, indexes connectivity between neural populations, and groups large-scale networks during cognition and rest. Additional work is required to better understand the mechanisms of endogenous frequency sliding and its implications for neural computation, but it stands out as one of the few data analysis approaches that can be used to translate findings across myriad spatiotemporal scales of brain function (Varela et al., 2001; Le Van Quyen, 2011; Cohen and Gulbinaite, 2014; Pesenson, 2013).

## Footnotes

This research was funded by a VIDI grant from the Netherlands Organization for Scientific Research. Thanks to Rasa Gulbinaite for critical discussions.

The author declares no competing financial interests.

- Correspondence should be addressed to Michael X Cohen, Psychology Department, University of Amsterdam, 4 Weesperplein, 1018 XA, The Netherlands. mikexcohen{at}gmail.com