Delineating single subject oscillatory brain networks with Spatio-Spectral Eigenmodes

The spatial and spectral structure of oscillatory networks in the brain provide a readout of the underlying neuronal function. Within and between subject variability in these networks can be highly informative but also poses a considerable analytic challenge. Here, we describe a method that simultaneously estimate spectral and spatial network structure without assumptions about either feature distorting estimation of the other. This enables analyses exploring how variability in the frequency and spatial structure of oscillatory networks might vary both across the brain and across individuals. The method performs a modal decomposition of an autoregressive model to describe the oscillatory signals present within a time-series based on their peak frequency and damping time. Moreover, an alternate mathematical formulation for the system transfer function can be written in terms of these oscillatory modes; describing the spatial topography and network structure of each component. We define a set of Spatio-Spectral Eigenmodes (SSEs) from these parameters to provide a parsimonious description of oscillatory networks. Crucially, the SSEs preserve the rich between-subject variability and are constructed without pre-averaging within specified frequency bands or limiting analyses to single channels or regions. After validating the method on simulated data, we 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 variability in peak frequency and network structure of alpha oscillations and reveal 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

1 Introduction 1 e synchronised activity of neuronal populations is observable via the wide variety of 2 oscillatory phenomena in electrophysiological recordings of brain function. Such oscillations 3 are thought to reflect coordinated activity in neuronal networks carrying out specific 4 functionality within the brain [1,2]. ese oscillatory signatures have a rich frequency 5 spectrum that is often simplified, in analysis or interpretation, to a set of discrete band-limited 6 components (typically delta, theta, alpha, beta and gamma). Whilst there is a strong basis for 7 the functional, spatial and spectral separation of these components, analyses assuming a 8 strong separation into a-priori defined bands risk overlooking spatial and spectral variability 9 both within and between subjects. 10 For example, the alpha oscillation is often characterised as an 7-13Hz signal originating 11 from occipital cortex [3,4] whose function has been associated with a wide range of cognitive 12 and clinical states [5,6]. Yet there is strong and growing evidence that alpha oscillations are 13 not homogeneous in across different frequencies, brain regions or individual participants. e 14 lower and higher edges of the 7-13Hz alpha band have distinct task responses indicating that 15 they shown to relate to different aspects of cognition [5,7]. Individual Alpha Frequency (IAF) 16 is variable across populations [7] and modulated by task demands within individuals [8]. 17 Moreover, IAF may be a valuable clinical marker; the slowing of alpha peak frequency is a 18 robust characteristic of both Alzheimer's Disease and Mild Cognitive Impairment [9][10][11][12][13][14][15]. e 19 spatial origin of the alpha rhythm is typically localised to the midline occipito-parietal and 20 occipital cortex [16,17] although individual subjects' alpha networks show a wide variability 21 in topography [18]. is network variability is heritable [19] and likely to reflect biologically 22 relevant subject differences. Some of this spatial variability may arise from functionally 23 distinct alpha generators in different brain regions [5,20]. For instance, distinctions have been 24 shown between occipito-parietal and occipito-temporal alpha [21], as well as alpha arising 25 from visual and parietal sources [20]. 26 ese lines of evidence emphasise the functional relevance of spatial and spectral 27 variability in neuronal oscillations whilst illustrating the difficulty of untangling the many 28 sources of within and between subject variability. It remains a substantial analytic challenge to 29 characterise variability in both the peak frequency, spatial distribution and network structure 30 of neuronal rhythms. We present a novel, data-driven approach which characterises both the 31 spatial and spectral structure of an oscillatory network. We use the modal 32 decomposition [22,23] of a multivariate autoregressive (MVAR) model to define a set of 33 Spatio-Spectral Eigenmodes (SSEs). Each mode contains a unimodal (single peak) frequency 34 response whose dynamical importance is represented by a damping time; rapidly damped 35 modes will be quickly extinguished and contribute less to the observed dynamics in the data. 36 Here, we show that the contribution of each mode to the system transfer function can be 37 computed from the parameters of each mode. Both the peak oscillatory frequency and spatial 38 representation are simultaneously estimated within each SSE, without needing to impose 39 arbitrary a priori frequency bands or spatial regions of interest. Each SSE is a property of the 40 whole system with a contribution to the whole network and whole power spectrum. 41 In this paper, we apply the Spatio-Spectral Eigenmodes to explore macro-scale spatial and characterised with an MVAR model of order p. 68 x where A k is an m × m array of regression parameters at lag k and is an m-variate white 69 noise process. is is a form of linear time-invariant (LTI) system in which future values of x(t) 70 are predicted from a linearly weighed combination of its past values. e parameter matrix A k 71 contains these linear dependencies between the past and future values of the time-series at a 72 given lag, k. e off-diagonal elements of A k describe the degree to which the different 73 channels within the system contain lagged interactions. A (without subscript) denotes the 74 3-dimensional parameter matrix containing A k for all fi ed values of k (from 1 to p). the ratio of the input to a system to the output of the system and is computed from the 78 z-transform of the A matrix. 79 where z = e ıω ≡ cos ω + i sin ω and ω = 2πf ∆t (4) ∆t denotes the sampling interval, ω the normalised frequency in radians and f the frequency causality, Directed Transfer Function and Partial Directed Coherence [28,29]. In practise, 92 analyses will typically compute these metrics across a range of frequencies before integrating 93 between specified frequency bands to isolate frequency specific structure. can then be used to simultaneously [23] explore the peak frequency and spatial structure of 102 brain networks. e modal decomposition of MVAR coefficients is closely related to linear 103 filter theory.

