Delineating between-subject heterogeneity in alpha networks with Spatio-Spectral Eigenmodes.

Between subject variability in the spatial and spectral structure of oscillatory networks can be highly informative but poses a considerable analytic challenge. Here, we describe a data-driven modal decomposition of a multivariate autoregressive model that simultaneously identifies oscillations by their peak frequency, damping time and network structure. We use this decomposition to define a set of Spatio-Spectral Eigenmodes (SSEs) providing a parsimonious description of oscillatory networks. We show that the multivariate system transfer function can be rewritten in these modal coordinates, and that the full transfer function is a linear superposition of all modes in the decomposition. The modal transfer function is a linear summation and therefore allows for single oscillatory signals to be isolated and analysed in terms of their spectral content, spatial distribution and network structure. We validate the method on simulated data and explore the structure of whole brain oscillatory networks in eyes-open resting state MEG data from the Human Connectome Project. We are able to show a wide between participant variability in peak frequency and network structure of alpha oscillations and show a distinction between occipital 'high-frequency alpha' and parietal 'low-frequency alpha'. The frequency difference between occipital and parietal alpha components is present within individual participants but is partially masked by larger between subject variability; a 10Hz oscillation may represent the high-frequency occipital component in one participant and the low-frequency parietal component in another. This rich characterisation of individual neural phenotypes has the potential to enhance analyses into the relationship between neural dynamics and a person's behavioural, cognitive or clinical state.

Introduction

The wide variety of oscillatory phenomena in electrophysiological recordings of brain function reflect the synchronised activity of underlying neuronal networks (Bressler, 1995, Fries, 2005). These oscillatory signatures have a rich frequency spectrum that shows meaningful between subject variability across cortex and between participants. In practice, this variability is often simplified in either the spatial or spectral domain by the choice of a priori frequency bands or regions of interest. There is a need for analytic approaches that simultaneously characterise both the spectral and spatial features of oscillatory signals within a single framework. Here, we present a data driven modal decomposition analysis approach which identifies oscillations in multivariate time-series and characterises their peak frequency, damping time, spatial topography and network structure. We apply this method to explore the macro-structure of alpha oscillations in Human neocortex and how this varies between subjects.

The alpha oscillation is often characterised as a 7–13Hz signal originating from occipital cortex (Berger, 1929, Lopes da Silva, 1991) whose function has been associated with a wide range of cognitive and clinical states (Jensen, Mazaheri, 2010, Klimesch, Sauseng, Hanslmayr, 2007). Yet there is strong and growing evidence that alpha oscillations are not homogeneous across different frequencies, brain regions or individual participants. The lower and higher edges of the 7–13Hz alpha band have distinct task responses indicating that they relate to different aspects of cognition (Klimesch, 1999, Klimesch, Sauseng, Hanslmayr, 2007). Individual Alpha Frequency (IAF) is variable across populations (Klimesch, 1999) and modulated by task demands within individuals (Haegens et al., 2014). Moreover, IAF may be a valuable clinical marker; the slowing of alpha peak frequency is a robust characteristic of both Alzheimer’s Disease and Mild Cognitive Impairment (Engels, Hillebrand, van der Flier, Stam, Scheltens, van Straaten, 2016, Garcés, Vicente, Wibral, Pineda-Pardo, López, Aurtenetxe, Marcos, de Andrés, Yus, Sancho, Maestú, Fernández, 2013, Hughes, Henson, Pereda, Bruña, López-Sanz, Quinn, Woolrich, Nobre, Rowe, Maestú, the BioFIND Working Group, 2019, López-Sanz, Bruña, Garcés, Camara, Serrano, Rodríguez-Rojo, Delgado, Montenegro, López-Higes, Yus, Maestú, 2016, Osipova, Ahveninen, Jensen, Ylikoski, Pekkonen, 2005, Peraza, Cromarty, Kobeleva, Firbank, Killen, Graziadio, Thomas, O’Brien, Taylor, 2018, Poza, Hornero, Abásolo, Fernández, García, 2007). Alpha power tends to peak in the midline occipito-parietal and occipital cortex in source reconstructions of resting state EEG and MEG recordings (Ciulla, Takeda, Endo, 1999, Hari, 1997). Yet, variants of the alpha rhythm are observed throughout the human brain in a range of contexts, likely arising from a combination of thalamo-cortical and cortico-cortical interactions (Foxe, Snyder, 2011, Hughes, Crunelli, 2005). Alpha band network connectivity can be widely variable between individual subjects (Wens et al., 2014); a part of this network variability is heritable (Colclough et al., 2017) and likely to reflect biologically relevant between subject heterogeneity. There is increasing evidence supporting functionally distinct alpha generators in different brain regions (Klimesch, Sauseng, Hanslmayr, 2007, Sokoliuk, Mayhew, Aquino, Wilson, Brookes, Francis, Hanslmayr, Mullinger, 2019) within individual participants. For instance, distinctions have been shown between occipito-parietal and occipito-temporal alpha (Barzegaran et al., 2017), as well as alpha arising from visual and parietal sources (Sokoliuk et al., 2019).

These lines of evidence emphasise the functional relevance of spatial and spectral variability in neuronal oscillations whilst illustrating the difficulty of untangling the many sources of within and between subject variability. It remains a substantial analytic challenge to characterise variability in both the peak frequency, spatial distribution and network structure of neuronal rhythms. We present a novel, data-driven approach which characterises both the spatial and spectral structure of an oscillatory network. We use the modal decomposition (Neumaier, Schneider, 2001, von Storch, Bürger, Schnur, von Storch, 1995) of a multivariate autoregressive (MVAR) model to define a set of Spatio-Spectral Eigenmodes (SSEs). Each mode contains a unimodal (single peak) frequency response whose dynamical importance is represented by a damping time; rapidly damped modes will be quickly extinguished and contribute less to the observed dynamics in the data. Both the peak oscillatory frequency and spatial representation are simultaneously estimated within each SSE, without needing to impose arbitrary a priori frequency bands or spatial regions of interest. Each SSE is a property of the whole system with a contribution to the whole network and whole power spectrum. Here, we show the system transfer function can be computed from a linear superposition of these SSEs. The standard autoregressive transfer function is exactly equivalent to the modal transfer function when all modes are included. Crucially, the SSE approach allows individual oscillatory signals to be isolated from the superposition. This property allow us to objectively detect oscillations using a non-parametric permutation scheme to identify SSEs from within the superposition. Once identified, the spatial and spectral parameters of the SSEs can be analysed at the group or individual subject level.

In this paper, we apply the Spatio-Spectral Eigenmodes to detect and characterise oscillations in Human MEG data, exploring macro-scale spatial and spectral variability in oscillatory resting state networks both within and between subjects. The method is first validated with simulations before being applied to resting state MEG data from the Human Connectome Project (Larson-Prior, Oostenveld, Della Penna, Michalareas, Prior, Babajani-Feremi, Schoffelen, Marzetti, de Pasquale, Di Pompeo, Stout, Woolrich, Luo, Bucholz, Fries, Pizzella, Romani, Corbetta, Snyder, 2013, Van Essen, Smith, Barch, Behrens, Yacoub, Ugurbil, 2013). Source time-courses are estimated from the pre-processed sensor data using a LCMV-beamformer (Van Veen et al., 1997) before voxels are combined within regions of the Automated Anatomical Labelling (AAL) atlas (Rolls et al., 2015). This source-parcellated data is described with an MVAR before the modal decomposition is used to describe the oscillatory features in the data. We estimate the source power distribution from the MVAR parameters and identify dynamically important modes based on an objective non-parametric permutation scheme. We describe the principal components of spatial variability across the dynamically important SSEs revealing large-scale patterns in network structure across the whole brain. Importantly, though the individual variability in peak frequency is of a similar magnitude to the spatial variability, this method is able to show that the frequency difference within different gradients are largely consistent within individuals, even though the overall IAF is highly variable between individuals.

Methods

Spatio-spectral eigenmodes from multivariate autoregressive models

Here, we give an overview of a standard approach to Multivariate Autoregressive (MVAR) modelling and spectrum estimation before outlining the modal decomposition and definition of Spatio-Spectral Eigenmodes. An illustrative summary of the analysis pipeline used is given in Fig. 3

Fig. 3

A graphical outline of the analysis procedures used in this paper. Upper: Summary of the procedures used to calculate the MVAR model and analyse the results using a Fourier (upper section) or Modal (lower section) approach. The Fourier-integration approach is used in Fig. 4 and the Modal decomposition is explored in example participants in 5 and at the group level in 6Lower: Outline of the procedures used to take the modal decomposition of the MVAR model and compute spatial principal components each of which can explain variability in different frequencies within and between-participants. The group results of the PCA analysis are presented in Fig. 7 and a summary of the results of individual participants in Fig. 8.

