Alpha blocking and 1/fβ spectral scaling in resting EEG can be accounted for by a sum of damped alpha band oscillatory processes

The dynamical and physiological basis of alpha band activity and 1/fβ noise in the EEG are the subject of continued speculation. Here we conjecture, on the basis of empirical data analysis, that both of these features may be economically accounted for through a single process if the resting EEG is conceived of being the sum of multiple stochastically perturbed alpha band damped linear oscillators with a distribution of dampings (relaxation rates). The modulation of alpha-band and 1/fβ noise activity by changes in damping is explored in eyes closed (EC) and eyes open (EO) resting state EEG. We aim to estimate the distribution of dampings by solving an inverse problem applied to EEG power spectra. The characteristics of the damping distribution are examined across subjects, sensors and recording condition (EC/EO). We find that there are robust changes in the damping distribution between EC and EO recording conditions across participants. The estimated damping distributions are found to be predominantly bimodal, with the number and position of the modes related to the sharpness of the alpha resonance and the scaling (β) of the power spectrum (1/fβ). The results suggest that there exists an intimate relationship between resting state alpha activity and 1/fβ noise with changes in both governed by changes to the damping of the underlying alpha oscillatory processes. In particular, alpha-blocking is observed to be the result of the most weakly damped distribution mode becoming more heavily damped. The results suggest a novel way of characterizing resting EEG power spectra and provides new insight into the central role that damped alpha-band activity may play in characterising the spatio-temporal features of resting state EEG.