104
To perform the modal decomposition, we first rewrite the order-p A matrix as an order 1 105 system in a square block form. e autoregressive model in equation 1 can be restructured into 106 a blocked form using a delay embedding of X(t) = {x(t), x(t − 1), . . . , x(t − p)} and the 107 companion form C of the MVAR parameter matrix [23].
C is a blocked mp × mp matrix with the sparse p − 1 rows at the bo om shifting the 109 corresponding rows in X(t − 1) down to create space for the x(t) in the prediction. e 110 simplified matrix form of this equation is of almost identical form to an order 1 autoregressive model in the standard formulation 112 in equation 1. e eigendecomposition of the square parameter matrix C then yields λ, V and 113 W as the eigenvalues, right eigenvectors and left eigenvectors respectively. e eigenvalues λ 114 are the roots of the characteristic equation of the matrix C and as such directly define the 115 frequency response of the pole. e characteristic frequency of each pole can be calculated as: 116 Oscillations are represented by complex conjugate pairs of poles within λ whilst single 117 poles lying on the real line represent non-oscillatory parts of the signal. e damping time of a 118 mode is also computed from its eigenvalue: is describes the rate at which the amplitude of an oscillation would drop to zero if the 120 system were energised with an impulse response. Longer damping times indicate less damped 121 modes which will oscillate for longer durations following a single input. Short damping times 122 indicate that the behaviour of the mode is quickly extinguished once the system is energised. 123 e complex valued eigenvector matrices W and V are the same mp × mp size as C. ey 124 have a specific Vandermonde structure in which a row contains p blocks of m values raised to 125 successive powers of their corresponding eigenvalue λ j .
Due to the repeating structure in rows of W , we reduce analysis of the eigenvector of a  139 where λ is the modal eigenvalue and R j is the mode residue matrix. e Fourier form comes 140 from substituting equation 2 into equation 5. In the modal form, the mode residue R is the 141 coefficient of each term in the expansion (and distinct from the residuals of the autoregressive 142 model fit) computed from the outer product of the first m terms in the left and right 143 eigenvectors R j = v j ⊗ outer w * j where * denotes the complex conjugate. R j is then an 144 m × m matrix whose elements are coefficients denoting the strength of the mode at each node 145

5/40
and connection in the system. In other words, it acts to project the oscillation defined by λ j in 146 the signal within each node and the connections between them.
derived from a given eigenmode j in the eigendecomposition above. e transfer function of 167 an SSE is evaluated only at the peak frequency of the mode in question. We will often split the 168 total set of SSEs to explore the properties of a subset defined by permutation, frequency range 169 or both. 170

171
To illustrate the MVAR modal decomposition and Spatio-Spectral Eigenmodes we explored a 172 single simulated dataset and a group simulation designed to exhibit realistic inter-run or 173 inter-participant variability. e simulation scheme is summarised in figure 1A and described 174 in detail in section 4.4. Briefly, the simulated activity in this network contained two resonances 175 with pre-specified spatial and spectral structures. Two modes with distinct spectral structures 176 were defined by direct pole placement and used to generate time-courses which were projected 177 into a 10 node network structure. 20 realisations of 300 seconds in duration were simulated 178 from this network structure.  e Fourier based transfer function and power spectra were computed from this model and the 187 spectrum of each node is shown in figure 1C. ese spectra show the contributions from the 188 two modes across the ten nodes. Some nodes contain signal from mode one (e.g. node 1), mode 189 A: Summary of the simulation. ten nodes are generated from linear combinations of two modes with different spectra. e 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. C: e power spectral density from ten nodes. ese capture the gradual drop in frequency from 0Hz and a peak around 9Hz split across the different nodes. D: e modal power spectrum for the node highlighted in black in C. e gradual slope from 0Hz and the 9Hz peak are clearly isolated as distinct resonances. e remaining modes have very low amplitudes with no clear peaks. E: z-plane representation of the modal eigenvalues (shown as crosses). e 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). e two high-power modes from D are clearly visible as the modes with the largest magnitude (red and blue). e remaining modes have small magnitudes and have evenly distributed angles (black crosses). F: Damping-time of each mode as a function of frequency. e red and blue modes have significantly longer damping times than the null distribution (99% threshold shown as the do ed line). ese relate to 0Hz and 9Hz resonances in the data. e 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). H: Modal PSD matrix computed from the eigenvectors associated with the red mode (via the transfer function).

