Spectral graph theory of brain oscillations

Abstract The relationship between the brain's structural wiring and the functional patterns of neural activity is of fundamental interest in computational neuroscience. We examine a hierarchical, linear graph spectral model of brain activity at mesoscopic and macroscopic scales. The model formulation yields an elegant closed‐form solution for the structure–function problem, specified by the graph spectrum of the structural connectome's Laplacian, with simple, universal rules of dynamics specified by a minimal set of global parameters. The resulting parsimonious and analytical solution stands in contrast to complex numerical simulations of high dimensional coupled nonlinear neural field models. This spectral graph model accurately predicts spatial and spectral features of neural oscillatory activity across the brain and was successful in simultaneously reproducing empirically observed spatial and spectral patterns of alpha‐band (8–12 Hz) and beta‐band (15–30 Hz) activity estimated from source localized magnetoencephalography (MEG). This spectral graph model demonstrates that certain brain oscillations are emergent properties of the graph structure of the structural connectome and provides important insights towards understanding the fundamental relationship between network topology and macroscopic whole‐brain dynamics. .


| The structure-function problem in neuroscience
It is considered paradigmatic in neuroscience that the brain's structure at various spatial scales is critical for determining its function. In particular, the relationship between the brain's structural wiring and the functional patterns of neural activity is of fundamental interest in computational neuroscience. Brain structure and function at the scale of macroscopic networks, that is, among identifiable gray matter (GM) regions and their long-range connections through white matter (WM) fiber bundles, can be adequately measured using current noninvasive measurement techniques. Fiber architecture can be measured from diffusion tensor imaging (DTI) followed by tractography algorithms (Hagmann et al., 2008;Iturria-Medina, 2013). Similarly, brain function manifested in neural oscillations can be measured noninvasively using magnetoencephalography (MEG) and reconstructed across whole-brain networks. Does the brain's white matter wiring structure constrain functional activity patterns that arise on the macroscopic Chang Cai and Xihe Xie contributed equally to this work. network or graph, whose nodes represent gray matter regions, and whose edges have weights given by the structural connectivity (SC) of white matter fibers between them? We address this critical open problem here, as the structural and functional networks estimated at various scales are not trivially predictable from each other (Honey et al., 2009).
A more detailed accounting of the structure-function relationship requires that we move beyond statistical descriptions to mathematical ones, informed by computational models of neural activity. Numerical simulations are available of mean field (Destexhe & Sejnowski, 2009;El Boustani & Destexhe, 2009;Wilson & Cowan, 1973) and neural mass (David & Friston, 2003;Deco et al., 2012) approximations of the dynamics of neuronal assemblies. By coupling many such neural field or mass models (NMMs) using anatomic connectivity information, it is possible to generate via large-scale stochastic simulations a rough picture of how the network modulates local activity at the global scale to allow the emergence of coherent functional networks (Deco et al., 2012).
However, simulations are unable to give an analytical (i.e., closed form) encapsulation of brain dynamics and present an interpretational challenge in that behavior is only deducible indirectly from thousands of trial runs of time-consuming simulations. Consequently, the essential minimal rules of organization and dynamics of the brain remain unknown. Furthermore, due to their nonlinear and stochastic nature, model parameter inference is ill posed, computationally demanding and manifest with inherent identifiability issues (Xie et al., 2018).
How then do stereotyped spatiotemporal patterns emerge from the structural substrate of the brain? How will disease processes perturb brain structure, thereby impacting its function? While stochastic simulations are powerful and useful tools, they provide limited neuroscientific insight, interpretability and predictive power, especially for the practical task of inferring macroscopic functional connectivity from long-range anatomic connectivity. Therefore, there is a need for more direct models of structural network-induced neural activity patterns-a task for which existing numerical modeling approaches, whether for single neurons, local assemblies, coupled neural masses or graph theory, are not ideally suited. Here we use a spectral graph model (SGM) to demonstrate that the spatial distribution of certain brain oscillations are emergent properties of the spectral graph structure of the structural connectome. Therefore, we also explore how the chosen connectome alters the functional activity patterns they sustain.
1.2 | A hierarchical, analytic, low-dimensional and linear spectral graph theoretic model of brain oscillations We present a linear graph model capable of reproducing empirical macroscopic spatial and spectral properties of neural activity. We are interested specifically in the transfer function (defined as the frequency-domain input-output relationship) induced by the macroscopic structural connectome, rather than in the behavior of local neural masses. Therefore, we seek an explicit formulation of the frequency spectra induced by the graph, using the eigendecomposition of the structural graph Laplacian, borrowing heavily from spectral graph theory used in diverse contexts including clustering, classification, and machine learning (Auffarth, 2007;Kondor, 2002;Larsen, Nielsen, Sporring, Zhang, & Hancock, 2006;Ng & M. Jordan YW., 2002). This theory conceptualizes brain oscillations as a linear superposition of eigenmodes. These eigen-relationships arise naturally from a biophysical abstraction of fine-scaled and complex brain activity into a simple linear model of how mutual dynamic influences or perturbations can spread within the underlying structural brain network, a notion that was advocated previously (Abdelnour et al., 2014;Galán, 2008;Goni et al., 2014). We had previously reported that the brain network Laplacian can be decomposed into its constituent "eigenmodes," which play an important role in both healthy brain function (Abdelnour et al., 2014;Abdelnour et al., 2018; Abdelnour, Dayan, Devinsky, Thesen, & Raj, 2015b; Atasoy, Donnelly, & Pearson, 2016;Wang, Owen, Mukherjee, & Raj, 2017) and pathophysiology of disease (Abdelnour, Mueller, & Raj, 2015c;Abdelnour, Raj, Devinsky, & Thesen, 2016;Raj, Kuceyeski, & Weiner, 2012;Wang et al., 2017).
We show here that a graph-spectral decomposition is possible at all frequencies, ignoring nonlinearities that are operating at the local (node) level. Like previous NMMs, we lump neural populations at each brain region into neural masses, but unlike them we use a linearized (but frequency-rich) local model-see Figure 1a. The macroscopic connectome imposes a linear and deterministic modulation of these local signals, which can be captured by a network transfer function. The sequestration of local oscillatory dynamics from the macroscopic F I G U R E 1 The linearized spectral graph model. (a) Conventional neural mass models typically instantiate a large assembly of excitatory and inhibitory neurons, which are modeled as fully connected internally. External inputs and outputs are gated through the excitatory neurons only, and inhibitory neurons are considered strictly local. The proposed linear model condenses these local assemblies into lumped linear systems f e (t) and f i (t), Gamma-shaped functions having time constants τ e and τ i -see panel (b). The recurrent architecture of the two pools within a local area is captured by the gain terms g ee , g ii , g ei , indicating the loops created by recurrents within excitatory, inhibitory and cross-populations. (c) The absolute value of eigenvalues of the complex Laplacian L ω ð Þ are plotted against the eigenvector index. Each dot represents one eigenvalue λ(ω); its color represents the frequency ω-low (blue) to high (yellow). Clearly, these eigenvalues change somewhat by frequency; small eigenvalues change more compared to large ones. (d) Frequency response of each eigenmode plotted on the complex plane with default choices of model parameters and a template structural connectome from the human connectome project (HCP). Each curve represents the transit in the complex plane of a single eigenmode's frequency response, starting at low frequencies in the bottom right quadrant, and moving characteristically to the upper left quadrant at high frequencies. The magnitude of the response, given by the distance from the origin, suggests that most eigenmodes have two prominent lobes, roughly corresponding to lower frequency alpha rhythms and higher frequency gamma rhythms, respectively. In contrast, the lowest few eigenmodes start off far from the origin, indicative of a low-pass response. (e) Magnitude of the frequency response of each eigenmode reinforces these impressions more clearly, with clear alpha peak, as well as slower rhythms of the lowest eigenmodes. The spectral profile of the eigenmodes, especially the peak frequencies, are sensitive to the choice of model parameters. (f) The spatial patterns of the top 5 eigenmodes of L ω ð Þ, evaluated at the alpha frequency, 10 Hz. The first 4 eigenmodes u 1 − u 4 , produce posterior and temporal spatial patterns, including many elements of the default mode network; u 4 resembles the sensorimotor network; and u 5 the structural core of the human connectome. However, these patterns are not exclusive and greatly depend on the frequency at which they are evaluated, as well as the model parameters and the connectome network in this way enables the characterization of whole brain dynamics deterministically in closed form in Fourier domain, via the eigen-basis expansion of the network Laplacian. As far as we know, this is the first closed-form analytical model of frequency-rich brain activity constrained by the structural connectome. We applied this model to and validated its construct against measured source-reconstructed MEG recordings in healthy subjects under rest and eyes-closed. The model closely matches empirical spatial and spectral MEG patterns. In particular, the model displays prominent alpha and beta peaks, and, intriguingly, the eigenmodes corresponding to the alpha oscillations have the same posterior-dominant spatial distribution that is repeatedly seen in eyes-closed alpha power distributions. In contrast to existing less parsimonious models in the literature that invoke spatially-varying parameters or local rhythm generators, to our knowledge, this is the first account of how the spectral graph structure of the structural connectome can parsimoniously explain the spatial power distribution of alpha and beta frequencies over the entire brain measurable on MEG.
2 | METHODS 2.1 | Spectral graph model development 2.1.1 | Notation In our notation, vectors and matrices are represented by boldface, and scalars by normal font. We denote frequency of a signal, in Hertz, by symbol f, and the corresponding angular frequency as ω = 2πf. The connectivity matrix is denoted by C = {c jk }, consisting of connectivity strength c ij between any two pair of regions j, k.

