Introduction

Analysis of macroscopic electromagnetic brain activity (e.g., by EEG and MEG) has long been focusing on the investigation of ‘rhythmic’ neural oscillations. In the frequency domain, neural oscillations appear as distinct spectral peaks, also referred to as the periodic part of the spectrum (Buzsáki & Draguhn, 2004; Engel et al., 2001; Schnitzler & Gross, 2005). The full spectrum, however, also consists of a continuous component whose analysis has, so far, seen less attention. This aperiodic or ‘arrhythmic’ part of the spectrum (Freeman & Zhai, 2009; Miller et al., 2009) has been related to the integration of underlying synaptic currents (Buzsáki et al., 2012). Since the time series of the aperiodic part is typically self-similar across many temporal scales, it is also referred to as “fractal” or “scale-free” activity. The power spectral density (PSD) of the aperiodic component follows a power law \(P\propto 1/{f}^{\beta }\) (Miller et al., 2009) and is sometimes called 1/f activity for that reason. In this text, we will refer to the scaling exponent \(\beta\) in this equation as 1/f exponent.

The investigation of neural oscillations has received much attention in electrophysiological studies (Buzsáki & Draguhn, 2004; Singer, 1999; Ward, 2003). However, the standard analysis of assessing periodic power through bandpass filtering is problematic because the pass-band comprises both periodic and aperiodic activity. If the power of aperiodic activity changes between two conditions, analyzing neural oscillations in bandpass filtered signals would hence be confounded by these changes in the aperiodic part of the spectra. For that reason, estimating the 1/f component before determining the power of periodic activity has recently been suggested (Donoghue et al., 2021; Wen & Liu, 2016).

Besides investigating neural oscillations, the investigation of the aperiodic component has recently gained considerable interest (He, 2014; Kello et al., 2010). For example, the 1/f exponent was shown to change with task (Ouyang et al., 2020; Podvalny et al., 2015; Waschke et al., 2021), age (Bódizs et al., 2021; Dave et al., 2018; Waschke et al., 2017; Cellier et al., 2021; He et al., 2019; Schaworonkow & Voytek, 2021; Voytek et al., 2015), psychoactive drug administration (Muthukumaraswamy & Liley, 2018; Stock et al., 2020; Timmermann et al., 2019), and disease (Molina et al., 2020; Robertson et al., 2019; Veerakumar et al., 2019; van Heumen et al., 2021; Ostlund et al., 2021, Karalunas et al. 2022) and it decreases with cortical depth (Halgren et al., 2021). Furthermore, using computational modeling, (Gao et al., 2017) suggested the 1/f exponent \(\beta\) as an estimator of excitation–inhibition (E–I) balance. Many studies comparing conscious states to unconscious states seem to support this concept. Conscious states are associated with increased excitation while unconscious states, such as NREM sleep and anesthesia, are associated with pronounced inhibitory processes. Experiments comparing the 1/f exponent in these two conditions did indeed find larger values for NREM sleep (Lendner et al., 2020; Miskovic et al., 2019) and anesthesia (Colombo et al., 2019; Muthukumaraswamy & Liley, 2018; Waschke et al., 2021; Zhou et al., 2021) compared to awake states.

But how to best estimate the 1/f exponent? This will be the main question discussed in this study. One option is to simply fit a straight line using (robust) linear regression. (Gao et al., 2017) used this method in the frequency ranges apart from pronounced oscillatory peaks in electrocorticography (ECoG) data and identified distinct 1/f exponents during wakefulness and anesthesia. However, in the presence of periodic components, this method is error-prone because larger periodic peaks will bias the linear regression fit.

Irregular-resampling autospectral analysis (IRASA) (Wen & Liu, 2016) aims to separate periodic components from the aperiodic part of the spectrum. Due to their fractal nature, aperiodic time series remain robust against resampling, whereas periodic components are strongly affected by this procedure. IRASA takes advantage of this dichotomy and ‘removes’ the periodic parts from a spectrum. The–ideally–pure aperiodic part of the spectrum obtained with this method can then be used for fitting the 1/f exponent.

Another method, ‘fitting oscillations & one over f’ (FOOOF) (Donoghue et al., 2020), aims at modeling the periodic components: It iteratively applies Gaussian fits to all periodic components and hereby obtains a model of the periodic part. This model of periodic activity is subtracted from the spectrum to obtain an–ideally–pure aperiodic component which can be used for fitting \(\beta\). In addition, the periodic model allows for analyzing the periodic components (e.g., regarding center frequencies, bandwidths, and power) without the bias from aperiodic activity.

This article highlights and discusses the general challenges of estimating 1/f exponents. In addition, we also discuss method-specific challenges of FOOOF and IRASA, the most commonly used algorithms for that purpose.

In the Methods section, we will introduce our simulations, our datasets, and both algorithms FOOOF and IRASA. We will analyze challenges by the example of FOOOF in section FOOOF and by the example of IRASA in the section IRASA. To aim for broad applicability of our assessment, we will apply these methods to simulations with known ground truth in addition to various electrophysiological signals obtained from empirical EEG, gradiometer MEG, magnetometer MEG, source-reconstructed voxel activity from MEG, and subthalamic nucleus-(STN-)LFP data acquired by three independent research groups. We will discuss these challenges in the section Discussion, and we will provide some guidance on how to use these methods in the Conclusion section.

Methods

Simulations

We simulate aperiodic 1/f activity by constructing a Fourier power spectrum following a preset \(1/{f}^{\beta }\) power-law. The corresponding phases of the Fourier spectrum are distributed uniformly randomly. To add oscillations, we add Gaussian-shaped peaks to the Fourier power spectrum with amplitudes \(A\) and a spectral extent given by center frequencies \({f}_{\text{center}}\) and variances \(\sigma_f^2\). The corresponding time series, consisting of both ‘neural’ oscillations and aperiodic activity, is then obtained by applying the inverse fast Fourier transform. The simulated time series either have a duration of 180 s at a sampling rate of \({f}_{\text{sample}}=2400\) Hz or are matched to the empirical data to which a simulation might be compared. If noted in the text, Gaussian white noise may be added to the time series afterward. Since most algorithms to generate 1/f activity lead to identical power spectra, the specific choice of the algorithm has no impact on the present analysis.

Empirical Data

We compare the results from our simulations to three empirical datasets.

Dataset 1

Dataset 1 was re-analyzed from (Litvak et al., 2010, 2011) and contains MEG and LFP data of 14 Parkinsonian patients after bilateral implantation of subthalamic nucleus (STN) stimulation electrodes (Medtronic, Minneapolis, MN, USA with four platinum-iridium cylindrical surfaces of diameter 1.27 mm, length 1.5 mm, and center-to-center separation 2 mm) for deep brain stimulation (DBS). The joint ethics committee of the National Hospital of Neurology and Neurosurgery and the University College London Institute of Neurology approved the study, and all patients gave their written informed consent. The patients were recorded three days after surgery when the electrode leads were still externalized. The recordings were obtained during a Parkinsonian state OFF medication (after overnight withdrawal) and an ON medication state. MEG (275 channels, CTF/VSM MedTech, Vancouver, Canada) and DBS-LFP were recorded simultaneously during three minutes of resting-state at a sampling rate of \({f}_{\text{sample}}=2400\) Hz. The LFP recordings were referenced to the right mastoid during recording and later re-referenced to a bipolar montage between adjacent electrode contacts. This results in 3 bipolar LFP channels per hemisphere. All data were bandpass filtered in hardware between 1–600 Hz. MEG source reconstruction was performed with varying regularization by Linearly Constrained Minimum Variance beamformer (Van Veen et al., 1997). Aside from the six LFP channels, the dataset contains three MEG channels per patient from voxels located in the supplementary motor area (SMA), left primary motor cortex (M1), and right M1. In this study, we draw examples from voxel data located in the supplementary motor area (SMA) of patients 5 and 6 and bipolarly recorded LFPs from the STN of patients 9 and 10. Details regarding the data recording, processing, and inverse modeling can be obtained from the original publications of this dataset (Litvak et al., 2010, 2011). In this study, we further process this dataset by applying a notch filter at 100, 150, …, 600 Hz power line noise for visualization purposes (multi-taper estimation of sinusoidal components “spectrum_fit” of MNE python (Gramfort et al., 2013). Note, that a notch filter should not be applied to frequency ranges used for FOOOF fitting. Therefore, we exclude 50 Hz from the notch filter (see SI Fig. 1).

Fig. 1
figure 1

Algorithms for 1/f estimation. IRASA: a) Simulated time series. b) PSDs of resampled time series on the y-axis and frequencies on the x-axis. In this figure, the time series is upsampled by the resampling factors \({h}_{i}\) of the \({h}_{\text{set}}=\left\{1.3, 1.6, 2\right\}\) and downsampled by \(1/{h}_{i}\). c) The geometric mean of all resampling pairs (\({h}_{i}\), \(1/{h}_{i}\)) is calculated. d) The aperiodic component (orange) is the median of the geometric means. A final fit (dashed-blue) estimates the y-intercept and the 1/f exponent \(\beta\). FOOOF: e) A PSD is calculated from the time series. f) FOOOF applies an initial linear fit (dashed-blue) to the PSD in log–log space and g) subtracts the obtained linear trend from the spectrum. h) A Gaussian model (dotted-green) is fitted to the largest peak exceeding the thresholds (dashed-grey) and removes it. The relative threshold is recalculated from the peak-removed flattened spectrum (pink). The procedure is repeated until no peak exceeds the relative threshold. d) Subtraction of all Gaussian models from the original PSD yields the aperiodic component, which is then finally re-fit

