A general, noise driven mechanism for the 1/f-like behavior of neural field spectra

1. Introduction

Transient oscillations are a prominent feature of macroscopic neural field activity (Buzsaki, 2011), linked to brain function (Buzsaki & Draguhn, 2004) and dysfunction (Gibbs & Lennox, 2002; Uhlhaas & Singer, 2012). Oscillations appear as narrowband increases in the spectrum above an aperiodic background in which the power P decreases proportional to the frequency f raised to an exponent $\beta :P\propto f^{\beta }$. Characterizing oscillations while accounting for the aperiodic background is important for understanding neural spectra (Donoghue et al., 2020). Sophisticated methods (Donoghue et al., 2020; Wen & Liu, 2016; L. E. Wilson et al., 2022) support estimation of the 1/f-like, scale-free (He, 2014) or power-law (Newman, 2005) behavior of neural field spectra. Changes in $\beta$, the aperiodic exponent, have been investigated in many domains, including sleep (Bódizs et al., 2021; Freeman & Zhai, 2009; G. Horváth et al., 2022; Lendner et al., 2020), aging (Cesnaite et al., 2023; Schaworonkow & Voytek, 2021; Voytek et al., 2015), and disease (Lanzone et al., 2022; Numan et al., 2022; Robertson et al., 2019). While many factors impact estimation of the aperiodic exponent (e.g., the frequency range analyzed (Gerster et al., 2022; Stumpf & Porter, 2012)), values of the exponent reported at higher frequencies (>20 Hz) typically range between −4 and −2 (Table 1).

The universal observation of 1/f-like neural field spectra with a restricted range of aperiodic exponents across different recording modalities, experiments, and neural systems suggests a common generative mechanism. However, many complex, specific mechanisms have been identified to generate this phenomenon. These include excitatory/inhibitory balance (Gao et al., 2017), low-pass frequency filtering by dendrites (Buzsáki et al., 2012) or the extracellular medium (Bédard et al., 2006), nonideal resistive components in the cell membrane (Bédard & Destexhe, 2008), stochastic firing of neurons convolved with an exponential relaxation process (Bédard et al., 2006; Miller et al., 2009; Milstein et al., 2009), stochastic synaptic conductances (Rudolph et al., 2005), stochastically driven damped oscillators with different relaxation rates (Evertz et al., 2022), local homogenous connectivity (Jirsa, 2009), combinations of many transient oscillations at different frequencies and amplitudes (He et al., 2010), or network mechanisms linking slower rhythms to broad neuronal recruitment and therefore larger amplitude field potentials (Buzsáki et al., 2012). Theoretically, scale-free phenomena – with 1/f-like behavior – have been linked to fractal properties (Pritchard, 1992), critical transitions (Newman, 2005; O’Byrne & Jerbi, 2022), and self-organized criticality (Bak et al., 1987; Cocchi et al., 2017). How these proposed biological and mathematical mechanisms contribute – separately or combined – to the range of aperiodic exponents observed across diverse neural field spectra remains unclear.

Here we demonstrate how, in general, a dynamical system with stochastic drive generates 1/f-like behavior at higher frequencies in neural field spectra. We show that two noise terms – representing correlated and uncorrelated noise inputs – produce the range of aperiodic exponents observed in vivo. We illustrate these general results in nonlinear models of neural and non-neural activity to demonstrate the ambiguity in determining the specific mechanisms given only the observed 1/f-like behavior in the spectrum. While more complex underlying mechanisms may exist, we instead illustrate how the range of aperiodic exponents observed in vivo occurs in general for dynamical systems with stochastic drive, without requiring a specific biological mechanism or tuning to collective critical behavior.

2. Results

As a general model of neural activity, we consider the n-dimensional dynamical system:

where, for each $k=\{\mathrm{1,2},\dots ,n\}$, $X_k$ is a 1-dimensional variable, $f_k$ is a nonlinear function, $B_{kj}$ is a constant, and ${ϵ}_{X,j}$ is an independent Gaussian white noise source with mean 0 and variance ${\sigma }_{X,j}^2$. We assume the variable $X_1$ (i.e., $k=1$) is an observable quantity (e.g., the voltage recorded in the EEG or LFP) and all other variables ($X_2,X_3,\dots ,X_n$) represent n-1 unobserved or latent variables impacting the observable dynamics. The unspecified model (1) is general and therefore consistent with diverse models of neural activity. To derive the main result from this general model requires no specific biophysical mechanism. In what follows, we illustrate these general results by making specific model choices for the variables (X_k) and functions (f_k).