Spectrum estimation from multivariate autoregressive models

We start with a vector time series, x(t) with channels , . Time-lagged linear dependencies within and between the channels can be characterised with an MVAR model of order .where is an array of regression parameters at lag and is an -variate white noise process. This is a form of linear time-invariant (LTI) system in which future values of are predicted from a linearly weighed combination of its past values. The parameter matrix contains these linear dependencies between the past and future values of the time-series at a given lag, . The off-diagonal elements of describe the degree to which the different channels within the system contain lagged interactions. (without subscript) denotes the 3-dimensional parameter matrix containing for all fitted values of (from 1 to ).

The interactions described by may be expressed in the frequency domain by computing the system transfer function as a function of frequency . The transfer function describes the ratio of the input to a system to the output of the system and is computed from the -transform of the matrix.whereand denotes the sampling interval, the normalised frequency in radians and the frequency in Hertz. Eq. (2) can be evaluated for any value of in the complex plane. Here, we only evaluate on the unit circle (where ) as the output of these points can be directly related to an oscillatory frequency and Eq. (2) is equivalent to the discrete time Fourier transform.

The power spectrum of can be computed from the frequency transform of the autoregressive parameters via the transfer function and the residual covariance matrix .

Where denotes the complex conjugate transform of . contains the power spectrum of each of the regions in the diagonal and the cross-spectrum between regions in the off-diagonal terms. The properties of , and form the basis of a range of spectral connectivity metrics including Magnitude Squared Coherence, Geweke-Granger causality, Directed Transfer Function and Partial Directed Coherence (Baccalá, Sameshima, 2001, Blinowska, 2011). In practise, analyses will typically compute these metrics across a range of frequencies before integrating between specified frequency bands to isolate frequency specific structure.

MVAR Modal decomposition

An autoregressive model is a form of Infinite Impulse Response (IIR) filter whose spectral characteristics are completely described by the polynomial roots of its parameters. These roots directly relate to resonances in and describe how the filter extracts an input at frequency to obtain the filter output. This is well characterised for univariate systems and can be generalised to multivariate systems to provide an intuitive description of the frequency information contained in an MVAR model. This modal representation of the transfer function can then be used to simultaneously explore the peak frequency and spatial structure of brain networks (Neumaier and Schneider, 2001). The modal decomposition of MVAR coefficients is closely related to linear filter theory.

To perform the modal decomposition, we first rewrite the order- matrix as an order 1 system in a square block form. The autoregressive model in Eq. (1) can be restructured into a blocked form using a delay embedding of and the companion form of the MVAR parameter matrix (Neumaier and Schneider, 2001).

is a blocked matrix with the sparse rows at the bottom shifting the corresponding rows in down to create space for the in the prediction. The simplified matrix form of this equation is of almost identical form to an order 1 autoregressive model in the standard formulation in Eq. (1). The eigendecomposition of the square parameter matrix then yields , and as the eigenvalues, right eigenvectors and left eigenvectors, respectively. The eigenvalues are the roots of the characteristic equation of the matrix and as such directly define the frequency response of the pole. The characteristic frequency of each pole can be calculated as:

Oscillations are represented by complex conjugate pairs of poles within whilst single poles lying on the real line represent non-oscillatory parts of the signal. The damping time of a mode is also computed from its eigenvalue:

This describes the rate at which the amplitude of an oscillation would drop to zero if the system were energised with an impulse response. Longer damping times indicate less damped modes which will oscillate for longer durations following a single input. Short damping times indicate that the behaviour of the mode is quickly extinguished once the system is energised.

The complex valued eigenvector matrices and are the same size as . They have a specific Vandermonde structure in which a row contains blocks of values raised to successive powers of their corresponding eigenvalue .

Due to the repeating structure in rows of , we reduce analysis of the eigenvector of a mode to a vector of the first values in each row ( or ). This reduced vector (also called mode shape) describes the structure of the resonance across the dimensions of the input.

The modal form of the transfer function

When using autoregressive models for spectrum analysis, the transfer function is typically estimated from the Fourier transform of the time-domain parameters (Eq. (2)). Here, we show that it may equivalently be computed from the modes from the eigenvalue decomposition. The eigenvalues and eigenvectors defined above form the parameters of a partial fraction expansion of the transfer function. This converts the transfer function from a ratio of two long polynomials to the sum across a set of fractions with simple denominators. A modal form of the transfer function can then be defined as a summation of a ratio of the properties of the modes.where is the modal eigenvalue and is the mode residue matrix. The Fourier form comes from substituting Eq. (2) into Eq. (5). A full derivation of the modal form is included in appendix Appendix B and details on estimating the modal parameters is presented in appendix Appendix C. In the modal form, the mode residue is the coefficient of each term in the expansion (and distinct from the residuals of the autoregressive model fit) computed from the outer product of the first terms in the left and right eigenvectors where denotes the complex conjugate. is then an matrix whose elements are coefficients denoting the strength of the mode at each node and connection in the system. In other words, it acts to project the oscillation defined by in the signal within each node and the connections between them. When all modes are included in the summation, the Fourier and Modal forms of are exactly equivalent. The modal form is related to Gilbert’s Realisation (Gilbert, 1963, Kailath, 1980) which expresses a rational transfer function as a partial fraction expansion.

This modal form of has several benefits. Firstly, the relation to the other modal parameters provides important context to . Though we can evaluate at any frequency up to the Nyquist limit, the resolution of the power spectrum is limited by the number of modes. A decomposition with a higher model order will have more modes in its decomposition and therefore a richer spectral structure. Secondly, as the modal form is a linear superposition (or summation) across modes, the contribution of a single resonance can be easily isolated or removed from altogether. The computation of reduced transfer functions provides a convenient way to summarise network state from a subset of modes. Selection of modes by peak frequency can be useful as an alternative to integrating across the spectrum within specified frequency bands. In cases where a spectral peak lies close to the edge of a specified band, mode selection will allow the full contribution of that mode to enter the average without cropping its width to fit the band. The mode selection scheme can be tuned to fit the priorities of the research question at hand.

Spatio-Spectral Eigenmodes

We define a Spatio-Spectral Eigenmode by the resonant frequency, damping time and transfer function of a single component in the modal decomposition of an MVAR model.derived from a given eigenmode in the eigendecomposition above. The transfer function of an SSE is evaluated only at the peak frequency of the mode in question. We will often split the total set of SSEs to explore the properties of a subset defined by permutation, frequency range or both.

Software

All simulation, MVAR modelling and model decomposition steps are computed in Python 3.7.3 the Spectral Analysis In Linear Systems toolbox (Quinn and Hymers, 2020), https://vcs.ynic.york.ac.uk/analysis/sails and https://sails.readthedocs.io). Underlying python dependencies are numpy (Harris et al., 2020) and scipy (SciPy 1.0 Contributors et al., 2020) for computation and matplotlib (Hunter, 2007) for visualisation. MEG data pre-processing and beamforming was performed using Fieldtrip and the OHBA Software Library (https://github.com/OHBA-analysis/osl-core) in Matlab version R2019a on a cluster of x86-64 systems. The Bayesian statistical analysis was carried out in R version 3.5.2 using the BRMS (version 2.11.1; Bürkner (2018); Colclough et al. (2015)) and loo (version 2.2.0; Vehtari et al. (2017)) packages. Full scripts for the preprocessing, data analysis and statistical assessment in simulated and HCP MEG data are available online (https://vcs.ynic.york.ac.uk/analysis/rs-mvar). The scripts include a tool for checking out the correct versions of the external toolboxes which are used.

Simulations

The MVAR Modal Decomposition is first explored with simulations. 20 realisations (representative of 20 participants) of 300 s of data from a 10 node network are generated with a sampling frequency of 128 Hz. Each dataset is built from two subnetworks with different spatial and spectral profiles, the first is defined by a real-valued pole at 0 Hz and the second by a complex-conjugate pair of poles between 8 and 12 Hz, jittered across participants. Oscillatory data for these networks are generated by placing the poles within the z-plane and transforming them back to their polynomial form. These polynomials are then used as coefficients to filter white noise to produce oscillatory time-series. Each of the two oscillations are then projected through the network using a set of weights defining the relative strength of the oscillation in each of the 10 nodes. Finally the two oscillatory networks are added together with white noise to create the final signal.