| Model summary
Details of the spectral graph model (SGM) is described in detail below.
There are very few model parameters, seven in total: τ e , τ i , τ G , v, g ii , g ei , α, which are all global and apply at every node. See Table 1 for their meaning, initial value and range. Note that the entire model is based on a single equation of graph dynamics, Equation (1), which is repeatedly applied to each level of the hierarchy. Here we used two levels: a mesoscopic level where connectivity is all-to-all, and a macroscopic level, where connectivity is measured from fiber architecture. In theory, this template could be refined into finer levels, where neural responses become increasingly nonlinear, and connectivity becomes sparser and structured.

| Canonical rate model over a graph
We use a canonical rate model to describe neural activity across two hierarchical levels-local cortical mesoscopic levels and long-range macroscopic levels. At each level of the hierarchy of brain circuits, we hypothesize a simple linear rate model of recurrent reverberatory activity given by where x e/i (t) is the mean signal of the excitatory/inhibitory populations and p e/i (t) is internal noise source reflecting local cortical column computations or input. The transit of signals, from presynaptic membranes, through dendritic arbors and axonal projections, is sought to be captured into ensemble average neural impulse response func- respectively. We disregard the nonlinearity of the neural response, hence the output at the terminal to a presynaptic input u(t) is the simple convolution   (Roland et al., 2014). Brief deflection of a single barrel whisker in the mouse evokes a somatotopically mapped cortical depolarization that remains localized to its C2 barrel column only for a few milliseconds, thence rapidly spreading to a large part of sensorimotor cortex within tens of milliseconds, a mechanism considered essential for the integration of sensory information (Ferezou et al., 2007;Polack & Contreras, 2012). Interestingly, the evoked response curve in S1 from the (Ferezou et al., 2007) study had a prominent Gamma shape. Of note, the duration of S1 response (~50 ms) was considerably longer than the time to first sensory response in C2 (7.2 ms) (Ferezou et al., 2007). Interestingly, feedback projections from higher to lower areas take~50 ms, hence have a much slower apparent propagation velocity (0.15-0.25 m/s) than what would be predicted by axonal conduction alone (1-3 m/s) (Roland et al., 2014).
Individual neural elements are connected to each other via connection strengths c jk . Let the cortico-cortical fiber conduction speed be v, which here is assumed to be a global constant independent of the pathway under question. For a given pathway connecting regions j and k of length d jk , the conduction delay of a signal propagating from region j to region k will be given by τ v jk = d jk v . Hence, signals from neighboring elements also participate in the same recurrent dynamics, giving the second term of Equation (1). Equation (1) will serve as our canonical rate model, and will be reproduced at all levels of the hierarchy, and only the connectivity strengths will vary depending on the level of hierarchy we are modeling, as explained below.

