Beyond rhythm – a framework for understanding the frequency spectrum of neural activity

Introduction

The Fourier transform, or analogous methods, is routinely used to analyze the frequency content (i.e., the spectrum) of neural activity (Penttonen and Buzsáki, 2003; Bruns, 2004). Frequency spectra are highly sensitive to changes in the processes driving a signal (Percival and Walden, 1993). They can thus be used to detect subtle variations in the dynamics of the brain’s electric fields during attention, memory formation or retrieval, locomotion, and motor responses (Gray et al., 1989; O’Keefe and Recce, 1993; Skaggs et al., 1996; Fries et al., 2001; Girardeau et al., 2009; Niell and Stryker, 2010; Bosman et al., 2012; Vinck et al., 2015; Chen et al., 2017; Veit et al., 2017; Uran et al., 2022). However, the behavior of frequency spectra is often counterintuitive. For instance, the finite nature of recordings induces distortion of their spectral content (Percival and Walden, 1993). Furthermore, relating frequency spectra to continuously evolving neural processes is a complex problem. As a result, the interpretation of the frequency spectrum of neural activity is not straightforward.

Increased spectral power within some frequency bands is often believed to reflect oscillations. In fact, the words “oscillation” and “rhythm” are generally used to describe the subfield of neuroscience dedicated to study of the brain’s electrical activity (Buzsáki and Vöröslakos, 2023). Some neural patterns do indeed show strong periodicity. This includes well-characterized rhythms such as hippocampal theta (Vanderwolf, 1969; Winson, 1974; O’Keefe and Recce, 1993; Skaggs et al., 1996), thalamic spindles (Steriade, 2006; Niethard et al., 2018) and slow waves observed in the cortex during sleep (Steriade et al., 1993; Contreras et al., 1996; Sanchez-Vives and McCormick, 2000). Nonetheless, neuronal processes can also display non-periodic dynamics (Burns et al., 2011; Xing et al., 2012; Weber and Pillow, 2017; Naud and Sprekeler, 2018; van Ede et al., 2018; Donoghue et al., 2020; Steinmetz et al., 2021; Williams et al., 2021; Spyropoulos et al., 2022). In addition, most signals, periodic or not, have spectra where some frequencies are enhanced. Gaussian functions are non-periodic, yet exhibit a spectrum with a strong representation of low frequencies (Starosielec and Hägele, 2014). Thus, the notion that oscillations underlie changes in the spectrum of neural activity appears generally restrictive.

Here, we provide a didactic discussion for the often-intimidating interpretation of the frequency spectrum of neural recordings. Using a well-developed and broadly applicable theoretical framework (Cox and Lewis, 1966; Verveen and DeFelice, 1974; Rice, 1977; Bédard et al., 2006; Díaz et al., 2023), we show how neural patterns are usefully conceptualized as processes where discrete events can occur with varying degrees of periodicity. Using simulations, we illustrate how multiple processes can sum up in the local field potential (LFP) and how the basic properties of the Fourier transform shape the frequency spectrum of neural recordings. We show that the spectral signature of neural patterns depends essentially on 2 aspects: (1) the distribution of events in time and (2) the waveform of the voltage induced by individual events. Finally, we show examples of how this conceptual framework can be applied to decipher processes inducing changes in the spectral profile of real-world data acquired in the primary visual cortex of behaving mice.

Materials and methods

All simulations and analysis were performed in Matlab 2018b (Mathworks). The scripts used to generate the figures are available online. Simulated time series had a sampling rate of 1 kHz and a duration of 1000 s. Event timing pulse trains had a value of 1 at the time of events and zero everywhere else. Except in the case of perfectly periodic processes, event intervals were drawn from a gamma distribution:

P(t)=1/(b*aΓ(a))*t*a-1exp(-t/b)

where P(t) is the probability density at time t, Γ (x) is the gamma function of x (i.e., the gamma function is a generalization of the factorial to real numbers), and a and b are parameters determining the shape of the distribution. The process is sub-Poissonian if a < 1, Poisson if a = 1 and super-Poissionian if a > 1. Parameter b was always set to 1/(a * f) where f is the overall event frequency of the process. Values of a and f for all simulations are detailed in Table 1.

Table 1. Parameters used for simulation.

Recurring neural events were simulated by convolving event pulse trains with waveform traces of 500 ms. Various waveforms (i.e., impulse response functions) were considered. Synaptic events were modeled with the alpha function (Markram et al., 1997):

A⁢(t)=H⁢(t)*t*exp⁡(-t/tau)

Where A(t) is the alpha function at time t, H(t) is Heaviside step function, and tau is a parameter determining the time course of the synaptic response. Parameter tau was always set to 5.6 ms.

Spikes and spindles were modeled with real valued Morlet wavelets (i.e., the multiplication of a sinusoid with a gaussian window):

W(t)=exp(-(t-μ)/2(2*σ))*sin((t-t)0*(2*π)/p)