7/40
two (e.g. node 4) or both modes (e.g. node 7). Next, the modal decomposition was computed   Figure 2A shows the 217 spectra of node 7 across the realisations of the simulation. e alpha peak frequency has a 218 uniform +/-2Hz variability across realisations (gray lines represent individual subject); the 219 group average can be seen as the solid black line. As in the single case, node 7 contains a 220 contribution from both oscillatory networks; showing a 1/f type slope and a peak at around 221 9Hz. e simulated variance in peak frequency can also be seen clearly in the Fourier spectra 222 shown in figure 2B; the frequencies-of-interest are highlighted in red and blue. e Fourier 223 spectra captures the average features well, but the use of pre-determined frequency bands 224 leads to clipping at the edges of some peaks. In addition, we can see contributions from the 225 0Hz peak influencing the shape and magnitude of the PSD around the 10Hz oscillation. Whilst 226 adapting the frequency band of interest to the individual peak frequency could reduce the 227 effect of peak clipping in the Fourier integration approach, it is harder to reduce interference 228 between oscillations with overlapping spectra. As an alternative, the modal PSDs are shown in 229 figure 2C. ese are computed from the reduced transfer functions using poles selected by 230 their damping time and driving frequency. In contrast to Fourier integration, this approach 231 extracts single-peaks which vary depending on specific frequency content of the data. A: e PSD of node 7 for all 20 realisations with variable spectra (gray) and the group average (black). B: Fourier-based PSD of node 7. e PSDs are split into pre-specified 'low' and 'high' frequency bands (in blue and red respectively). ough 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. e 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). e modal spectrum contains a single peak per mode and allows for variability in peak frequency between participants. D: Original network structure matrices. is figure shows the ground truth for the simulations generated in figure 1. E: Fourier-based network structure reconstruction. e 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. F: Modal-based network structure reconstruction. e network structure estimated from the modal bands seen in 2C. Here, the two spectrally distinct networks are properly resolved and there is li le interference between the two. e diagonal structure in 2E is contained within the noise modes that did not survive the thresholding. G: e 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. e modal matrices show a much larger correlation with the ground truth than the Fourier-integration derived matrices. H: e correlation between the noise modes and the ground truth and Fourier network structures. e Fourier-integration matrices have a large correlation with the diagonal structure which is not explicitly associated with either of the simulated structures. 9/40 frequencies and tune itself to variance in individual peak frequency. Both methods achieve a 241 high (r >.9) correlation between true and estimated network structure across with the whole 242 network and all realisations, however the modal networks are much more tightly clustered 243 close to 1 ( figure 2G). In addition, the noise network estimated from the residual modes 244 correlates between r=.4 and r=.5 in the case of the Fourier-integration estimated networks 245 whilst the same correlation in the modal networks is much lower (around .2) (figure 2H). 246 247 We next explored the frequency structure of oscillatory networks in the Human Connectome 248 Project MEG data. We demonstrate that an autoregressive model can capture the frequency 249 specific content of a whole-brain functional network using standard Fourier band-integration 250 before moving to explore the SSEs. A summary of the whole analysis pipeline for the HCP is 251 given in figure 3. Pre-processed MEG data were source localised to a 5mm grid throughout the 252 brain using an LCMV beamformer before groups of voxels were combined into parcels based 253 on the 78 cortical regions in the AAL atlas. e parcel time-courses were then orthogonalised 254 to reduce leakage (details on the MEG processing are included in section 4.5). MVAR models 255 were fi ed to each recording session before their Power and Cross Spectral Densities were 256 estimated using the Fourier-integration approach (Model fi ing and validation is described in 257 detail in section 4.6 and 4.7). 258 e topography of MEG functional networks vary as a function of frequency [32][33][34][35][36]. As a 259 result the spectrum is typically split into a set of independently analysed frequency bands. Our 260 autoregressive model fits were able to capture frequency specific power distributions and 261 network structure within these a priori defined bands. e average Fourier derived PSD (across 262 participants) from each node of the AAL can be seen in figure 4A. Overall, each node shows a   e 3-7Hz theta band has strongest power in medial prefrontal regions with connectivity 273 including connections with the parietal and occipital cortex. In contrast, alpha (7-13Hz) power 274 is strongly localised to occipital cortex with strong connections within the occipital region and 275 between the occipital and temporal regions, with a smaller number of parietal connections.

