Parameterizing neural power spectra into periodic and aperiodic components

Abstract

Electrophysiological signals exhibit both periodic and aperiodic properties. Periodic oscillations have been linked to numerous physiological, cognitive, behavioral and disease states. Emerging evidence demonstrates that the aperiodic component has putative physiological interpretations and that it dynamically changes with age, task demands and cognitive states. Electrophysiological neural activity is typically analyzed using canonically defined frequency bands, without consideration of the aperiodic (1/f-like) component. We show that standard analytic approaches can conflate periodic parameters (center frequency, power, bandwidth) with aperiodic ones (offset, exponent), compromising physiological interpretations. To overcome these limitations, we introduce an algorithm to parameterize neural power spectra as a combination of an aperiodic component and putative periodic oscillatory peaks. This algorithm requires no a priori specification of frequency bands. We validate this algorithm on simulated data, and demonstrate how it can be used in applications ranging from analyzing age-related changes in working memory to large-scale data exploration and analysis.

Extended Data Figure 1 |. False oscillatory power changes and illusory oscillations.

(A) Here we took a real neural PSD (blue) and artificially introduced a change in the aperiodic exponent, similar to what is seen in healthy aging. This PSD was then inverted back to the time domain (right panels). The exponent change manifests as amplitude differences in the time domain. This affects apparent narrowband power when an a priori filter is applied. This is despite the fact that the true oscillatory power relative to the aperiodic component is unaffected. (B) Even when no oscillation is present, such as the case with the white and pink (1/f) noise examples here (blue and green, respectively), narrowband filtering gives rise to illusory oscillations where no periodic feature exists in the actual signal, by definition.

Extended Data Figure 2 |. Algorithm performance on simulated data across a broader frequency range.

(A-C) Power spectra were simulated across the frequency range (1–100 Hz), with two peaks, one in a low range, and one in a high range (see Methods), across five distinct noise levels (1000 spectra per noise level). (A) Example power spectra with simulation parameters as aperiodic [offset, knee, exponent] and periodic [center frequency, power, bandwidth] (B) Absolute error of algorithmically identified peak center frequency, separated for the low (3–35 Hz) and high range (50–90 Hz) peaks. (C) Absolute error of algorithmically identified aperiodic parameters, offset, knee, and exponent. All violin plots show full distributions, where small white dots represent median values and small box plots show median, first and third quartiles, and ranges. Note that the error axis is log-scaled in B,C.

Extended Data Figure 3 |. Algorithm performance on simulated data that violate model assumptions.

(A-C) Power spectra were simulated across a broader frequency range (1–100 Hz), with two peaks, one in a low range, and one in a high range (see Methods), across five distinct knees values (1000 spectra per knee value), with a fixed noise level (0.01). Power spectra were parameterized in the ‘fixed’ aperiodic mode (without a knee) to evaluate how sensitive performance is to aperiodic mode. (A) Example power spectra with simulation parameters as aperiodic [offset, knee, exponent] and periodic [center frequency, power, bandwidth], showing spectra with knee values of 0 and 150, both fit in the ‘fixed’ aperiodic mode. (B) Absolute error of algorithmically identified aperiodic exponent, across spectra with different knee values. Notably, exponent reconstruction is high when spectra with knees are fit without a knee parameter. (C) The number of peaks fit by the model, across knee values. Note that all spectra in this group have two peaks, indicating here that the presence of knee’s in ‘fixed’ mode leads to overfitting peaks. (D-F) A distinct set of simulations were created in which power spectra were created with asymmetric or skewed peaks (see Methods), across five distinct skew levels (1000 spectra per skew level). (D) Example simulated spectra, showing two different skew levels. (E) Absolute error of algorithmically identified peak center frequency, across peak skewness values. (F) The number of peaks fit by the model, across peak skewness. Note that all spectra in this set have one peak. (G-I) A distinct set of simulations, in which time series were generated with asymmetric oscillations in the time domain, from which power spectra were calculated (see Methods), across five distinct levels of oscillation asymmetry (1000 spectra per asymmetry value). (G) Example simulation of an asymmetric oscillation, simulated in the time domain, and the associated power spectrum. Note that the power spectrum displays harmonic peaks. (H) Absolute error of algorithmically identified peak center frequency, across oscillation asymmetry values. (I) The number of peaks fit by the model, compared across oscillation asymmetry values. Note that these simulations all contained one oscillation in the time domain. All violin plots show full distributions, where small white dots represent median values and small box plots show median, first and third quartiles, and ranges. Note that the error axis is log-scaled in B,E,H.