| Local neural assemblies
The local connectivities c local jk are assumed to be all-to-all, giving a complete graph. Further, the axonal delays τ v jk associated with purely local connections were already incorporated in the lumped impulse responses f e/i (t). Hence, we assert: From spectral graph theory, a complete graph has all equal eigenvalues, which allows the local network to be lumped into gain constants, and the summation removed. Indeed, rewriting x e/i (t) as the mean signal of all the excitatory/inhibitory cells and setting the gains g ee = 1 − c e N e and g ii = 1 − c i N i we get Given the Fourier transform pairs d dt $ jω , f e=i t ð Þ $ we take the Fourier transform of Equation (1) and obtain the local assembly's frequency spectrum: Writing this in terms of transfer functions X e (ω) = H e (ω)P e (ω), we get the lumped local system illustrated in Figure 1a. Finally, we must also account for signals that alternate between the two populations, which is given by the transfer function We fix g ee = 1 without loss of generality, and let the other terms g ii , g ei be model parameters to be fitted. Finally, the total cortical transfer function is the sum and X local (ω) = H local (ω)P(ω) represents all neural activity in this region, whether from excitatory or inhibitory cells. The canonical local activity is therefore defined by the Fourier transform pair: x local (t) $ X local (ω).

| Macroscopic scale: Signal evolution on the entire graph
For the macroscopic level, we use the same canonical network dynamics as Equation (1), but now the inter-regional connectivity c jk is nonzero and given by the structural connectome. Similarly, axonal conductance delays are determined by fiber length and conductance speed τ v jk = d jk =v. Further, the external driving signals at each node is the local neural activity x local (t) defined above rather than a noise process p(t). In the interest of parsimony we set each node of the macroscopic graph to have the same internal power spectrum X local (ω), that is, all regions are experiencing the same transfer function, driven by identically distributed (but of course not identical) noise. At this scale, activity measured at graph nodes is no longer excitatory or inhibitory, but mixed, and the corticocortical connections are all between long, pyramidal excitatory-only cells. Thus, for the kth node Here we have introduced a global coupling constant α, similar to most connectivity-coupled neural mass models, that seeks to control the relative weight given to long-range afferents compared to local signals. We have also introduced a new time constant, τ G , which is an excitatory time constant and it may be the same as the previously used constant τ e . However, we allow the possibility that the longrange projection neurons might display a different capacitance and morphology compared to local circuits, hence we have introduced τ G (subscript G is for "graph" or "global").
Stacking all equations from all nodes and using vector valued sig- where the braces {Á} represent all elements of a matrix indexed by j, k.
We wish to evaluate the frequency spectrum of the above. In Fourier space, delays become phases; hence we use the transform connectivity matrix" at any given angular frequency ω as We then define a normalized complex connectivity matrix at frequency ω as where the degree vector deg is defined as deg k = P j c jk . Taking the Fourier transform of Equation (5), we get where we assumed identically distributed noise signals driving both the excitatory and inhibitory local populations at each node, such that P e, k (ω) = P i, k (ω) = P k (ω) at the k-th node. We then collected all nodes' driving inputs in the vector P(ω) = {P k (ω), 8k}. Here, we define the complex Laplacian matrix where I is the identity matrix of size N × N. This complex Laplacian will be evaluated via the eigen-decomposition where is a diagonal matrix consisting of the eigenvalues of the complex Laplacian matrix of the connectivity graph C ω ð Þ, at the angular frequency ω.
where we invoke the eigen-decomposition of L ω ð Þ, and that U(ω)U (ω) H = I. It can then be shown easily that This is the steady state frequency response of the whole brain dynamics. In steady state, we assume that each cortical region is driven by internal noise that spans all frequencies, that is, white noise. Hence, we assume that the driving function p(t) is an uncorrelared Gaussian noise process, such that P ω ð Þ= , where  is a vector of ones. This asserts identical cortical responses at each brain region.