Each network is described with an order 5 MVAR model. The Fourier-based cross and power spectral density (CPSD) is computed and the averaged within two frequency bands of interest 0–4 Hz and 8–12 Hz reflecting the two simulated oscillations. The modal decomposition is then computed and a modal form of split into three reduced models, two models for the poles which survive the permutation thresholding procedure and a residual model. The poles-of-interest for the simulation are taken as those which are identified as surviving the thresholding procedure. The surviving poles are then assigned to the low or high frequency band of interest based on having a characteristic frequency lying within 4Hz of the relevant frequency (the same ranges as used for the Fourier analysis).

An supplemental simulation with an additional oscillatory mode at 42 Hz was run to explore whether the method is able to describe oscillatory systems with a wide range of peak frequency values. The analysis parameters are the same as in the main text simulation with the exception that the model order in increased to 7 to accommodate the additional oscillation. A variant of this simulation was run with a variable level of noise to establish that the SSEs are able to handle low SNR signals. Both the standard Fourier power spectrum and modal parameters were robust to increasing noise levels, with all three modes in the simulation resolvable with noise levels up to five times the standard-deviation of the signal. Details are presented in supplemental section Appendix E. A detailed analysis was then run on the simulation at a moderate noise level. The results show that the method is straightforwardly able to resolve both the frequency content and network structure of the third mode. Details are presented in supplemental section Appendix F.

Ethics statement