Full size image

Dataset 2

Dataset 2 contains EEG data from a 12-year-old boy with absence epilepsy recorded at the Department of Epileptology at the University of Bonn (Gerster et al., 2020). The university’s ethics committee approved the study, and a parent gave written informed consent that the clinical data might be used and published for research purposes. EEG data were acquired at a sampling rate of \({f}_{\text{sample}}=256\) Hz (16-bit A/D conversion) within a bandwidth of 0.3–70 Hz from 19 electrodes in a bipolar montage. The locations and nomenclature of these electrodes are standardized by the American Electroencephalographic Society (Sharbrough, 1991). The EEG was recorded over several hours and contains 5 absence seizures. In this study, we present 40 s of the bipolar EEG channel “F3 − C3” during one absence seizure.

Dataset 3

Dataset 3 contains MEG recorded with gradiometers and magnetometers and LFP data from a Parkinsonian patient recorded at the Universitätsklinikum Düsseldorf. The data were acquired using a whole-head MEG system with 306 channels (Elekta Vectorview, Elekta Neuromag, Finland), and segmented “1–3–3–1” electrode DBS-LFP (Abbott St. Jude Medical model 6172, contact height: 1.5 mm with 0.5 mm vertical spacing) during the ON- and OFF-medication state (after overnight withdrawal). The patient was recorded 1 day after surgery when the electrode leads were still externalized. The resting-state was recorded for 10 min at a sample rate of \({f}_{\text{sample}}=2400\) Hz. The LFP recordings were referenced to the right mastoid during recording and later re-referenced to a bipolar montage between adjacent electrodes. The data were offline band-pass filtered between 0.3 Hz and 600 Hz and notch-filtered at 50, 100, …, 600 Hz power line noise (with a second-order IIR filter of bandwidth 1 Hz). Note that notch filtering in the fitting range at 50 Hz is unproblematic with using IRASA. The patient gave written consent to participate in the study, which was approved by the Ethics committee of the Universitätsklinikum Düsseldorf. In this study, we analyze data from one gradiometer channel, one magnetometer channel, and one LFP channel of the subject.

Power Spectral Densities (PSDs)

We calculate the PSDs from the simulated and recorded time series using the Welch algorithm. We use a segment length of 1 s which corresponds to a frequency resolution of 1 Hz, and the Hann-windowed segments overlap by 50%. Please note that other segment lengths can be used depending on the properties of the data. However, for FOOOF, the PSDs should be sufficiently smooth to avoid fitting noise peaks. IRASA receives time series as input and calculates the PSDs internally. For IRASA, the PSD resolution should be sufficiently high. We, therefore, use a segment length of 4 s (corresponding to a resolution of 0.25 Hz), Hann windows, and 50% overlap.

Irregular-Resampling Auto-Spectral Analysis (IRASA)

Irregular-resampling auto-spectral analysis (IRASA) aims at separating periodic components from the aperiodic part of the signal (Wen & Liu, 2016). In contrast to FOOOF, the algorithm requires time series as input (Fig. 1a) and does not explicitly model the signals’ spectra. The input time series is upsampled by a set of predefined resampling factors \({h}_{i}{\in h}_{\text{set}}\). By default, \({h}_{\text{set}}\) ranges from 1.1 to 1.9 with increments of 0.05, yielding 17 resampling factors \({h}_{\text{set}}=\left\{1.1, 1.15, ..., 1.9\right\}\). In addition, the time series is downsampled by all inverse resampling factors \(1/{h}_{i}\), with \({h}_{i}{\in h}_{\text{set}}\). For each of the 17 pairs of up- und downsampled spectra (Fig. 1b), the geometric mean of the PSD is calculated (Fig. 1c). For illustration purposes in Fig. 1, we use a very small \({h}_{\text{set}}=\left\{1.3, 1.6, 2\right\}\). Finally, the median is calculated from all 17 geometric means, yielding the aperiodic component (Fig. 1d). The compound oscillatory part of the spectrum is obtained by subtracting the aperiodic component from the original PSD. After applying IRASA, the slope \(\beta\) can be obtained by fitting the aperiodic component in double logarithmic space in the predefined fitting range.

As parameters, IRASA requires the fitting range, the resampling factors \({h}_{\text{set}}\), and the segment length for the PSD calculation. In this study, we vary the fitting range and the \({h}_{\text{set}}\) but keep the segment length at 4 s. IRASA’s Python implementation used for this article was adapted from the YASA toolbox (Vallat, 2019) and is published along with the complete code for this study on GitHub at https://github.com/moritz-gerster/oscillation_and_1-f_separation.

Fitting-Oscillations-&-One-Over-F (FOOOF)

FOOOF was introduced to parameterize neural power spectra as a combination of an aperiodic component and peaks representing oscillatory processes (Donoghue et al., 2020). The Python-based toolbox works as outlined in Fig. 1. First, the PSD of the time-series of interest (Fig. 1a) is calculated and input into the algorithm, Fig. 1e. Next, FOOOF calculates an initial robust linear fit of the spectrum in double logarithmic space, Fig. 1f, and subtracts the result from the spectrum, Fig. 1g. In this flattened spectrum, a relative threshold is calculated based on the standard deviation (SD) of the spectrum, Fig. 1h. The relative threshold is set to two times the SD, by default \(\text{thresh}_{\text{rel}}=2 \cdot \text{SD}\). Optionally, FOOOF also allows setting an additional absolute threshold for the peak heights, but it is set to 0 by default (\(\text{thresh}_{\text{abs}}=0\)). A Gaussian function is fitted to the largest peak of the flattened PSD exceeding both thresholds and then subtracted from the spectrum. Note, that this fit is not applied to negative peaks in the spectrum subceeding both thresholds. Therefore, spectral dips caused by notch filtering should be avoided. This procedure is iterated for the next largest peak after subtracting the previous peak until no peaks are exceeding the thresholds. The oscillatory components are finally obtained by fitting a multivariate Gaussian to all extracted peaks simultaneously. After the iterations, the initial fit is added back to the flattened peak-free PSD, which results in the aperiodic component of the PSD, Fig. 1d. Afterward, this aperiodic component is fitted again, leading to the final fit with y-intercept and slope \(\beta\) as parameters. The fitted Gaussian functions are parameterized by center frequency \({f}_{\text{center}}\)(“CF”), amplitude \({A}_{f}\) (“PW”), and bandwidth \(2\cdot \sigma^2_f\) (“BW”).

The algorithm can be used with the following parameters: “peak_width_limits” allows setting the minimum and maximum peak width limits of the Gaussian fits. The default is set to 0.5 Hz and 12 Hz, respectively. “max_n_peaks'' determines the maximum number of peak fitting iterations. The default is set to infinity. “min_peak_height” is the absolute threshold and corresponds to the smallest peaks that will be fitted in the units of the input data. The absolute threshold is set to 0 by default. “peak_threshold” is the relative threshold in SD multiples that peaks must exceed to be fitted and defaults to “peak_threshold” = 2 SD. The peak fitting stops when all remaining peaks are below either of these two thresholds. “aperiodic mode” allows for two modes of modeling: “fixed” and “knee” which allows modeling a bend (i.e., a “knee”) in the PSD of the aperiodic component. The default is “fixed.” Finally, FOOOF accepts a fitting range for which the algorithm performs the given steps. For a detailed description of this algorithm, we refer the reader to the methods section of the original publication (Donoghue et al., 2020).