Where W(t) is the Morlet wavelet at time t, μ and σ are the mean and standard deviation of the gaussian component, t0 is the offset of the sine function and p is its period. Parameter μ was always set to zero (i.e., centered). Values of σ, t0, and p for all simulations are detailed in Table 1.

To study the effect of waveform width, waveforms were modeled as Hann window function (i.e., one cycle of a sinusoid) with varying periods p (Table 1).

Background synaptic noise was modeled by convolving white noise with a 500 ms waveform corresponding to the impulse response function of a passive resistive capacitive (RC):

F(t)RC=(1/tau)*exp(-t/tau)

Where FRC(t) is the value of the impulse response function of the RC filter at time t, and tau is a parameter determining the time course of the response (note the similarity to the alpha synaptic response). Parameter tau was set to 35 ms based on values inferred from recording of the membrane resistance and capacitance of excitatory neurons in the rat barrel cortex (Karagiannis et al., 2009).

The spectrum of waveforms was computed by tapering with a Hann window and taking the squared amplitude of the Fourier transform. The spectrum of recurring event time series was estimated with Welch’s method, that is by averaging the squared amplitude of the Fourier transform within non-overlapping 500 ms segments tapered with a Hann window (Percival and Walden, 1993).

The experiment described in Figure 4 is part of a data set used in Perrenoud et al. (2022). Briefly, a mouse was chronically implanted in the primary visual cortex (∼2.5 mm left and ∼4mm posterior from Bregma) with A16 probes (Neuronexus) having 16, 50 μm spaced, recording sites. The mouse was habituated to run head-fixed on a wheel and recorded in this configuration using a Digital Lynx SX system (Neuralynx). The method used for the extraction of gamma events is called CBASS and is described and freely available at (https://github.com/cardin-higley-lab/CBASS). CBASS ties increased power within a frequency band of the spectrum [here: gamma (30–80 Hz)] during a particular state (here: running) to discrete events in time. The LFP is band-pass filtered at the frequency band of interest. Candidate events are then selected at the trough of the filtered LFP in an arbitrary channel. The dynamics of candidate events across channels are parameterized and a variant of the K-means algorithm is used to find a cluster of events whose dynamics are more prevalent during the state of interest.

Results

To understand the frequency spectrum of neural recordings, it is useful to first consider the origin of the brain’s electric field. Neurons encode signals through variations of their membrane potential which are induced by the diffusion of Na+, K+, Cl− and Ca2+ ions through tightly controlled ion channels (Buzsáki et al., 2012). This movement causes variations in the electric field whose behavior is described by Maxwell’s equations. With a perfect knowledge of the brain’s electric (and magnetic) field, Maxwell’s equations might in principle, permit to infer movements of charges at the neuronal level. However, the sensitivity of our recording methods and the complexity of ion movement through the several billions of neurons that constitute the brain make this intractable in practice.

While the brain’s electric activity can give the impression of a continuously evolving chaos, it is made out of many repeating elements. These elements are the activation of synapses, the firing of neurons or, at a larger scale, coordinated firing and synaptic barrages within one brain region, or from one area to another (Womelsdorf et al., 2014). Each such event has a defined impact on the brain’s electric field and may thus induce a voltage deflection in a recording electrode. A recurring event, such as the firing of one neuron, can have of a strong influence on a signal depending on its magnitude and its proximity (Schmitzer-Torbert et al., 2005). We will thus start by considering how one recurring event affects the spectrum.

Impact of a recurring event’s temporal distribution on the Fourier spectrum

Let us first examine the influence on the spectrum of how often and regularly an event occurs. To assess this, we will idealize recurring events and treat them as point processes (Cox and Lewis, 1966; Perkel et al., 1967). Point processes are completely defined by their inter-event interval distribution. The simplest and most commonly occurring type of point-process is the Poisson process. Poisson processes are characterized by a decaying exponential inter-event interval distribution (Figure 1B) and have 2 interesting properties: (1) the probability of an event at any given time is invariant and entirely unaffected by the process’s history; (2) the variance of the number of events per unit of time is equal to the event rate.

Figure 1. Super-Poissonian, Poisson, sub-Poissonian and periodic point-processes have recognizable Fourier spectra. (1) Event interval distribution, (2) excerpt of an event train and (3) Fourier spectrum of panel (A) super-Poissonian, (B) Poisson (C) sub-Poissonian and (D) perfectly periodic point processes. Distinct inter-event interval distributions translate into specific energy distribution in the Fourier spectrum. All processes have a mean frequency of 40 Hz.

Point processes can be generally categorized by how regularly events occur compared to a Poisson process. Processes happening more regularly have a less variable inter-event interval. The variance of the number of events per unit of time in thus lower than the event rate. These processes are called sub-Poissonian (Figure 1C). As the variance of the event interval drops, the process begins to resemble a periodic process (Supplementary Figure 1). In the most extreme case, the variance of the inter-event interval is null, and the process is perfectly periodic (Figure 1D). Conversely, the variance of the number of events per unit of time can be superior to the event rate. In this case, the process is called super-Poissonian (Figure 1A). The firing of bursting neurons is an example of super-Poissonian process (Williams et al., 2021).

The gamma distribution (among others) can be parameterized to generate point processes (called gamma processes) having Poisson, super-Poissonian or sub-Poissonian inter-event intervals. To illustrate how the variance of point-processes affect their frequency spectrum, we parameterized gamma processes to generate super-Poissonian, Poisson and sub-Poissionian event time series. A perfectly periodic event times series was also constructed. All four time-series had a simulated sampling rate of 1000 Hz and a matched event rate of 40 Hz (Figure 1). The Fourier spectrum of each time series was then estimated using Welch’s methods with a 500 ms Hann window (Materials and Methods), a method classically used to estimate Fourier spectra in real data.

As can be seen in Figure 1B3, the spectrum of a Poisson train of event is perfectly flat. While the rate of a Poisson process affects overall power (energy per unit of time), energy remains evenly distributed across frequencies. The spectrum of a super-Poissonian gamma process is a decaying exponential and has thus more energy at low frequencies (Figure 1A3). Conversely, the spectrum of sub-poissonian gamma process shows reduced power at low frequency, a pronounced peak around the rate of the process and a flat energy distribution for higher frequencies (Figure 1C3). As the process becomes more regular secondary peaks start appearing at the harmonics of the processes rate (i.e., integer multiples of the core frequency, in our case 80 Hz, 120 Hz and so on; Supplementary Figure 1). For a perfectly periodic point process, the power is entirely concentrated and evenly distributed across the process’s rate and its harmonic (Figure 1D3). It thus appears that the event distribution of a process has a substantial and recognizable influence on the shape of its spectrum.

Periodic, Poisson and sub or super-Poissionian gamma processes are idealized processes that might not be representative of what happens in the brain. To illustrate how the timing of events affect the Fourier spectrum in more arbitrary cases, we built a compound inter-event interval distribution by adding the distributions a Poisson process and 2 sub-Poissonian gamma processes having distinct rate and variance (Supplementary Figure 2). Accordingly, the inter-event interval distribution showed two prominent modes on top of an exponential decay. We generated an event times-series following this distribution and for each summed process for comparison. The Fourier spectrum of each time-series was computed as described above.

The Fourier spectrum of the compound process (Supplementary Figure 2D3) closely resembled the sum of the Fourier spectra of the summed process (Supplementary Figures 2A3, B3, C3). Two main peaks were observed respectively at the frequencies of each sub-Poissonian process. These results illustrate that, for a point process, peaks in the spectrum are broadly indicative of modes (i.e., preferred intervals) in the inter-event interval distribution.

Impact of a recurring event’s waveform on the Fourier spectrum

So far, we have only considered idealized time series where the impulse response function (i.e., the waveform of the voltage induced by one event) of single events is a perfect pulse (i.e., having a value of 1 at the time of the event and zero everywhere else). In actual recordings of the brain’s field activity, repeating events will have more complex waveforms. We will now consider how different waveforms affects the Fourier spectrum in conjunction with the event interval distribution. For the sake of simplicity, we will only consider that the waveform of all events in a process does not vary in shape or amplitude.

A train of recurring events can be seen as the convolution of its waveform with a signal representing the events timing such as the point-processes considered in the previous section. Such processes have a well-developed theory and are often referred to as generalized shot noise (Schottky, 1918; Rice, 1977; Díaz et al., 2023). The Fourier transform has the remarkable property that the spectrum of the convolution of two time series is the product of their spectrum. In other words, convolution in the time domain translate into a simple multiplication in the frequency domain (Percival and Walden, 1993). Within the accuracy of estimation methods, the Fourier spectrum of a time series of recurring events is thus equivalent to the product of the spectra of the event’s waveform and that of its event timing pulse train.

To illustrate how this interaction shapes the Fourier spectrum, we considered three basic waveforms: 1) the alpha function, a popular choice to model synaptic responses (Figure 2B1); and 2 real-valued Gabor wavelet (the product of a sinusoidal and a gaussian window) shaped respectively to resemble 2) a spike (Figure 2C1) and 3) a spindle (Figure 2D1). All impulse response functions have specific Fourier spectra (Figures 2B2, C2, D2). The alpha function has a spectrum where energy is concentrated in lower frequencies. The spectrum of Gabor wavelets is centered on the frequency of their sinusoid component and has a bandwidth inversely proportional to the width of their gaussian window.