Introduction Electroencephalography (EEG) is a non-invasive method used to measure the electrical activity of the brain at the surface of the scalp, with the recorded voltage fluctuations being generated by ionic current flows resulting from neural synaptic activity across the cortex [1]. EEG time series recordings display a rich array of dynamic activity that is distributed spatio-temporally across the head, providing a unique window into the inner workings of the human brain. Since its discovery the EEG has been widely used as a sensitive measure of brain state in disease and health [2,3]. We aim to explore two prominent features of the EEG that are well observed but of which the mechanistic origins are the subject of continued speculation; the alpha rhythm and its changes between eyes-closed (EC) and eyes-open (EO) resting states, and the 1/f β frequency scaling of the power spectrum.
First discovered and recorded in the early 20th century by Hans Berger [4], the alpha rhythm (8)(9)(10)(11)(12)(13), is arguably the most dominant feature of the resting EEG. Classically, the alpha rhythm is observed as waxing and waning oscillations in EEG time series recordings of the EEG, which in the frequency domain appear as a spectral peak between 8 and 13 Hz [1]. Alpha band activity, defined as the 8-13 Hz spectral power, will include such alpha rhythmic activity but may also include activity from mechanistically unrelated cortical and extra-cortical non-oscillatory sources.
There exists significant heterogeneity of alpha band activity across the population in the terms of its the magnitude of its activity, its peak/central frequency and its reactivity to changes in behavioural and physiological state [5]. Topographically alpha band activity can be recorded across the scalp with it being particularly prominent occipitally [6,7]. Despite nearly a century of detailed empirical investigation, the physiological mechanisms responsible for the genesis of alpha rhythm remain essentially unresolved. The prevailing view is that the thalamus is central to the generation of the alpha rhythm through reverberant feed forward and feedback interactions between thalamic and cortical neuron populations [8,9]. Such a view was born out of earlier conceptions where alpha oscillations, intrinsically generated in the thalamus were thought to directly drive activity in overlying cortex [10]. However, other attempts to explain the genesis of the alpha rhythm also exist that depend on a mean field or neuronal population framework to model reverberant activity solely between excitatory and inhibitory cortical neuronal populations [11]. Complicating efforts to develop mechanistically and physiologically coherent accounts of the genesis of alpha band activity is the well described reduction in alpha band power between EC and EO conditions and in response to the exertion of mental effort, known commonly as alpha blocking [9,12]. alternative mechanisms to criticality for the genesis of scale free activity: '1/f ' noise could arise simply from a mixture of alpha band relaxation processes that having a distribution of dampings. Such a hypothesis lies in contrast to the more common view that alpha band activity and '1/f ' noise are functionally and mechanistically distinct and thus electroencephalographically independent [36,37].
To date efforts to model the dynamic activity of the EEG fall into two broad categories. The first and most intuitive is the network modeling approach, which attempts to construct the EEG dynamics from the ground up by focusing on modelling individual cortical neurons and their vast networks of interactions. Such methods are limited by the computational uncertainty of simulating the interaction of hundreds of millions of neurons that lie under even a single EEG sensor. The ambiguity of how much physiological detail to include in these network approaches and how to meaningfully parameterize the model is a critical challenge [38]. An alternative approach are neural population models or mean-field theories. Motivated by the fact that a single EEG electrode records the aggregate activity of millions of cortical neurons, mean-field models account for EEG activity by modelling the average activity of interconnected populations of excitatory and inhibitory neurons in cortex and thalamus [39]. Mean-field models treat cortical and thalamic tissue as spatially continuous and model mesoscopic-scale cortical dynamics, and thus provide a major advantage in that they are constructed at spatial scales that are more commensurate with the physical coverage of an individual EEG channel (mm 2 -cm 2 ) [38]. Given that most mean field approaches generally divide the surface of the brain into populations of interconnected cortical neuron populations and investigates bulk dynamic properties, the number of modelled elements is drastically fewer than a neural network approach.
Mean field models can be mathematically quite compact and are typically formulated as multiple coupled nonlinear ordinary or partial differential equations [38]. Under suitable assumptions the defining non-linear equations can be linearised and in doing so reveal a rich range of electroencepahlographically plausible features including alpha band activity [11,40,41]. In such linearisations alpha band oscillatory activity is generally accounted through a linear time invariant transfer function driven by broad-band white noise, which pragmatically accounts for the spatio-temporal complexity of extra-cortical neuronal input. This approach is generally agnostic to the specific mean-field model used and enables the derivation of an electrocortical transfer function that describes the spectral response of the system to white noise driving.
Given that most mean-field models of the EEG attempt to produce the alpha rhythm while simultaneously requiring a stable linearisation, the general form of the electrocortical impulse response embodies two essential features: i) it will produce oscillations that fall within the alpha band, and ii) will decay away to zero after perturbation, ensuring its stability. Regardless of the particular mean field theory chosen, qualitatively the modelled spectral behaviour of the resting EEG will be similar in form; a sharp alpha resonance followed by decline in spectral power with frequency.
The essential dynamical premise for this linearisation, that resting EEG arises from the random driving of a stable quasi-linear system having an oscillatory alpha resonance, is consistent with the broader literature regarding efforts to characterise the dynamical architecture of resting EEG. In a range of states (resting state, sleep, anaesthesia) the spontaneous scalp-recorded EEG is assessed as being indistinguishable from a linear random process [29,[42][43][44][45][46][47][48]. Thus mean field models, and the approaches derived from them such as ours, are to be seen as physiologically plausible given the principles of their formulation and the dynamics of the modelled electroencephalographic activity.
The current uncertainty associated with the mechanisms responsible for alpha blocking and '1/f ' noise in resting state EEG leaves many open questions. For this reason we explore whether an alternative mechanism, that is theoretically motivated by mean-field models of the EEG and a general model for 1/f noise [18], can account for both of these phenomena within a single mechanistic framework.
On this basis we construct a physiologically well-motivated mathematical description of the resting EEG power spectrum arising from two well established empirical predicates: resting EC/EO EEG is dynamically equivalent to a random linear process, the power spectral density of EC/EO resting EEG has a clear alpha band spectral peak. Specifically, this is achieved this by constructing an electrocortical impulse response, that when driven by a broadband (i.e. spectrally flat) noise process, approximates the main spectral features of the recorded EEG signal. On this basis the simplest electrocortical impulse response for the system under investigation is a damped linear oscillation (specifically a decaying cosinusoid) that is parameterised by its oscillation frequency (f α = 8 -13 Hz) and damping (γ). In the frequency domain the power spectral density of such an impulse response will appear approximately Lorentzian in shape, in that it will have a sharp alpha peak and scale inversely with the square of the frequency (see Fig 1a). By extending this description to the more physiologically plausible case where the EEG signal recorded at a single location is composed of a sum of damped linear oscillations arising Resting EEG as a sum of damped independent alpha oscillatory processes. (a) Power spectral density of the model electrocortical impulse response cos(2πf α t) exp(−γt)Θ(t) for γ = {0.1, 1, 10, 50}s −1 , where Θ(t) is the Heaviside step function. Resting EEG is modeled as a continuous sum of such processes over some suitable interval of dampings [γ l , γ h ]. In general we aim to find the appropriate weights {w i } that account for the shape of the resting EEG power spectrum. (b) A 1/f spectrum (dashed blue) can be approximated as a sum of 1/f 2 relaxation processes (solid blue) that are uniformly distributed in damping over some interval which for the displayed case is [γ l = 10 −2 , γ h = 10 2 ] s −1 . For clarity we have assumed f α = 0 Hz. https://doi.org/10.1371/journal.pcbi.1010012.g001 from multiple neuronal populations, we can take advantage of a general approach in which 1/f power spectral profiles can arise as the result of a sum of exponentially decaying processes having a uniform distribution of decay rates [18] (Fig 1b). Such an assumption is the natural extension of the well known idea that the EEG arises from the superposition of multiple underlying oscillatory processes, such as excitatory postsynaptic potentials on the dendrites of excitatory [1], weakly connected alpha band limit cycle oscillators having a distribution of phases [1,49,50], or damped oscillatory processes having different frequencies [42,[51][52][53][54][55].
Thus, our central, physiologically well-motivated, mathematical hypothesis is that the resting EEG power spectrum can be well described as a superposition of populations of uncorrelated damped alpha band linear oscillations that have a distribution of dampings. On this basis we seek to determine whether changes in the underlying distribution of dampings of alpha band relaxation processes are able to sufficiently account for '1/f ' noise and the spectral changes associated with alpha blocking, using two freely available EEG data sets (combined total of 136 participants) containing EC and EO resting state time series data. The current work follows directly from the approach outlined in Muthukumaraswamy & Liley (2018) [16]. However, in contrast we significantly generalise this approach by making no a priori assumption on the functional form of the damping distribution. The empirical estimation of the distribution of dampings are found by solving an inverse problem by applying regularization methods to a defining Fredholm integral equation of the first kind. The numerical estimation of the subject level damping distribution for alpha band oscillatory activity provides an alternative methodology to derive information about the dynamic activity of the resting EEG and a novel way of characterising changes in recorded power spectral densities.
The paper is organized as follows. The method section begins with a detailed mathematical specification of the inverse problem to be solved in order to estimate the distribution of alpha relaxation oscillation dampings from recorded EEG power spectral densities. This is then followed by a description of the Tikhonov regularisation method used to numerically solve the the corresponding Fredholm integral equation of the first kind. Finally in the methods, we detail the empirical EEG data sets used and the procedures used to ensure the numerical fidelity, and the properties, of our inverse solutions. Detailed results are then presented, followed by a discussion regarding the essential results.

Numerical validation of Tikhonov regularised inverse method
To ensure that our inverse method was able to meaningfully recover damping probability distributions we conducted a simulation study. Six distinct prior damping distributions (Fig 2a-2f) together with 100 regularised parameters spanning the interval [10 −6 � λ � 1] were generated. The damping distributions for a given regularisation parameter that resulted in the minimum RSS error was deemed the best solutions. The regularisation scheme performed well for a unimodal distributions (Fig 2a). For bimodal distributions, the regularisation method correctly recovered the number and position of modes for both cases (equal and different peak magnitudes) with small discrepancies evident in the second mode where the peak amplitude was less in the regularised damping distribution than in the test distribution (Fig 2b and 2c). The trimodal distribution with constant peak amplitude (Fig 2d) was recovered well with only minor variations recorded in the width of the second and third mode. The recovery for the trimodal distributions with variable mode amplitudes (Fig 2e and 2f) performed in a similar manner to the bimodal case, where the second and third mode amplitudes of the recovered distributions were different to the test distribution. In both trimodal examples where the peak magnitudes are not equal, we observe the second distribution mode peak position has minor differences compared to the test distribution. Despite some overt variations between the recovered and test distributions, the key features (number of peaks, peak positions, peak width and peak magnitudes) are largely replicated and were well recovered using our regularised inverse method.

Empirically estimated damping distributions
The empirically estimated damping distributions, across both EC and EO conditions were largely found to be bimodal, at 82% and 70% respectively. Examples of empirically recovered distributions are shown in Fig 3 along with the respective EC and EO condition power spectral densities. To determine whether the bimodal structure was a numerical artefact, we tested whether the removal of the second mode still permitted good spectral fits in the cases where multiple modes were present. The second distribution mode was removed by setting all relevant weightings to 0. The altered damping distributions were all renormalized to appropriately account for the removal of the second mode weightings and were then used in the forward model Eq 14 to generate spectral fits. We found that eliminating the second mode resulted in an inability to fit the high frequency scaling of the power spectrum, indicating that its presence was necessary to obtain optimal fits. As expected we found that power spectra with more pronounced alpha peaks were generally associated with damping distributions that were more weakly damped using the weakly damped measure described in Methods. The area of the unnormalized empirically estimated damping distributions were not significantly different between EC and EO conditions, suggesting that state induced spectral power changes are predominantly the results of changes in empirically estimated probability distribution shape (See S1 Fig for details). Further, this shape did not depend significantly on the exact form of the integral kernel used. Using the kernel of Eq (5) instead of Eq (6) made little difference to the shape of the empirically estimated damping distributions and no difference to the number of inferred modes (results not shown). It is entirely possible though, that using alternative kernels (Gaussian, Gamma) may result in distributions of differing structure, that may not present with the same bimodal distribution in dampings. However, these functional forms are not immediately compatible with the theoretical and experimental framework that we have used. As such, it would be more difficult to motivate their use over the approximate Lorentzian form employed here.

Spectral fit
By using the empirically estimated damping distributions in the forward problem (i.e Eq 14) to generate model spectra, we are able to demonstrate that we are able to obtain good spectral fits that replicate the power-law like scaling measured in the original experimental data. All the model spectra were generated using the maximum entropy damping distributions obtained by evaluating a range of regularisation parameters with fit quality evaluated using a RSS between the original EEG spectra and the forward model spectra.  (a-f) Damping distributions obtained from EC (blue) and EO (red) power spectra for six individuals (three per data set) who presented with strong alpha blocking. Plotted alongside the experimental power spectra are the respective model fits (EC-blue, EO-red) generated using the EC and EO damping distributions in the forward problem Eq 14. The inferred damping distributions are clearly bimodal in shape, a feature common across both data sets. The first mode of the EC damping distribution is generally peaked at a smaller damping value and has a larger respective area when compared to the EO distributions. We analysed the distribution of RSS to gain insight into the fit quality most prevalent in our model. The median of the distribution (EC and EO) corresponds to model fits that are of high quality. Even in instances of poor fit qualities (bottom 5% of fit quality) the experimental spectra are reproduced reasonably well. The median value for the RSS distribution in EC (0.0079) was slightly lower than EO (0.0090) but was statistically significant (Δ median = −0.0011, p = 10 −6 ).
The power spectrum scaling exponents, β, were computed across all subjects and channels for both the experimental spectra and the model by fitting a 1/f β profile over the frequency range 15-40 Hz. Pooled data (N = 8677 β values for experimental and model spectra across each state) analysis revealed a range of experimental spectral scaling across EC and EO conditions (Fig 5a) with the most probable values of β E � 1.8 and β E � 1.7 respectively. The broad range of scaling behaviour recorded in the experimental power spectra was reproduced in the forward model spectral fits (Fig 5b). Across the pooled data the model and experimental spectral exponents are well correlated in EC (r = 0.64, p = 10 −6 ) and EO conditions (r = 0.81, p = 10 −6 ) (Fig 5c). In general the shape of (in particular alpha peak location, amplitude and high frequency scaling) most power spectral densities could be well accounted for as a sum of damped alpha band oscillatory processes. However, the alpha peak being heavily attenuated (damped) and/or the presence of significant resonant beta band (13-30 Hz), which were not explicitly accounted for in our approach, were major contributing factors to poor model fit. We discuss this issue in more detail in the Discussion.

Alpha blocking in resting state EEG
To quantify alpha blocking we used resting state data recorded in EC and EO conditions from two separate EEG data sets having a combined total of 136 participants, and computed channel based power spectral densities using Welch's method. A total of 8677 (109 subjects × 64 channels + 27 subjects × 63 channels) power spectra were estimated for each recording condition. The Jensen-Shannon divergence was then computed between each EC/EO spectral pair to measure the degree of alpha blocking across the respective data sets. A range of alpha blocking was observed across the both data sets; some participants evinced little blocking whereas others exhibited strong alpha band attenuation between states. Example cases are plotted in Fig 6a-6f where power spectra from an occipital electrode (O1) are plotted together in descending order of the degree of alpha blocking. Fig 6g and 6h shows the topography for two strong blockers exhibiting classical occipital dominance. Such an occipital dominance was clearly seen in the respective group averages (Fig 6i and 6j).

Changes in damping distributions can account for alpha blocking
Our estimated damping distributions suggest that alterations in the alpha peak power (see Fig  3) were predominantly driven by changes in the weakly damped mode. We therefor quantified this relation further by investigating the spatial distribution of both the JSD and the grouplevel differences in damping. Changes in the damping distribution between EC and EO conditions were quantified by calculating the difference in weakly damped measure (Eq 25) between EC and EO states (i.e. WDM EC − WDM EO ), on an electrode by electrode basis for both data sets. Group level averages were computed and plotted topographically such that direct comparisons can be made between the topography of alpha blocking quantified using the Jensen-Shannon Divergence (see Fig 6i and 6j).
As expected group-level weakly damped measure difference values were largest in the occipital regions for both data sets (Fig 7a and 7b). However, in addition to a clear occipital dominance, statistically significant differences in the weakly damped measure were topographically widespread, particularly in the second data set. The similarity of the topographic variation of the difference in the weakly damped measure between EC and EO states and the Jensen-Shannon Divergence (Fig 6i and 6j) are clearly apparent. Averaging across subjects in each data set and plotting electrode-wise we see (Fig 7c) that these respective measures are well correlated across data set 1 (r = 0.838, p = 10 −6 ) and data set 2 (r = 0.901, p = 10 −6 ) On this basis we can reasonably conclude that alpha blocking is being driven by increases in the damping of a weakly damped population of alpha band oscillatory processes.

Discussion
Our aim was to explore whether the resting EEG could be accounted for in terms of a sum of damped alpha oscillatory processes and in particular whether changes in the power spectrum between EC and EO states could be reasonably accounted for solely by changes in the damping of such oscillatory processes. Through the use of regularisation methods, we estimated EC and EO damping distributions across a total of 136 participants, available from two independently collected EEG data sets, and characterised their features. We found the model performed well when applied to real EEG spectra, providing high quality fits to both the alpha band and the high frequency power-law-like scaling. The attenuation of the peak alpha power observed in Red asterisks indicate significant (p < 0.05) differences in the weakly damped measure between EC and EO states calculated using the Max/Min-t permutation test (see Methods) and corrected for multiple comparisons. (c) Scatter plot of electrode-wise Jensen-Shannon divergence and the difference in the weakly damped measure between EC and EO, averaged across all participants in each data set (red = data set 1, blue = data set 2). In general the two measures are well correlated with each other: data set 1, r = 0.838, p = 10 −6 ; data set 2, r = 0.901, p = 10 −6 . Note that the grouplevel quantities differ between the data sets, in particular the <JSD>. The reason for this is not clear. It is possible that it is due to differing experimental configurations of the two studies, or more likely, it is due to a larger prevalence of strong alpha blockers in data set 2 (Fig 6j). Statistical significance was validated using a nonparametric permutation testing.
https://doi.org/10.1371/journal.pcbi.1010012.g007 PLOS COMPUTATIONAL BIOLOGY heavy alpha blockers seems to be well explained by an increase in damping between EC and EO states as we found that the distribution of damping became more heavily damped in EO states than in EC when alpha blocking was present. Our results therefore suggest an economical approach to describing the resting EEG spectrum in terms of the summed activity of stochastically driven alpha band damped linear oscillations and represents an alternative way to characterising the spectral properties of resting EEG.
Our findings suggest both alpha blocking and '1/f ' noise in resting EEG, in contrast to the dominant view, might be accounted in terms of a single underlying mechanism. This is significant as the mechanistic genesis for the alpha rhythm and its EC-EO induced changes remains uncertain. The received view for the origin of the alpha rhythm is that it arises principally via reverberant thalmo-cortical neuronal population interactions [8]. In contrast the predominant explanation for alterations in alpha band power, in response to various tasks and changes in behavioural state, are changes in the dynamical synchronization of microscopic neural population activity [15]. Using a model that is theoretically founded in mean-field model approaches, which assumes the EEG signal is generated by the aggregate activity of multiple coupled excitatory and inhibitory cortical neural populations, we have shown that it is possible to account for alpha blocking in terms of changes in the underlying damping of stochastically driven alpha band damped linear oscillations. Our conclusion regarding the role that damping plays in regulating alpha band activity is consonant with the conclusions arrived at using time series modelling methods [7], methods quite different to the frequency domain method employed here. The detailed approach discussed here inherently incorporates a description of '1/f ' noise in EEG by including a distribution in the dampings of the alpha oscillatory processes. On this basis viewing the resting EEG as being composed of rhythmic processes occurring on a background of arrhythmic activity, suggested to be scale-free in nature and indicative of self-organised criticality [30], is not necessarily required. The heterogeneity of spectral scaling seen in resting state EEG (see Fig 5) can be readily accounted for without the need to invoke the presence of fractal temporal dynamics or other non-linear processes.
From a practical standpoint, our approach follows naturally from what is known about the dynamical activity of scalp based EEG. Dynamically the emergence of alpha band oscillatory activity in the EEG can arise either through limit cycle oscillatory activity (non-linear mechanisms) or as a fluctuation spectrum in response to the random forcing/driving of a near equilibrium (linear) system [1,38]. Therefore, given that non-linear time series analyses of recorded EEG have shown that resting EEG cannot be reliably distinguished from a pseudorandom linear process [43][44][45], the latter mechanism is to be preferred as the explanatory basis for the dynamical features of resting EEG. Significantly, our approach does not result in the differentiation of rhythmic and arrhythmic EEG components, as is the goal of a number of other approaches and methods aimed at characterising the dynamical structure of resting EEG [36,37]. Rhythm implies repetition, or as is often meant in the context of neuronal oscillations, a limit cycle (non-linear) oscillation. A randomly forced (or driven) linear system, of the type we consider, is by definition not rhythmic-the absence of any determinism means phase is not defined and thus it cannot be said to be rhythmic even though it contains oscillatory activity. Our approach implies to first approximation that the resting EEG spectrum can be accounted for in terms of randomly driven alpha relaxation oscillatory activity-with is both arrhythmic and oscillatory.
The most surprising result we obtained was the predominant bimodal structure in the distribution of dampings. Interpreting this from a model fit perspective though is in hindsight fairly simple i.e. two modes are required to achieve a superior fit to both the alpha peak and the '1/f ' tail. From a neurophysiological perspective however, interpretation of these two persistent distribution modes is more challenging. One interpretation would be that large changes in damping could be indicative of a more responsive neural system. We found that in instances where a large degree of alpha blocking was present, the mode separation was largest in the EC state which was subsequently reduced in EO. In some cases, distributions went from bimodal to unimodal in EC and EO states for particularly strong alpha blockers. This could mean that participants who are subject to a large difference in damping between states, are more efficient at corralling the underlying neural population activity into task engagement and then subsequently dispersing it when no longer needed. When the eyes are closed the brain is in a more restful idling state where a larger proportion of the cortical neural populations are weakly damped. Upon opening of the eyes, most neural population activity is suppressed in response to incoming stimuli while a small undamped proportion is engaged in the task, as reflected in systematic large scale increases in the damping. Closing the eyes once more engages neuronal population activity and sets them to idling as the demand for cognitive resources reduces. Such a view lends support to the notion that highly reactive alpha band activity is correlated with better memory engagement and learning/intelligence [56,57,57,58]. Indeed, research has shown that alpha power and peak frequency respond in a variety of ways to mental effort and memory retention and/or retrieval [12]. While this interpretation is highly speculative, it does provide a clear motivation for further explorations of the functional significance of our approach. For example we could assume that our damping distribution can be well represented as a sum (or mixture) of one or more parametric distributions with our goal being the optimal estimation of the corresponding parameters. In contrast we assumed no prior functional form for our damping distributions.

Comparison with existing theoretical accounts
Our results suggest that the resting alpha rhythm corresponds to a state of increased excitability in that increases in alpha band power imply a reduction in damping, and thus a state closer to the imaginary axis. Such a conclusion receives independent corroboration in that time series modelling of resting EEG, using physiologically specified fixed-order auto-regressive movingaverage models, clearly indicates that eyes-closed (EC) alpha band activity is more weakly damped than eyes-open (EO) alpha band activity [7,16]. Indeed, such a result is consistent with the theoretical results of a range of cortical and thalamo-cortical mean field models [11,59,60], in which the alpha rhythm arises as a consequence of a stable linearised system (having an alpha resonance) being driven by 'extra-cortical' white noise. Our results therefore seem to be at odds with recent, experimentally supported, hypotheses in which the alpha rhythm is conceived as reflecting a mechanism of functional inhibition [61][62][63][64]. Numerous studies seem to indicate that strong ongoing alpha oscillatory activity is related to a state of low excitation, as inferred by reductions in measures of neuronal activity, that include, at the microscopic level, neuronal firing [65][66][67][68][69] and, at the meso-/macroscopic level, broadband high-frequency activity (BHA; [61]). BHA, typically characterized as the power between 70-150 Hz (also known as 'high gamma' and a key analytic signal recordable using the bandwidth available through invasive intracranial recordings), is thought to be a reliable measure of local neuronal excitation [70,71].
However, dynamically speaking, excitation is not the same as excitability, and as such reductions in excitation (e.g., mean population neuronal firing rate) will not necessarily reflect reductions in excitability (as characterized by increases in damping). For example, fitting a well-known neuronal population model of the alpha rhythm to empirical EC and EO EEG spectra (Hartoyo et al 2020 in particular the datasets and code available at https://github.com/ cds-swinburne/Hartoyo-et-al-2020-DATA-n-CODE) reveals that, consistent with the numerous experimental studies (see above), high amplitude EC alpha activity (which is theoretically conceived as weakly damped stochastically perturbed alpha oscillatory activity) is associated with lower mean excitatory and inhibitory neuronal firing rates than the lower amplitude EO alpha band activity, which is more heavily damped and thus less excitable. In other words, low excitation is associated with high excitability, whereas high excitation is associated with low excitability. Parametrically it is found that changes in excitatory input (subcortically and from other cortical areas) to inhibitory cortical neurons (p ei ) most sensitively modulates alpha band power, with increases in p ei hypothesised to be the principal driver of reductions in alpha band power between EC and EO states. Thus, increases in afferent cortical excitatory input are predicted to attenuate alpha band activity (through reductions in excitability) while being associated with increases in neuronal population excitation, and thus changes in neuronal population activity are driving changes in alpha band activity rather than the other way around. On this basis alpha can be functionally viewed as an inhibited rhythm rather than an inhibitory rhythm. From this physiologically plausible theoretical perspective, our inferred changes in the excitability of alpha between EC and EO states can be unproblematically reconciled with the corresponding extensive experimental literature.

Physiological plausibility of the decomposition
The usefulness of any parametrisation or decomposition of the resting EEG ultimately depends upon its connection with theory and thus its physiological plausibility. Physiological plausibility ensures that these parametrisations/decompositions can be used as genuine tools of discovery and not simply arbitrary lenses by which to view and partition seemingly complex phenomena. For example, EEG decompositions that separate the resting power spectrum into 'periodic' and 'aperiodic' components (e.g. FOOOF- [37], IRASA- [36], BOSC- [72]) appear to be predicated upon the idea that well defined neuronal level circuits are ontogenetically/ phylogenetically configured to produce well defined 'periodic' (oscillatory) activity on a background of 'aperiodic' (or 'scale-free') activity thought to reflect the complex self-organized temporal dynamics of neural activity at many different physical and temporal scales [20]. How physiologically plausible then is our decomposition?
Even structurally relatively simple systems, such as a metal plate or a drum skin, when vibrated at different frequencies, will produce a range of characteristic spatio-temporal patterns or eigenmodes, such that any arbitrary perturbation can be decomposed into a weighted sum of these eigenmodes. Indeed, such a decomposition can be applied to more complex systems having a range of local and non-local interactions, by defining an eigen-decomposition involving a structural graph Laplacian defined on the basis of these local and non-local interactions [73][74][75]. Recently this and similar approaches [76] have been applied to investigating the macroscopic emergence of spatial patterns of oscillatory activity in linearised neural field EEG models constrained by the estimated cortico-cortical connectivity of the human brain [73,74,[77][78][79][80]. However, neuronal population connectivity is not restricted to these long-range excitatory connections but will also include the many 'intra-columnar' feedforward and feedback connections between inhibitory and excitatory neurons that traverse distances of the order of micrometres and millimetres [81]. By applying the previously developed graph theoretic approaches to spatially more fined grained linearised neural field models, we could reasonably expect the emergence of eigenmodes having much smaller characteristic physical scales. On this basis, and ignoring volume conduction, the EEG recorded at a single electrode would represent the superposition of these fine grained spatio-temporal modes. Our hypothesis would be that these modes are spatially distinct, of relatively uniform frequency (alpha) and be parameterizable by their rate of decay following excitation. At this stage it is an open question under what conditions neural field models [40,81,82], as currently configured, are able to support such activity. Empirically the determination of such 'intra-columnar' modes could be achieved using methods such as dynamic mode decomposition [83,84] applied to microelectrode array [85] or two photon calcium imaging recordings [86] of neuronal activity, though we are unaware of any implementation of this or similar methods to the analysis of such data. Further, we speculate that during a perceptual act the majority of these modes would be suppressed (for example through increases in p ei as discussed previously), with the small number remaining unsuppressed potentially available to be excited to autonomous oscillation, through some form of dynamical bifurcation ('cell assembly ignition' e.g., [87], [88]). Such a bifurcation could be expected to result in the production of harmonics of alpha, and thus contribute to the beta band power that currently we do not explicitly account for (as discussed further in 'Model limitations and extensions'). Viewed from this perspective alpha is hypothesized to be a functionally selectively inhibited rhythm rather than a functionally inhibitory rhythm.

Model limitations and extensions
While the simplicity of our model is one of the features that makes it useful, there exists some clear limitations that could not be currently be addressed within this work. These include the low frequency scaling of the power spectrum and the appearance of significant beta band activity (13-30 Hz) manifesting as a clear peak in the power spectrum.
Resting EEG is often characterised to exhibit two frequency scaling regions, one at lower frequencies ranging over delta (0-4 Hz) and theta (4-7 Hz) band activity, and one at higher frequencies, with a transitional region between the two occurring somewhere in the alpha band [16,20], that are speculated to be functionally distinct [16,35,89]. Given such putative functional independence, and in the interests of simplicity, we neglected the low frequency scaling and chose to investigate a broadband region starting in the alpha band. Future work needs to extend our model to include either a single low frequency relaxation oscillation or a distribution of such processes. This inclusion would then posit that the resting EEG is almost entirely explained by the joint activity of low frequency and alpha band processes that have a distribution of relaxation rates. In support of this idea Chaudhuri et. al. (2018) [35] demonstrated that broadband ECoG power spectrum could be modeled by the sum of two Lorentzian processes, suggesting that the underlying neural activity had two distinct time courses with dominant low frequency and high frequency components.
Our exclusion of beta band activity could be rectified in two ways. The simplest would be to simply filter out the beta band activity and attempt to maintain the scaling present in the spectrum. This would correct any errors introduced when modelling the spectrum when a large beta peak is present, like those seen in the EC state which resulted in a weaker correlation between experiment and model. However, beta band activity is typically highly correlated with the alpha band, with some hypothesizing an explicit harmonic relationship [62,90]. On this basis it would be better to modify the underlying model assumptions to incorporate beta band activity. However, it is not immediately clear how this would be achieved and any such alteration would undoubtedly complicate the model by suggesting, possibly non-linear, interactions between damped alpha and beta band processes.
Median model fit was found to be slightly, though significantly, better in the EC compared to the EO state. This difference arose from the automatic identification of peak alpha frequency (PAF) in the respective EEG spectra. Our current method identifies the PAF by fitting a single peaked 'Lorentzian' function to experimental power spectra (see Methods section). This approach will inevitably poorly resolve the PAF when the alpha peak is heavily attenuated. Because EO state alpha band activity is usually heavily attenuated, often to the point where no discernable peak is present in the power spectral density, there is increased uncertainty in determining the PAF and thus in the associated model fit. In future work we will address this by employing the use of fixed order autoregressive-moving-average time series models, which has been empirically demonstrated to provide successful estimation of the PAF in EO states [7] where the alpha peak is heavily attenuated. Further contributing to differences in fit quality we assumed that EC-EO changes were solely driven by changes in damping. However, systematic changes in the PAF between EC and EO states are often observed [12,16]. In the current study the PAF, estimated separately for each state, was passed to the regularisation scheme as a free parameter, and thus any systematic dependency between f α and damping in EC and EO states was ignored. A general solution would be to extend the model to include a distribution of frequencies in addition to dampings, resulting in a two dimensional probability density function over both alpha frequency and damping. However, such inclusions complicate the problem substantially. Nevertheless, preliminary analyses, not reported here, indicate that it may be possible to solve such a problem using gradient descent methods.
One final limitation worth mentioning is volume conduction, which is a persistent problem in EEG recordings [1]. Generally speaking the principal effect of volume conduction will be to reduce the spatial resolution of recorded EEG activity, and thus we can reasonably expect that the observed topographic variations of spectral power, spectral scaling and damping have been blurred and thus underestimated. However as we see clear topographical heterogeneity in our estimates of damping and given we are principally interested in differences between EC and EO states there is no real reason to believe that our empirically estimated damping distributions have been qualitatively altered by volume conduction effects.

Conclusion
We have demonstrated that alpha blocking and '1/f ' noise in resting EEG can be accounted for by a simple singular mechanism consisting of damped alpha band oscillatory activity. Through the use of inverse methods, the distribution of damping was typically revealed to be bimodal in both EC and EO resting states, with the degree of alpha blocking measured between EC and EO shown to be driven by systematic increases in the damping of a weakly damped mode. The topographic distribution of these changes in damping paralleled well the corresponding topographic changes in alpha blocking. The power-law scaling of the generated spectral fits were shown to be well correlated with the corresponding experimental power spectra being able to reproduce the 1/f β tail. While our assumption that the EEG signal is composed of the bulk collective activity of many uncorrelated stochastic alpha band damped linear oscillatory processes is simplistic it does provide a novel way of interpreting the spectral behaviour of resting EEG.

EEG as a sum of damped relaxation oscillatory processes
The fundamental assumption of our model is that the resting EEG signal (measured at the sensor level) is generated by a vast number of uncorrelated stochastically perturbed alpha band damped linear oscillations that arise from multiple neural populations, which have a distribution of dampings (relaxation rates), the form of which is unknown to us. Our assumed model is theoretically motivated by a general model for generating 1/f noise [18] and by mean-field approaches to modelling the resting EEG [38].
We begin by noting that most linearised mean-field models generally account for the EEG signal, eeg(x, t) in a quantitatively similar manner i.e. a linear time invariant transfer function driven by broad-band noise eegðx; tÞ ¼ ffi ffi ffi k p g e ðx; tÞ � nðx; tÞ ð1Þ where g e (x, t) is the (excitatory) electrocortical impulse response [81,91] for cortical location x, n(x, t) is the driving broad-band noise process, � is the convolution operator and ffi ffi ffi k p is a constant taking into account the (bio)physical processes associated with EEG recordings. In the frequency domain this becomes where the conjugate variables are EEG (x, ω), G e (x, ω) and N(x, ω). Stable linearisation requires an impulse response that has a decaying envelope. Therefore, we choose the simplest parametric electrocortical impulse response function as a damped cosinusoid assuming it has an instantaneous rise time in response to perturbation (rise time is much shorter than the characteristic decay rate). On this basis the functional form is where the spatial dependency on x has been removed for clarity, ω α = 2πf α is the oscillating frequency with f α the (parametric) centre alpha frequency, γ is the parametric damping and Θ(t) is the Heaviside step function (it is entirely expected that γ and f α will depend on a set of population model biophysical parameters [γ � γ(p), f α � f α (p)], however the current parameterisation is ignorant to this). From Eq (2), the power spectral density is then By assuming the electrocortical impulse response of Eq (3), |G e (ω)| 2 can be calculated as: where � is the complex conjugate. Rewriting in terms of f α , Eq (6) (given that there is uncertainty about the exact form of g e (t) and our choice is merely the simplest, we choose the numerically simpler reduction of Eq (6) instead of Eq (5)) is then We now assume that the signal measured at a single location is the sum of many such processes arising from multiple neural populations, where the centre alpha frequency, f α is assumed constant and remains physiologically fixed, while the damping γ varies: The estimated power spectral density for this model process, S( f ), is then, We now assume • the white noise driving processes N i (f) are all uncorrelated such that their cross-spectrum vanishes, • the N i (f) are broadband and of equal power (RMS) i.e. N i (f) ! N (i.e. is a constant), reducing our estimated model power spectral density to We now make one final assumption that the total number of damped linear oscillators, M, is sufficiently large that we can define a distribution of dampings, p(γ), such that the number of oscillators with damping between γ and γ + dγ, is Mp(γ)dγ, where R p(γ)dγ = 1. On this basis Eq (12) can be rewritten in the continuum limit as where c = kNM and S e (f; f α , γ) � |G e (f; f α , γ i )| 2 . For our purposes S(f) will be empirically estimated directly from EEG data and S e (f; f α , γ) is our theoretically specified spectral response (Eq 7). Substituting this into our equation we arrive at our final model description which is a Fredholm integral equation of the first kind [92]. In practice, the peak alpha frequency is selected as described in the Methods section and the upper (γ h ) and lower (γ l ) limits of the damping distribution are chosen such that a broad range of spectral behaviour is covered i.e. both large and sharp alpha peaks and those which are heavily damped.

Numerical solution of the model Fredholm integral equation using Tikhonov regularization
Fredholm integral equations of the first kind are inherently ill-posed and difficult to solve [93].
For this reason to obtain useful results, regularization methods are required. We employ the use of Tikhonov regularization, a well known and studied method for solving ill-posed problems [93]. The process used to estimate the damping distributions using Tikhonov regularization is described as follows. Beginning with our model formulation, Eq (14), we define our Lorentzian-shaped frequency domain electrocortical transfer function as Our general Fredholm integral equation of the first kind is then where we note that we have incorporated our constant c into Eq (15). Making one further simplification by setting k i (γ) = k(f i , γ), we then discretize the integral Finally, we can arrange the above summation as a set of linear equations in matrix notation where S is an N × 1 data vector (power spectrum), K is an N × M coefficient matrix with entries consisting of the kernel evaluated at each damping value and P is the unknown distribution vector M × 1 in dimensions. Simply inverting the equation where Γ = λI is the Tikhonov Matrix, which is generally chosen to be identity matrix, I, multiplied by a regularization parameter, λ and k�k 2 is the Euclidean (or L 2 ) norm. In practice, we solve Eq (20) by setting it up as a constrained general least squares minimization problem of the form and use MATLAB's lsqnonneg algorithm. With all problems requiring regularization the difficulty in obtaining an appropriate solution is further compounded by the determining the 'optimal' regularization parameter. The degree to which a solution is regularized is ultimately a trade off between over or underfitting. When the magnitude of the regularization parameter is large, the determined solutions will be overly smooth. When the regularization is small, the solution will still be noisy and susceptible to the effects of numerical noise. Methods exists that allow an optimal regularization parameter to be determined, but these are largely heuristic with optimal defined as the value which satisfies the constraints of the experiment in question [93]. For the current study, we explore a regularization parameter range of 10 −5 � λ � 10 0 divided into 100 logarithmically spaced steps, with the optimal regularization parameter chosen such that it finds the maximum entropy solution that falls within an error bound of approximately 2.5% of the minimum residual sum of squares error (RRS) between the generated power spectrum and the data power spectrum. The maximum entropy distribution within that error bound was chosen as it favours distributions that are smoother and broader. See S1 Appendix for a graphical example of the optimal regularization parameter selection process.

Simulated tests
A simulation study was conducted to validate the Tikhonov regularisation method. Damping distributions of known form that were uni/bi/tri-modal with differing relative mode areas were used in the forward problem to generate simulated power spectra. The simulated power spectra were then used as the input to the regularised inversion method of the previous section to recover the known distributions. The simulated damping distributions were Gaussian in formp with r p the ratio between the peak magnitudes, μ γ the mean and σ γ the standard deviation. For bi/trimodal structures, the damping distributions were produced by summing multiple Gaussians. The relevant model parameters tested are presented below: • Unimodal: Damping interval [0 � γ � 10] s −1 , μ p = 5 s −1 , σ = 5 s −1 .
• A successful recovery of a known distribution was measured using a RSS cost function between the known distribution and the recovered distribution. A regularisation parameter range of [10 − 6 � λ � 1] was used and divided into 100 logarithmically spaced steps.

EEG data
Two freely available EEG data sets collected from two independent experiments were chosen because both explicitly included EC and EO resting condition recordings. The first set (subsequently referred to as Dataset 1) was made available by [94] and consisted of data from 30 participants recorded across multiple drug/task conditions during EC and EO resting state across 64 channels using a modified 10-20 layout. We used the data of 27 participants in the EC and EO pre-drug conditions. Each EEG recording was approximately 2 minutes in duration and was recorded at 500 Hz. Three participant recordings were not used due to significant artefact. The second set (subsequently referred to as Dataset 2) is an open source EEG data made available by [95] (https://physionet.org/content/eegmmidb/) and includes 109 participants who had their EEG recorded, using the BCI2000 system, across 14 experimental tasks. We make use of the 2 minute EC and EO recordings that were collected from 64 electrodes using a modified 10-20 layout at a sampling frequency of 160 Hz. We chose not to group data sets together for analysis as each dataset used slightly different channel layouts.

Power spectral analysis
All spectral analysis was performed at the sensor level. Power spectra for the EEG was calculated by computing the Welch periodogram across all channels using the measured time series with 2 second Hamming windows with 50% overlap. Power spectra broadband activity between 7 Hz and 40 Hz is used for the EEG data as it encompasses alpha band activity and the high frequency scaling which are our main focus. The choice of the specific upper and lower frequency limits that were used was motivated by two main reasons: • The majority (95%) of the spectral power in resting EEG lies below 25 Hz [96]. Therefore, extending the choice of upper limit out to arbitrarily higher frequencies has diminishing returns. Additionally, the 40 Hz upper limit enabled us to avoid us to avoid the 50/60 Hz mains line noise profile.
• Empirical characterisation of the resting EEG has shown the presence of at least two distinct scaling regions; one low frequency and one high frequency with a 'knee' frequency generally occurring in the alpha band [16,35,89,[97][98][99][100][101][102]. The specific intervals in each scaling regime vary between studies. However, the consistent feature is that the boundary between these two scaling regions is located within the alpha band and that the two are regarded as arising from dynamically distinct processes. Therefore, we exclude the low frequency portion of the spectrum as it is not mechanistically relevant to the genesis of alpha band oscillatory activity.
The peak alpha frequency, f α , for each experimental spectrum was calculated by fitting a single function of the form 1/[(f − f α ) 2 + b 2 ] within the alpha band using MATLAB curvefit toolbox functions. Because a comparison between the experimental and model spectral slope for frequencies is explored, a spectral scaling parameter β was calculated by fitting a 1/f β profile to the power spectrum over a frequency interval from 15-40 Hz. This procedure was performed on all the experimental spectra and their respective model spectra.

Quantifying alpha blocking using the Jensen-Shannon divergence
The degree of alpha blocking present between EC and EO conditions varies between participants: some presenting with significant alpha attenuation and others showing minimal blocking. To explore and characterize the presence of alpha blocking within the EEG data we require a measure that can quantify the degree to which an individual attenuates alpha band activity between EC and EO states. On this basis we employ the same methodology used by [41], where the Jensen-Shannon divergence (JSD), D JS , between the respective EC and EO power spectra is used to quantify the level of alpha blocking. The JSD is a symmetric and nonnegative measure of the distance between two probability distributions. By normalizing the power spectra to unity, we can treat each as a probability distribution allowing us to compute the distance between the EC and EO spectra in probability space. The JSD provides a useful way to measure the degree of alpha blocking between resting states. The expectation is that the larger the attenuation of alpha band activity between EC and EO states, the larger the JSD will be between the two distributions. The JSD is based upon the Kullback-Liebler divergence [103], D KL , which is defined as follows; where p(x) and q(x) are two probability distributions. The JSD for two spectra is then Weakly damped distribution measure To explore changes in the damping distributions between EC and EO states requires a measure that can quantify how 'weakly damped' a given distribution is. An obvious choice would be the peak position of the distribution. However, in instances when multiple distribution peaks are present, this simple method does not accurately characterise the dominant decay rate. Furthermore, given that the respective weighting of each damping plays an essential role in the outcome of the spectral profile, the peak position alone is not sufficient to quantify which distribution is more weakly damped when an overlap in the distribution peak may occur across conditions. Therefore, we constructed a measure that takes into account both the peak position and the respective area that lies under each mode in each recording condition with respect to the other state. The measure can be computed for any given pair of distributions and assigns a specific value to each resting state. A distribution (i.e p(γ) in Eq 13) is considered more weakly damped if it contains a larger area that falls to the left side of the other distributions peak. This measure is applied to the distribution peak or first mode if the distribution is multimodal. In practice, calculating the weakly damped measures (WDM) for two distributions requires calculating two specific values for each resting state damping distribution, achieved as follows; where A1 P1 is the area in distribution 1 that falls before the peak of distribution 1, A1 P2 is the area of distribution 1 that falls before the peak of distribution 2, with A2 P2 and A2 P1 being similarly defined. The weakly damped measures are then the sum of the two quantities for each state. The two values can then be compared and that which is larger is deemed more weakly damped i.e. WDM 1 > WDM 2 implies distribution 1 is more weakly damped than distribution 2. The inclusion of A1 P1 and A2 P2 terms are to account for situations comparing bimodal distributions where the first mode in one condition is much smaller than in the other (specifically with respect to weighting/area of the first mode) but happens to be positioned to the left of the other distributions first mode peak. In these instances, without the inclusion of these additional terms, the distribution with the smaller mode would have a weakly damped measure that was non-zero and the other distribution, which may have a much larger first distribution mode, would have a measure that is equal to 0, due to the way the A1 P2 and A1 P2 are defined. A graphical example of this and a further discussion is provided in S2 Appendix.

Characterisation of numerically calculated damping distributions
For numerically estimated damping distributions the following features were quantified: • Number of modes-the number of peaks, calculated using MATLABs peakfind function, in the respective EC and EO damping distributions.
• Mode positions-the peak position of each mode.
• Difference in damping between EC and EO states-calculated subject and electrode-wise as the difference between the respective weakly damped measures.
These quantities were calculated for both resting states across all subjects and channels providing a total of 8677 damping distributions [(109 Subs × 64 Channels) + (27 Subs × 63 Channels)] for EO and EC conditions. Given the differing channel layouts of the respective data sets separate group level topographic plots were calculated with any summary statistics being computed separately across the data sets.

Permutation testing
We employed permutation methods to test whether there exists statistically significant differences between EC and EO resting state data. We use EEG channel based data averaged across participants and compute t-values from the difference in means between the two states. For each EEG sensor a null hypothesis, H 0,. . .,i is assumed in which there is no significant difference between EC and EO resting states. Given that multiple hypothesis tests are being conducted in parallel, the permutation test was corrected for the multiple comparison problem. Here we use the Max-t/Min-t method detailed in [104] which corrects for the inflated Type I error rate introduced by multiple comparisons. For our experimental permutation testing we permuted the data 5000 times. For the sake of clarity we provide a brief description of the Max-t/Min-t algorithm: Algorithm 1: Max-t/Min-t Multiple Comparison Correction Permutation Test