We assume an equilibrium exists in the noise-free model so that $\left(\frac{d{X}_1}{dt},\frac{d{X}_2}{dt},\dots ,\frac{d{X}_n}{dt}\right)=(\mathrm{0,0},\dots ,0)$ at $\left(X_1,X_2,\dots ,X_n\right)\equiv {\overrightarrow{X}}^0$. Near this equilibrium, the dynamics for one variable (X_k) of the nonlinear system (1) can be approximated by the corresponding linear system,

where $x_k$ represent small deviations of the variable $k$ from the equilibrium $X_k^0$, we evaluate the partial derivatives of the nonlinear function $f_k$ at the equilibrium ${\overrightarrow{X}}^0$ (Guckenheimer, John & Holmes, Philip, 1983; Izhikevich, E, 2007), and we include the same stochastic perturbations as in the nonlinear system (1). We express the system (2) as

where $\overrightarrow{x}$ is the $n$-by-1 vector of deviations from the equilibrium, $\overrightarrow{ϵ}$ is the m-by-1 vector of independent $\delta$-correlated Gaussian white noise sources, A is the n-by-n Jacobian matrix of the nonlinear system (1) evaluated at the equilibrium, and B is the n-by-m noise matrix,

where the parameter $B_{ij}$ determines the contribution of noise source j to variable i. We assume that A satisfies the conditions required for the linear system (3) to accurately approximate the dynamics of the noise-free nonlinear system (1) near the equilibrium (i.e., we assume A has no eigenvalues with zero real part and the equilibrium is therefore hyperbolic (Guckenheimer, John & Holmes, Philip, 1983; Izhikevich, E, 2007)).

For the linear system (3), the cross-spectral matrix $\mathit{S}[\omega ]$ can be obtained from the expression (Gardiner, C, 2004; Kleeman, 2011; Thomas & Lindner, 2019),

where A is real, I is the identity matrix, $i=\sqrt{-1}$, and $\omega =2\pi f$ is the frequency. Evaluating the asymptotic behavior of the cross-spectral matrix at high frequencies ($\omega$ larger than any frequency associated with a natural rhythm of the linear system in (3)), the spectrum of the observable variable X₁ is

where $𝒪\left({\omega }^k\right)$ indicates the limiting behavior of the spectrum as a function of the k^th power of $\omega$ as $\omega \to \mathrm{\infty }$; see Appendix.

The result in (5) shows that, at high frequencies, the aperiodic exponent $\beta$ depends on the relative noise to the observable variable. Without stochastic drive to the observable variable X₁ (i.e., with B_1k = 0 for all k in {1,2, … m}),

so that the aperiodic exponent $\beta =-4$ at high frequencies (i.e., as $\omega \to \mathrm{\infty }$). Alternatively, with stochastic drive to the observable variable X₁ (i.e., with $B_{1k}\ne 0$ for any k in {1,2, … m}),

so that the aperiodic exponent $\beta =-2$ at high frequencies (i.e., as $\omega \to \mathrm{\infty }$). We note that, in this case, random walk dynamics dominates the spectrum at high frequencies, with well-known power-law behavior (e.g., see (Milotti, 2002)).

To summarize, we consider a general, n-dimensional nonlinear dynamical system with stochastic drive (1). We assume an equilibrium exists in this system with dynamics well-approximated by the linearized system (2). Near this equilibrium, the dynamics produce aperiodic exponents between −4 and −2 at high frequencies (i.e., frequencies beyond the natural frequencies or spectral peaks of the system) consistent with the range of values observed in vivo. The value of the aperiodic exponent depends on the noise in the observed and latent variables; when noise in the observable variable X₁ dominates, the aperiodic exponent ≈ −2, while when the noise in the latent variables dominates, the aperiodic exponent ≈ −4.

The main result (5) and implications for the aperiodic exponent ($-4\le \beta \le -2$) are for the general model (1). These general results do not require a specific biophysical model of neural activity. In what follows, we illustrate the generality of these results in four example models, in which we choose the nonlinear functions f_k in model (1) and assume a square diagonal noise matrix B to simplify the presentation of numerical simulation results. In doing so, we show that each model produces aperiodic exponents consistent with in vivo data ($-4\le \beta \le -2$) but with different physical interpretations.