Figure 2. The Fourier spectrum of recurring event trains depends on the inter-event interval distribution and the waveform of single events. (1) (A) a Poisson process was convolved with 3 waveforms (top) mimicking the shape of panel (B) a synaptic event, (C) a spike and (D) a spindle, resulting in three distinct recurring event time series having the same event timing (bottom). (2) Convolution in the time domain translates into a simple multiplication into the frequency domain. Thus, the Fourier spectrum of the recurring event time series in B1, C1, and D1 (bottom) is simply the product of the spectrum of their waveform (top) and the spectrum of the Poisson pulse train (A2). Here, as Poisson processes have a flat spectrum, the spectrum of event times series is essentially determined by that of their waveform. The timing train (A) is identical for all processes and is not meant to be realistic.

Each impulse response function was then convolved with a Poissonian (Figure 2) or a periodic timing pulse train (Supplementary Figure 3). The Fourier spectrum of resulting time-series was estimated as before. As predicted by theory, all constructed event time-series had spectra nearly equal to the product of that of their timing pulse train and that of their waveform. Due to the flat shape of their timing pulse train’s spectrum, Poisson distributed event time-series have a Fourier spectrum that is essentially that of their waveform (Figures 2B2, C2, D2). Thus, as exemplified in the Poissonian spikes train constructed in Figure 2C, a pronounced peak in the spectrum can be inherited solely from the shape of a recurring event’s voltage response without periodicity in its timing.