In this study, we keep FOOOF at the default parameters if not stated otherwise and input PSDs with a spectral resolution of 1 Hz.

FOOOF

Challenge 1: The Spectral Plateau Disrupts the 1/f Power Law

The 1/f power law is sometimes called “scale-free” because log-transformed power is typically approaching a linear function across an extended range of log-transformed frequencies. For electrophysiological PSDs, however, this concept should be exercised carefully. For example, (He et al., 2010) measured different values for the slope β for frequency ranges 0.01–0.1 Hz and 1–100 Hz and found a small plateau in the range from 0.1 Hz to 1 Hz. In many studies, these low-frequency ranges < 0.1 Hz are eliminated by a hardware high-pass filter. However, this finding underlines the importance of selecting a representative frequency range to fit the 1/f slope.

In addition to the aforementioned low-frequency plateau, one regularly encounters a high-frequency spectral plateau (or flattening) in spectra of electrophysiological data. Such plateaus might be due to the presence of Gaussian noise which appears as a horizontal line with a slope \(\beta =0\) in double-logarithmic space and disrupts the 1/f power law. The origin of such white noise is often due to EMG artifacts and electronic noise of the recording system (Waterstraat et al., 2015b). It has been shown in EEG (Scheer et al., 2006; Waterstraat et al., 2015a) and MEG (Waterstraat et al., 2021) that extremely low-noise recording devices can shift this high-frequency plateau into the kHz range – leaving a wider unaffected frequency range for fitting the spectra. In conventional data, however, spectral plateaus are regularly present and will be discussed in this section because this can pose a severe challenge for estimating the aperiodic exponent: it shrinks the frequency range at which the 1/f exponent may be examined.

In Fig. 2a, a simulation of an aperiodic PSD with an exponent of \(\beta =2\) is shown. By adding white noise, a plateau can be observed starting at 100 Hz in the high-frequency range. Here, we define the onset of the plateau as the lowest frequency of a 50 Hz frequency interval with a vanishing exponent (i.e., flattening of the spectrum). Specifically, we apply FOOOF without periodic peak fitting to measure the slope from 1–50 Hz. We then gradually shift this interval by 1 Hz towards higher frequencies and fit the slope again. We repeat this procedure until the estimated slope reaches a value below \({\beta }_{\text{thresh}}=0.05\).

Fig. 2
figure 2

The spectral plateau disrupts the 1/f power law. The x-axis and the y-axis indicate frequency and PSD, respectively. a) Simulation of an aperiodic PSD (black) with a plateau starting at \(100\) Hz (grey). The spectrum starts to deviate from the ground truth (dashed line) after around \(10\) Hz. Applying FOOOF yields smaller 1/f exponent estimates with larger upper fitting range borders. b) A Parkinsonian LFP spectrum from the subthalamic nucleus shows large oscillations that hinder the plateau onset’s precise detection. c) Adding oscillations of various powers and widths on top of different aperiodic ground truths yields the same 1/f estimation of \(\beta \approx 0.77\) in FOOOF. The ground truths are \(\beta =1\) (blue), \(\beta =1.5\) (green), and \(\beta =2\) (orange)

Full size image

We apply FOOOF in the frequency intervals 1-10 Hz, 1-50 Hz, 1-100 Hz, and 1-200 Hz which yields estimated 1/f exponents of \(\beta =1.97\), \(\beta =1.64\), \(\beta =1.17\), and \(\beta =0.70\), respectively. The spectral plateau gradually biases the estimated 1/f exponents towards smaller values starting already at 10 Hz. This challenge might be to some extent alleviated if the analysis aims to study differences between groups or experimental conditions such that relative changes of the exponent are most important. However, when the precise onset of the spectral plateau varies across conditions, the upper fitting range border should be chosen as low as possible to minimize this unequal bias. Even if the plateaus seem to be similar across conditions, a lower upper fitting range border will increase the signal-to-noise ratio of the exponent estimates if we define the 1/f ground truth as signal and the impact of the plateau as (possibly Gaussian) noise. We, therefore, recommend choosing low upper fitting range borders and determining the precise onset of the plateau across conditions in order to estimate the exponents under equivalent conditions.

Yet, measuring the onset of the plateau can be difficult in practice if oscillatory peaks mask it. For example, the spectrum in Fig. 2b appears to have a spectral plateau onset at 120 Hz. The presence of the high-frequency oscillation peaking at 360 Hz, however, produces a positive slope between 160 and 320 Hz, potentially masking a continued 1/f trend of the spectrum. Accordingly, one cannot exclude the possibility that the actual onset of the plateau is at a higher frequency value (e.g., 200 Hz). On the other hand, the large oscillation ranging from 6 to 100 Hz, peaking in the beta range at 25 Hz, counteracts this effect: In theory, one also cannot exclude that the actual flattening occurs already at 20 Hz.

To demonstrate this effect, for Fig. 2c we simulate three power spectra with three different 1/f exponent ground truths of \(\beta =1\) (blue), \(\beta =1.5\) (green), and \(\beta =2\) (orange). Next, we add eight oscillations at \(3\) Hz, \(5\) Hz, \(10.5\) Hz, \(16\) Hz, \(23\) Hz, \(42\) Hz, \(50\) Hz, and \(360\) Hz and tune the oscillation amplitude and width parameters in all three examples to match the recording of Fig. 2b (purple). Finally, using FOOOF, we estimate the 1/f exponent in the frequency range from 1-95 Hz in the three simulated and the real PSD. Despite strongly diverging 1/f exponent ground truths, FOOOF estimates an 1/f exponent of about \(\beta \approx 0.77\) in all four cases. The diverging ground truths are apparent in Fig. 2c because the true aperiodic component (which is shown in light grey and is invisible to FOOOF) has a plateau onset at high frequencies in the blue curve, at intermediate frequencies in the green, and at low frequencies in the orange curve. However, neither for FOOOF nor for the experimental observer, it is possible to know which of these three scenarios best reflects the real spectrum in Fig. 2b. Therefore it is difficult to determine at which frequency scale an 1/f estimate might be valid.

Note that this challenge applies not only to cases where the goal is to estimate the 1/f exponent but also when such an estimate is used in order to remove the aperiodic component from the spectrum. While fitting a “shoulder” allows for modeling such a plateau, the strongly varying “shoulder” onsets in Fig. 2c cannot be captured, given that the three power spectra share the same appearance. Specifically, the oscillation power estimates based on FOOOF would be almost identical in all three spectra. However, in the simulation, the oscillation power increases considerably from the blue to the green to the orange curve.

Recommendations

Scenario A - the power spectra have a plateau onset at higher frequencies, and oscillations do not mask it:

Challenge 1 does not Apply

Scenario B - the power spectra have a plateau onset at lower frequencies, and the onset is easily discernible (because no or just small oscillations are present):

Determine the precise plateau onsets across conditions. Choose the upper fitting range border as low as possible to increase SNR.

Scenario C - the power spectra potentially have a plateau onset at lower frequencies, but oscillations mask the exact onset:

The upper fitting range border must be lower than the onset of the masking oscillation. If the remaining frequency range is too small (as in Fig. 2b), aperiodic fitting should be avoided.

Challenge 2: Avoiding Oscillations Crossing Fitting Range Borders

When choosing the fitting range to model the aperiodic component, oscillations crossing the fitting range borders must be avoided for all investigated power spectra. FOOOF assumes all oscillation peaks lying within the fitting range because it does not fit partial Gaussian peaks. Consequently, the estimated 1/f exponent error becomes large if the lower or upper fitting range border overlaps with a spectral peak.