276
Beta power (13-30Hz) is predominantly seen in the bilateral motor cortices with a broad range 277 of connections. Overall, these results suggest that MVAR model based Fourier    A: e Fourier Power Spectral Density averaged across participants for each node in the AAL parcellation. e model captures a clear 1/f type trend across the spectrum as well as a distinct alpha peak. e 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. 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 bo om 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). A: e standard Fourier power spectrum, coloured lines indicate brain regions following the colour code in figure 4A. B: e 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: e 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.

13/40
from chance in this dataset (see section 4.9 for details). Figure 5 shows the Fourier power 290 spectrum and SSEs for four individual participants (in rows). e first participant has a strong 291 alpha peak at around 11Hz (Fourier power spectrum shown in column A) which is well 292 represented by the SSEs with long damping times (sca er plot of SSE damping times by peak 293 frequency in column B). e SSEs in red are included in further analyses having survived the 294 non-parametric permutations. e power spectrum from full and included set of SSE (black 295 and red lines in column C respectively) again indicate that the included SSEs do capture the 296 prominent oscillations in the full power spectrum. Furthermore, the included SSEs have a 297 bilateral spatial distribution and network structure around the occipital pole with a bias toward 298 the the right hemisphere (included SSE network connectivity matrix and surface plot in 299 columns D and E). e second participant has a similar alpha peak in the power spectrum with 300 a slightly lower peak frequency. In contrast to the first participant, this participant's alpha is 301 broadly distributed around bilateral occipital and parietal regions. e third participant has a 302 small alpha peak corresponding to significant SSEs with relatively short damping times 303 compared to participants 1 and 2. A single SSE in the beta band survives the permutation 304 scheme and contributes to a diffuse power and network structure between occipital, parietal 305 and motor cortex. Finally, example participant 4 shows two separate alpha peaks in two 306 different brain regions as shown by the 8 and 10Hz peaks in the power spectrum and SSE 307 damping time plots. e average spectrum shows a prominent, relatively low frequency alpha 308 peak which is will described by the significant SSEs. Similar to participants 2 and 3, this 309 participant has a relatively diffuse alpha power distribution across occipital and parietal cortex. 310

Group variability in alpha peak frequency 311
To describe the oscillatory frequency content across the whole brain and group, the damping 312 times of all modes across the full HCP dataset are plo ed as a function of peak alpha frequency 313 in figure 6A. e included sets of SSEs (as identified by non-parametric permutations) are 314 indicated in red with the remaining SSEs in black. As in the individual cases, the modes with 315 the longest damping times occur around the strongest peaks in the Fourier spectrum 316 ( figure 4A). e majority of these (for the present eyes-open resting state data) lie within the 317 traditional alpha range with a smaller number in the delta, theta and beta ranges and a single 318 mode above 30Hz. In contrast, the excluded SSEs are relatively distributed across the whole 319 frequency range. Figure 6B shows the average power spectrum across all regions and 320 participants (black) and the spectrum reconstructed from only the significant SSEs (red). e 321 relatively small number of SSEs preserve the prominent oscillations in the signal at the group 322 level suggesting that the permutation scheme is successful in extracting the SSEs related to the 323 largest resonances in the system. e predominance of alpha in the surviving SSEs reflects the 324 prominence of the alpha rhythm in power spectra across eyes-open resting state MEG scans. 325 is is a static power spectrum estimate across the whole duration of each scan. As such it is 326 possible that individual variability in these alpha peaks are driven by temporal dynamics as 15/40 therefore well represented by the SSEs, the standard alpha band does not contain the full width 341 of the oscillatory peak for these participants. As such, an analysis that imposed a strict cut-off 342 would likely clip the edges of these alpha peaks leading to potential distortions or 343 misrepresentation of the spectral content. A relatively small number of oscillatory modes in 344 some participants occur within the delta, theta and beta bands. Each SSE is a property of the whole brain rather than a single region or ROI, allowing it to 348 represent the distribution of an oscillation across space and network connections. e spatial 349 variability in the SSEs PSD networks whose peak lies within the 7-13Hz alpha range was Crucially, this PCA was computed on the spatial network structure of alpha SSE without 358 knowledge of the corresponding resonant frequencies. We next quantified the correspondence 359 between the spatial content and the peak frequencies of the SSEs across participants. A