Contrasting with Poisson trains, the energy of periodic timing pulse trains is concentrated at the process’s core frequency and its harmonics. Thus, in this case, the spectrum is more generally determined by the spectrum of the timing pulse train than by that of the waveform. However, as exemplified in Supplementary Figure 3B2, some waveforms will greatly reduce the power of a periodic timing pulse train when the energy distribution of their spectrum does not overlap.

The examples above illustrate that the spectrum of a train of recurring events is shaped both by the event timing (i.e., its timing pulse train) and the impact of single events on the signal (i.e., its impulse response function or waveform). How much each shapes the spectrum varies depending on the case. One general parameter will however influence the relative contribution of one over the other, namely the duration of the waveform relative to an events overall rate of occurrence.

To illustrate this, we considered a moderately periodic sub-Poissonian pulse train (Supplementary Figure 4). The train was then convolved with Hann functions (i.e., one cycle of a sinusoid and a popular window function for spectral estimation) of increasing width relative to the processes core frequency (40 Hz). For short Hann pulses (5 times shorter the average inter-event interval, Supplementary Figure 4B), the spectrum of the resulting event times series closely resembles that of the timing pulse train. However, for longer waveforms (5 times longer the average inter-event interval, Supplementary Figure 4D), the shape of the spectrum was mostly determined by that of the waveform.

These results indicate that neural processes with short waveforms, such as action potentials, have a spectrum that tends to be influenced by their timing (i.e., timing pulse train). Conversely, neural process with a longer waveform relative to their rate of occurrence, such as synaptic events, will have a spectral imprint that tends to be determined by the waveform of single events (i.e., impulse response function).

How do multiple neural processes sum up in the Fourier spectrum?

We have examined the Fourier spectrum of single recurring events occurring in isolation. Neural activity is composed of many such events. We will thus now consider how multiple recurring events contribute to the voltage recorded by an electrode and how they collectively shape the Fourier spectrum of neural recordings. Let us illustrate how this may happen in a simple recording configuration.

As can be derived from Gauss’s law (the first of Maxwell’s equation), the voltage deflection caused by multiple sources at an electrode is the sum of the voltage induced by each source. Thus, the voltage induced by neural events simply sums up in the local field potential. The Fourier transform of the sum of two signals is quite simply the sum of their Fourier transforms (Percival and Walden, 1993). As this summation occurs in the complex plane, signals may at times cancel each other when locked in opposite phases. However, when signals are independent, the power of their sum will approach the sum of their power (the power is the squared norm of the Fourier transform over time). Most often, multiple neural signals thus tend to add their energy distributions to the spectrum of neural recordings. The voltage deflection induced by a single source is inversely proportional to its distance to the recording electrode. Sources that are small and relatively far away from the recording site, such as synaptic events, tend to be indistinguishable and sum into a background signal. We will thus first simulate this background signal.

In awake cortical recordings, background synaptic events often occur at a sustained high rate, a regime that is called desynchronized (Steriade and Deschenes, 1984; Destexhe and Paré, 1999; Chance et al., 2002; Destexhe and Contreras, 2006; Haider and McCormick, 2009; Petersen and Crochet, 2013). As synaptic responses tend to be slower, the spectrum has properties that are mostly shaped by the waveform of synaptic events (i.e., impulse response function, previous section). Synaptic events display a characteristic 1/f2 power law energy distribution (Figure 2B2). This energy distribution is likely inherited in part from the biophysical properties of neuronal membranes which tend to behave as passive resistor-capacitor (RC) filters (Koch, 1984; Bédard et al., 2006; Freeman and Zhai, 2009), though it has received numerous interpretations (Bédard et al., 2006; Miller et al., 2009; Touboul and Destexhe, 2010; Gao et al., 2017; Donoghue et al., 2020; Schaworonkow and Voytek, 2021; O’Byrne and Jerbi, 2022). Here, we simply simulated the background LFP signal as white noise passing through an RC filter (Figure 3A1, Materials and methods). The spectrum of that times series was computed as before (Figure 3A2).