In the upper panel of Fig. 3a, we simulate a PSD with a slope of \(\beta =2\) and oscillation peaks at \(5\) Hz, \(15\) Hz, and \(35\) Hz (black graph). We fix the upper fitting range border at \(100\) Hz and measure the 1/f exponent for all lower fitting ranges from 1–100 Hz up to 80–100 Hz. The lower panel in Fig. 3a indicates the absolute error of the estimated slope as a function of the lower fitting range border (red). The error is the absolute deviation from the ground truth \(\left|{\beta }_{\text{truth}}-{\beta }_{\text{FOOOF}}\right|\). Note that the error function resembles the oscillatory peaks, with the greatest errors occurring approximately at the peak center frequencies. The FOOOF parameters are kept at the default setting.

Fig. 3
figure 3

Oscillations must not cross fitting range borders. a) Upper panel: PSD of a simulated spectrum with β = 2 and oscillations at 5 Hz, 15 Hz, and 35 Hz (black). The x-axis and the y-axis indicate frequency and PSD, respectively. Lower panel: The exponent \(\beta\) is measured using FOOOF for all 80 frequency ranges from 1–100 Hz to 80–100 Hz (red). The x-axis indicates the lower fitting range border, while the y-axis shows the absolute deviation from the ground truth. b) Various frequency ranges commonly used for E–I estimation are applied to an STN-LFP PSD of a Parkinsonian patient (purple). Since many of the chosen ranges overlap with spectral peaks, the estimated exponents \(\beta\) are strongly differing. FOOOF parameters: max_n_peaks = 0 (for 30–45 Hz); max_n_peaks = 1 (for 40–60 Hz); peak_width_limits = (1, 100) (for 1–45 Hz and 1–95 Hz). c) The simulated PSD in the middle panel (green) was tuned to match the empirical PSD in b) (purple). FOOOF estimates a similar aperiodic exponent for the simulated and the real spectrum (β = 0.61). When decreasing the power of the 2 Hz delta oscillation (blue), the estimated aperiodic exponent decreases (β = 0.50) despite a constant exponent for the simulated spectrum. When increasing the power of the delta oscillation (orange), the estimated aperiodic exponent increases (β = 0.72)

Full size image

If the peaks do not completely lie within the fitting range, very error-prone fits are obtained, as shown in another exemplary STN-LFP recording from a Parkinsonian patient (purple) in Fig. 3b. The fit from 30–45 Hz (turquoise), a frequency range commonly used for estimation of E–I balance, measures the slope of the beta-to-gamma peak, not of the aperiodic component. For this fit, we set the FOOOF parameter for the maximum number of allowed peaks (max_n_peaks) to 0 since this frequency range is usually chosen to avoid oscillations altogether. Therefore, FOOOF just fits a straight line without peak modeling. The 40–60 Hz range (green) lies on top of the beta-to-gamma peak, too. Here, we set (max_n_peaks) to 1 to account for the power line noise. Further, the 1–45 Hz range (orange) is inappropriate because its upper fitting range border at 45 Hz lies in the middle of the gamma peak. While these are obvious examples of ill-chosen fitting ranges, in practice more subtle (but similar) errors might occur. Therefore the presence of oscillations at fitting range borders must be carefully checked for every single PSD of interest.

The 1–95 Hz range (purple dotted) seems to be the only acceptable range for this spectrum: The upper fitting range border extends beyond the beta-to-gamma peak but ends before the onset of the spectral plateau. The estimated exponent has a value of \({\beta }_{\text{FOOOF}}=0.61\). For these two frequency ranges (1–45 Hz and 1–95 Hz), we increased the peak width limits from 0.5–12 Hz (default) to 1–100 Hz to account for the chosen spectral resolution (1 Hz) and enable modeling of the very broad (> 12 Hz) beta-to-gamma peak. The corresponding FOOOF fits of Fig. 3b are shown in SI Fig. 2.

While the 1–95 Hz range seems best, it appears almost impossible to avoid low-frequency oscillations crossing the lower fitting range border. If some delta oscillations are present, they lead to a steepening of the spectrum which impacts the estimation of the 1/f exponent. We visualize this effect by reproducing the empirical LFP spectrum in three simulations in Fig. 3c. We set the oscillation frequencies to 2 Hz, 12 Hz, 18 Hz, 27 Hz, 50 Hz (gamma), 50 Hz (power line), and 360 Hz. In the panels in Fig. 3c from left to right, we only vary the delta power at 2 Hz while keeping the aperiodic component and all other oscillations’ amplitudes and widths fixed. Since the delta oscillation has a bandwidth crossing the lower fitting range border of 1 Hz, FOOOF-estimates of the 1/f exponent diverge strongly between the three scenarios (same FOOOF parameters as for the 1–95 Hz range in Fig. 3b). While the aperiodic (white noise-free) ground truth remains unchanged at \(\beta =1.5\) for all three simulations, FOOOF estimates an \(18\%\) lower 1/f exponent (blue, \(\beta =0.50\)) if the delta oscillation from the middle panel (green, \(\beta =0.61\)) is removed. On the other hand, it estimates an \(18\%\) larger 1/f exponent (orange, \(\beta =0.72\)) if we double the power of the delta oscillations. The power of the true delta oscillations in the purple curve is, of course, unknown.

Overall, fitting and removing delta oscillation peaks seems unfeasible since they rarely occur as a single distinguishable peak in the double logarithmic representation. Furthermore, FOOOF requires smooth input spectra to reduce the impact of noise which at the same time hinders fitting sharp peaks. Therefore, we recommend 1/f estimation for a higher lower-border of the fitting range to avoid the impact of these low-frequency oscillations. For high lower borders of the fitting range, oscillations can be better avoided, and if they are present, they likely have less impact on the estimation.

Estimating the power of low-frequency oscillations by removing the aperiodic part of the spectrum poses a special challenge in this regard. Many studies (Donoghue et al., 2020; El Boustani et al., 2009; Fransson et al., 2013; Freeman & Zhai, 2009; Miller et al., 2009; Wen & Liu, 2016) have conceptualized the aperiodic part of the spectrum as self-similar, or fractal, across a wide range of frequencies, such that the estimation of the 1/f exponent should be independent of the chosen fitting range. If this assumption does not hold, however, the aperiodic component must be fitted in the given frequency range of interest. If this frequency range of interest coincides with low-frequency oscillations, this challenge cannot be avoided.

The impact of (sub-)delta oscillations should therefore be kept in mind as a limitation. If one finds a difference of the 1/f exponent between groups of investigation, one should check whether the delta power of the FOOOF-fits varies across conditions. If delta power is similar across conditions but the slope varies, it seems likely that indeed the aperiodic component causes these differences in the estimated slopes and not a distortion by delta oscillations. If delta power does change across conditions (without a global offset of the PSD across all frequencies), the change of slopes could either be caused by a change of oscillatory delta activity (as shown in SI Fig. 1) or by a change in the aperiodic component itself, and these two possibilities cannot be differentiated with full certainty.

Recommendation

Scenario A - the 1/f exponent needs to be estimated:

Use a fitting range at higher frequencies (for example 40–60 Hz) to avoid distortion by low-frequency oscillations.\newline

Scenario B - the aperiodic component needs to be removed from the PSD:

If the assumption of self-similarity across a wide range of frequencies holds for the aperiodic part of the spectrum, both slope and intercept of its linear fit could theoretically be obtained from any frequency range. In reality, different exponents could be present in different frequency ranges. In that case, the exponent should be estimated in the broadband range starting at very low frequencies. For this lower fitting range border (starting often at around 1 Hz), the challenge cannot be avoided and should be kept in mind as a potential limitation of the results.

Challenge 3: FOOOF Cannot Characterize Oscillation Peaks that are not Clearly Distinguishable

As illustrated in Fig. 1, FOOOF models oscillations as Gaussian functions fitted to peaks in the flattened PSD. While this does not impose a severe challenge for clearly isolated peaks, the modeling becomes complicated when peaks overlap partially. If many different peaks overlap, the resulting PSD can be caused by various combinations of oscillations with different frequencies and powers that are impossible to disentangle on a single spectrum. Furthermore, whereas spectral leakage from oscillations at neighboring frequencies but same Fourier phase can add up in different combinations to yield similar power spectra, oscillations at different phases can also subtract power from other peaks.