Resting-state MEG datasets recorded from 79 participants in the Human Connectome Project (http://www.humanconnectome.org) (Larson-Prior, Oostenveld, Della Penna, Michalareas, Prior, Babajani-Feremi, Schoffelen, Marzetti, de Pasquale, Di Pompeo, Stout, Woolrich, Luo, Bucholz, Fries, Pizzella, Romani, Corbetta, Snyder, 2013, Van Essen, Smith, Barch, Behrens, Yacoub, Ugurbil, 2013) are analysed in this manuscript. Participant recruitment and data collection were carried out by Washington University and the University of Minnesota. All participants provided written informed consent prior to data collection (Van Essen et al., 2013). The experimental procedures were approved by the Institutional Review Board (IRB) at Washington University (IRB number 201204036; “Mapping the Human Connectome: Structure, Function, and Heritability”). For the analysis in this study, the preprocessed dataset was downloaded from ConnectomeDB (https://db.humanconnectome.org).

Resting state MEG data

Each participant underwent three separate runs of a 6 min eyes-open resting state protocol MEG data were collected using a 4D Neuroimaging WH-3600 scanner, equipped with 248 magnetometer sensor channels and 23 reference channels, and were sampled at 2034.51 Hz. Participant headshapes were digitised using a polhemus tracker system prior to MEG data collection.

The HCP pre-processed resting state MEG datasets were used along with the room noise recordings for the relevant session and information regarding the ICA components from the de-noising process. Co-registrations for the MEG and MRI data for each participant were taken from the models provided by the HCP project.

Seventy-eight areas from the AAL2 atlas (Rolls, Joliot, Tzourio-Mazoyer, 2015, Tzourio-Mazoyer, Landeau, Papathanassiou, Crivello, Etard, Delcroix, Mazoyer, Joliot, 2002) were used as target regions of interest. Beamformer weights were calculated for locations on an 8mm-spaced grid spaced inside each of the regions of interest. A Linearly-Constrained Minimum Variance beamformer (Van Veen, Van Drongelen, Yuchtman, Suzuki, 1997, Woolrich, Hunt, Groves, Barnes, 2011) at the orientation which showed maximum power. The source virtual electrode time-series were then resampled to 240Hz. The individual time-series from the grid locations within each region were then reduced to a single time-series per region by taking the first principal component across the voxels within the region. The time-series across all regions were then orthogonalised to reduce the impact of spatial leakage (Colclough et al., 2015). The main analyses were explored with and without this orthogonalisation step on a subset of the HCP dataset (a single run from each participant) in supplemental section Appendix M. This additional analysis indicates that the orthogonalisation has the greatest impact on the relatively small, highly damped modes. The larger oscillatory modes contain relatively consistent information whether orthogonalisation is applied or not. Finally, the beamformed time-series were downsampled 2-fold using a windowed fourier-domain method, giving a final sampling rate of 120 Hz.

Model order selection

Prior to analysis of an autoregressive model of any dataset, the model order must be selected. This choice can be informed by metrics such as Akaike’s Information Criterion (AIC: Akaike (1974), however this often gives a monotonically decreasing profile with no clear optimal model (Brovelli, Ding, Ledberg, Chen, Nakamura, Bressler, 2004, Ding, Bressler, Yang, Liang, 2000, Jansen, 1991, Schlögl, Supp, 2006). Even if there is a local minimum in the AIC time-domain metric, this does not guarantee that the resulting power spectrum will provide a good representation of the data. In addition, the choice of sampling frequency can have an impact on autoregressive power spectra (Quirk and Bede Liu, 1983). The effect of changing these parameters is relatively predictable and is analogous to the choice of window length and sample rate in a more conventional measure such as Welch’s Periodogram. For example, changing model order increases the number of modes that the MVAR model can represent, with lower model orders having fewer modes and therefore smoother spectra. In contrast, changing sample rate affects the Nyquist frequency of the power spectrum and therefore the frequency range that the available modes represent.

To demonstrate these effects and to directly assess which combination of model order and sampling frequency to use in the main analysis, an exploration of the power spectra and modal parameters was carried out on a subset of the HCP data (a single run per participant). The results are presented in detail in supplemental section Appendix H. In short, the we find that model orders below 8 and sampling rates above 120Hz produce results which do not well represent the approximately 1–40 Hz frequency range of interest in this paper. The final analysis uses a model order of 12 and a sampling rate of 120Hz producing a clear physiological power spectrum whose oscillatory modes are evenly distributed across frequency. This choice has the consequence of removing relatively high frequency oscillations (greater than 50 Hz) from the analysis. A different parametrisation would be more sensitive to these frequencies at the small expense of resolution at lower frequencies. Due to the prevalence of alpha power in resting state scans, we decided to optimise analysis for the relatively low frequency bands. Critically, the power spectrum and modal decomposition are relatively robust to moderate changes to both model order and sampling rate around our chosen values.

Model fitting and validation

MVAR Models were fitted with order 12 across all 78 parcels in the HCP data. 12 was chosen by a combination of the AIC and manual inspection of the model spectra. An order of 12 produced good spectra and was not before an inflection point in the AIC.

After fitting, the models were checked for stability (using the Stability Index (SI): (Lütkepohl, 2007, pages 15–16)), residual autocorrelation (using the Durbin-Watson index: Durbin and Watson (1950)) and variance explained. The models were able to fit between 21 and 29% of variance (mean=24.990%, SD=4.003%) within each recording session. We consider this to be a good proportion of variance to explain with a single stationary and linear model of a whole brain functional parcellation. All models were stable, having SI values below 1 (mean=.959, SD=.016) and no substantial autocorrelation could be found in the residuals according to the Durbin-Watson test (mean=2.003, SD.001).

Once the MVAR models ( matrix) were fitted for each scan session, the transfer function and spectral matrix were computed between 0 and 48Hz using the Fourier method. The matrices were averaged within the set of specified frequency bands to summarise the frequency-specific spatial topologies captured by the MVAR models.

Fourier and SSE network connectivity estimation

We validate that a single MVAR model is able to describe the spatial and spectral content of a whole brain functional connectome estimated from MEG data using a standard Fourier-based approach. The system transfer function is estimated using the Fourier Eq. (2) for all frequencies between 0-60Hz in 100 steps. Subsequently the spectral matrix is computed for each frequency using the and the residual covariance matrix . Finally, we integrate within a set of pre-specified frequency bands to summarise how the network structure of oscillatory brain networks changes across frequency.

Modal decomposition and non-parametric permutation

The modal decomposition of each MVAR model was computed using the methods described above and the peak frequency, damping time, and were computed for each mode. The modal decomposition of a system returns modes which could number in the hundreds or thousand for a large system. Many of these modes are likely to be modelling noisy characteristics of the system or its measurement rather than physiologically interesting oscillatory activity. In order to select the most dynamically relevant modes, a non-parametric permutation testing method was used on the damping times of the modes. Each individual time series was split into non-overlapping temporal epochs resulting in a 3d data array [channels x samples x epochs]. Permutations are carried out by randomising both the channels and epochs in order to construct null datasets in which the relationships between nodes have been destroyed whilst maintaining the overall spectral nature of the data. At each permutation, a MVAR model is estimated on the surrogate dataset and the modal decomposition computed. A maximum statistic method was then used (Nichols and Holmes, 2002) in which the maximum damping time of all of the modes within the given model was entered into the null distribution. This was repeated for each permutation, resulting in a null distribution of damping times for each participant, for each run. A threshold which represented the 1% tail of the null distribution was then established in this way for each run, for each participant. Modes were then selected from the un-permuted data using these individual damping time thresholds.

Spatial principal components analysis of SSE networks

Patterns of spatial and network variation in the SSE surviving the permutation scheme was performed using a principal components analysis (PCA). The [nnodes x nnodes] PSD matrices for the significant SSEs were vectorised and concatenated a [modes x nnodes*nnodes] matrix and demeaned before a PCA was used to identify the principal axes of variation across the connections within the network across modes. The components of the PCA then show patterns in the spatial distribution of oscillatory power across a number of modes regardless of the characteristic frequency of the modes which significantly contribute to the component. Whilst each network is computed at its peak resonant frequency, these resonances are free to vary (within the specified alpha range) across networks both between and within individuals.

The PCA was computed for subsets of SSE whose peak frequency lies within each of three frequency bands. Theta (1–7Hz), alpha (7–13Hz) and beta (13–30Hz). Crucially, the inclusion of an SSE in a band depends only on its peak frequency. Once included, all information on the network structure within that SSE is included in the analysis, even if part of the spectral peak goes outside the specified band. Reproducibility of the components arising from the PCA were assessed using a split-half correlation. 500 split halves of the SSEs included in a PCA were computed and the PCA computed on each half independently before the spatial components of each half are then correlated and stored. Both the split-half correlation and proportion of variance explained by each component was used in determining whether the component was included in further analyses. The components describe the pattern of variability across space captured by that PC whilst the PC-score indicates the extent to which that shape is expressed in each individual SSE. An example spatial map is computed for maximum and minimum observed score in each PC by projecting that score back into the original data-space.

Relatively few SSEs survived the permutation scheme in the theta and beta bands. Potentially as a result, the PCA components from these bands also showed relatively low reproducibility. For completeness, the SSE and the first 4 PCs for these bands are included in section Appendix K. In contrast, over 1200 SSEs were included in the alpha band and the first two PCA components showed high variance explained and split-half reliability. These are interpreted in the main text and carried forward for further statistical analysis.

Relationship between mode frequency and PC projection score

In order to examine whether there was a relationship between the frequency of each mode and the score with which it projected onto a given component, we performed a Bayesian linear regression using the BRMS package (Bürkner, 2017, Bürkner, 2018). For each PC, we scaled the scores by its standard deviation and fit a model of , allowing an overall change in mean frequency per-participant. Model inference was performed using the standard NUTS sampler used by STAN through BRMS.

The prior for the Frequency parameter was chosen to be normally distributed with a mean of 0 and a standard deviation of 1; reflecting our default position that there was no a-prior reason to expect frequency to vary with score. Altering the standard deviation of the prior to other plausible ranges had no significant effect on the overall results. Examination of diagnostic plots showed that the parameters have converged in all cases.

To assess whether there was a relationship for each principal component, we fit an intercept-only model and a model with frequency as an additional linear regressor and compared evidence for the models using a Leave-One-Out (LOO) cross-validation methods (Vehtari et al., 2017). Our criteria for determining that a model with the additional frequency regressor has more evidence is that the difference in LOO should be more than twice the estimate of the LOO standard standard error. For models where the model with frequency was assessed as having more evidence, we then report and assess the magnitude of the frequency parameter in the full model along with is 95% Credible Interval (CI). Due to the scaling of the scores, the frequency parameters is expressed in terms of the standard deviation of the score.

Results

Validation in simulated data

To illustrate the MVAR modal decomposition and Spatio-Spectral Eigenmodes we explored a single simulated dataset and a group simulation designed to exhibit realistic inter-run or inter-participant variability. The simulation scheme is summarised in Fig. 1A and described in detail in Section 2.3. Briefly, the simulated activity in this network contained two resonances with pre-specified spatial and spectral structures. Two modes with distinct spectral structures were defined by direct pole placement and used to generate time-courses which were projected into a 10 node network structure. 20 realisations of 300 s in duration were simulated from this network structure. Independent realisations of white noise were added to the simulated data prior to MVAR and SSE analyses.

Fig. 1

The simulations and modal decomposition for a single realisation of the simulated data. A: Summary of the simulation. ten nodes are generated from linear combinations of two modes with different spectra. The modes are shown in blue and red with the nodes in black, horizontal bars indicate the weighting of the two modes into each of the 10 nodes. B: Summary of the true network structure of the two modes. Darker colours indicate higher power in each node or connection. C: The power spectral density from ten nodes. These capture the gradual drop in frequency from 0Hz and a peak around 9Hz split across the different nodes. D: The modal power spectrum for the node highlighted in black in C. The gradual slope from 0Hz and the 9Hz peak are clearly isolated as distinct resonances. The remaining modes have very low amplitudes with no clear peaks. E:-plane representation of the modal eigenvalues (shown as crosses). The frequency of each mode relates to its angle as it increases counter-clockwise from (1,0) to (-1,0). Negative frequencies correspond to angles increasing clockwise from (1,0) to (-1,0). The two high-power modes from D are clearly visible as the modes with the largest magnitude (red and blue). The remaining modes have small magnitudes and have evenly distributed angles (black crosses). F: Damping-time of each mode as a function of frequency. The red and blue modes have significantly longer damping times than the null distribution (99% threshold shown as the dotted line). These relate to 0Hz and 9Hz resonances in the data. The remaining modes have short damping times indicating that the influence of these modes is very short-lived. G: Modal PSD matrix computed from the eigenvectors associated with the blue mode (via the transfer function). Darker colours indicate higher power in each node or connection. H: Modal PSD matrix computed from the eigenvectors associated with the red mode (via the transfer function). Darker colours indicate higher power in each node or connection.

Modal decomposition of a single session

An example segment of simulated data with its generating modes is shown in Fig. 1A. A non-oscillatory (blue) and an oscillatory (red) source time-course is created and projected across a network to create 10 node time courses (black). The red and blue horizontal bars indicate the weighting of each model time-course into each of the 10 nodes. The “true” network matrix containing the structure of each mode is shown in Fig. 1B. The time-series were described with an order-5 MVAR model fitted across the whole 300 s simulation. The Fourier based transfer function and power spectra were computed from this model and the spectrum of each node is shown in Fig. 1C. These spectra show the contributions from the two modes across the ten nodes. Some nodes contain signal from mode one (e.g. node 1), mode two (e.g. node 4) or both modes (e.g. node 7). Next, the modal decomposition was computed and the PSD con for node 7 (shown in black in Fig. 1C) is shown in Fig. 1D. Whilst node 7 contains contributions from both modes which are mixed in the Fourier-based analysis (Fig. 1C), these are clearly split into separate peaks (blue and red) in the modal power-spectrum (Fig. 1D).

The frequency of each mode was computed directly from the eigenvalue of the fitted MVAR model. A -plane plot of the eigenvalues of the decomposition (Fig. 1E) reveals that the -type mode is represented by a single real-valued mode (blue cross), in contrast the 9Hz mode is modelled by a complex conjugate pair of modes (red crosses). In the z-plane, frequency is represented by the angle of the complex eigenvalue, whilst the magnitude of the mode is its distance from the origin. As a more intuitive alternative, we show a scatter plot with individual modes with peak frequency on the x-axis and damping time on the y-axis (Fig. 1F). The damping time indicates how quickly an oscillation in that mode would be extinguished, longer damping times indicate that energy in the oscillation will dissipate more slowly. The damping time plots emphasise the dynamically important modes with long damping times whilst the frequency can be directly read out from the x-axis. In addition, the calculated damping time threshold for the simulation run is shown as a dotted line, demonstrating that the two relevant modes are easily separable from the background (Supplementary figure D.10 contains this data for all 20 realisations). Finally, the network structure of each of the two modes can be reconstructed from their modal transfer functions, computed from the relevant eigenvectors. Fig. 1G and H shows the Modal PSD of modes 1 and 2 respectively, each evaluated at its peak frequency (as determined from the respective eigenvalue). These reproduce the ground-truth structure shown in Fig. 1B.

Modal decomposition of group-level networks

Next we examined how the Fourier band-integration and SSE approaches can describe oscillations with between subject variability in peak frequency. We computed 20 realisations (representative of 20 participants) of the simulation in Fig. 2 with varying peak frequency and amplitude in peak 2 whilst keeping the network structure itself static. Fig. 2A shows the spectra of node 7 across the realisations of the simulation. The alpha peak frequency has a uniform +/-2Hz variability across realisations (gray lines represent individual subject); the group average can be seen as the solid black line. As in the single case, node 7 contains a contribution from both oscillatory networks; showing a type slope and a peak at around 9Hz. The simulated variance in peak frequency can also be seen clearly in the Fourier spectra shown in Fig. 2B; the frequencies-of-interest are highlighted in red and blue. The Fourier spectra captures the average features well, but the use of pre-determined frequency bands leads to clipping at the edges of some peaks. In addition, we can see contributions from the 0Hz peak influencing the shape and magnitude of the PSD around the 10Hz oscillation. Whilst adapting the frequency band of interest to the individual peak frequency could reduce the effect of peak clipping in the Fourier integration approach, it is harder to reduce interference between oscillations with overlapping spectra. As an alternative, the modal PSDs are shown in Fig. 2C. These are computed from the reduced transfer functions using poles selected by their damping time and driving frequency. In contrast to Fourier integration, this approach extracts single-peaks which vary depending on specific frequency content of the data.

Fig. 2

Power spectra and networks from the group simulation. A: The PSD of node 7 for all 20 realisations with variable spectra (gray) and the group average (black). B: Fourier-based PSD of node 7. The PSDs are split into pre-specified ‘low’ and ‘high’ frequency bands (in blue and red respectively). Though these capture the features around each frequency, they do not account for either individual variance in peak frequency or overlap between adjacent frequencies. C: Modal-based PSD of node 7. The modal spectra identified by thresholding the damping times of each mode of the modal decomposition and assigning each mode to its closest band (low in blue or high in red). The modal spectrum contains a single peak per mode and allows for variability in peak frequency between participants. D: Original network structure matrices. This figure shows the ground truth for the simulations generated in Fig. 1. Darker colours indicate higher power in each node or connection. E: Fourier-based network structure reconstruction. The network structure estimated from the Fourier-integration approach based on the bands in 2B, this captures the main structure with some interference from the adjacent frequency band. Darker colours indicate higher power in each node or connection. F: Modal-based network structure reconstruction. The network structure estimated from the modal bands seen in 2C. Here, the two spectrally distinct networks are properly resolved and there is little interference between the two. The diagonal structure in 2E is contained within the noise modes that did not survive the thresholding. Darker colours indicate higher power in each node or connection. G: The level of correlation (across the twenty realisations) between the ground-truth network structure and the network structure extracted for each individual run and both the Fourier and Modal analyses. The modal matrices show a much larger correlation with the ground truth than the Fourier-integration derived matrices. H: The correlation between the noise modes and the ground truth and Fourier network structures. The Fourier-integration matrices have a large correlation with the diagonal structure which is not explicitly associated with either of the simulated structures.

Fig. 4

Fourier-based frequency-specific networks extracted from fitted MVAR parameters. A: The Fourier Power Spectral Density averaged across participants for each node in the AAL parcellation. The model captures a clear 1/f type trend across the spectrum as well as a distinct alpha peak. The diagram below shows the colour code for each region of the AAL atlas. B: Surface plots showing the average PSD for in each cortical parcel within the theta, alpha and beta bands. C: Network matrices showing the average CSD between parcels within the theta, alpha and beta bands. Colourscales correspond to those in sub-panel B. D: Circular network plot showing the average CSD between parcels within the theta, alpha and beta bands. Region colouring is shown in the legend at the bottom of the figure: red: frontal, yellow: medial, purple: temporal, green: parietal, blue: occipital. Lighter colours refer to the left hemisphere (and are on the left of the network matrices and circular plots) and darker colours refer to the right hemisphere (and are on the right of the network matrices and circular plots).

Fig. 5

Spatial and spectral variability illustrated by data from four example participants. Each row contains the Fourier power spectrum and modal decomposition for that participant. A: The standard Fourier power spectrum, coloured lines indicate brain regions following the colour code in Fig. 4A. B: The damping times of the modal decomposition as a function of mode peak frequency. All modes are shown with dynamically important modes shown in red. C: The power spectrum averaged across all brain regions computed using the Fourier method (black) and the reduced spectrum computed from the modes surviving the permutation tests is shown in red. D: Network connectivity matrices computed from the reduced modal transfer function (corresponding to the red line in C). E: Spatial distribution of power from the reduced modal transfer function corresponding to the diagonals in D.

Fig. 6

Assessment of modal poles in individual participants. A: Plots of damping time against frequency for all poles, for all participants. Poles coloured in red survived the non-parametric individual subject thresholding and were carried through to later analyses. The left hand plot shows the normalised distributions of poles which survive (red) and do not survive (black) the thresholding procedure against damping time. B: The power spectrum of all modes (black) and modes surviving the non-parametric permutation scheme (red) across all datasets. Each region in the AAL atlas is shown as an individual line. C: Plot of the modes for all participants sorted by peak alpha frequency. Each row on the y-axis shows an individual participant and each dot in the graph an individual mode (all three runs for each individual are combined in one graph). Below threshold modes are shown as small black dots whilst modes which survived thresholding are shown as larger dots in colours. The x-axis shows the frequency of the mode and a set of canonical frequency bands are shown in gray boxes.

Fig. 7

Results of the PCA decomposition of the modal analysis for PCs 1 and 2. For each component we show: A: a surface plot of the component structure across nodes in the AAL atlas. B: histograms showing the distribution of SSE frequencies for SSEs with positive (purple) and negative (green) scores. C: scatter plot showing the relationship between SSE peak frequency and PC-score. For PC2, this also contains the regression line quantifying the modelled relationship. The frequency component of the model was not significant in PC1. D: network matrix showing the component structure across network connections. E: circle plot showing the component structure across network connections (top 15% of connections shown).

Fig. 8

Correspondence between the PC scores and peak frequency of the surviving modes for PC1 and PC2. For each component we show A: a surface plot of the component structure across nodes in the AAL atlas. B: Histograms of the SSE peak frequency for each mode split by positive (purple) and negative (green) score for the component. The left hand histogram shows the absolute mean peak frequency for each individual for the positive and negative scores. The right hand histogram shows the frequency of the components relative to the participant mean. C: Per-participant scatter plot of the SSE peak frequency for each mode split by positive (purple) and negative (green) score for the component. The left hand scatter plot shows the absolute mean peak frequency for each individual for the positive and negative scores. The right hand scatter plot shows the frequency of the components relative to the participant mean. The order of participants is sorted in the same manner as in Fig. 6 C.

The group simulation used the network structure defined in Fig. 2D, with the network connectivity pattern driven by the -type signal on the left and the simulated alpha oscillation on the right. The network structure estimated by the Fourier frequency band integration approach captures the core features of the ground-truth simulations, but show spectral ‘leakage’ between the two underlying network patterns (Fig. 2E). The high and low resonances overlap in the frequency axis leading to low frequency content leaking into the high frequency integration window and vice-versa. In contrast, this mixing is absent in the modal estimation (Fig. 2F) which is able to separate the contribution from each pole to all frequencies and tune itself to variance in individual peak frequency. Both methods achieve a high (r >0.9) correlation between true and estimated network structure across with the whole network and all realisations, however the modal networks are much more tightly clustered close to 1 (Fig. 2G). In addition, the noise network estimated from the residual modes correlates between r=.4 and r=.5 in the case of the Fourier-integration estimated networks whilst the same correlation in the modal networks is much lower (around 0.2) (Fig. 2H).

Additional supplemental simulations were able to show that the MVAR model and SSE decomposition is robust to increasing the level of noise in the simulation (See supplemental section Appendix E and to the addition of a third oscillatory mode with a distinct spatial structure (See supplemental section Appendix F). Finally, a simulation with a non-sinusoidal shows that the MVAR-SSE approach introduces harmonic components to represent non-linear waveform shapes. This limitation is shared with all Fourier based methods and can lead to ambiguity in the interpretation of modes in the presence of strong distortions in waveform shape Huang et al. (1998); Quinn et al. (2021); See supplemental section Appendix G).

Oscillatory networks in MEG data

We next explored the frequency structure of oscillatory networks in the Human Connectome Project MEG data. We demonstrate that an autoregressive model can capture the frequency specific content of a whole-brain functional network using standard Fourier band-integration before moving to explore the SSEs. A summary of the whole analysis pipeline for the HCP is given in Fig. 3. Pre-processed MEG data were source localised to a 5mm grid throughout the brain using an LCMV beamformer before groups of voxels were combined into parcels based on the 78 cortical regions in the AAL atlas. The parcel time-courses were then orthogonalised to reduce leakage (details on the MEG processing are included in Section 2.5). MVAR models were fitted to each recording session before their Power and Cross Spectral Densities were estimated using the Fourier-integration approach (Model fitting and validation is described in detail in Sections 2.6 and 2.7).

The topography of MEG functional networks vary as a function of frequency (Brookes, Woolrich, Luckhoo, Price, Hale, Stephenson, Barnes, Smith, Morris, 2011, Colclough, Brookes, Smith, Woolrich, 2015, Hipp, Hawellek, Corbetta, Siegel, Engel, 2012, Marzetti, Della Penna, Snyder, Pizzella, Nolte, de Pasquale, Romani, Corbetta, 2013, Wens, Bourguignon, Goldman, Marty, Op de Beeck, Clumeck, Mary, Peigneux, Van Bogaert, Brookes, De Tiège, 2014). As a result the spectrum is typically split into a set of independently analysed frequency bands. Our autoregressive model fits were able to capture frequency specific power distributions and network structure within these a priori defined bands. The average Fourier derived PSD (across participants) from each node of the AAL can be seen in Fig. 4A. Overall, each node shows a trend with the strongest oscillations visible in the alpha band. Frequency-specific source topographies are shown in the remaining columns of Fig. 4. These maps were computed by averaging the MVAR PSD estimates across participants within a priori frequency bands. Each panel includes source-space images containing the diagonal of the PSD matrix within the frequency band as well as a network matrix showing the off-diagonal Cross-Spectral Densities (CSD) and a circular connectivity plot showing the network connectivity based on the CSD. The circular connectivity plots show the connections whose magnitude falls within the top 5% of the off-diagonal CSD distribution for each frequency band. Details of the labelling and colouring of each cortical region can be found in the figure caption.

The 3–7Hz theta band has strongest power in medial prefrontal regions with connectivity including connections with the parietal and occipital cortex. In contrast, alpha (7–13Hz) power is strongly localised to occipital cortex with strong connections within the occipital region and between the occipital and temporal regions, with a smaller number of parietal connections. The detailed structure of the alpha network is explored by applying an eigenvalue decomposition to the network matrix. The results reproduced the main occipital structure in the first component before revealing more detailed occipito-parietal, parieto-temporal and lateralised structures (see supplemental section Appendix I for details). Beta power (13–30Hz) is predominantly seen in the bilateral motor cortices with a broad range of connections. Overall, these results suggest that MVAR model based Fourier frequency-domain power and connectivity estimates are able to represent whole-brain functional connectivity patterns in line with expectations from the literature (Hillebrand et al., 2012).

Spatio-Spectral Eigenmodes capture individual variability in oscillatory networks

The Spatio-Spectral Eigenmodes (SSEs) provide an alternate description of network power spectra based on their characteristic frequency and network structure. The modal decomposition was computed for the MVAR model of each recording session yielding SSEs (in this case an MVAR model over 78 regions with model order 12 yields 938 SSEs), though only a minority of these reflect dynamically important structure in the data. Non-parametric permutations were used to split the full set of SSEs into an included set of dynamically relevant functional modes and an excluded set of modes whose damping times are not distinguishable from chance in this dataset (see Section 2.9 for details). Fig. 5 shows the Fourier power spectrum and SSEs for four individual participants (in rows). The first participant has a strong alpha peak at around 11 Hz (Fourier power spectrum shown in column A) which is well represented by the SSEs with long damping times (scatter plot of SSE damping times by peak frequency in column B). The SSEs in red are included in further analyses having survived the non-parametric permutations. The power spectrum from full and included set of SSE (black and red lines in column C respectively) again indicate that the included SSEs do capture the prominent oscillations in the full power spectrum. Furthermore, the included SSEs have a bilateral spatial distribution and network structure around the occipital pole with a bias toward the right hemisphere (included SSE network connectivity matrix and surface plot in columns D and E). The second participant has a similar alpha peak in the power spectrum with a slightly lower peak frequency. In contrast to the first participant, this participant’s alpha is broadly distributed around bilateral occipital and parietal regions. The third participant has a small alpha peak corresponding to significant SSEs with relatively short damping times compared to participants 1 and 2. A single SSE in the beta band survives the permutation scheme and contributes to a diffuse power and network structure between occipital, parietal and motor cortex. Finally, example participant 4 shows two separate alpha peaks in two different brain regions as shown by the 8 and 10Hz peaks in the power spectrum and SSE damping time plots. The average spectrum shows a prominent, relatively low frequency alpha peak which is will described by the significant SSEs. Similar to participants 2 and 3, this participant has a relatively diffuse alpha power distribution across occipital and parietal cortex.

Group variability in alpha peak frequency

To describe the oscillatory frequency content across the whole brain and group, the damping times of all modes across the full HCP dataset are plotted as a function of peak alpha frequency in Fig. 6A. The included sets of SSEs (as identified by non-parametric permutations) are indicated in red with the remaining SSEs in black. As in the individual cases, the modes with the longest damping times occur around the strongest peaks in the Fourier spectrum (Fig. 4A). The majority of these (for the present eyes-open resting state data) lie within the traditional alpha range with a smaller number in the delta, theta and beta ranges and a single mode above 30Hz. In contrast, the excluded SSEs are relatively distributed across the whole frequency range. Fig. 6B shows the average power spectrum across all regions and participants (black) and the spectrum reconstructed from only the significant SSEs (red). The relatively small number of SSEs preserve the prominent oscillations in the signal at the group level suggesting that the permutation scheme is successful in extracting the SSEs related to the largest resonances in the system. The predominance of alpha in the surviving SSEs reflects the prominence of the alpha rhythm in power spectra across eyes-open resting state MEG scans. This is a static power spectrum estimate across the whole duration of each scan. As such it is possible that individual variability in these alpha peaks are driven by temporal dynamics as well as oscillatory amplitude, similarly there may be transient bursts in other frequencies that are not detected in the average spectrum (Quinn et al., 2019).

Across all participants the frequency distribution of SSEs provides a straightforward summary of the spectral variability in the HCP dataset. (Fig. 6C; surviving modes are shown as blue-green dots in rows) though the three runs within each participant are quite consistent. The rows of Fig. 6C are sorted by Individual Alpha Frequency (IAF; derived from an average frequency of SSEs between 6 and 14Hz, weighted by damping time) showing the variability across participants. The majority of surviving modes fall within the traditional 7–13Hz alpha band. 76/79 participants have at least one significant mode within alpha although peak frequencies are widely variable across participant (median peak frequency of alpha modes across participants range from 7.4Hz to 12.9Hz; mean equivalent from 7.4Hz to 12.9Hz). The participants with the lowest and highest IAFs lie very close to the boundaries of the standard 7–13Hz frequency range. Though the peak frequency lies within these bounds and are therefore well represented by the SSEs, the standard alpha band does not contain the full width of the oscillatory peak for these participants. As such, an analysis that imposed a strict cut-off would likely clip the edges of these alpha peaks leading to potential distortions or misrepresentation of the spectral content. A relatively small number of oscillatory modes in some participants occur within the delta, theta and beta bands.

Alpha peak frequency varies between occipital and parietal cortex

Each SSE is a property of the whole brain rather than a single region or ROI, allowing it to represent the distribution of an oscillation across space and network connections. The spatial variability in the SSEs PSD networks whose peak lies within the 7–13Hz alpha range was described by a small number of components in a Principal Components Analysis (PCA). The components of PCA analysis describe the axes of spatial and network variability across modes whilst the Principal Component (PC) scores indicate the extent to which a particular PC component is expressed within a given SSE. The reproducibility of each PC was evaluated by split-half correlations, this indicated that the first two components were highly replicable across halves of the data (the validations for the PCA analysis are described in detail in Section 2.10 and the results are shown in the Supplementary Material in section Appendix J)

Crucially, this PCA was computed on the spatial network structure of alpha SSE without knowledge of the corresponding resonant frequencies. We next quantified the correspondence between the spatial content and the peak frequencies of the SSEs across participants. A Bayesian regression was used to assess the extent to which peak frequency can be used to predict PC score across SSEs. There is large between subjects variability in alpha peak frequency (Fig. 6C) so we include each participant as a random effect term in the model. This allows us to directly quantify the between subject variability in peak frequency and explore whether there is a consistent relationship between network structure and peak frequency across the dataset even where the absolute peak frequency itself is variable. Further, to assess whether a term in the model provides good out of sample predictions, we used a Leave-One-Out (LOO) cross-validation procedure. Details of the Bayesian model inference and validation are described in Section 2.11.

This PCA analysis was repeated for SSEs lying within the theta and beta bands. SSEs within these bands reflected the expected spatial structure of theta and beta activity (details in section Appendix K). As relatively few theta and beta SSEs survived thresholding by non-parametric permutation, we did not go on to perform a detailed investigation of the their spatio-spectral covariation. The small number of theta and beta modes surviving permutation reflects the predominance of alpha oscillations in resting state recordings. The distribution of modes reflects the content of the signal, so we would expect to find more theta and beta modes in the task evoked data.

The first principal component in the alpha band (PC1: 23.0% variance explained, average split-half correlation) relates to the average power across an occipital network similar to the standard alpha network. The component values for PC1 have the same sign in all brain regions indicating that changing PC-scores will act to increase or decrease power across this whole network. In other words, SSEs with a positive score in PC1 will have a high power across this distribution, whilst SSEs with a negative score will have low power across the whole brain (Fig. 7 PC1). For the first principal component, the difference in LOO scores between the intercept-only and full models was -9.7 (SE: 5.3), indicating that there is insufficient evidence to conclude that PC-score (corresponding to overall power) is related to peak frequency in this component.

The second alpha component (PC2: 12.7% variance explained, average split-half reproducibility) contains a spatial gradient with the occipital pole at one end and parietal lobes at the other. Power at the two ends of this gradient are in counterpoint, a positive score for PC2 indicates high power in occipital lobes with suppressed power in parietal lobes and vice versa for negative scores. The Bayesian model was used to assess whether SSE peak frequency can be used to predict PC-score for PC2. In this case, the difference in LOO scores between the intercept-only and full modes was -102.6 (SE: 13.0) indicating that there is sufficient evidence to warrant assessing the full model. The frequency parameter in the full model had a central parameter estimate of 0.46, with a 95% CI of 0.40–0.52. This indicates that an increase of peak frequency of 1Hz would correspond to an increase in PC-score of around 0.46 of a standard deviation of the distribution of scores. In other words, increasing peak frequency across SSEs corresponds toward increased power in occipital cortex and decreased power in parietal cortex. The overall distribution of SSEs between the two principal components is shown in figure L.19.

Within and between subject variability in alpha frequency

The distribution of SSE damping times as a function of frequency (Fig. 5) and the random effects term in the Bayesian linear model indicate that there is very wide individual variability in alpha peak frequency. Next, we visualise how this between subject variability interacts with occipto-parietal gradient in PC2. Fig. 8 shows the distribution of SSEs with positive and negative PCs scores as a function of frequency for PCs 1 and 2. PC1 shows a wide distribution of peak frequencies between 7 and 13Hz (Fig. 8A). The distribution of frequency differences between SSEs with positive and negative scores are nearly completely overlapping (Fig. 8B).

In contrast, the distribution of PC2 frequencies shows a mean shift between positive and negative PC scores in both the absolute and relative distributions. Parietal alpha has, on average, a lower peak frequency distribution than occipital alpha. The two between subject distributions begin to overlap around 9Hz. Using the SSE methods to un-mix spatial and spectral variability we can see that parietal alpha SSEs occur between 7–11Hz across participants. The higher end of this distribution would otherwise be masked by the stronger occipital power at frequencies above 9Hz. A key source of this mixing is that variability in overall alpha peak frequency between subject is larger than the frequency difference between parietal and occipital alpha. Specifically, the overall alpha peak distribution ranges between 7 and 13Hz (range of 6Hz), though the relative difference between the two ends of occipito-parietal gradient is around 1Hz.

The spatial maps in Figs. 7 and 8 show how PC scores are distributed at the group level. We can project this mapping back to individual subjects to characterise how the gradient structure in PC2 varies across participants. Next, we show the spatial variability in four example participants with SSEs that have both positive and negative scores for PC2. The left hand column of Fig. 9 shows the average power in alpha across cortex based on PC1. We see that all four participants show the strongest power in occipital regions though there is variability in the more anterior end of the distribution towards parietal cortex. The right hand column of Fig. 9 shows the score of PC 2 across cortex for the same four participants. The higher frequency ‘occipital’ end of this PC is consistently focused around the occipital pole whilst the lower frequency ’parietal’ end is more anterior. Crucially, the gradient structure seen at the group level is visible in individual subjects but shows substantial spatial variability.

Fig. 9

Projections of PCA scores for four example datasets from the HCP MEG data. Each row shows the spatial distribution of alpha SSEs for that subject with the two columns colour coding the spatial distribution of PCA scores for PC 1 and PC 2.

Discussion

Spatio-Spectral Eigenmodes defined from the properties of an autoregressive model provide a flexible representation of oscillatory brain networks with minimal pre-specification of regions of frequencies of interest. We introduce the theory behind this approach and demonstrate its application in simulations and resting state MEG data from the Human Connectome Project. Firstly, we established that an autoregressive model is able to describe frequency specific functional networks in whole brain MEG data. The modal decomposition is then computed on these models to identify Spatio-Spectral Eigenmodes (SSEs). The resonant frequencies and damping times of these SSEs are shown to provide a simple summary of the oscillatory content of whole functional network. The spatial distribution and network structure of each SSE can then be explored through its contribution to the system transfer function and subsequently, power and cross spectral density. We utilised these properties to explore spatial and spectral variability in alpha oscillations. The SSEs expressed the between subject variation in individual alpha frequency and showed that, on average, alpha oscillations in the parietal lobe are lower frequency than those in the occipital lobe. Though this is a robust group effect, wide between subject variability means that the definition of ’low’ or ’high’ frequency are overlapping across individuals. A 10Hz oscillation could load onto the parietal lobe in one participant and the occipital lobe in another. The SSE approach separates these sources of variability by computing network structure and peak frequency simultaneously in each dataset. The SSE parameters then provide a convenient and intuitive representation of the spectral shape, spatial topography and network connectivity of neuronal oscillations.

Spatio-Spectral organisation of alpha networks

The oscillatory signatures of brain function which are measured during eyes open resting state are dominated by the alpha rhythm. The first spatial component of alpha power identified by our analysis describes variations in mean power in a network centred, on average, around medial-occipital cortex (Ciulla, Takeda, Endo, 1999, Hari, 1997). The individual resonant frequencies of these SSEs support a wide previous literature showing that IAF estimates between subjects vary widely within and around the traditional 7–13Hz range (Haegens, Cousijn, Wallis, Harrison, Nobre, 2014, Klimesch, 1999). Crucially, this variability is functionally relevant and has been linked with a wide range of cognitive and clinical markers (Clayton et al., 2018). This presents a practical problem in that the estimation of spatial maps or networks in participants whose alpha peak lies close to these bounds. If the whole width of the alpha peak is not within the specified range parts of it will be cut off, leading to possible distortions in the estimated maps and networks. One solution to this is to tune the centre frequency or width of the frequency bands to the peak of each individual subject (Haegens, Cousijn, Wallis, Harrison, Nobre, 2014, Klimesch, 1999), however this depends on the accurate quantification of the peak. Here, we show that the parameters of MVAR-derived SSEs can overcome some of these limitations by characterising individual spectral peaks without pre-filtering data into frequency bands or locations of interest.

The second PC of spatial variability in alpha SSEs shows a distinction between networks loaded onto occipital or parietal cortex. On average, SSEs with high scores towards the occipital end of this component tend to have higher frequencies than those centered at the parietal end ( 10Hz compared to  8–9Hz). This split between low and high frequency alpha is similar to the low and high bands (typically defined around 8.5 and 10Hz in the literature) are suggested to reflect separate cognitive functionality (Klimesch, 1999). Though this difference is strong on the group level, the frequency distributions of occipital and parietal SSEs are overlapping, suggesting that an oscillation of 10Hz could correspond to one participant’s low frequency parietal alpha and another participants high frequency occipital alpha. A small minority of participants in this study show the reverse effect, with higher frequency oscillations in parietal rather than occipital cortex. The difference between occipital and parietal alpha when both are present within an individual (around  1Hz) is substantially smaller than the range of alpha peak frequencies across participants ( 7–13Hz). As such, the full range of frequency variability in these regions is only visible with an analysis approach that can simultaneously deconvolve spatial and spectral variability.

The paper presents an exploratory analysis of a publicly available, eyes-open resting state dataset with the aim of characterising the structure and variability in oscillatory networks. Whilst we can show that alpha oscillations have these spectral and spatial distributions across a large dataset (from the Human Connectome Project), without an experimental task or prior hypothesis we do not make strong claims about its functional interpretation based on these analyses. To guide future research, we propose three potential interpretations of the distinction between occipital high-alpha and parietal low-alpha found in PC2. Firstly, these rhythms may reflect spatially and functionally distinct generators of alpha (Sokoliuk et al., 2019). Occipital alpha is thought to represent the locus of visual attention (Jensen and Mazaheri, 2010) whilst, parietal alpha has been linked with attentional processing and is suggested to exert top-down control of visual alpha depending on attentional state (van Dijk, Schoffelen, Oostenveld, Jensen, 2008, Sokoliuk, Mayhew, Aquino, Wilson, Brookes, Francis, Hanslmayr, Mullinger, 2019). Our results show this distinction between occipital and parietal alpha may be present in resting-state data and that these alpha sources are additionally separated by peak oscillatory frequency. A second possibility is that the second PC identified in our analysis represents a continuous gradient of oscillatory behaviour between occipital and parietal cortex. Similar gradients in structural and functional MRI data have been proposed as an organising principle of the brain (Huntenburg et al., 2018), PC2 may then represent a occipito-parietal gradient organising alpha oscillations. A related idea is that PC2 could reflect an aspect of the posterior to anterior alpha travelling waves (Zhang et al., 2018). Finally, the parietal end of PC2 may represent the sensori-motor Mu rhythm rather than a distinct parietal alpha source. The Mu rhythm peaks over sensorimotor cortex and has a similar frequency but distinct waveform shape to occipital alpha (Pineda, 2005) Future research in this area using task-related data could distinguish between these hypotheses.

MVAR models: parameterisation & limitations

The Spatio-Spectral Eigenmode decomposition method is dependant on a good estimation of the power spectra of the system via the underlying MVAR model. In turn, the estimation of the PSD is dependent on adequate selection of the hyper-parameter of the MVAR model: the model order (p) and the sample rate of the data (Quirk and Bede Liu, 1983). In the current work we downsample the source time-courses to 120Hz and use a model order of 12. This provides a good trade off between the high spectral resolution arising from high model order and straightforward model estimation from low model order (see SI Appendix H for further details). Further, as autoregressive models will always fit the entire spectrum from zero to Nyquist, the low sample rate ensures that the spectrum fit focuses on the physiological range of interest. Though these parameters, give a good fit in this instance, it is not guaranteed that they will generalise to novel datasets and appropriate diagnostics must be performed in these cases.

The modal-form of the transfer function has a spatial constraint; a single SSE is associated with a rank-1 network structure. More complex network structure is described through a combination of SSEs. This is mathematically straightforward as the transfer function can be summed across modes, yet the method for identifying which modes to combine must be tuned to the application in hand. The linear summation of modes is only equal to the full Fourier model at the level of the transfer function. Though properties such as the PSD matrix can be defined from a single SSE, the summation of these modal-PSD matrices will not necessarily equal the Fourier equivalent. Here, we explore the spatial and spectral properties of PSD matrices across many SSEs without directly summing them. Other applications may wish to combine these SSEs at the level of the transfer function for each data recording prior to group analyses. In addition, the modal cross-spectral densities used in the network analyses represent both the shared power and phase-locking between each pair of nodes. A richer representation of connectivity could be gained by using coherence or directed transfer function based metrics rather than the CSD, however the normalisation of these measures is difficult with the rank-1 matrix structure limitation in SSE analysis. We are continuing work into these issues and the wider picture of how the SSE decomposition and modal transfer function relates to standard power spectrum and connectivity measures.

In this analysis we have elected to use a permutation scheme to restrict analysis to modes with long damping times when compared to a null distribution computed from spatially and temporally shuffled data. Though this has been effective in identifying the relatively large alpha oscillations in these resting state scans, this permutation scheme might not be optimal for all analyses and we do not mean to imply that the remaining modes are functionally irrelevant. The choice of SSEs to focus on in an analysis is flexible choice and is customisable for the analysis in hand. For example, an analysis may use all SSEs within a given frequency range irrespective of damping times and spatial structure. Alternatively, the permutation scheme could be adapted to identify sets of SSEs with consistent spatial structure irrespective of peak frequency and damping time. Such mode selection schemes might reveal different aspects of the dynamics in a dataset to the alpha-optimised damping time thresholding used in this manuscript.

Relation to other decompositions

The decomposition of autoregressive models of univariate EEG time-series into their natural frequencies, damping times and transfer function contributions has a long history (Franaszczuk, Blinowska, 1985, Gersch, Yonemoto, 1977, Wright, Kydd, Sergejew, 1990). Recently, the computation of natural frequencies and damping times has been generalised to multivariate autoregressive models (Neumaier and Schneider, 2001). We link these multivariate parameters to the system transfer function via Gilbert’s Realisation (Gilbert, 1963, Kailath, 1980) leading to the definition of the Spatio-Spectral Eigenmodes.

There are several mathematically related approaches in the literature. In particular the method in this paper are closely related to techniques for modal analysis which have widespread use in engineering. Firstly, SSEs are closely related to the Principal Oscillatory Patterns and Principal Interaction Pattern analyses of autoregressive models originally developed for analyses of climate systems (von Storch et al., 1995). The peak frequencies and damping times from the eigenvalues of these analysis have previously been used to investigate EEG recorded during epileptic seizures (Mullen et al., 2012). Next, a Hankel matrix can be used to identify a state-space parameters and permits a modal decomposition to identify mode frequencies and damping times (for example the Eigensystem Realization Algorithm; Juang and Pappa (1985)). Decompositions of the Hankel matrix have been previously applied to explore the frequency modes of epileptic seizures (Hunyadi et al., 2014). Finally, the Dynamic Mode Decomposition (DMD) represents oscillatory dynamics via Koopman modes (Schmid, 2010). It is optimised for image-type datasets where there are more regions or channels than time-points in a dataset and has previously been applied to fMRI (Casorso et al., 2019) and ECoG (Brunton, Johnson, Ojemann, Kutz, 2016, Shiraishi, Kawahara, Yamashita, Fukuma, Yamamoto, Saitoh, Kishima, Yanagisawa, 2020) recordings. The application of these methods and their deeper mathematical relationship is a point of active research in the Neuroscience and the wider dynamical systems literature.

A range of conceptually related methods look to isolate oscillatory activity in electrophysiology data using linear spatial filters (see Cohen (2017) for a review). Unlike the approaches above, these typically involve computing a frequency or time-frequency spectrum across the dataset and before carrying out the decomposition, often using using PCA, ICA or related techniques. The spectrum estimation and decomposition stages may be carried out and optimised separately. These decompositions tend to be relatively unconstrained in the frequency domain, the resulting component power spectra can be comprised of complex, multi-modal shapes which may be challenging to interpret as clear oscillatory signals.

Conclusion

We have shown that a modal decomposition of MVAR parameters can be used to simultaneously estimate spatial and frequency structure within human resting state MEG data. In the SSE framework, brain networks are decomposed into oscillatory signals on an individual whole-brain basis with minimal pre-specification and averaging. Using this method, we have demonstrated that multiple, spatially overlapping, sub-networks exist within the normal alpha band activity. Detailed within-subject networks can be identified despite large between-participant variance. These structure captured by the SSEs can be used enhance investigation into how individual oscillatory phenotypes relate to individual difference in cognitive and clinical states.

Funding

This work was supported by an ESRC PhD Studentship from the White Rose Doctoral Training Centre. The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z). This work was also supported by the NIHR Oxford Health Biomedical Research Centre and the Medical Research Council grant (RG94383/RG89702).

Data & code availability

All scripts for data simulation, processing, analysis and visualisation in this paper are available online at https://vcs.ynic.york.ac.uk/analysis/rs-mvar-scripts. The HCP data used for the human MEG component of the paper is available from the Human Connectome Project.

MEG data from this study are available to download through the Human Connectome Project (HCP; www.humanconnectome.org). Prior to downloading data, users must register with the HCP and agreed to the data use terms (https://www.humanconnectome.org/study/hcp-young-adult/data-use-terms).

Symbol	Definition	Dimension
Time-Series
x	A vector time-series	m×T
x	A single channel time-series	T
m	Number of channels in a recording
T	Number of time samples in a recording
Δt	The data sampling interval
Autoregressive Models
A	Autoregressive parameters	m×m×p
p	Autoregressive model order
ϵ	A white noise process
Σ	Autoregressive residual covariance matrix	m×m
C	Autoregressive parameters A in the companion form	mp×mp
Frequency transforms
z	A z transform (Me^iω)
M	Magnitude of a point in the z-plane
ω	The complex argument of a point in the z-plane (2πfΔt)
f	A frequency in Hertz
A(f)	Frequency transformed autoregressive parameters	m×m
H(f)	A system transfer function at a given frequency	m×m
S(f)	A spectral matrix at a given frequency	m×m
Modal Decomposition
λ	Eigenvalues of C	mp
V	Right Eigenvectors of C	mp×mp
v_j	First m values of row/mode j of V	m
W	Left Eigenvectors of C	mp×mp
w_j	First m values of row/mode j of W	m
ρ	Characteristic frequency of an AR mode
δ	Damping time of an AR mode
R	Residue matrix of an AR mode	m×m