Figure 3. Neural processes sum their energy in the Fourier spectrum. (1) Excerpt, and (2) Fourier spectrum of panel (A) white noise filtered with a passive Resistive-Capacitive (RC) filter mimicking desynchronized background synaptic activity, (B) a Poisson spike train and (C) the sum of the time series in panels (A,B). The spectrum of background synaptic activity in A has a characteristic 1/f2 shape which is inherited from the RC filter. The spectrum of the Poisson spike train in B has a bell shape peaking at 75 Hz which is inherited from the spike waveform. The spectrum of the sum of two signals approaches the sum of their spectra.

Figure 4. The spectrum of gamma activity in the mouse visual cortex is shaped by the timing of gamma events. (A) Schematic of the recording configuration. A linear electrode array is used to record the local field potential (LFP) across the layers of the primary visual cortex (V1) of a mouse freely running on a wheel. (B) Locomotion (purple) induces an increase in LFP power in the gamma range (30–80 Hz) having a narrow band peak around 55–60 Hz and broad band component. (C) Excerpt of the recording around locomotion onset. Locomotion modulated gamma activity can be tied to discrete network event [orange; further detail on the method can be found in Perrenoud et al. (2022)]. (D) (1) Average LFP around gamma events and (2) its Fourier spectrum. The average LFP shows a broad band and a narrow band peak in its energy distribution. (E) (1) Inter-event interval distribution of gamma events and (2) Fourier spectrum of the event timing train of gamma events. The gray bar in E2 represents the 95% confidence interval of the spectra of matched Poisson trains (1,000 randomizations). The broad-band and a narrow-band components of locomotion modulated gamma activity are visible in the Spectrum of the event’s timing.

Recordings of the local field potential will also often be impacted by action potentials originating from one or several neurons (i.e., units) located in the vicinity (within ∼30–50 μm) of the recording electrode (Schmitzer-Torbert et al., 2005). The voltage deflection induced by these action potentials can sometimes have a significant impact on the spectrum. To illustrate this, we simulated neuronal firing by convolving a Poisson timing pulse train with a spike shaped Morlet wavelet (Figure 3B1). As the timing of the train is Poissonian, the distribution of energy of the spectrum is purely inherited from the waveform of the Morlet spike (i.e., the impulse response function). The spectrum of the process showed a bell-shaped energy distribution with a peak at 75 Hz (Figure 3B2).

The two time-series were then summed to simulate a common recording scenario where neuronal firing at the vicinity of the electrode adds to a background synaptic signal (Figure 3C1). The spectrum of the resulting times series was estimated as before (Figure 3C2). The spectrum of the sum of the two signal is nearly equivalent to the sum of their spectra. The Poisson spike train added to the spectrum of the background LFP an induced an upward inflection of the energy distribution peaking at 75 Hz. This simulation illustrates how the contributions of multiple neuronal processes can be added up in the Fourier Spectrum. These results also show that a relatively complex spectrum can arise from the simple addition of non-periodic processes.

Quasi-periodic processes can be modeled as recurring events

As illustrated in the example above, background synaptic activity has considerable influence on the shape of the LFP spectrum, especially at frequencies under ∼50 Hz (Buzsáki et al., 2012). We have also seen that under a desynchronized regime, the spectral energy of this background activity tends toward a 1/f2 distribution (Donoghue et al., 2020; Schaworonkow and Voytek, 2021). However, background synaptic activity will at times include patterns of synchronized synaptic barrages (Buzsáki et al., 2012; Buzsáki and Vöröslakos, 2023). These barrages tend to target specific subregions of the somato-dendritic compartment and thus, to induce return currents as neurons equilibrate their membrane potential. In oriented structures, such as the cortex and hippocampus, return currents add up and greatly amplify the voltage induce by barrages in the LFP (Buzsáki et al., 2012). Thus, localized synaptic barrage often have a strong influence on the spectrum.

Some well-characterized patterns of synaptic barrages show strong periodicity. These include thalamocortical spindles (Steriade and Deschenes, 1984; Timofeev et al., 1996; Steriade, 2006), the theta rhythm of the hippocampus (Vanderwolf, 1969; O’Keefe and Recce, 1993; Skaggs et al., 1996) or cortical activity during slow wave sleep (Steriade et al., 1993; Contreras et al., 1996). At times neural patterns having distinct frequencies may also interact (Canolty and Knight, 2010; Cohen, 2017). It is often natural to treat quasi-periodic or transiently periodic neural pattern as oscillations. However, we will illustrate how it is also useful and generally applicable to conceptualize these patterns as a train of recurring events such as those described in the previous sections.