In the right panel of Fig. 4a, we present a real PSD that might exemplify a spectrum containing many strongly overlapping oscillation peaks. The underlying time series was recorded from a subject with epilepsy during an absence seizure using a bipolar montage of EEG electrodes F3-C3. Pre-seizure activity is highlighted in turquoise, seizure activity in red, and post-seizure activity in yellow. Absence seizures are proposed to be related to cortico-thalamic E–I dysbalance (Onat et al., 2013) caused by reduced cortical inhibition (Tan et al., 2007), hyperexcitable somatosensory neurons (Karpova et al., 2005), GABAB receptor dysfunctions (Inaba et al., 2009; Merlo et al., 2007), changes in NMDA (D’Arcangelo et al., 2002; Pumain et al., 1992), or mGLU2/3 receptors (Ngomba et al., 2005). It would be interesting to complement these molecular rodent studies by non-invasive human electrophysiological recordings. Specifically, using the 1/f exponent as a biomarker of E–I balance before, during, and after the seizure might help to gain new insights into absence seizures. However, the non-sinusoidal 3 Hz spike-wave discharges might create many harmonic peaks throughout the spectrum. Applying FOOOF (default parameters) in a frequency range of 1–100 Hz yields estimated 1/f exponents of \({\beta }_{\text{pre}}=1.52\), \({\beta }_{\text{seiz}}=2.31\), and \({\beta }_{\text{post}}=1.52\). One could interpret this finding as an increase of the aperiodic 1/f exponent during the seizure, indicating (quite counterintuitively) stronger neural inhibition. However, even though FOOOF subtracts a substantial part of the harmonic peaks by modeling them as four broad peaks with center frequencies at 11 Hz, 22 Hz, 37 Hz, and 50 Hz (see SI Fig. 3), it is not clear whether it can correctly estimate the peak heights. A peak height is the power of an oscillation on top of the aperiodic component. However, there is no reference point for the aperiodic component from which the height could be measured in the scenario of many overlapping oscillations. Hence, it might be that the aperiodic exponent does not change during the absence seizure–instead, the inaccurately removed 3 Hz-harmonics likely caused the increased 1/f exponent value.

Fig. 4
figure 4

FOOOF cannot characterize oscillation peaks that are not clearly distinguishable. a) Left: Time series of an absence seizure measured using EEG. Turquoise: Pre-seizure, red: seizure, yellow: post-seizure activity. Right: Corresponding PSDs and aperiodic FOOOF fits. Note the increase of the 1/f exponent during the seizure. b) Left: Simulated 1/f noise and temporarily (red) added 3 Hz saw-tooth signal. Right: Aperiodic FOOOF fits. Note the increase of the 1/f exponent despite constant ground truth of \({\beta }_{\text{truth}}=1.8\)

Full size image

Figure 4b shows a time series of simulated 1/f noise with an exponent of \({\beta }_{\text{sim}}=1.8\). During the same time interval in which the absence seizure in Fig. 4a occurs, we add a saw-tooth oscillation of 3 Hz to the signal. As in the example of the real seizure in Fig. 4b, FOOOF estimates a strongly increased 1/f exponent even though the ground truth exponent remains constant. The corresponding model (default parameters) is shown in SI Fig. 3.

Note that it is possible to enable FOOOF fitting of the many harmonious peaks by reducing the maximum peak width limits to 1 Hz. While it is not feasible to tune the parameters across conditions (there is an alpha peak with a peak width larger than 1 Hz in the pre-and post-condition), even with the specifically tuned parameters, FOOOF returns increased 1/f exponents \(\beta =2.28\) and \({\beta }_{\text{sim (tuned)}}=1.91\) for real and simulated data, respectively SI Fig. 3.

Recommendation

Scenario A - the PSD appears as a straight line with well-distinguishable peaks on top of this line:

Challenge 3 does not apply.

Scenario B - the PSD might contain overlapping peaks:

The more peaks overlap, the less accurate the model results will be. The lower and upper fitting range borders must vastly extend the overlapping oscillation peaks (challenge 2) to enable peak removal. Estimating the power of overlapping peaks will be difficult.

Scenario C - almost the full PSD seems to consist of overlapping peaks (as in Fig. 4):

Avoid fitting the aperiodic component.

IRASA

Challenge 1: The Evaluated Frequency Range is Larger than the Fitting Range

While FOOOF tries to iteratively fit all oscillatory peaks to obtain a periodic model, IRASA takes the median of spectra after up- and down-sampling to eliminate the peaks, as shown in Fig. 1. As a result, it aims to obtain the pure aperiodic component that is assumed to be invariant to resampling. As an advantage over FOOOF, IRASA can overcome challenge 2: Even if a peak crosses the fitting range border (at the original sampling rate), it can be removed successfully due to the resampling procedure.

In Fig. 5a, we replot the spectrum of Fig. 3a and estimate the 1/f exponent for all frequency ranges between 1–100 Hz and 80–100 Hz. In contrast to FOOOF, IRASA has minimal errors for all frequency ranges. The reason is that the fitting range of FOOOF has well-defined borders: If the lower border is set to 5 Hz (the center frequency of the first peak), it cannot identify and model the 5 Hz peak correctly. On the other hand, for IRASA, the fitting range is blurry: By up- and down-sampling the spectrum, the peaks are shifted towards lower and higher frequencies. Therefore, the evaluated frequency range of IRASA is much more extensive than the actual fitting range.

Fig. 5
figure 5

IRASA’s evaluated frequency range is larger than the fitting range. a) Upper panel: Same simulation as in Fig. 3a. Lower panel: The lower fitting range border is shown on the x-axis, the absolute deviation from the ground truth on the y-axis. IRASA correctly estimates the 1/f exponent for all used fitting ranges. b) Simulated aperiodic PSD with a ground truth of \(\beta =2\). A 1 Hz highpass filter disrupts the 1/f power law. IRASA’s fitting range for the maximum resampling factor \({h}_{\text{max}}\in \{2, 8, 15\}\) is indicated as bright-colored lines upon the fitted aperiodic components, with the evaluated frequency ranges after up- and down-sampling indicated in corresponding transparent colors. IRASA’s error of the 1/f estimation increases with larger resampling rates \({h}_{\text{max}}\) (and lower resampling rates 1/\({h}_{\text{max}}\), respectively). c) Same as b) with a spectral plateau disrupting the 1/f power law. d) FOOOF 1/f estimate within 1–30 Hz for a spectrum obtained from voxel data after MEG source reconstruction. e) IRASA 1/f estimates for an evaluated frequency range of 1–30 Hz (green) and an evaluated frequency range of 0.3–90 Hz (green-dashed, corresponding to a fitting range of 1–30 Hz at hmax = 3). f) FOOOF (blue) and IRASA (green) estimates of the 1/f exponent for the same fitting range of 1–30 Hz

Full size image

For example, if we chose only two resampling factors \({h}_{\text{set}}=\{2, 3\}\), the spectrum would be up-sampled by \({h}_{\text{up1}}=2\) and \({h}_{\text{up2}}=3\) and down-sampled by \({h}_{\text{down1}}=1/2\) and \({h}_{\text{down2}}=1/3\). As a result, a fitting range of 10–100 Hz would correspond to four evaluated frequency ranges of 20–200 Hz, 30–300 Hz, 5–50 Hz, and 3.3–33 Hz. Of these four resampled spectra, IRASA takes the median. The lower border of the evaluated frequency range \({f}_{\text{eval. min}}\) can be calculated from the minimum fitting range border \({f}_{\text{fit min}}\) divided by the maximum resampling factor \({h}_{\text{max}}\) according to Eq. (1). The upper evaluated frequency border \({f}_{\text{eval. max}}\) corresponds to the upper fitting range border \({f}_{\text{fit max}}\) multiplied by \({h}_{\text{max}}\) according to Eq. (2).

$${f}_{\text{eval. min}} = {f}_{\text{fit min}}/{h}_{\text{max}}$$
(1)
$${f}_{\text{eval. max}} = {f}_{\text{fit max}}\cdot {h}_{\text{max}}$$
(2)

While evaluating a larger frequency range than the actual fitting range can be advantageous, as shown in Fig. 5a, it can also lead to severe challenges, as shown in Fig. 5b. Here, we simulate an aperiodic PSD with \(\beta =2\), which is highpass filtered at 1 Hz. We then apply IRASA in a fitting range of 2–30 Hz for three different h-sets with maximum resampling factors \({h}_{\text{max}}=2\), \({h}_{\text{max}}=8\), and \({h}_{\text{max}}=15\), respectively. The fitting ranges, indicated in green, orange, and red, are the same, but the evaluated frequency ranges, shown in the corresponding transparent colors, increase with increasing \({h}_{\text{max}}\).