3. Discussion

A universal feature of neural field potential spectra is 1/f-like scaling at high frequencies. To produce this power-law, many generative mechanisms have been proposed with diverse biological implementations and interpretations. Here, we do not propose a specific biological mechanism produces 1/f-like scaling. Instead, to understand how different neural mechanisms can produce similar 1/f-like scaling, we consider a general nonlinear dynamical system with stochastic drive. We show that dynamics near a (hyperbolic) equilibrium in this general model of neural activity produce aperiodic exponents between −4 and −2, consistent with the range of values reported in vivo for higher frequencies (e.g., >20 Hz). We illustrate these results in neural and non-neural models. We propose that the range of aperiodic exponents observed across recording modalities, experiments, and neural systems is a natural consequence of a noise driven dynamical system.

We considered here a single statistic – the aperiodic exponent – reflecting the 1/f-like feature of the neural field spectrum. We note that many different models can explain the same observed statistic. For example, observations from neural systems can produce spectra with broadband peaks in the gamma band (approximately 30-80 Hz) (Fries et al., 2007). Many approaches exist to explain these observed spectral peaks, including statistical approaches (e.g., an autoregressive model of order two (Spyropoulos et al., 2020)), mechanical approaches (e.g., a damped driven oscillator (Spyropoulos et al., 2020)), or biophysical approaches (e.g., the interneuron network model, or the pyramidal-interneuron network model (Whittington et al., 2000)). Our understanding of the gamma rhythm in a particular experiment depends on the model choice. In the same way, observations from neural systems produce spectra with 1/f-like scaling at high frequencies, and many models exists to explain this scaling. Our results show mathematically why many models can produce the range of aperiodic exponents observed in vivo ($-4\le \beta \le -2$). Due to the general nature of this mathematical result, it is not surprising that many different proposals exist to explain the 1/f-like spectrum. In practice, we expect biological mechanisms must exist to create the scaling observed in vivo. Identifying these biological mechanisms, and their expression in the diverse observations of 1/f-like scaling reported, remains an important challenge.

Our results are consistent with previous work showing that power-law scaling occurs in simple stochastic or physical systems. For example, in (Gao et al., 2017) the authors relate estimates of the aperiodic exponent to changes in the balance between excitation and inhibition. To do so, the authors simulate the local field potential as the summed synaptic current generated by independent stochastic spiking excitatory and inhibitory cells. In (Bédard et al., 2006) the authors propose that 1/f-like scaling in the spectrum does not rely on critical states, but instead depends on the filtering properties of the extracellular medium (although this mechanism remains debated (Logothetis et al., 2007)). In (Touboul & Destexhe, 2010) the authors show that simple models of stochastic processes (high-frequency shot-noise processes or Ornstein-Uhlenbeck processes) produce peak-amplitude distributions consistent with power-law distributions. In (Priesemann & Shriki, 2018) the authors show that inhomogeneous Poisson processes can produce approximate power law distributions in the size and duration distributions of avalanches (i.e., activity cascades). Consistent with our results, these examples generate power-law scaling without requiring a sophisticated biological or mathematical mechanism. Distinct from these previous works, we consider an (unspecified) n-dimensional dynamical system and show that, near a (hyperbolic) equilibrium, stochastic drive produces aperiodic exponents consistent with values observed in vivo.

Under the general framework considered here, the range of aperiodic exponents reflects different types of noise. The observable dynamics (e.g., variable X₁ in (1)) depend directly on the stochastic drives to the observable variable, and indirectly on the stochastic drives to the latent variable(s). The latent dynamics introduce correlations in the uncorrelated latent noise process before this noise reaches the observable dynamics. Therefore, the observable dynamics depend on both uncorrelated and correlated noise inputs, and we may interpret the aperiodic exponent in terms of these different types of noise; an aperiodic exponent near $\beta =-4$ indicates correlated noise inputs dominate the high frequency observable dynamics, while $\beta =-2$ indicates uncorrelated noise dominates the high frequency observable dynamics. These general results are independent of a specific biophysical mechanism, and a biophysical interpretation of the relationship between different types of noise and the aperiodic exponent depends on the specific model choice.