The theta rhythm of the hippocampus (∼7 Hz) is remarkably regular and is generally observed during locomotion (Vanderwolf, 1969; O’Keefe and Recce, 1993; Skaggs et al., 1996). Here we modeled hippocampal theta by convolving a quasi-periodic sub-poissonian timing train having a mean frequency of 7 Hz with an inverted Hann function having a matching period (Table 1 and Supplementary Figure 5A1). The spectrum of this process shows a peak at 7 Hz and closely resemble a sine wave (Supplementary Figure 5A). Hippocampal theta activity is often accompanied by high frequency activity in the gamma range (30–80 Hz) (Bragin et al., 1995; Csicsvari et al., 2003). Gamma activity generally occurs during the downward phase of theta activity (Lisman and Jensen, 2013). This modulation of high frequencies by lower frequency activity is known as cross frequency coupling (Canolty and Knight, 2010). We model theta-nested gamma activity as a smaller inverted Hann function driven by a quasi-periodic timing train having a mean frequency of 80 Hz. Events occurring during the upward phase of theta were dropped (Supplementary Figure 5B1). The spectrum of theta-modulated gamma showed prominent peaks in the gamma range (∼80 Hz), as well as theta (7 Hz) and its first harmonic (15 Hz) (Supplementary Figure 5B2). This illustrates how high frequency processes modulated by a lower frequency display increased energy at that low frequency simply because of the longer latencies that this modulation induces in their event timing. As for the simulation displayed in Figure 3, the spectrum of the sum of theta and theta-modulated gamma closely approaches the sum of their spectra (Supplementary Figure 5C).

Identifying recurring patterns of neural activity in vivo

The temporal dynamics of synaptic barrages and their influence on the spectrum tend to vary from region to region and involve specific pathways. Characterizing these patterns and their influence on neural processing is a subject of active research (Girardeau et al., 2009; Cei et al., 2014; Veit et al., 2017, 2023; Uran et al., 2022; Buzsáki and Vöröslakos, 2023). Conceptualizing neural patterns as train of recurring events can also be a powerful way to explore their properties in vivo. To illustrate this, we will use data that we have acquired in a recent study focusing on activity induced in the gamma range (30–80 Hz) in the visual cortex of mice (Perrenoud et al., 2022).

Gamma activity arises in many brain regions under varying behavioral contingencies (Bosman et al., 2012; Cardin, 2016; Chen et al., 2017; Uran et al., 2022; Fernandez-Ruiz et al., 2023). In the mouse visual cortex, a specific increase of gamma activity is observed during locomotion (Niell and Stryker, 2010; Vinck et al., 2015). This activity comprises a narrow band component around 60 Hz (Saleem et al., 2017; Shin et al., 2023) and a broad band component starting from ∼30 Hz and expending to high frequencies (Figure 4B). In a recent study (Perrenoud et al., 2022), we took advantage of multichannel electrode arrays to show that this spectral increase can be tied to recurring events having a specific pattern of propagation across cortical layers (Figures 4A, C, Materials and methods). Gamma events occurred more frequently during locomotion but also happened at a sustained rate during quiescence (Figure 4C). Decomposing gamma activity as an event train allows examination of how it entrains cortical neurons in relation with behavior with a high temporal resolution (Perrenoud et al., 2022) and provides some insight into how its spectral signature arises and contributes to the LFP.

Gamma events showed a complex inter-event distribution with a prominent mode around 17.8 ms, showing that gamma events tend to occur in rapid sequences or bursts (Figure 4E1). Averaging the field potential around events shows a rhythmic pattern of activity (Figure 4D1), whose spectral energy distribution resembles that observed during locomotion (taking away the 1/f2 RC filtering induced by neural membrane, Figure 4D2). As gamma events shift phase across recording locations (Figure 4D1), it is difficult to precisely isolate the waveform of single events. However, we can now estimate how much of the shape of cortical gamma activity is determined by event timing. To do this we computed an idealized perfect pulse times series having a value of one at the time of events and zero everywhere else. The spectrum of that time series was computed as described above. We found that the characteristic broad-band and narrow-band components of the distribution of energy of gamma activity can be derived solely by considering the events timing (Figure 4E2). To estimate the significance of this spectral distribution, we compared it to the spectrum of Poisson trains having a matching number of events (Figure 4E2, 1,000 randomizations). The broad band and narrow band gamma components of the spectrum were both well outside the 95% confidence interval of randomized trains. Gamma activity in mouse V1, is thus a real-world example of how a complex spectrum can be shaped by the temporal distribution of recurring neural events.

Discussion

Neural activity has highly complex dynamics (Buzsaki, 2006; Buzsáki et al., 2012). The Fourier transform is a remarkably powerful tool allowing the detection of subtle changes in neural processes (Percival and Walden, 1993; Bruns, 2004). However, relating variation in the Fourier spectrum to actual neural processes is a complex problem. Thus, the interpretation of the spectrum of neural activity is rarely straightforward.