360
Bayesian regression was used to assess the extent to which peak frequency can be used to 361 predict PC score across SSEs. ere is large between subjects variability in alpha peak 362 frequency (figure 6C) so we include each participant as a random effect term in the model.

363
is allows us to directly quantify the between subject variability in peak frequency and 364 explore whether there is a consistent relationship between network structure and peak 365 frequency across the dataset even where the absolute peak frequency itself is variable. Further, 366 to assess whether a term in the model provides good out of sample predictions, we used a 367 Leave-One-Out (LOO) cross-validation procedure. Details of the Bayesian model inference and 368 validation are described in section 4.11.

369
is PCA analysis was repeated for SSEs lying within the theta and beta bands. SSEs 370 within these bands reflected the expected spatial structure of theta and beta activity (details in 371 section 9.5). As relatively few theta and beta SSEs survived thresholding by non-parametric 372 permutation, we did not go on to perform a detailed investigation of the their spatio-spectral 373 covariation. e small number of theta and beta modes surviving permutation reflects the 374 predominance of alpha oscillations in resting state recordings. e distribution of modes 375 reflects the content of the signal, so we would expect to find more theta and beta modes in the 376 task evoked data.  participants. e higher end of this distribution would otherwise be masked by the stronger 414 occipital power at frequencies above 9Hz. A key source of this mixing is that variability in 415 overall alpha peak frequency between subject is larger than the frequency difference between 416 parietal and occipital alpha. Specifically, the overall alpha peak distribution ranges between 7 417 and 13Hz (range of 6Hz), though the relative difference between the two ends of 418 occipito-parietal gradient is around 1Hz. Spatio-Spectral Eigenmodes defined from the properties of an autoregressive model provide a 421 flexible representation of oscillatory brain networks with minimal prespecification of regions 422 of frequencies of interest. We introduce the theory behind this approach and demonstrate its 423 application in simulations and resting state MEG data from the Human Connectome Project. 424 Firstly, we established that an autoregressive model is able to describe frequency specific 425 functional networks in whole brain MEG data. e modal decomposition is then computed on 426 these models to identify Spatio-Spectral Eigenmodes (SSEs). e resonant frequencies and 427 damping times of these SSEs are shown to provide a simple summary of the oscillatory content 428 of whole functional network. e spatial distribution and network structure of each SSE can 429 then be explored through its contribution to the system transfer function and subsequently, 430 power and cross spectral density. We utilised these properties to explore spatial and spectral  wide previous literature showing that IAF estimates between subjects vary widely within and 446 around the traditional 7-13Hz range [7,8]. Crucially, this variability is functionally relevant 447 and has been linked with a wide range of cognitive and clinical markers [39]. is presents a 448 practical problem in that the estimation of spatial maps or networks in participants whose 449 alpha peak lies close to these bounds. If the whole width of the alpha peak is not within the 450 specified range parts of it will be cut off, leading to possible distortions in the estimated maps 451 and networks. One solution to this is to tune the centre frequency or width of the frequency 452 bands to the peak of each individual subject [7,8], however this depends on the accurate  . is split between low and high frequency alpha is 460 similar to the low and high bands (typically defined around 8.5 and 10Hz in the literature) are 461 suggested to reflect separate cognitive functionality [7]. ough this difference is strong on 462 the group level, the frequency distributions of occipital and parietal SSEs are overlapping, 463 suggesting that an oscillation of 10Hz could correspond to one participant's low frequency 464 parietal alpha and another participants high frequency occipital alpha. A small minority of 465 participants in this study show the reverse effect, with higher frequency oscillations in parietal 466 rather than occipital cortex. e difference between occipital and parietal alpha when both are 467 present within an individual (around 1Hz) is substantially smaller than the range of alpha 468 peak frequencies across participants ( 7-13Hz). As such, the full range of frequency variability 469 20/40 in these regions is only visible with an analysis approach that can simultaneously deconvolve 470 spatial and spectral variability. 471 e paper presents an exploratory analysis of a publicly available, eyes-open resting state 472 dataset with the aim of characterising the structure and variability in oscillatory networks.