Note that the highpass filter disrupts the 1/f power law for low frequencies and violates IRASA’s assumption of a resampling-invariant aperiodic component. With increasing \({h}_{\text{max}}\), IRASA evaluates substantially larger parts of the low-frequency stopband which increasingly biases its 1/f estimates towards smaller values. A good agreement with the ground truth of \(\beta =2\) is only obtained for \({h}_{\text{max}}=2,\) which corresponds to an evaluated frequency range of 1–60 Hz, avoiding the stopband of the highpass-filtered spectrum.

Apart from low-frequency fitting artifacts due to highpass filtering, care must also be taken to avoid fitting artifacts at high frequencies. For example, in Fig. 5c, the high-frequency spectral plateau disrupts the 1/f power law. Even though the upper fitting range border of IRASA is set well below the plateau onset to \({f}_{\text{fit max}}=30\) Hz, IRASA does nevertheless evaluate the plateau due to the upsampling step. Therefore, with growing \({h}_{\text{max}}\), IRASA biases the 1/f estimates towards smaller values again.

Even in the absence of a spectral plateau, care must be taken to avoid the resampled Nyquist frequency. For example, for a sampling rate of \({f}_{\text{sample}}=2400\) Hz and \({h}_{\text{max}}=10\), the resampled Nyquist frequency reduces from \({f}_{\text{Nyquist}}=1200\) Hz to \({f}_{\text{Nyquist resampled}}=120\) Hz. Accordingly, the upper fitting range border must not exceed this value. The same holds true for a potentially applied lowpass filter. In general, to avoid accidentally fitting spectra above Nyquist frequency or in the stopbands of lowpass or highpass filters, it is advisable to choose \({h}_{\text{max}}\) as small as possible. Furthermore, the evaluated frequency range should always be checked by calculation from the fitting range and \({h}_{\text{max}}\).

Given that IRASA evaluates a more extensive frequency range than the fitting range, the meaning of the fitting range becomes imprecise. For example, if we are interested in fitting the 1/f exponent from 1–30 Hz and use \({h}_{\text{max}}=2\), we should choose 2–15 Hz as a fitting range for the IRASA algorithm. However, since only the minimum and maximum resampled spectra contain the 1 Hz and 30 Hz borders of interest, IRASA emphasizes the estimation of the 1/f exponent from intermediate frequency values above 1 Hz and below 30 Hz. Therefore, 1/f exponents estimated by IRASA cannot be directly compared to 1/f exponents estimated by FOOOF.

We visualize this effect for a spectrum of voxel data obtained by MEG source reconstruction in the lower panels d) – f) of Fig. 5. In d), FOOOF estimates an 1/f exponent of \({\beta }_{\text{FOOOF}}=1.41\) in the fitting range of 1–30 Hz. Due to the highpass filter, IRASA obtains a lower value of \({\beta }_{\text{IRASA}}=1.09\) for the same fitting range which, however, actually corresponds to an evaluated frequency range of 0.33–90 Hz at \({h}_{\text{max}}=3\) (Fig. 5e). Hence, setting the evaluated frequency range to 1–30 Hz (by setting the fitting range to 3–10 Hz) yields \({\beta }_{\text{IRASA}}=1.38\) which is similar to the FOOOF estimate.

Matching the evaluated frequency range of IRASA to the fitting range of FOOOF is not always possible, though. Consider, for example, the fitting range of 30–45 Hz shown in Fig. 5e. At \({h}_{\text{max}}=3\), the evaluated frequency range of IRASA is 10–135 Hz. Due to the spectral plateau, IRASA estimates a much smaller exponent of \({\beta }_{\text{IRASA}}=1.22\) compared to \({\beta }_{\text{FOOOF}}=2.11\). This time, we cannot shrink IRASA’s fitting range to match its evaluated frequency range with FOOOF’s fitting range. At \({h}_{\text{max}}=3\), the lower fitting range border of IRASA must be \(3\cdot 30\) Hz \(=90\) Hz to match the lower fitting range border of FOOOF at 30 Hz. However, the upper fitting range border needs to be \(45\) Hz \(/ 3=15\) Hz to match the upper fitting range border of FOOOF. This would lead to an inverse fitting range of 90–15 Hz. Here, it cannot be avoided that IRASA evaluates a much more extensive frequency range than 30–45 Hz. As a result, FOOOF and IRASA cannot yield comparable 1/f estimates for this frequency range.

Recommendations

Always calculate the evaluated frequency range from the fitting range and \({h}_{\text{max}}\) according to Eqs. (1) and (2). Choose the maximum resampling factor \({h}_{\text{max}}\) as small as possible in order to 1) avoid fitting artifacts, 2) to improve comparability with other methods, and 3) to improve the interpretability of the investigated frequency range.

Set the evaluated frequency range–and not the fitting range–to the frequency range of interest.

Challenge 2: Broad Peak Widths Require Large Resampling Factors

In challenge 1, we recommend choosing the maximum resampling factor \({h}_{\text{max}}\) as small as possible. However, for IRASA to work correctly, the resampling factors must be sufficiently large. This is because IRASA shifts the peaks in the frequency scale up and down through up- and downsampling. Therefore, a single peak appears multiple times on the frequency scale (Fig. 1b). For a range of sufficiently large (and small) resampling factors, the resampled peaks are completely separated and, by taking the median of their geometric mean, subsequently eliminated. However, if the range of resampling factors is too small or the peaks too broad, the resampled peaks overlap. In that case, peak removal by taking the median will not be successful.

In Fig. 6a, we replot Fig. 5a. However, by increasing the peak widths from the left to the right panels, the 1/f estimation error of IRASA increases strongly. This is because the peaks cannot be fully separated. As a result, IRASA’s calculated aperiodic component, shown in grey, still contains the up and downsampled peaks after taking the median. Note that not the peak width \(\Delta f\) itself must be sufficiently small to get separated, but instead, \(\Delta {f}_{log}\), as it appears in the logarithmic frequency scale, the logarithmic peak width needs to be sufficiently small. For this reason, a peak width of 4 Hz at a center frequency of 5 Hz has a similar effect as a peak width of 12 Hz at a center frequency of 35 Hz.

Fig. 6
figure 6

Broad peak widths require large resampling factors. a) Upper panel: Similar as in Fig. 5a) but with increasing peak widths from left to right. Note that removal of peaks from the aperiodic component (grey) worsens with broader peak widths. Lower panel: The lower fitting range border is on the x-axis, the absolute deviation from the ground truth on the y-axis. The 1/f exponent estimation error increases with larger peak widths. b) Simulation of a 30 Hz and 300 Hz peak with increasing peak widths from left to right. Larger peak widths require larger resampling factors. Note that not the absolute peak width but rather the logarithmic peak width \(\Delta {f}_{log}\) determines the minimum resampling factors

Full size image

We visualize this effect in panel Fig. 6b by simulating a PSD with two oscillations at \({f}_{1}=30\) Hz and \({f}_{2}=300\) Hz. The peak width of the second peak is 70 Hz and therefore 10 times as large as the peak width of the first peak. However, on the logarithmic frequency axis, they appear with the same width. A maximum resampling factor of \({h}_{\text{max}}=2\) is sufficient to remove the peaks correctly. Thus, they are fully eliminated from the aperiodic component shown in turquoise. However, when the peak widths are increased to a logarithmic value of 0.2 log(Hz), \({h}_{\text{max}}=2\) is not sufficient anymore: The up- and down-sampled peaks remain visible in the estimate of the aperiodic component. If we increase \({h}_{\text{max}}\) to a value of 8, however, peak removal works well. For a further increase of the logarithmic peak width to 0.3 log(Hz), however, \({h}_{\text{max}}=35\) is necessary. We visualize this challenge on empirical data of MEG and LFP data of dataset 3 in SI Fig. 4.

We calculated the logarithmic peak width as Eq. (1) \(\Delta {f}_{log}=lo{g}_{10}{(f}_{2}/{f}_{1})\) where \({f}_{1}\) corresponds to the lower bound of the peak and \({f}_{2}\) to the upper bound of the peak. The bounds were found by calculating the first bin of the PSD, which deviates above a threshold of 0.001 from the aperiodic ground truth. Note that there is no exact equation/heuristic to calculate the minimum \({h}_{\text{max}}\) as a function of peak width because always many resampling factors h are calculated, which will lead to a gradual peak removal depending on the degree of peak separation.