| MRI
A 3 Tesla TIM Trio MR scanner (Siemens, Erlangen, Germany) was used to perform MRI using a 32-channel phased-array radiofrequency head coil. High-resolution MRI of each subject's brain was collected using an axial 3D magnetization prepared rapid-acquisition gradient-  The T1-weighted images were parcellated into 68 cortical regions and 18 subcortical regions using the using the Desikan-Killiany atlas available in the FreeSurfer software (Fischl et al., 2002). To do this, the subject specific T1-weighted images were back-projected to the atlas using affine registration, as described in our previous studies

| Structural connectivity networks
We constructed different structural connectivity networks with the same Desikan-Killiany parcellations to access the capabilities of our proposed model. Firstly, we obtained openly available diffusion MRI data from the MGH-USC Human Connectome Project to create an average template connectome. As in our previous studies (Abdelnour et al., 2014;Owen et al., 2013), subject specific structural connectivity was computed on diffusion MRI data: Bedpostx was used to determine the orientation of brain fibers in conjunction with FLIRT, as implemented in the FSL software (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012). In order to determine the elements of the adjacency matrix, we performed tractography using probtrackx2. We initiated 4,000 streamlines from each seed voxel corresponding to a cortical or subcortical gray matter structure and tracked how many of these streamlines reached a target gray matter structure. The

| Alternative benchmark model for comparison
In order to put the proposed model in context, we also implemented for comparison a Wilson-Cowan neural mass model (Destexhe & Sejnowski, 2009;Muldoon, Pasqualetti, Gu, et al., 2016;Wilson & Cowan, 1973;Xie et al., 2018) with similar dimensionality. Although NMMs like this can and have been implemented with regionally varying local parameters, here we enforced uniform, regionally nonvarying local parameters, meaning all parcellated brain regions shared the same local and global parameters. This is a fair comparison since the proposed model is also regionally nonvarying. The purpose of this exercise is to ascertain whether a nonregional NMM can also predict spatial power variations purely as a consequence of network transmission, like the proposed model, using the same model optimization procedure (see below). This NMM incorporates a transmission velocity parameter that introduces a delay based on fiber tract lengths extracted from diffusion MRI, but, unlike our model, does not seek to explicitly evaluate a frequency response based on these delays.

| Model optimization
We computed maximum a posteriori estimates for parameters under a flat noninformative prior. A simulated annealing optimization algorithm was used for estimation and provided a set of optimized parameters {τ e , τ i , τ c , g ei , g ii , α, υ}. We defined a data likelihood or goodness of fit (GOF) as the Pearson correlation between empirical source localized MEG power spectra and simulated model power spectra, averaged over all 68 regions of a subject's brain. The proposed model has only seven global parameters as compared to neural mass models with hundreds of parameters, and is available in closedform. To improve the odds that we capture the global minimum, we chose to implement a probabilistic approach of simulated annealing (Kirkpatrick, Gelatt, & Vecchi, 1983). The algorithm samples a set of parameters within a set of boundaries by generating an initial trial solution and choosing the next solution from the current point by a probability distribution with a scale depending on the current "temperature" parameter. While the algorithm always accepts new trial points that map to cost-function values lower than the previous costfunction evaluations, it will also accept solutions that have costfunction evaluations greater than the previous one to move out of local minima. The acceptance probability function is 1= 1 + Δ where T is the current temperature and Δ is the difference of the new minus old cost-function evaluations. The initial parameter values and boundary constraints for each parameter are given in Table S1. All simulated annealing runs were allowed to iterate over the parameter space for a maximum of N p × 3000 iterations, where N p is the number of parameters in the model. As a comparison, we performed the same optimization procedure to a regionally nonvarying Wilson-Cowan neural mass model (Muldoon et al., 2016;Wilson & Cowan, 1973).
We have recently reported a similar simulated annealing optimization procedure on this model (Xie et al., 2018).