Extended Data Figure 4 |. Oscillation band occurrences and relative powers in the MEG dataset.

(A) The proportion of participants for whom an oscillation peak was fit, at each vertex, per band. (B) The group level relative power, per band. For each participant, the oscillation power within the band was normalized between 0 and 1, and then averaged across all participants, such that a maximal relative power of 1 would indicate that all participants have the same location of maximal band-specific power. Note that alpha and beta have maximal values approaching 1, reflecting a high level of consistency in location of maximal power, whereas in theta the values are lower, reflecting more variability. The “oscillation score” metric, as presented in Fig. 7, is the result of multiplying the occurrence probability map with the power maps.

Extended Data Figure 5 |. Methods comparison of measures of the aperiodic exponent.

(A-C) Comparisons of spectral parameterization to a linear fit (in log-log space), as used in BOSC, on simulated power spectra (1000 per comparison). (A) Comparison of a linear fit and spectral parameterization on a low frequency range (2–40 Hz) with one peak. (B) Comparison across the same range with multiple (3) peaks (C) Comparison across a broader frequency range (1–150 Hz) with two peaks and an aperiodic knee. In all cases (A-C), spectral parameterization outperforms the linear fit. (D-F) Comparisons of spectral parameterization to IRASA. Groups of simulations mirror those used in (A-C). Note that for these simulations, the data were simulated as time series (1000 per comparison) (see Methods). IRASA and spectral parameterization are comparable for the one peak cases (D), but spectral parameterization is significantly better in the other cases (E,F). Note that as IRASA has both greater absolute error and a systematic estimation bias (see Supplemental Modeling Note). (G-I) Example of IRASA and spectral parameterization applied to real data. (G) The IRASA-decomposed aperiodic component, using default settings, in orange, is compared to the original spectrum, in blue. There are still visible non-aperiodic peaks, meaning IRASA did not fully separate out the periodic and aperiodic components. (H) The IRASA-decomposed aperiodic component, with increased resampling. This helps remove the peaks, but also increasingly distorts the aperiodic component, especially at the higher frequencies, due to its multi-fractal properties (the presence of a knee). (I) The isolated aperiodic component from spectral parameterization (computed as the peak-removed spectrum), showing parameterization can account for concomitant large peaks and knees, providing a better fit to the data. All violin plots show full distributions, where small white dots represent median values and small box plots show median, first and third quartiles, and ranges. Note that the error axis is log scaled in (A-F). * indicates a significant difference in between the distributions of errors between methods (paired samples t-tests). Discussion of these methods and results are reported in the Supplemental Modeling Note.

Figure 1 |. Overlapping nature of periodic and aperiodic spectral features.

(A) Example neural power spectrum with a strong alpha peak in the canonical frequency range (8–12 Hz, blue shaded region) and secondary beta peak (not marked). (B) Same as A, but with the alpha peak removed. (C-D) Apparent changes in a narrowband range (blue shaded region) can reflect several different physiological processes. Total power (green bars in the inset) reflects the total power in the range, and relative power (purple bars in the insets) reflect relative power of the peak, over and above the aperiodic component. (C) Measured changes, with a peak present, including: (i) oscillatory power reduction; (ii) oscillation center frequency shift; (iii) broadband power shift, or; (iv) aperiodic exponent change. In each simulated case, total measured narrowband power is similarly changed (inset, green bar), while only in the true power reduction case (i) has the 8–12 Hz oscillatory power relative to the aperiodic component actually changed (inset, purple bar). (D) Measured changes, with no peak present. This demonstrates how changes in the aperiodic component can be erroneously interpreted as changes in oscillation power when only focusing on a narrow band of interest.

Figure 2 |. Algorithm schematic on real data.