In the present study, we use a simple yet generally applicable conceptual framework to guide the interpretation of the Fourier spectrum of neural recordings (Cox and Lewis, 1966; Rice, 1977; Bédard et al., 2006; Díaz et al., 2023). We note that neural activity is made of recurring events such as synaptic currents, action potentials and, at larger scales, coordinated synaptic barrages within a region or from one region to another. This allows to model neural processes straightforwardly as the convolution of (1) a timing train and (2) a waveform modeling the voltage deflection induced by single events (Rice, 1977; Bédard et al., 2006). Using simulations, we show how the spectrum of recurring neural events depends critically on these two aspects. We then illustrate how multiple recurring events tend to sum up in the LFP and how this shapes the spectrum in a typical recording scenario. Finally, we illustrate how this framework can be used to describe the dynamics of gamma activity (30–80 Hz) in real data obtained in the mouse visual cortex (Perrenoud et al., 2022).

Increased energy in the frequency spectrum within some frequency bands is often interpreted as the reflection of oscillatory neural activity. Accordingly, neural activity is frequently studied with concepts related to oscillations such as phase, amplitude, and coherence (Vinck et al., 2023). Some neural patterns do show strong periodicity (Buzsaki, 2006; Buzsáki and Vöröslakos, 2023), including slow wave sleep cortical activity (Steriade et al., 1993; Contreras et al., 1996; Sanchez-Vives and McCormick, 2000), theta rhythm in the hippocampus (O’Keefe and Recce, 1993; Skaggs et al., 1996) or thalamic spindles (Steriade, 2006; Niethard et al., 2018). Treating these patterns as oscillatory is thus a powerful way to quantify their properties (Bruns, 2004). However, these notions will not appropriately capture the properties of non-periodic signals (Donoghue et al., 2020). When then, might it be appropriate to treat neural patterns as oscillations?

Conceptualizing neural process as recurring events can a powerful way to address this question. Recurring events can be modeled simply and straightforwardly as the convolution of a timing train and a waveform (also called impulse response function) (Bédard et al., 2006; Díaz et al., 2023). The theory underlying such models is well developed (Cox and Lewis, 1966; Perkel et al., 1967; Rice, 1977) and has found a wide variety of applications starting with the description of shot noise by Schottky (1918). Accordingly, such processes are often referred to as generalized shot noise (Rice, 1977). Here, this framework allowed us to investigate how the spectrum is shaped by processes having arbitrary degrees of periodicity. Our work indicates that a tendency toward periodicity will indeed translate into increased power within a defined frequency band. However, we also provide simple examples of how band-specific power increase can arise from process that are entirely non-periodic. This indicates that some caution is warranted when concluding that peaks in the spectrum are the signature of oscillations or rhythmicity. One must indeed make sure that there is some energy in the spectrum within a specific frequency band. However, it is important that verification be performed in the time domain. A simple check for consistency in the amplitude and phase of the filtered signal can be an appropriate way to address these concerns.

What to do then when neural signals do not show strong signs of periodicity? Here, with the example of one of our recent studies (Perrenoud et al., 2022), we illustrate how treating of neural process as recurring events can also yield key insights into their properties (Figure 4). Thinking of neural patterns as recurring events is a straightforward way to link them to tangible neural processes such as synaptic responses, action potentials, or synaptic barrages. It also permits direct quantification of how the distribution of events in time might shape the spectrum. When the waveform of the voltage induced by single events can be estimated with precision (as for action potentials) it is potentially possible to estimate the full contribution of a recurring event to the spectrum.

Identifying repeating events in the local field potential and relating them to neural processes is not a simple problem. We and others have had success in identifying recurring events by taking advantage of multichannel recordings to resolve regularity in the propagation of neural activity in space (Perrenoud et al., 2022; Sibille et al., 2022). The detection methods that we used in our case, functions on the premise that events manifest as single cycles in the filtered LFP (Perrenoud et al., 2022). Thus, it may not be broadly applicable to detect multicycle events such as spindles (Timofeev et al., 1996; Olbrich et al., 2011; Niethard et al., 2018). Further methodological development is needed in this direction. Recent progress in increasing the spatial resolution of multielectrode recordings may make this approach more feasible in the future (Jun et al., 2017; Steinmetz et al., 2021). Another important limitation of the framework introduced here is that, for the sake of clarity, we have only considered waveforms that are invariant in shape and amplitude. In a real-world scenario there might be variability in the intensity and shape of the waveform of recurring neural patterns over time. However, the simple conceptual framework introduced here may serve as a useful starting point to guide developments addressing these questions in the future.

As a final note, for the sake of brevity, we have not discussed the relationship between the spectrum of a process and its autocorrelation function (Lamarre and Raynauld, 1965; Perkel et al., 1967; Lang et al., 1999; Lévesque et al., 2020). Readers wishing to expend on the theory of spectral estimation can find extended discussion in the book ‘Spectral Analysis for Physical Applications’ (Percival and Walden, 1993).

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.5061/dryad.crjdfn394.

Ethics statement

The animal study was approved by the Institutional Animal Care and Use Committee (IACUC). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

QP performed the simulations and wrote the manuscript. JC acquired the funding, supervised the project, and critically reviewed the manuscript. Both authors designed the research, contributed to the article, and approved the submitted version.