Recommendations

Choose \({h}_{\text{max}}\) as small as possible (challenge 1) while keeping it large enough to obtain peak-free estimates of the aperiodic component (challenge 2). If the peaks are very broad and \({h}_{\text{max}}\) cannot be chosen sufficiently large without avoiding challenge 2, IRASA cannot be applied.

Challenge 3: IRASA Cannot Characterize Oscillation Peaks that are not Clearly Distinguishable

Similar to FOOOF, IRASA cannot separate strongly overlapping peaks. However, as shown in Fig. 7b, IRASA performs quite well for dataset 2 because the harmonic peaks do not strongly overlap above 10 Hz. Instead, many local power minima in between the harmonic peaks are very close to the power of the aperiodic ground truth. As a consequence, adding the 3 Hz sawtooth signal only slightly increases the estimated 1/f exponent from \({\beta }_{\text{pre/post}}=2.24\) to \({\beta }_{\text{seiz}}=2.46\). In the middle panel of SI Fig. 5b, the extracted oscillatory component of IRASA is shown in orange, indicating a good extraction of harmonic peaks at multiple integers of 3 Hz.

Fig. 7
figure 7

IRASA cannot characterize oscillation peaks that are not clearly distinguishable. a) and b) left panel: Same as Fig. 4a) and b) Right panel: Same as Fig. 4 but showing the 1/f fits by IRASA. c) IRASA’s performance on the simulation drops significantly if two strongly overlapping peaks in the alpha (10 Hz) and beta range (25 Hz) are added. Ground truth: \({\beta }_{\text{truth}}=1.8\)

Full size image

If, however, we now add two strongly overlapping oscillations at 10 Hz and 25 Hz, IRASA is no longer capable of successfully removing the peaks. As a result, it now estimates an exponent of \({\beta }_{\text{seiz}}=3.05\) – much larger than the ground truth at \({\beta }_{\text{truth}}=1.8\).

Discussion

Both periodic and aperiodic components of power spectra are frequent targets of investigation in electrophysiological studies. The separation of both components before analysis helps to disentangle their relative contribution to the spectrum. FOOOF and IRASA are commonly used for this purpose. It should be highlighted, though, that the methods follow different concepts: Whereas FOOOF models periodic components and a single aperiodic component and outputs the corresponding parameters, IRASA only separates them, allowing further independent processing. Other methods to separate periodic and aperiodic PSD components exist too, for example, eBOSC (Kosciessa et al., 2020). However, they cannot overcome the method-unspecific challenges in electrophysiological PSDs such as 1) spectral plateau onsets at relevant frequencies, 2) hidden low-frequency oscillations, and 3) overlapping peaks.

Here, we evaluated common challenges of the separation procedure based on two popular methods and summarized both general and method-specific challenges. These challenges apply to EEG, MEG, and LFP data obtained by independent research groups, indicating the general applicability of the results.

Neurophysiological Interpretation

Spectral Plateau

The spectral plateau can hinder a correct separation of the PSD components. If the fitting range of FOOOF or the evaluated frequency range of IRASA includes a spectral plateau, the 1/f exponent will be estimated too low. In the presence of periodic components at the spectral flattening, a faulty aperiodic power estimation will lead to a faulty periodic power estimation. In addition, large periodic components could hide the onset of the spectral plateau. This hinders a proper decision on the choice of the upper fitting range border. If the 1/f exponent is estimated too low due to the spectral flattening, this could be misinterpreted as an increased E–I ratio since typically flatter spectra are associated with more pronounced excitability (Gao et al., 2017).

E–I balance estimation is usually applied to estimate relative 1/f differences between conditions. If the spectral plateau onset were to occur at exactly the same frequency for all spectra, a relative 1/f comparison would still be viable. A random fluctuation of the onset would introduce noise to the estimates, and a systematic difference of the plateau onset between conditions would lead to type 1 errors.

The origin of the spectral plateau at high frequency is likely rooted in Gaussian properties of amplifier noise and impedance noise (Scheer et al., 2006; Waterstraat et al., 2015a, b). (Waterstraat et al., 2021) showed that the onset of the plateau starts at higher frequencies if recordings are done with a low-noise MEG system. In addition to system noise, biological high-frequency noise caused by electromyography (EMG) from the head muscles can contribute to spectral flattening, although EMG does not necessarily have a flat spectrum (Muthukumaraswamy, 2013). However, even when recording LFPs from the subthalamic nucleus using a low-noise amplifier, which can be considered as hardly affected by EMG activity, a spectral plateau could be observed one order of magnitude above the system's noise level (unpublished data). Neuronal population spiking activity probably contributes to this spectral plateau (Belluscio et al., 2012; Buzsáki et al., 2012; Zanos et al., 2011). Better understanding the origins of the spectral plateau is of major interest and requires further research. If it is caused by noise such as systems noise or EMG, an identification of the origin could help to clean the data from this high-frequency plateau. If it has a neurophysiological origin, a thorough analysis of the plateau might yield novel neurophysiological insights.

Fitting Ranges

The choice of the fitting range depends on the goal of the study and the properties of the data. In the literature, the 1/f exponent was investigated for different frequency ranges such as 0.01–0.1 Hz (He et al., 2010), 0.5–35 Hz (Miskovic et al., 2019), 1–10 Hz (Schaworonkow & Voytek, 2021), 1–20 Hz (Bédard et al., 2006), 1–30 Hz (Wen & Liu, 2016), 1–40 Hz (Colombo et al., 2019), 1–20 and 20–40 Hz (Colombo et al., 2019), 1-15 and 15-125 Hz (Chaoul & Siegel, 2021), 1–100 Hz (He et al., 2010), 2–24 Hz (Voytek et al., 2015), 3–30 Hz (Pereda et al., 1998), 3–55 Hz (Waschke et al., 2021), 10–100 Hz (Freeman & Zhai, 2009), 20–65 Hz (Bédard et al., 2006), 30–50 Hz (Gao et al., 2017; Lendner et al., 2020; Stolk et al., 2019) and 40–60 Hz (Gao et al., 2017).

If a study using FOOOF aims to generally estimate an 1/f exponent and is free to choose any fitting range for that purpose, we generally recommend avoiding low lower fitting range borders. It is unknown how hidden low-frequency oscillations at the lower fitting range border might impact the 1/f estimate. For example, if the 1/f exponent is compared between two conditions and in one condition there are larger delta oscillations, a fitting range starting from 1 Hz could have a larger y-intercept due to the presence of low-frequency oscillations. This could lead to a larger 1/f exponent and could be misinterpreted as stronger neural inhibition. The same holds for any other lower fitting range border. Therefore, the lower border should be chosen to best avoid known oscillation frequencies, depending on the study.

The same problem applies to IRASA but less severely since this method is not based on a single frequency range but rather, due to up- and downsampling, to a set of different frequency ranges. This comes at the cost of only vaguely defined upper and lower fitting range borders, hindering an easy comparison with fitting ranges used in other studies.

In general, there is no one-range-fits-all fitting range applicable to all kinds of PSDs. Therefore, we recommend examining the PSDs of interest carefully and choosing the fitting range that best avoids the challenges discussed so far in addition to further possible data-specific or goal-specific challenges.

Finally, if the purpose of the 1/f estimation is not to obtain the 1/f exponent but rather the removal of the aperiodic component for better periodic power assessment, a broadband range (such as 1–100 Hz) should be chosen.

Overlapping Peaks

The stronger periodic components overlap, the more difficult estimating their power becomes. In the shown exemplary data, overlapping peaks occurred mainly in STN data and in dataset 2. MEG and EEG cortical data of healthy participants typically have peaks in the ranges 8–13 Hz and 18–25 Hz which does not impose considerable challenges for the estimation of the aperiodic part of the spectrum. Assessing the 1/f exponent is still feasible if the overlapping periodic components make up only a minor part of the fitting range. On the other hand, if the overlapping peaks make up a majority of the frequency range to investigate, as in Figs. 4, 7, and 8b, a separation of the periodic and aperiodic components is not recommended and will likely lead to imprecise results. In the case of the absence seizure shown in Figs. 4 and 7, neural inhibition during the seizure is likely overestimated due to overlapping peaks leading to a false 1/f estimation.