(A) The power spectral density (PSD) is first fit with an estimated aperiodic component (blue). (B) The estimated aperiodic portion of the signal is subtracted from the raw PSD, the residuals of which are assumed to be a mix of periodic oscillatory peaks and noise. (C) The maximum (peak) of the residuals is found (orange). If this peak is above the noise threshold (dashed red line), calculated from the standard deviation of the residuals, then a Gaussian (solid green line) is fit around this peak based on the peak’s frequency, power, and estimated bandwidth (see Methods). The fitted Gaussian is then subtracted, and the process is iterated until the next identified point falls below a noise threshold or the maximum number of peaks is reached. The peak-finding at this step is only used for seeding the multi-Gaussian in D, and, as such, the output in D can be different from the peaks detected at this step. (D) Having identified the number of putative oscillations, based on the number of peaks above the noise threshold, multi-Gaussian fitting is then performed on the aperiodic-adjusted signal from B to account for the joint power contributed by all the putative oscillations, together. In this example, two Gaussians are fit with slightly shifted peaks (orange dots) from the peaks identified in C. (E) This multi-Gaussian model is then subtracted from the original PSD from A. (F) A new fit for the aperiodic component is estimated—one that is less corrupted by the large oscillations present in the original PSD (blue). (G) This re-fit aperiodic component is combined with the multi-Gaussian model to give the final fit. (H) The final fit (red)—here parameterized as an aperiodic component and two Gaussians (putative oscillations)—captures >99% of the variance of the original PSD. In this example, the extracted parameters for the aperiodic component are: broadband offset = −21.4 au; exponent = 1.12 au/Hz. Two Gaussians were found, with the parameters: (1) frequency = 10.0 Hz, power = 0.69 au, bandwidth = 3.18 Hz; (2) frequency = 16.3 Hz, power = 0.14 au, bandwidth = 7.03 Hz.

Figure 3 |. Algorithm performance on simulated data.

(A-C) Power spectra were simulated with one peak (see Methods), at five distinct noise levels (1000 spectra per noise level). (A) Example spectra with simulation parameters are shown (black), as aperiodic [offset, exponent] and periodic [center frequency, power, bandwidth]. Spectral fits (red), for the one-peak simulations in a low- and high-noise scenario. Simulation parameters for plotted example spectra are noted. (B) Median absolute error (MAE) of the algorithmically identified aperiodic offset and exponent, across noise levels, as compared to ground truth. (C) MAE of the algorithmically identified peak parameters—center frequency, power, and bandwidth—across noise levels. In all cases, MAE increases monotonically with noise, but remains low. (D-F) A distinct set of power spectra were simulated to have different numbers of peaks (0–4, 1000 spectra per number of peaks) at a fixed noise level (0.01; see Methods). (D) Example simulated spectra, with fits, for the multi-peak simulations. Conventions as in A. (E) Absolute model fit error for simulated spectra, across number of simulated peaks. (F) The number of peaks present in simulated spectra compared to the number of fitted peaks. All violin plots show full distributions, where small white dots represent median values and small box plots show median, first and third quartiles, and ranges. The algorithm imposes a 6.0 Hz maximum bandwidth limit in its fit, giving rise to the truncated errors for bandwidth in C. Note that the error axis is log-scaled in B,C,E.

Figure 4 |. Algorithm performance compared to human raters on real EEG and LFP data.

(A,B) Examples of two different EEG spectra labeled by expert human raters, highlighting cases of strong (A) and weak (B) consensus amongst raters. The black line is the PSD of real data against which center frequency estimates were made. The red line is the algorithm fit; the red stars are the center frequencies identified by the algorithm. The dots are each individual expert’s center frequency rating(s). Note that even when human consensus was low, with many identifying no peaks, as in B, the algorithm still provides an accurate fit (in terms of the R fit and error). Nevertheless, the identified center frequencies in B would all be marked as false positives for the algorithm as compared to human majority rule, penalizing the algorithm. (C) Human raters show a strong precision/recall tradeoff, with some variability amongst raters. Inset is Spearman correlation between precision and recall. (D) Despite the penalty against the algorithm for potential overfitting, as in B, it performs comparably to the human majority rule. n.s.: algorithm not significantly different from human raters.