Funding

This work was supported by funding from the NIH (EY022951 and MH113852 to JC, EY026878 to the Yale Vision Core), a McKnight Scholar Award (to JC), an award from the Kavli Institute of Neuroscience (to JC), and a BBRF Young Investigator Grant (to QP).

Acknowledgments

We are grateful to Garrett Neske and Tom Morse for helpful discussion and criticism.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2023.1217170/full#supplementary-material

Supplementary Figure 1 | Transition from a Poisson to a periodic point process. (1) Event interval distribution, (2) excerpt of an event train and (3) Fourier spectrum of panel (A) a Poisson process, (B–D) sub-Poissonian processes having diminishing variance in their inter-event interval distribution and (E) a perfectly periodic point process. As the variance of the inter-event interval diminishes, secondary peak start to appear at the harmonics (i.e., integer multiples) of the mean frequency. All processes have a mean frequency of 40 Hz.

Supplementary Figure 2 | Peaks in the spectrum of a point process reflect peaks in the event interval distribution. (1) Event interval distribution, (2) excerpt of an event train and (3) Fourier spectrum of panel (A) a Poisson process (mean frequency 50 Hz), (B,C) two sub-Poissonian processes (mean frequency, 80 and 20 Hz, respectively) and (D) a mixture of the processes in panels (A–C). Peaks in the Fourier spectra of point processes indicate the presence of some preferred event interval.

Supplementary Figure 3 | The Fourier spectrum of recurring event trains depends on the inter-event interval distribution and the waveform of single events. (1) (A) a perfectly periodic process was convolved with 3 waveforms (top) mimicking the shape of panel (B) a synaptic event, (C) a spike and (D) a spindle, resulting in three distinct recurring event time series having the same event timing (bottom). (2) Convolution in the time domain translates into a simple multiplication into the frequency domain. Thus, the Fourier spectrum of the recurring event time series in B1, C1 and D1 (bottom) is simply the product of the spectrum of their waveform (top) and the spectrum of the periodic pulse train (A2). The timing train (A) is identical for all processes and is not meant to be realistic. Note that the spindle waveform (D) and the regular time train in panel (A) have none overlapping spectra and cancel each other.

Supplementary Figure 4 | The width of the event waveform relative to the process rate influences how much each shapes the Fourier spectrum. (1) (A) a quasi-periodic sub-Poissonian process was convolved with 3 Hann waveform (top) having widths (B) 5 time shorter, (C) equal, and (D) five time longer than the mean period of the process, resulting in three distinct recurring event time series having the same event timing (bottom). (2) Convolution in the time domain translates into a simple multiplication into the frequency domain. Thus, the Fourier spectrum of the recurring event time series in B1, C1, and D1 (bottom) is simply the product of the spectrum of their impulse response functions (top) and the spectrum of the pulse train (A2). Shorter event waveforms make the spectrum look more like that of the event timing whereas longer waveforms tend to dominate the spectrum.

Supplementary Figure 5 | Complex spectra arise from cross-frequency coupling. (1) Excerpt, and (2) Fourier spectrum of panel (A) a 7 Hz quasi-periodic sub-poissonian event train convolved with a matched Hann window simulating hippocampal theta, (B) a 80 Hz train constructed similarly but where event occurrence is modulated by theta phase, and (C) the sum of the time series in panels (A,B). Cross frequency modulation induces a clear peak in the spectrum of gamma at theta (7 Hz) and its first harmonic (15 Hz). The spectrum of the sum of two signals approaches the sum of their spectra.

Footnotes References

Bédard, C., Kröger, H., and Destexhe, A. (2006). Does the 1/f frequency scaling of brain signals reflect self-organized critical states? Phys. Rev. Lett. 97:118102. doi: 10.1103/PhysRevLett.97.118102

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosman, C. A., Schoffelen, J. M., Brunet, N., Oostenveld, R., Bastos, A. M., Womelsdorf, T., et al. (2012). Attentional stimulus selection through selective synchronization between monkey visual areas. Neuron 75, 875–888. doi: 10.1016/j.neuron.2012.06.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Bragin, A., Jandó, G., Nádasdy, Z., Hetke, J., Wise, K., and Buzsáki, G. (1995). Gamma (40-100 Hz) oscillation in the hippocampus of the behaving rat. J. Neurosci. 15(1 Pt. 1), 47–60. doi: 10.1523/JNEUROSCI.15-01-00047.1995

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruns, A. (2004). Fourier-, Hilbert- and wavelet-based signal analysis: Are they really different approaches? J. Neurosci. Methods 137, 321–332. doi: 10.1016/j.jneumeth.2004.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Burns, S. P., Xing, D., and Shapley, R. M. (2011). Is gamma-band activity in the local field potential of V1 cortex a “clock”; or filtered noise? J. Neurosci. 31, 9658–9664. doi: 10.1523/JNEUROSCI.0660-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the brain. Oxford: Oxford University Press.