Fig. 8
figure 8

“Easy” and “hard” PSDs. a) Left: Voxel MEG PSD of a Parkinsonian patient on a semilogarithmic scale. Right: Same PSD on a double logarithmic scale. FOOOF, IRASA, and simply connecting the PSD value at 1 Hz to the PSD value at 95 Hz as a straight line (“straight”) yield similar 1/f exponents. We regard such a PSD as “easy” because it avoids all the discussed challenges. b) LFP data of a Parkinsonian patient on a semilogarithmic scale. Right: Same PSD on a double logarithmic scale. FOOOF, IRASA, and “straight” yield different 1/f exponents. We regard such a spectrum as “hard” because it contains many challenges

Full size image

Broad Peak Widths

In contrast to FOOOF, IRASA cannot handle very broad peaks well. This limitation is especially severe for the analysis of LFPs of Parkinsonian patients (datasets 1 and 3). In the original article (Wen & Liu, 2016), IRASA was only evaluated on pure sine oscillations, for which the method works very well. We, therefore, do not recommend using IRASA if the peaks seem to have broad logarithmic peak widths \(\Delta {f}_{log}\). It is not possible to give a threshold value for a maximum logarithmic peak width as IRASA is, in theory, able to fit any peaks if the h-values are chosen sufficiently large. However, in practice, the h-values must not exceed a certain range to avoid too low (highpass, Fig. 5b) or too high (spectral plateau, Fig. 5c) fitting ranges, calculated using Eq. (1) and (2).

Estimating E–I Balance

(Gao et al., 2017) proposed to use the 1/f exponent as an indicator of E–I balance. Subsequent studies indicated the usefulness of this idea also for non-invasive EEG/MEG data (Colombo et al., 2019; Gao et al., 2017; Lendner et al., 2020; Miskovic et al., 2019; Waschke et al., 2021). While we outlined that 1/f exponent estimation is affected by many possible error sources, we do not argue that it should be avoided altogether. While proper 1/f estimation seems to be beyond reach for some PSDs (for example, the one shown in Fig. 8b), it seems to be a promising measure for others (Fig. 8a). Thus we suggest that existing methods could be enhanced by more elaborate data cleaning, such as spatio-spectral decomposition (SSD) (Nikulin et al., 2011), independent component analysis (ICA), or inverse modeling. Moreover, it might be possible to develop new methods that measure the 1/f exponent more reliably than the ones discussed in the present study. For example, if the periodic and aperiodic components are assumed to vary over time independently, it could be possible to disentangle them using machine learning algorithms such as non-negative matrix factorization (Lee & Seung, 1999). And it might become possible to measure E–I balance through other electrophysiological measures thus further validating the 1/f exponent of the PSD. We elaborate on this below.

(Bruining et al., 2020), for example, proposed to measure E–I balance based on the alpha band amplitude envelope and its detrended fluctuation analysis (DFA) exponent (Peng et al., 1995). (Stephani et al., 2020) related the N20 of somatosensory evoked potentials to cortical excitability. Other researchers related spontaneous fluctuations of alpha-band power to E–I balance (Romei et al., 2008) and (Iemi et al., 2019) found alpha- and beta-band power to predict suppression of ERP-components, which was interpreted as increased inhibition. This relationship held true even after controlling for fluctuations in the 1/f exponent, which correlates with alpha power (Muthukumaraswamy & Liley, 2018). It might be possible to estimate E–I balance by measuring transcranial magnetic stimulation (TMS) evoked potentials using EEG. By combining these two methods, (Massimini et al., 2005) showed a breakdown of effective cortical connectivity during non-rapid eye movement (REM) sleep. Effective connectivity was also related to the 1/f exponent by (El Boustani et al., 2009). The perturbational complexity index (Casali et al., 2013) follows these lines to separate unconscious states of low excitability (non-REM sleep, anesthesia) from conscious states of high excitability (wakefulness, REM sleep). Indeed, (Colombo et al., 2019) could link this index to the 1/f exponent during wakefulness and anesthesia yielding similar results with both methods. However, it should be noted that REM sleep (a conscious state of mind) is associated with a larger 1/f exponent compared to NREM sleep (unconscious) while NREM sleep is associated with a larger 1/f exponent compared to wakefulness (conscious) (Lendner et al., 2020). These findings agree with in vivo calcium imaging measurements of E–I balance in mice during wakefulness, NREM sleep, and REM sleep (Niethard et al., 2016).

In the best scenario, different methods used for E–I estimation will lead to similar results and might be used in conjunction. So far, the relationship between 1/f exponent and E–I balance remains a hypothesis to be further validated.

Computational Cost and Parameter Tuning

From a computational perspective, FOOOF is much faster than IRASA. When applied to 9 time series of dataset 1 (ca. 180 s at \({f}_{sample}=2400\) Hz corresponding to \(\approx 9\times \mathrm{440,000}\) data points), parameterization with FOOOF was about 50 times faster than separation with IRASA when the PSD calculation time was included. FOOOF was 100 times faster if the PSDs were precalculated. For the comparison, we used 7 runs and fitted a frequency range from 1–30 Hz. For FOOOF, the default parameters were chosen, and for IRASA a window length of 4 s and a set of 17 resampling factors \({h}_{set}=\left\{1.1, 1.15, ..., 1.9\right\}\) was used. FOOOF computation slows down when the PSDs have a very high resolution leading to many iterations of fitting noise peaks. IRASA computation slows down when the number of resampling factors is increased and when their values are increased.

Increasing IRASA’s resampling values can help with very broad peak widths (challenge 2) but simultaneously enlarges the evaluated frequency range (challenge 1). Increasing the number of resampling factors beyond 17 or changing the window length does not help with the challenges presented in this article. FOOOF requires extensive parameter tuning for optimum results, but the posed challenges cannot be resolved by parameter selection. In general, the fitting range of FOOOF and the evaluated frequency range of IRASA are the most critical parameters for each method.

Conclusion

To study either periodic or aperiodic PSD components, it is useful to disentangle both components. As there are theoretically infinite solutions to this inverse problem, it is probably neither possible to perfectly separate them nor to evaluate and verify a performed separation since the ground truth remains unknown. Some PSDs seem to be particularly easy to separate because they avoid most of the discussed challenges. For those PSDs, we generally recommend performing a separation to study the periodic or aperiodic components in a more isolated manner. We give an example of such an “easy” PSD in Fig. 8a.

These “easy” PSDs appear as an almost straight line in double logarithmic space with some well-distinguishable, narrow periodic peaks on top of it. There is no spectral plateau disrupting the 1/f power law and the y-axis of the PSD extends over 4 orders of magnitude. When applying FOOOF and IRASA from 1–95 Hz or simply connecting values at 1 Hz and 95 Hz to a straight line in double logarithmic space, similar values (\({\beta }_{\text{FOOOF}}=1.65\), \({\beta }_{\text{IRASA}}=1.71\), \({\beta }_{\text{straight}}=^{1.53}\)) are obtained for the 1/f exponent.

Other PSDs seem to be very difficult to separate. For such, we recommend avoiding the separation since the results will be arbitrary and might lead to ill-informed interpretations. An example of such a “hard” PSD is shown in Fig. 8b. These spectra do not appear as a straight line. They have very broad and overlapping peaks and a spectral plateau onset at lower frequencies. As a result of this plateau, the y-axis spans only one order of magnitude. When applying FOOOF, IRASA, or a straight-line connection between 1 and 95 Hz, strongly diverging 1/f exponent values (\({\beta }_{\text{FOOOF}}=0.82\), \({\beta }_{\text{IRASA}}=1.10\), \({\beta }_{\text{straight}}=0.62\)) are obtained.

Checking PSDs for the challenges discussed in this work will help to decide whether a technique to separate neural oscillations from aperiodic 1/f activity should be applied, which algorithm to use, and which parameters to choose.

Information Sharing Statement

The datasets analyzed in the present study are not publicly available due to privacy regulations of patient health information but are available from the corresponding author upon reasonable request. The code of the entire study is available at https://github.com/moritz-gerster/separating_periodic_from_aperiodic_PSDs.