Some observations report aperiodic exponents greater than −2 (i.e., $\mathrm{\beta }>-2$), beyond the range of aperiodic exponents derived here. Experimental factors, such as measurement noise, which flattens the spectrum and shifts $\beta$ towards 0, might contribute to these observations. In addition, analysis factors may impact reported results. For example, the frequency range considered varies widely across studies (see discussion of fitting ranges in (Gerster et al., 2022)). In general, larger aperiodic exponents (more negative $\mathrm{\beta }$) occur at higher frequencies (Bédard et al., 2006; Colombo et al., 2019; Ibarra Chaoul & Siegel, 2021; Racz et al., 2021), although not always (Chaudhuri et al., 2018). In the framework considered here, the result $-4\le \mathrm{\beta }\le -2$ holds in the high frequency limit (when the measured frequency exceeds any natural frequency of the neural population) and near an equilibrium of the dynamical system. In lower frequency bands (e.g., below the natural frequency), our simplifying asymptotic arguments no longer apply, and the relationship between power and frequency will depend on the model parameters. In addition, we note that the aperiodic exponent reported in lower frequency bands is more difficult to interpret. To assess low frequency rhythms requires long durations of data, which increases the chance of nonstationarity. Artifacts (e.g., slow drifts) and analysis choices (e.g., whether to subtract the signal mean) also impact the low frequency power. Finally, while power-law features at low frequencies may reflect the same power-law features at high frequencies, these different phenomena unlikely reflect the same neural mechanisms. We also expect analyzing dynamics away from an equilibrium of the dynamical system will increase the aperiodic exponent. Away from an equilibrium, nonlinear terms in the model have a greater impact on the dynamics and resulting spectrum. These nonlinearities may increase power at high frequencies (e.g., (Kramer, Tort, et al., 2008)) and therefore increase the aperiodic exponent beyond the range derived for the linear dynamics near the equilibrium.

Here we considered additive Gaussian noise as the stochastic inputs to the dynamical models. However, alternative noise sources would provide more accurate approximations of biological systems. For example, for conductance based neural models, more biophysically plausible noise sources would include noise in the membrane current dynamics with variance proportional to the opening and closing rates of ion channels (Pu & Thomas, 2020, 2021). For models of neural population activity, stochastic perturbations may be scaled by the square root of the population size (Benayoun et al., 2010; Candia et al., 2021). For models of ecological population dynamics, noise may be scaled proportional to the birth and death rates (Barendregt & Thomas, 2023; Huynh et al., 2023; Strang et al., 2019). Understanding the impact of more biologically realistic noise (e.g., state-dependent noise) on the power-law behavior remains an important topic for future investigation.

Many well-supported observations of power-laws appear in neuroscience (e.g., avalanches of population voltage discharges (Beggs & Plenz, 2003), amplitudes of narrowband oscillations (Linkenkaer-Hansen et al., 2001)). Here, we consider one type of power-law: the 1/f-like neural field spectrum and a general, noise-driven dynamical system. Under this general model, the aperiodic exponent represents the impact of noise in the observable and latent dynamics, without requiring a sophisticated biological or dynamical mechanism. We propose that the range of aperiodic exponents $-4\le \beta \le -2$ observed in vivo represents the expected dynamics near an equilibrium in a nonlinear dynamical system driven by noise. The generality of the model is consistent with the universality of 1/f-like field spectra, reflecting a basic dynamical feature present in many different neural systems. However, this simplicity may also limit the computational utility of this mechanism and the role of the aperiodic exponent in measuring neural computations.

Materials And Methods

To estimate the aperiodic exponent $\beta$, we first compute the spectrum in the standard way. To a simulated voltage time series $\left(V_t\right)$ with sampling interval $\mathrm{\Delta }$ and duration $T$, we subtract the mean, apply a Hanning taper, compute the Fourier transform ($V_f$), and multiply by the complex conjugate $\left(V_f^{\mathrm{*}}\right)$:

We note that the square of the Fourier coefficients is essential for consistent interpretation of the aperiodic exponent across studies; omitting the square is a common mistake identified in previous work (see discussion in (Milotti, 2002)). For frequencies f, we fit a linear model to the logarithm base 10 of the spectrum (${\mathrm{l}\mathrm{o}\mathrm{g}}_{10}P$) with predictor logarithm base 10 of the frequency $\left({\mathrm{l}\mathrm{o}\mathrm{g}}_{10}f\right)$,

where $\beta$ is the estimate of the aperiodic exponent. Code to compute the spectrum and estimate the aperiodic exponent is available at https://github.com/Mark-Kramer/Aperiodic-Exponent-Model.