| Graph Laplacian eigenmodes mediate a diversity of frequency responses
First, we demonstrate the spectra produced by graph eigenmodes as per our theory using default choices of model parameters. Figure 1c shows the eigen-spectrum of the complex Laplacian, with eigenvalue magnitude ranging from 0 to 1. The absolute value of eigenvalues of the complex Laplacian L ω ð Þ are plotted against the eigenvector index.
Each dot represents one eigenvalue λ(ω); its color represents the frequency ω-low (blue) to high (yellow). Clearly, these eigenvalues change somewhat by frequency. Small eigenvalues undergo a larger shift due to frequency, while the large ones stay more stable and tightly clustered around the nominal eigenvalue (i.e., at ω = 0). Each eigenmode produces a frequency response based on its frequencydependent eigenvalue (Figure 1d,e). Figure  The spatial patterns of the first 5 eigenmodes of L ω ð Þ, evaluated at the alpha peak of 10 Hz, are shown in Figure 1f. Eigenmodes u 1 − 4 produce posterior and temporal spatial patterns, including many elements of the default mode network; u 4 resembles the sensorimotor network; and u 5 the structural core of the human connectome. However, these patterns are not exclusive and greatly depend on the frequency at which they are evaluated, as well as the model parameters. Higher eigenmodes are especially sensitive to axonal velocity and frequency (not shown here).
Since the spectral graph model (SGM) relies on connectome topology, we demonstrate in Figure 2  for the HCP template were used. We can observe the spectral profile of the eigenmodes, especially the peak frequencies, are sensitive to the choice of the connectome and the model parameters. All modeled power spectra show a broad alpha peak at around 10 Hz and a narrower beta peak at around 20 Hz. This is expected, since these general spectral properties are governed by the local linearized neural mass model. It is important to note that different eigenmodes accommodate a diversity of frequency responses; for instance, the lowest eigenmodes show a lowfrequency response with no alpha peak whatsoever. In the frequency responses from biologically realistic individual and HCP template connectomes, there is a diversity of spectral responses among eigenmodes that is lacking in the response produced by the unrealistic uniform and randomized matrices. As we will see below, graph topology is critical to the power spectrum it induces, hence we explored whether and how sparsity of random graphs mediates spectral power (Figure 2d-f). At incrementally increasing sparsity levels, the diversity of spectral responses of different eigenmodes increases and approaches that of realistic connectomes. Therefore, graph eigenmodes induce unique and diverse frequency responses that depend on the topology of the graph.

| Spectral distribution of MEG power depends on model parameters but not connectivity
Network eigenmodes exhibit strong spatial patterning in their frequency responses, even with identical model parameters (Figure 3). We evaluated the model spectral response using the subject-specific C ind matrices of four representative subjects (Figure 3a). The model power spectra strikingly resemble empirical MEG spectra, displaying both the alpha and beta peaks on average, and similar regional variability as in real data.
Regional averages of empirical and modeled power spectra of the entire group after full parameter optimization over individual subjects are shown in Figure 3b. The model closely replicates the observed power spectrum (red circles) equally well with both C ind (black triangles) and C HCP (purple triangles). Thus, in most cases we can safely replace the subject-specific connectome with the template connectome. In contrast, when nonoptimized average parameters were used (golden green triangles), it resulted in a worse fit, especially at high frequencies, suggesting that individualized parameter optimization is essential to produce realistic spectra. We also examined the model behavior for a random connectomes with 80% sparsity (bright green triangles), or a distance-based connectome (blue triangles) was chosen with identical sparsity (80%) to the actual connectome, and found that even with optimized parameters the average spectra could be accounted for by these connectomes. sparsity. (f) Model using randomized connectivity matrix with 95% sparsity. In all cases the connectome modulates the spectral response in delta-beta range, leaving the higher gamma frequencies unchanged. In general, the low eigenmodes (u 1 − u 20 ) appear to modulate the lower frequency range, up to beta, and may be considered responsible for the diversity of spectra observed in the model across subjects. The time constants τ e , τ i showed tight clustering but the rest of the parameters showed high variability across subjects.
The optimal parameters are in a biologically plausible range, similar to values reported in numerous neural mass models. The optimization algorithm aimed to maximize a cost function proportional to the posterior likelihood of the model, and was quantified by the Pearson's correlation between MEG and modeled spectra ("Spectral correlation"). The convergence plots shown in Figure 4b, -20). Therefore, we conclude that with the graph spectral model, the overall regional spectra appear to be dependent on global model parameters rather than on the actual structural connectome.

| Spectral graph model recapitulates the spatial distribution of MEG power
Next, we establish that the model is able to reproduce region-specific spectra, even though it uses identical local oscillations. We integrated the spectral area in the range 8-12 Hz for alpha and 13-25 Hz for beta, of each brain region separately. We define "spatial correlation" (as compared to spectral correlation above) as Pearson's R between the regional distribution of empirical MEG and model-predicted power within a given frequency band.