473
Whilst we can show that alpha oscillations have these spectral and spatial distributions across 474 a large dataset (from the Human Connectome Project), without an experimental task or prior 475 hypothesis we do not make strong claims about its functional interpretation based on these 476 analyses. To guide future research, we propose three potential interpretations of the 477 distinction between occipital high-alpha and parietal low-alpha found in PC2. Firstly, these 478 rhythms may reflect spatially and functionally distinct generators of alpha [20]. Occipital 479 alpha is thought to represent the locus of visual a ention [40] whilst, parietal alpha has been 480 linked with a entional processing and is suggested to exert top-down control of visual alpha 481 depending on a entional state [20,41]. Our results show this distinction between occipital and 482 parietal alpha may be present in resting-state data and that these alpha sources are 483 additionally separated by peak oscillatory frequency. A second possibility is that the second 484 PC identified in our analysis represents a continuous gradient of oscillatory behaviour 485 between occipital and parietal cortex. Similar gradients in structural and functional MRI data 486 have been proposed as an organising principle of the brain [42], PC2 may then represent a 487 occipito-parietal gradient organising alpha oscillations. A related idea is that PC2 could reflect 488 an aspect of the posterior to anterior alpha travelling waves [43]. Finally, the parietal end of 489 PC2 may represent the sensori-motor Mu rhythm rather than a distinct parietal alpha source. 490 e Mu rhythm peaks over sensorimotor cortex and has a similar frequency but distinct 491 waveform shape to occipital alpha [44] Future research in this area using task-related data 492 could distinguish between these hypotheses. between the high spectral resolution arising from high model order and straightforward model 500 estimation from low model order (see SI 9.3 for further details). Further, as autoregressive 501 models will always fit the entire spectrum from zero to Nyquist, the low sample rate ensures 502 that the spectrum fit focuses on the physiological range of interest. ough these parameters, 503 give a good fit in this instance, it is not guaranteed that they will generalise to novel datasets 504 and appropriate diagnostics must be performed in these cases. 505 e modal-form of the transfer function has a spatial constraint; a single SSE is associated 506 with a rank-1 network structure. More complex network structure is described through a 507 combination of SSEs. is is mathematically straightforward as the transfer function can be 508 summed across modes, yet the method for identifying which modes to combine must be tuned 509 to the application in hand. e linear summation of modes is only equal to the full Fourier 510 model at the level of the transfer function. ough properties such as the PSD matrix can be 511 defined from a single SSE, the summation of these modal-PSD matrices will not necessarily 512 equal the Fourier equivalent. Here, we explore the spatial and spectral properties of PSD 513 matrices across many SSEs without directly summing them. Other applications may wish to 514 combine these SSEs at the level of the transfer function for each data recording prior to group 515 analyses. Finally, the modal cross-spectral densities used in the network analyses represent 516 both the shared power and phase-locking between each pair of nodes. A richer representation 517 of connectivity could be gained by using coherence or directed transfer function based metrics 518 rather than the CSD, however the normalisation of these measures is difficult with the rank-1 519 matrix structure limitation in SSE analysis. We are continuing work into these issues and the 520 21/40 wider picture of how the SSE decomposition and modal transfer function relates to standard 521 power spectrum and connectivity measures. frequencies, damping times and transfer function contributions has a long history [46][47][48]. 525 Recently, the computation of natural frequencies and damping times has been generalised to 526 multivariate autoregressive models [23]. We link these multivariate parameters to the system 527 transfer function via Gilbert's Realisation [30,31] leading to the definition of the developed for analyses of climate systems [22]. e peak frequencies and damping times from 534 the eigenvalues of these analysis have previously been used to investigate EEG recorded during 535 epileptic seizures [49]. Next, a Hankel matrix can be used to identify a state-space parameters 536 and permits a modal decomposition to identify mode frequencies and damping times (for 537 example the Eigensystem Realization Algorithm; [50]). Decompositions of the Hankel matrix 538 have been previously applied to explore the frequency modes of epileptic seizures [51]. Finally, 539 the Dynamic Mode Decomposition (DMD) represents oscillatory dynamics via Koopman 540 modes [52]. It is optimised for image-type datasets where there are more regions or channels 541 than time-points in a dataset and has previously been applied to fMRI [53] and ECoG [54, 55] 542 recordings. e application of these methods and their deeper mathematical relationship is a 543 point of active research in the Neuroscience and the wider dynamical systems literature.