Figure 5 |. Age-related shifts in spectral EEG parameters during resting state.

(A) Visualization of individualized oscillations as parameterized by the algorithm, selected as highest power oscillation in the alpha (7–14 Hz) range from visual cortical EEG channel Oz for each participant. There are clear differences in oscillatory properties between age groups that are quantified in C. (B) A comparison of alpha captured by a canonical 10 ± 2 Hz band, as compared to the average deviation of the center frequency of the parameterized alpha, for the younger group (left, blue) and older group (right, green). In this comparison, the canonical band approach captures 84% of the parameterized alpha in the younger adult group, and only 71% in the older adult group, as quantified in the middle panel, reflecting a significant different (p=0.031). Red represents the alpha power missed by canonical analysis, which disproportionately reflects more missed power in the older group. (C) Comparison of parameterized alpha center frequency (p=0.036), aperiodic-adjusted power (p=0.018), and bandwidth (p=0.632), split by age group. (D) Comparison of aperiodic components at channel Cz, per group. For this visualization, the aperiodic offset and exponent, per participant, were used to reconstruct an “aperiodic only” spectrum (removing the putative oscillations). Red shaded regions reflect areas where there are significant power differences at each frequency between groups (p < 0.05 uncorrected t-tests). In this comparison, significant group differences in both aperiodic offset (p<0.0001) and exponent (p=0.0001) (E) drive group-wise differences that otherwise appear to be band-specific in both low (< ~10 Hz) and high (> ~40 Hz) frequencies, when analyzed in a more traditional manner. Bars in B, C, E represent mean values while stars indicate statistically significant difference, (two-sided paired samples t-test, uncorrected) at p < 0.05. n.s. not significant.

Figure 6 |. Event-related spectral parameterization of working memory in aging.

(A) Contralateral electrodes (filled blue dots on the electrode localization map) were analyzed in a working memory task, with spectral fits to delay period activity per channel, per trial, as well as to the pre-trial baseline period (see Methods). Task-related measures of each spectral parameter were computed by subtracting the baseline parameters from the delay period parameters. (B) Parameters were collapsed across channels, to provide a measure per trial, and collapsed across trials, to provide a measure per working memory load (condition). (C) For this analysis, condition average spectral parameters were used to predict behavioral performance, measured as d’, per condition. (D) The average evoked difference in spectral parameters, between baseline and delay periods, for each group, presented as spectra reconstructed from the spectral fits, including aperiodic and oscillatory alpha parameters. The inset tables present the changes in each parameter, shaded if significant (one sample t-test, p < 0.05; green: positive weight; red: negative weight). For further details of the statistical comparisons, see Supplementary Table 1. (E) Parameters for regression models predicting behavioral performance for each group. All behavioral models use the evoked spectral parameters from A. Note that all models also include an intercept term, and a covariate for load. Models were fit using ordinary least squares, and the model r-squared and results of F-test for model significance are reported (see Methods).

Figure 7 |. Large-scale analysis of MEG resting state data uncovers cortical spectral features.

Spectral parameterization was applied to a large MEG dataset (n=80 participants, 600,080 spectra). (A) Oscillation topographies reflecting the oscillation score: the probability of observing an oscillation in the particular frequency band, weighted by relative band power, after adjusting for the aperiodic component (see Methods). These topographies quantify the known qualitative spatial distribution for canonical oscillation bands theta (3–7 Hz), alpha (7–14 Hz), and beta (15–30 Hz). (Right) The topography of resting state aperiodic exponent fit values across the cortex. Group exponent—calculated as the average (mean) exponent value, per vertex, across all participants—shows that the aperiodic exponent is lower (flatter) for more anterior cortical regions. (B) The distribution of all center frequencies, extracted across all participants (n = 80), roughly captures canonical bands, though there is substantial heterogeneity. Peaks in the probability distribution are labeled for approximate canonical bands. (C) The distribution of aperiodic exponent values, fit across all vertices and all participants. (D) Correlations between the oscillation topographies and the exponent topography (as plotted in A) show that theta is spatially anti-correlated with the other parameters.