| A small number of eigenmodes capture spatial distributions of alpha and beta band activity
We noticed during our experimentation that only a few eigenmodes appear to contribute substantially to observed MEG alpha and beta patterns. Hence, we hypothesized that spatial correlations could be improved by selecting a small subset of eigenmodes. Therefore, we developed a sorting strategy whereby we first rank the eigenmodes in descending order of spatial correlation for a given subject and given frequency band. Then we perform summation over only these eigenmodes according to Equation (10), each time incrementally adding a new eigenmode to the sum. The spatial correlation of these "sortedsummed" eigenmodes against empirical alpha power are plotted in Figure 5c as a function of increasing number of eigenmodes; one curve for each subject. The thick black curve represents the average over all subjects. The spatial correlation initially increases as we add more well-fitting eigenmodes, but peaks around, and begins declining thereafter. Addition of the remaining eigenmodes only serves to reduce the spatial correlation. This behavior is observed in almost all subjects we studied.
Examples of predicted alpha patterns: Figure 5 shows brain surface renderings of the spatially distributed patterns of alpha band power for two representative subjects. Regions are color coded as a heatmap of regional power scaled by mean power over all regions.
The observed MEG spatial distribution pattern of alpha band shows higher power in posterior regions of the brain, as expected, with strong effect size in temporal, occipital and medial posterior areas.
This pattern is matched by one of the eigenmodes (#10, shown in middle panel, giving R = 0.65), and slightly better by a weighted combination of two eigenmodes (R = 0.69). However, the model did not reproduce parietal and parieto-occipital components seen in real data.
The other subject produced similar results, but with six eigenmodes.
F I G U R E 5 Alpha power spatial distribution depicted by specific spectral graph model eigenmodes. (a, b) The spatially distributed patterns of alpha band power for two representative subjects are displayed in brain surface renderings. For each four brain panels shown, the medial surface is rendered on the left column while the lateral surface is rendered on the right, the left hemisphere rendering is shown on top while the right hemisphere rendering is shown in the bottom row. Left column: The observed MEG spatial distribution pattern for alpha band power showing higher power in posterior regions of the brain. Middle column: Spatial distribution of the best matching eigenmode from the spectral graph model. The spatial correlation values are shown on top. Right column: Spatial distribution of the best cumulative combination of eigenmodes from the spectral graph model. Spatial correlation values and the number of eigenmodes are shown on top. (c) Across subject distribution of the alpha band spatial correlation values from spectral graph model simulations for the best fit eigenmodes and the cumulative combination of an increasing number of eigenmodes. Individual subject specific alpha band spatial correlation curves are shown in cyan (n = 36). Panels A and B correspond to the subjects indicated by red and blue curves respectively. Black curve is the average performance across all subjects In this instance, the parietal component seen in real data were reasonably reproduced by the model.
Examples of predicted beta patterns: Empirical beta power (Figure 6, left) is spread throughout the cortex, especially frontal and premotor cortex. A combination of four and six best matching eigenmodes produced the best model match to the source localized pattern of two representative subjects, respectively, with R = 0.55 and 0.48. Figure 6c shows how the spatial correlation changes as more eigenmodes are used in the "sorted summed" algorithm, analogous to that of alpha pattern. Here too, a peak is achieved for a small number of eigenmodes, typically under 10.
3.5 | Spatial correlation achieved by the spectral graph model is significantly higher than alternative models The distribution of peak spatial correlations in the alpha band, using optimized parameters and individual connectomes of all subjects, is plotted in Figure 7a. For comparison we show results for four models:  Figure 7b. For each connectome model, a selection of the cumulative best set of eigenvectors were separately obtained for spatial correlation calculations. Across all subjects the proposed model, SGM on C Ind , gives excellent spatial correlations in alpha band (R distribution centered at 0.6) as well as in the beta band (R distribution centered at 0.5). For both alpha and beta spatial distribution patterns, paired t-tests between SGM with C Ind and all other models show that, the SGM with C Ind significantly outperforms all other models, as determined by a paired t-test; p < .012 in each case, denoted by asterisk.

| Alternate nonlinear model
The Wilson-Cowan neural mass model also did not succeed in predicting the spatial patterns of alpha or beta power, with poor F I G U R E 6 Beta power spatial distributions depicted by specific spectral graph model eigenmodes. Legend is identical to Figure 5 but shown for beta power spatial distributions correlations (r centered at 0). This could be because in our implementation we enforced uniform local parameters with no regional variability. However, this is the appropriate comparison, since our proposed model also does not require regionally varying parameters. Interestingly, the random connectomes and geodesic distance based connectome also appear to have some ability to capture these spatial patterns (r centered at 0.4 and 0.2 respectively), perhaps due to the implicit search for best performing eigenmodes, which on average will give at least a few eigenmodes that look like MEG power purely by chance.
Collectively, we conclude that the graph model is able to fit both the spectral and spatial features of empirical source localized MEG data, and that the optimal fits performed on individual subjects occurs at widely varying subject-specific parameter choices.

| DISCUSSION
The proposed hierarchical graph spectral model of neural oscillatory activity is a step towards understanding the fundamental relationship between network topology and the macroscopic whole-brain dynamics. The objective is not just to model brain activity phenomenologically, but to analytically derive the mesoscopic laws that drive F I G U R E 7 Spatial correlation performance analysis of the spectral graph model. Distribution of the best fit spatial correlations of the spectral graph model across all subjects. (a) Alpha band spatial correlations. (b) Beta band spatial correlations. For both panels, spatial correlations are shown for spectral graph model (SGM) with subject specific individual connectomes (C Ind , black), random connectomes with 80% sparsity comparable to individual connectomes and 95% sparsity where the SGM model eigenmodes show spectral diversity (C Rdm , blue), geodesic distance based connectomes (C Dist ,green) and for a neural mass model (NMM) with subject specific individual connectome (C Ind , pink). For each model, a selection of the cumulative best set of eigenvectors were separately obtained for spatial correlation calculations. For both alpha and beta spatial distribution patterns, paired t-tests between SGM with C Ind and all other models show that, the SGM with C Ind significantly outperforms all other models, as determined by a paired t-test; p < .001 for all alpha spatial correlations and p < .012 for all beta spatial correlations, denoted by asterisk We designed a comprehensive parameter optimization algorithm on individual subjects' MEG data of a suitably defined cost function based on Pearson R statistic as a way to capture all relevant spectral features. Using this fitting procedure, we were able to obtain the range of optimally fitted parameters across the entire study cohort. As shown in Figure 4a, the range is broad in most cases, implying that there is significant inter-subject variability of model parameters, even if a template connectome is used for all. We tested the possibility that a group-averaged parameter set might also succeed in matching real spectral data on individuals. But as shown in Figures 3b and 4c, this was found to be a poor choice, supporting the key role of individual variability of model parameters (but not variability in the connectome).
However, no model is capable of reproducing higher frequencies in the higher beta and gamma range seen in MEG, since by design and by biophysical intuition, these frequencies arise from local neural assemblies rather than from modulation by macroscopic networks.

| Revealing sources of heterogeneity in spatial patterns of brain activity
The spatial match between model and data is strongest when the model uses empirical macroscopic connectomes obtained from healthy subjects' diffusion weighted MRI scans, followed by tractography. The use of "null" connectomes-randomized connectivity matrices of varying levels of sparsity and distance-based connectivity matrices, respectively, did far worse than actual human connectomes (Figure 7), supporting the fact that the latter is the key mediator of spatial patterns of real brain activity. The match was also significantly different when using a template HCP connectome versus the individual subject's own connectomes (Figure 7), and when compared to spatial patterns predicted by an NMM. In conclusion, for the purpose of predicting the spatial topography of brain activity, it is important to use individual connectomes and optimized model parameters.

| Macroscopic brain rhythms are governed by the connectome
A predominant view assumes that different brain rhythms are produced by groups of neurons with similar characteristic frequencies, which might synchronize and act as "pacemakers." How could this view explain why alpha and beta power are spatially stereotyped across subjects, and why the alpha signal is especially prominent in posterior areas? Although practically any computer model of cortical activity can be tuned, with suitable parameter choice, to oscillate at alpha frequency, for example, (David & Friston, 2003;Deco et al., 2012;Liley et al., 1999;Nakagawa et al., 2014;Nunez & Srinivasan, 2006;Robinson et al., 2005;Vijayan, Ching, Purdon, Brown, & Kopell, 2013), none of them were able to parsimoniously recapitulate the posterior origin of alpha. Thus, the prominence of posterior alpha might be explained by the hypothesized existence of alpha generators in posterior areas. Indeed, most oscillator models of local dynamics are capable of producing these rhythms at any desired frequency (David & Friston, 2003;Liley et al., 1999;Liley, Cadusch, & Dafilis, 2002;Spiegler & Jirsa, 2013;van Rotterdam, Lopes da Silva, van den Ende, Viergever, & Hermans, 1982), and therefore it is common to tweak their parameters to reproduce alpha rhythm. Local networks of simulated multicompartmental neurons can produce oscillations in the range 8-20 Hz 5 , and, in a nonlinear continuum theory, peaks at various frequencies in the range 2-16 Hz were obtained depending on the parameters (Liley et al., 2002). Specifically, the role of thalamus as pacemaker has motivated thalamocortical models (Izhikevich & Edelman, 2008;Robinson et al., 2005) that are capable of resonances in various ranges. Neural field models of the thalamocortical loop (Robinson et al., 2005) can also predict slow-wave and spindle oscillations in sleep, and alpha, beta, and higher-frequency oscillations in the waking state. In these thalamocortical models, the posterior alpha can arise by postulating a differential effect in weights of the posterior versus anterior thalamic projections, for example, (Vijayan et al., 2013). Ultimately, hypotheses requiring local rhythm generators suffer from lack of parsimony and specificity: a separate pacemaker must be postulated for each spectral peak at just the right location (Nunez, 1981).
An alternative view emerges from our results that macroscopic brain rhythms are governed by the structural connectome. Even with global model parameters, using the exact same local cortical dynamics captured by the local transfer function H local (ω), driven by identically distributed random noise P(ω), our model is capable of predicting prominent spectral (Figures 3 and 4) and spatial ( Figures 5 and 6) patterning that is quite realistic. This is especially true in the lower frequency range: indeed, the model is able to predict not just the frequency spectra in alpha and beta ranges, but also their spatial patterns-that is, posterior alpha and distributed but roughly frontal beta. Although this is not definitive proof, it raises the intriguing possibility that the macroscopic spatial distribution of the spectra of brain signals does not require spatial heterogeneity of local signal sources, nor regionally variable parameters. Rather, it implies that the most prominent patterning of brain activity (especially alpha) may be governed by the topology of the macroscopic network rather than by local, regionally varying drivers. Nevertheless, a deeper exploration is required of the topography of the dominant eigenmodes of our linear model, in order to understand the spatial gradients postulated previously (Robinson et al., 2005;Vijayan et al., 2013).

| Emergence of linearity from chaotic brain dynamics
The nonlinear and chaotic dynamics of brain signals may at first appear to preclude deterministic or analytic modeling of any kind. Yet, vast swathes of neuroscientific terrain are surprisingly deterministic, reproducible and conserved across individuals and even species. Brain rhythms generally fall within identical frequency bands and spatial maps (Freeman & Zhai, 2008;He et al., 2010;Robinson et al., 2005). Based on the hypothesis that the emergent behavior of long-range interactions can be independent of detailed local dynamics of individual neurons (Abdelnour et al., 2014;Destexhe & Sejnowski, 2009;Miši c et al., 2014;Miši c et al., 2015;Robinson et al., 2005;Shimizu & Haken, 1983), and may be largely governed by long-range connectivity (Abdelnour, Raj, et al., 2015a;Deco et al., 2012;Jirsa et al., 2002;Nakagawa et al., 2014), we have reported here a minimal linear model of how the brain connectome serves as a spatial-spectral filter that modulates the underlying nonlinear signals emanating from local circuits. Nevertheless, we recognize the limitations of a linear model and its inability to capture inherent nonlinearities across all levels in the system.
First, the goal of many prior efforts using DCMs is to examine effective connectivity in EEG, LFP and fMRI functional connectivity data, typically for smaller networks (Daunizeau, Kiebel, & Friston, 2009;Pinotsis et al., 2017), or dynamic effective connectivity (Park, Friston, Pae, Park, & Razi, 2018;Preti, Bolton, & Van De Ville, 2017; Van de Steen, Almgren, Razi, Friston, & Marinazzo, 2019). Hence, they address the second order covariance structures of brain activity. In particular, recent spectral DCM and regression DCM models (Frässle et al., 2017;Frässle et al., 2018;Razi et al., 2017) with local neural masses are formulated in the steady-state frequency-domain, and the resulting whole-brain cross-spectra are evaluated. The goals of these models are to derive model cross-spectra that define the effective connectivity in the frequency domain and are compared with empirical cross-spectra.
Based on second-order sufficient statistics, these models attempt to derive effective connectivity from functional connectivity data. These DCMs have so far only been applied to small networks or to BOLD fMRI regime. In contrast, our goal is to examine the role of the eigenmodes of the structural connectome and their influence on power spectral distributions in the full MEG frequency range, and over the entire whole brain. In subsequent work, we intend to extend our efforts to examining effective connectivity but such an effort currently remains outside the scope of the work in this paper. Here, we focus on models that directly estimate the first order effects of observed power spectra and its spatial distributions and compare them with empirical MEG source reconstructions. Our primary motivation is to examine whether spatial distribution of observed power spectra can arise from graph structure of the connectome, hence our focus on the effects of model behavior as a function of the underlying structural connectomewhether it is individualized, template-based, uniform, random or distance based. DCM methods have not reported first order regional power spectra as we do here, nor have they explored how the structural connectome influences model spectral distributions. Second, our model is more parsimonious compared to most of these above-mentioned models, which have many more degrees of freedom because they often allow for regions and their interactions to have different parameters. Our model parameterization, with only a few global parameters, lends itself to efficient computations over fine-scale whole-brain parcellations, whereas most DCMs (with the exception of recent spectral and regression DCMs (Frässle et al., 2017;Frässle et al., 2018;Razi et al., 2017)) are suited for examining smaller networks but involve large effective connectivity matrices and region-specific parameters. Furthermore, parameters of our model remain grounded and interpretable in terms of the underlying biophysics, that is, time constants and conductivities. In contrast, spectral and regression DCM models of cross-spectra have parameters that are abstract and do not have immediate biophysical interpretation.
The third major difference is in the emphasis placed on Variational Bayesian inference in DCM. Since our focus was on exploring model behavior over a small number of global parameters and a set of structural connectomes (whether anatomic or random) of identical sparsity and complexity, it was sufficient to use a maximum a posteriori (MAP) estimation procedure for Bayesian inference of our global model parameters with flat noninformative priors with predetermined ranges based on biophysics. Like most DCM efforts our model can be easily be extended to Variational Empirical Bayesian inference for parameter estimation, for instance to compute a full posterior of the structural connectivity matrix. In such a formulation, we can assume that the observed structural connectome will serve as the prior mean of the connectivity matrix. We reserve such extensions to our future work with this spectral graph model.

| Other limitations and extensions
The model currently examines resting-state activity, but future extensions will include prediction of functional connectivity, task-induced modulations of neural oscillations and causal modeling of external stimuli, for example, transcranial magnetic and direct current stimulation. The current implementation does not incorporate complex local dynamics, but future work will explore using nonwhite internal noise and chaotic dynamics for local assemblies. This may allow us to examine higher gamma frequencies. Although our model incorporates latency information derived from path distances, we plan to explore path-specific propagation velocities derived from white matter microstructural metrics such as axon diameter distributions and myelin thickness. Future work will also examine the specific topographic features of the structural connectome that may best describe canonical neural activity spectra. Finally, we plan to examine the ability of the model to predict time-varying structure-function relationships.

| Potential applications
Mathematical encapsulation of the structure-function relationship can potentiate novel approaches for mapping and monitoring brain diseases such as autism, schizophrenia, epilepsy and dementia, since early functional changes are more readily and sensitively measured using fMRI and MEG, compared to structural changes. Because of the complementary sensitivity, temporal and spatial resolutions of diffusion MRI, MEG, EEG and fMRI, combining these modalities may be able to reveal fine spatiotemporal structures of neuronal activity that would otherwise remain undetected if using only one modality. Current efforts at fusing multimodalities are interpretive, phenomenological or statistical, with limited cognizance of underlying neuronal processes. Thus, the ability of the presented model to quantitatively and parsimoniously capture the structure-function relationship may be key to achieving true multimodality integration.

DATA AVAILABILITY STATEMENT
The data and code that support the findings of this study are available from the GitHub repository at https://github.com/Raj-Lab-UCSF/ spectrome (Xie, Stanley, & Damasceno, 2019). The code used to produce basic figures can be run as interactive Jupyter notebooks via Binder (Jupyter, Bussonnier, Forde, et al., 2018). Some raw imaging data, for example, MRI scans and MEG recordings are not appropriate for public sharing and are too large to be saved in an online repository.
However, they could be made available by corresponding author upon reasonable request.