544
A range of conceptually related methods look to isolate oscillatory activity in 545 electrophysiology data using linear spatial filters (see [56] for a review). Unlike the approaches 546 above, these typically involve computing a frequency or time-frequency spectrum across the 547 dataset and before carrying out the decomposition, often using using PCA, ICA or related We have shown that a modal decomposition of MVAR parameters can be used to 554 simultaneously estimate spatial and frequency structure within human resting state MEG data. 555 In the SSE framework, brain networks are decomposed into oscillatory signals on an individual 556 whole-brain basis with minimal pre-specification and averaging. Using this method, we have 557 demonstrated that multiple, spatially overlapping, sub-networks exist within the normal alpha 558 band activity. Detailed within-subject networks can be identified despite large 559 between-participant variance. ese structure captured by the SSEs can be used enhance   e transfer function describes the filtering carried out by a linear system as it transforms an 565 input into an output. In this case our input is a noise process defined by the residual 566 covariance of a fi ed autoregressive model and the output is the frequency transformed data. 567 In the main text, we limit values of z to the unit circle (where |z| = 1) as these can be Here we see that the transfer function of a system is the fraction of two polynomials 583 defining how the system transforms inputs into outputs. We can evaluate H(z) around the e value of H(z) decreases as the distance of z from each pole increases. Where there are 592 multiple poles in a system, the value of H(z) at a given point is dependant on all of the poles, 593 23/40 though as the influence of a pole falls off rapidly with increasing distance the closest poles to a 594 point z have the greatest influence. As such, the value of the transfer function at a given point 595 z is linearly dependant on the distance of z from this set of poles, providing an alternate 596 formulation of H(z) (assuming no repeat poles).

597
Here, R j is the is the coefficient of the term for the jth pole λ j in the partial fraction the roots of the denominator of H(z)) and R the outer product of the eigenvectors for a given 601 mode.
R j is a m × m matrix of rank one containing the coefficients for term j of the partial fraction 603 expansion. is scales the mode into each node and connection in the system.
we can remove the powers on z and further simplify equation 21 to Equation 24 is also known as Gilbert's realisation and it re-parametrises the system 610 transfer function as a partial fraction expansion [30]. is realisation is valid when there are 611 no repeated poles in λ. is approach has wide applications in the field of engineering controls 612 systems and a modal analysis [57]. Equation 24 shows that H(z) from a linear system may be 613 wri en entirely in terms of the roots of its time domain parameters. As such, the number of 614 degree of freedom in the frequency domain representation of a linear system is completely 615 determined by the number of roots of a. It follows that the configuration of λ in the z-plane 616 forms a natural basis for a set of oscillatory modes as derived by [23]. Moreover, each mode 617 (linked to a real valued or complex conjugate pair of pole) has a direct physical interpretation 618 as an oscillator with a peak frequency and network structure throughout the system. does not permit an eigenvalue decomposition. erefore we replace the order-p parameter 625 matrix A with its companion form C as defined in equation 7.
Once computed from the eigenvalue decomposition of C, the parameters λ, W and V can be 645 plugged into the equations in sections 2.1 and 4.1 to define the spatio-spectral eigenmodes.  (representative of 20 participants) of 300 seconds of data from a 10 node network are generated 662 with a sampling frequency of 128Hz. Each dataset is built from two subnetworks with different 663 spatial and spectral profiles, the first is defined by a real-valued pole at 0Hz and the second by 664 a complex-conjugate pair of poles between 8Hz and 12Hz, ji ered across participants.

665
Oscillatory data for these networks are generated by placing the poles within the z-plane and 666 transforming them back to their polynomial form. ese polynomials are then used as 667 coefficients to filter white noise to produce oscillatory time-series. Each of the two oscillations 668 are then projected through the network using a set of weights defining the relative strength of 669 the oscillation in each of the 10 nodes. Finally the two oscillatory networks are added together 670 with white noise to create the final signal.   Seventy-eight areas from the AAL2 atlas [27,64] were used as target regions of interest.

692
Beamformer weights were calculated for locations on an 8mm-spaced grid spaced inside each 693 of the regions of interest. A Linearly-Constrained Minimum Variance beamformer [26,65]  Once the MVAR models (A matrix) were fi ed for each scan session, the transfer function 721 H and spectral matrix S were computed between 0 and 48Hz using the Fourier method. e S 722 matrices were averaged within the set of specified frequency bands to summarise the 723 frequency-specific spatial topologies captured by the MVAR models. We validate that a single MVAR model is able to describe the spatial and spectral content of a 726 whole brain functional connectome estimated from MEG data using a standard Fourier-based 727 approach. e system transfer function is estimated using the Fourier equation 2 for all 728 frequencies between 0-60Hz in 100 steps. Subsequently the spectral matrix is computed for 729 each frequency using the H(f ) and the residual covariance matrix Σ. Finally, we integrate 730 within a set of pre-specified frequency bands to summarise how the network structure of 731 oscillatory brain networks changes across frequency. modal decomposition of a system returns m * p modes which could number in the hundreds or 736 thousand for a large system. Many of these modes are likely to be modelling noisy 737 characteristics of the system or its measurement rather than physiologically interesting 738 oscillatory activity. In order to select the most dynamically relevant modes, a non-parametric 739 permutation testing method was used on the damping times of the modes. Each individual 740 timeseries was split into non-overlapping temporal epochs resulting in a 3d data array 741 [channels x samples x epochs]. Permutations are carried out by randomising both the channels 742 and epochs in order to construct null datasets in which the relationships between nodes have 743 been destroyed whilst maintaining the overall spectral nature of the data. At each permutation, 744 a MVAR model is estimated on the surrogate dataset and the modal decomposition computed. 745 A maximum statistic method was then used [73] in which the maximum damping time of all of 746 the modes within the given model was entered into the null distribution. is was repeated for 747 each permutation, resulting in a null distribution of damping times for each participant, for 748 each run. A threshold which represented the 1% tail of the null distribution was then 749 established in this way for each run, for each participant. Modes were then selected from the 750 un-permuted data using these individual damping time thresholds.  Pa erns of spatial and network variation in the SSE surviving the permutation scheme was 753 performed using a principal components analysis (PCA). e [nnodes x nnodes] PSD matrices 754 for the significant SSEs were vectorised and concatenated a [modes x nnodes*nnodes] matrix 755 and demeaned before a PCA was used to identify the principal axes of variation across the 756 connections within the network across modes. e components of the PCA then show 757 pa erns in the spatial distribution of oscillatory power across a number of modes regardless of 758 the characteristic frequency of the modes which significantly contribute to the component.

759
Whilst each network is computed at its peak resonant frequency, these resonances are free to 760 vary (within the specified alpha range) across networks both between and within individuals. 761 e PCA was computed for subsets of SSE whose peak frequency lies within each of three 762 frequency bands. eta (1-7Hz), alpha (7-13Hz) and beta . Crucially, the inclusion of 763 an SSE in a band depends only on its peak frequency. Once included, all information on the 764 network structure within that SSE is included in the analysis, even if part of the spectral peak 765 goes outside the specified band. Reproducibility of the components arising from the PCA were 766 assessed using a split-half correlation. 500 split halves of the SSEs included in a PCA were 767 computed and the PCA computed on each half independently before the spatial components of 768 each half are then correlated and stored. Both the split-half correlation and proportion of 769 variance explained by each component was used in determining whether the component was 770 included in further analyses. e components describe the pa ern of variability across space 771 captured by that PC whilst the PC-score indicates the extent to which that shape is expressed 772 in each individual SSE. An example spatial map is computed for maximum and minimum 773 observed score in each PC by projecting that score back into the original data-space. 774 Relatively few SSEs survived the permutation scheme in the theta and beta bands.

775
Potentially as a result, the PCA components from these bands also showed relatively low 776 reproducibility. For completeness, the SSE and the first 4 PCs for these bands are included in 777 section 9.5. In contrast, over 1,200 SSEs were included in the alpha band and the first two PCA 778 components showed high variance explained and split-half reliability. ese are interpreted in 779 the main text and carried forward for further statistical analysis. In order to examine whether there was a relationship between the frequency of each mode and 782 the score with which it projected onto a given component, we performed a Bayesian linear 783 regression using the BRMS package [61,62]. For each PC, we scaled the scores by its standard 784 deviation and fit a model of Score Scaled F requency + (1|P articipant), allowing an overall 785 change in mean frequency per-participant. Model inference was performed using the standard 786 NUTS sampler used by STAN through BRMS. 787 e prior for the Frequency parameter was chosen to be normally distributed with a mean 788 of 0 and a standard deviation of 1; reflecting our default position that there was no a-prior 789 reason to expect frequency to vary with score. Altering the standard deviation of the prior to 790 other plausible ranges had no significant effect on the overall results. Examination of 791 diagnostic plots showed that the parameters have converged in all cases. regressor has more evidence is that the difference in LOO should be more than twice the 797 estimate of the LOO standard standard error. For models where the model with frequency was 798 assessed as having more evidence, we then report and assess the magnitude of the frequency 799 parameter in the full model along with is 95% Credible Interval (CI). Due to the scaling of the 800 scores, the frequency parameters is expressed in terms of the standard deviation of the score. 801     single SSE and colour indicates peak frequency. e x-axis contains the scores for PC-1 which correspond to overall alpha power. e PC scores are orthogonal and have no linear correlation. Low PC-1 scores indicate SSEs with low overall power whilst high scores indicate SSEs with strong alpha networks which, in turn, have greater variability in PC-2. e frequency correlation with PC2 can be seen as a greater density of high frequencies (red colours) above zero in the y-axis.