Spectral graph theory of brain oscillations.

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. .

INTRODUCTION

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 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).

Although numerical models of single neurons and local microscopic neuronal assemblies, ranging from simple integrate‐and‐fire neurons to detailed multi‐compartment and multi‐channel models (Freeman & Zhai, 2008; Liley, Alexander, Wright, & Aldous, 1999; Markounikau, Igel, Grinvald, & Jancke, 2010; Roland, Hilgetag, & Deco, 2014; Schaffer, Ostojic, & Abbott, 2013) have been proposed, it is unclear if these models can explain structure–function coupling at meso‐ or macroscopic scales. At one extreme, the Blue Brain Project (Markram, n.d.; Markram et al., 2015) seeks to model in detail all 10 neurons and all their connections in the brain. Indeed spiking models linked up via specified synaptic connectivity and spike timing dependent plasticity rules were found to produce regionally and spectrally organized self‐sustaining dynamics, as well as wave‐like propagation similar to real fMRI data (Izhikevich & Edelman, 2008). However, it is unclear whether such efforts will succeed in providing interpretable models at whole‐brain scale (Potjans & Diesmann, 2014).

Therefore, the traditional computational neuroscience paradigm at the microscopic scale does not easily extend to whole‐brain macroscopic phenomena, as large neuronal ensembles exhibit emergent properties that can be unrelated to individual neuronal behavior (Abdelnour, Voss, & Raj, 2014; Destexhe & Sejnowski, 2009; Mišić et al., 2015; Mišić, Sporns, & McIntosh, 2014; Robinson, Rennie, Rowe, O'Connor, & Gordon, 2005; Shimizu & Haken, 1983), and are instead largely governed by long‐range connectivity (Abdelnour, Raj, Dayan, & Devinsky, 2015a; Deco, Senden, & Jirsa, 2012; Jirsa, Jantzen, Fuchs, & Kelso, 2002; Nakagawa et al., 2014). At this scale, graph theory involving network statistics can phenomenologically capture structure–function relationships (Achard, Salvador, Whitcher, Suckling, & Bullmore, 2006; Bullmore, Bullmore, Sporns, & Sporns, 2009; Strogatz, 2001), but do not explicitly embody any details about neural physiology (Mišić et al., 2014; Mišić et al., 2015). Strong correlations between functional and structural connections have also been observed at this scale (Abdelnour, Dayan, Devinsky, Thesen, & Raj, 2018; Abdelnour, Voss, & Raj, 2014; Ghosh, Rho, McIntosh, Kötter, & Jirsa, 2008; Hermundstad et al., 2013; Honey et al., 2009; Park & Friston, 2013; Rubinov, Sporns, van Leeuwen, & Breakspear, 2009; van den Heuvel, Mandl, Kahn, & Hulshoff Pol, 2009), and important graph properties are shared by both SC and functional connectivity (FC) networks, such as small worldness, power‐law degree distribution, hierarchy, modularity, and highly connected hubs (Bullmore et al., 2009; He, Zempel, Snyder, & Raichle, 2010).

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.

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 eigen‐decomposition 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 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.

Figure 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 and , Gamma‐shaped functions having time constants and —see panel (b). The recurrent architecture of the two pools within a local area is captured by the gain terms , indicating the loops created by recurrents within excitatory, inhibitory and cross‐populations. (c) The absolute value of eigenvalues of the complex Laplacian 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 , evaluated at the alpha frequency, 10 Hz. The first 4 eigenmodes , produce posterior and temporal spatial patterns, including many elements of the default mode network; resembles the sensorimotor network; and 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

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.

METHODS

Spectral graph model development

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 , and the corresponding angular frequency as . The connectivity matrix is denoted by , consisting of connectivity strength between any two pair of regions .

Model summary

Details of the spectral graph model (SGM) is described in detail below. There are very few model parameters, seven in total: , 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.

Table 1

SGM parameters values and limits

Name	Symbol	Initial/default value	Lower/upper bound for optimization
Excitatory time constant	τe	12 ms	[5 ms, 20 ms]
Inhibitory time constant	τi	3 ms	[5 ms, 20 ms]
Graph time constant	τG	6 ms	[5 ms, 20 ms]
Excitatory gain	gee	1	n/a
Inhibitory gain	gii	1	[0.5, 5]
Excitatory gain	gei	4	[0.5, 5]
Transmission velocity	v	5 m/s	[5 m/s, 20 m/s]
Long‐range connectivity coupling constant	α	1	[0.1, 1]

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 bywhere is the mean signal of the excitatory/inhibitory populations and 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 functions and respectively. We disregard the nonlinearity of the neural response, hence the output at the terminal to a presynaptic input is the simple convolution . The neural responses are Gamma‐shaped responses (Figure 1b) parameterized by time constants that here represent the end result of both synaptic membrane capacitance and the distribution of dendritic/axonal delays introduced by the arborization. NMMs typically use a single classical exponential decay term for membrane capacitance only, since NMMs model highly local cell assemblies where multisynaptic connections are infrequent and axonal and dendritic transport delays are usually incorporated explicitly via connectivity weights and delays. Since our lumped model was designed for relatively large cortical regions, we employ the Gamma‐shaped to capture not just classical membrane capacitance but also the expected diversity of dendritic transport delays. The dynamics of the entire assembly modeled via a self‐decaying term , typically used in most rate or NMM models, but the difference here is that we chose to apply convolution with neural response within the decay process. We believe this is necessary to ensure that the dynamics of the population cannot participate in the internal recurrent dynamics of the region until the signal has passed through one instance of the neuronal response. Since this neural response is meant to capture a distribution of local circuit delays, its time constants are purposefully far longer (up to 20 ms) than expected from membrane capacitance alone. Studies of cortical lag times using paired electrode recordings between primary and higher cortices demonstrate this. A short visual stimulus causes a neural response in the ferret V1 within 20 ms post‐stimulus, in the primary barrel field within 16–36 ms, and the entire visual cortex becomes engaged 48–70 ms after stimulus (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 . Let the cortico‐cortical fiber conduction speed be , 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 , the conduction delay of a signal propagating from region j to region k will be given by . 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 are assumed to be all‐to‐all, giving a complete graph. Further, the axonal delays associated with purely local connections were already incorporated in the lumped impulse responses . 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 as the mean signal of all the excitatory/inhibitory cells and setting the gains and we get

Given the Fourier transform pairs , , we take the Fourier transform of Equation (1) and obtain the local assembly's frequency spectrum:

Writing this in terms of transfer functions 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 without loss of generality, and let the other terms be model parameters to be fitted. Finally, the total cortical transfer function is the sumand 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: .

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 is nonzero and given by the structural connectome. Similarly, axonal conductance delays are determined by fiber length and conductance speed . Further, the external driving signals at each node is the local neural activity defined above rather than a noise process . In the interest of parsimony we set each node of the macroscopic graph to have the same internal power spectrum , 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 k‐th 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, , which is an excitatory time constant and it may be the same as the previously used constant . However, we allow the possibility that the long‐range projection neurons might display a different capacitance and morphology compared to local circuits, hence we have introduced (subscript G is for “graph” or “global”).

Stacking all equations from all nodes and using vector valued signals , we can writewhere the braces {·} represent all elements of a matrix indexed by .

We wish to evaluate the frequency spectrum of the above. In Fourier space, delays become phases; hence we use the transform pairs and . Therefore, define a “complex connectivity matrix” at any given angular frequency as . We then define a normalized complex connectivity matrix at frequency as

where the degree vector is defined as . Taking the Fourier transform of Equation (5), we getwhere we assumed identically distributed noise signals driving both the excitatory and inhibitory local populations at each node, such that at the k‐th node. We then collected all nodes' driving inputs in the vector . Here, we define the complex Laplacian matrixwhere is the identity matrix of size This complex Laplacian will be evaluated via the eigen‐decompositionwhere is a diagonal matrix consisting of the eigenvalues of the complex Laplacian matrix of the connectivity graph , at the angular frequency .

Hencewhere we invoke the eigen‐decomposition of , and that 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 is an uncorrelared Gaussian noise process, such that where is a vector of ones. This asserts identical cortical responses at each brain region.

Experimental procedures

Study cohort

We acquired MEG, anatomical MRI, and diffusion MRI for 36 healthy adult subjects (23 males, 13 females; 26 left‐handed, 10 right‐handed; mean age 21.75 years (range: 7–51 years). All study procedures were approved by the institutional review board at the University of California at San Francisco (UCSF) and are in accordance with the ethics standards of the Helsinki Declaration of 1975 as revised in 2008.

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‐echo (MPRAGE) T1‐weighted sequence (echo time [TE] = 1.64 ms, repetition time [TR] = 2,530 ms, TI = 1,200 ms, flip angle of 7°) with a 256‐mm field of view (FOV), and 160 1.0‐mm contiguous partitions at a 256 × 256 matrix. Whole‐brain diffusion weighted images were collected at b = 1000s/mm with 30 directions using 2‐mm voxel resolution in‐plane and through‐plane.

Magnetoencephalography data

MEG recordings were acquired at UCSF using a 275‐channel CTF Omega 2000 whole‐head MEG system from VSM MedTech (Coquitlam, BC, Canada). All subjects were instructed to keep their eyes closed for 5 min while their MEGs were recorded at a sampling frequency of 1,200 Hz.

Data processing

Region Parcellations

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 (Abdelnour et al., 2014; Owen et al., 2013).

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 weighted connection between the two structures , was defined as the number of streamlines initiated by voxels in region that reach any voxel within region , normalized by the sum of the source and target region volumes (). This normalization prevents large brain regions from having high connectivity simply due to having initiated or received many streamlines. Afterwards, connection strengths are averaged between both directions ( and ) to form undirected edges. It is common in neuroimaging literature to threshold connectivity to remove weakly connected edges, as this can greatly influence the implied topology of the graph. In our work, we chose not to apply further thresholding, as unlike conventional graph theoretic metrics, linear models of spread and consequently network eigenmodes are relatively insensitive to implied topology induced by presence (or lack) of weak nonzero connections. However, to determine the geographic location of an edge, the top 95% of nonzero voxels by streamline count were computed for both edge directions. The consensus edge was defined as the union between both post‐threshold sets.

MEG processing and source reconstruction

MEG recordings were down‐sampled from 1,200 Hz to 600 Hz, then digitally filtered to remove DC offset and any other noisy artifact outside of the 1 to 160 Hz bandpass range. Since MEG data are in sensor space, meaning they represent the signal observable from sensors placed outside the head, this data needs to be “inverted” in order to infer the neuronal activity that has generated the observed signal by solving the so‐called inverse problem. Several effective methods exist for performing source localization (Jerbi et al., 2004; Wipf, Owen, Attias, Sekihara, & Nagarajan, 2010; Zumer, Attias, Sekihara, & Nagarajan, 2008). Here we eschew the common technique of solving for a small number of discrete dipole sources which is not fully appropriate in the context of inferring resting state activity, since the latter is neither spatially sparse not localized. Instead, we used adaptive spatial filtering algorithms from the NUTMEG software tool written in house (Dalal et al., 2004) in MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States). To prepare for source localization, all MEG sensor locations were co‐registered to each subject's anatomical MRI scans. The lead field (forward model) for each subject was calculated in NUTMEG using a multiple local‐spheres head model (three‐orientation lead field) and an 8 mm voxel grid which generated more than 5,000 dipole sources, all sources were normalized to have a norm of 1. Finally, the MEG recordings were projected into source space using a beamformer spatial filter. Source estimates tend to have a bias towards superficial currents and the estimates are more error‐prone when we approach subcortical regions, therefore, only the sources belonging to the 68 cortical regions were selected to be averaged around the centroid. Specifically, all dipole sources were labeled based on the Desikan–Killiany parcellations, then sources within a 20 mm radial distance to the centroid of each brain region were extracted, the average time course of each region's extracted sources served as empirical resting‐state data for our proposed model.

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 {. 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 closed‐form. 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 cost‐function evaluations, it will also accept solutions that have cost‐function evaluations greater than the previous one to move out of local minima. The acceptance probability function is , 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 iterations, where 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).

RESULTS

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 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 ). Each eigenmode produces a frequency response based on its frequency‐dependent eigenvalue (Figure 1d,e). Figure 1d shows the transit in the complex plane of a single eigenmode's frequency response, starting at low frequencies in the bottom right quadrant, and moving to the upper left quadrant at high frequencies. The magnitude, given by distance from origin, suggests that most eigenmodes have two prominent lobes, one roughly corresponding to lower frequency alpha rhythm and another corresponding to higher frequency beta or gamma rhythms, respectively. In contrast, the lowest few eigenmodes start off far from the origin, indicative of a low‐pass response. The magnitude of these complex‐valued curves shown in Figure 1e reinforces these impressions, 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 as demonstrated below.

The spatial patterns of the first 5 eigenmodes of , evaluated at the alpha peak of 10 Hz, are shown in Figure 1f. Eigenmodes produce posterior and temporal spatial patterns, including many elements of the default mode network; resembles the sensorimotor network; and 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 that different connectivity matrices produce different frequency responses: (a) the individual's structural connectivity matrix, (b) HCP average template connectivity matrix, (c) uniform connectivity matrix of ones, (d) a randomly generated matrix, (e) and (f) are randomly generated matrices with 75% and 95% sparsity respectively. For Figure 2a, optimized parameters for the individual subject's connectome were used. For Figures 2b–f, parameters optimized 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 low‐frequency 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.

Figure 2

Spectral graph model predictions of MEG spectra for one representative subject. Top—Observed MEG power spectrum for each of the 68 parcellated brain regions. Average spectra for each brain region are shown in blue, and the average spectrum across all brain regions is shown in thick black curve. The subsequent rows show each eigenmode's spectral magnitude response with model parameters optimized to match the observed spectrum (, , , , and ). Left column shows each eigenmode's frequency response in a differently‐colored curve, while the right column shows the same information as a heatmap. (a) Model using subject's individual structural connectivity matrix. (b) Model using a template structural connectivity matrix obtained by averaging structural connectivity from 80 HCP subjects. (c) Model using uniform connectivity matrix of ones. (d) Model using randomized connectivity matrix with no sparsity. (e) Model using randomized connectivity matrix with 75% 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 () appear to modulate the lower frequency range, up to beta, and may be considered responsible for the diversity of spectra observed in the model

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 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.

Figure 3

Spectral graph model depicts MEG spectra across subjects. (a) The observed spectra and spectral graph model's simulated spectra for four representative subjects. Red and cyan curves illustrate source localized empirical average spectra and region‐wise spectra respectively, while black and blue curves illustrate modeled average spectra and region‐wise spectra respectively. (b) Averaged observed spectrum across subjects is shown in red. The average simulated model spectra summing the first two‐third eigenmodes with optimized parameters for individual subject's connectome is shown in black. Model spectrum with optimized parameters and the HCP template connectome is shown in purple. Model spectrum with average parameter values and individual subject's connectome is shown in golden green. Model spectrum with optimized parameters and a distance connectome with 80% sparsity is shown in blue. Model spectrum with optimized parameters and symmetric random connectomes with 80% sparsity is shown in green. Finally, model power spectrum estimated by a neural mass model (NMM) with each subject's optimized global parameters and a HCP template connectome is shown in pink

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 (black triangles) and (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.

As another benchmark for comparison, a nonlinear neural mass model (Muldoon et al., 2016; Wilson & Cowan, 1973) using our in‐house MATLAB implementation (Xie et al., 2018), was generally able to produce characteristic alpha and beta frequency peaks, but this model does not resemble empirical wideband spectra. Note that no regionally varying NMM parameters were used in order to achieve a proper comparison with our model, but both models were optimized with the same algorithm.

Figure 4a shows violin plots of the optimized values, indicating that there is a large range of individually optimal model parameters across subjects. The time constants 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, one curve for each subject, indicates substantial improvement in cost function from default choice as optimization proceeds. The distribution of optimized spectral correlations is shown in 4C. Other model choices were evaluated for comparison: SGM on random connectomes with 80 and 95% sparsity, with and without a distance effect described in Methods, and SGM applied with average optimal model parameters instead of individually optimized ones. In order to test for significance, Fisher's R to z transform was applied, followed by a paired t‐test across all subjects between the optimal SGM and other models. The spectral fits for an SGM model with individual connectomes were significantly better than SGM models with average parameters, no matter what connectomes were chosen (p < 0.001). Interestingly, spectral fits for SGM model were comparable across all connectomes (p > 0.05). Furthermore, spectral fits for the SGM model were significantly better than that for NMM models with optimized parameters and individual connectome (p < 1e‐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.

Figure 4

Spectral graph model parameter optimization improves spectral fits. (a) Distribution of optimized model parameter values across all 36 subjects for the set of parameters { are shown in violin plots with each dot representing one subject. (b) Performance of optimization algorithm. Spectral Pearson correlation between model and source localized MEG spectra at each iteration. Each curve shows the spectral correlation achieved by the model optimized for a single subject, averaged over all regions, with increasing mean correlation values until the algorithm convergence to a set of optimized parameters. (c) Distribution of spectral correlations between optimized model and observed spectra across subjects. Correlations with optimized parameters are shown in the left three columns with individual connectomes (black), symmetric random connectomes with 80% and 95% sparsity (green) and geodesic distance‐based connectomes (blue). Correlation with average parameter values and individual connectomes are shown in golden green. Spectral correlations are highest for the SGM model with optimized parameters and the individual subject specific connectome when compared to SGM model with average parameters, regardless of the connectome and with an optimized NMM model, as denoted by asterisk (p < .001)

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 “sorted‐summed” 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.

Figure 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

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. 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.

Figure 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

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: (a) spectral graph model (SGM) on subject specific individual connectomes (CInd, black); (b) SGM with the HCP template connectome (CHCP, purple); (c) SGM on random connectomes with 80% sparsity comparable to individual connectomes or with 95% sparsity where the model shows spectral diversity (CRdm, green); (d) SGM on geodesic distance based connectomes (CDist, blue); and (e) a Wilson–Cowan neural mass model (NMM) with subject specific individual connectome (CInd, pink). Analogous results for beta band spatial correlations are contained in 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 CInd, 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 CInd and all other models show that, the SGM with CInd significantly outperforms all other models, as determined by a paired t‐test; p < .012 in each case, denoted by asterisk.

Figure 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 (CInd, black), random connectomes with 80% sparsity comparable to individual connectomes and 95% sparsity where the SGM model eigenmodes show spectral diversity (CRdm, blue), geodesic distance based connectomes (CDist,green) and for a neural mass model (NMM) with subject specific individual connectome (CInd, 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 CInd and all other models show that, the SGM with CInd 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

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 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 macroscopic dynamics. This model of the structure–function relationship has the following key distinguishing features. (a) Hierarchical: the model's complexity depends on the level of hierarchy being modeled: complex, nonlinear and chaotic dynamics can be accommodated at the local level, but linear graph model is sufficient at the macro‐scale. (b) Graph‐based: Macroscopic dynamics is mainly governed by the connectome, hence linear approximations allow the steady‐state frequency response to be specified by the graph Laplacian eigen‐decomposition, borrowing heavily from spectral graph theory (Auffarth, 2007; Kondor, 2002; Larsen et al., 2006; Ng & M. Jordan YW., 2002). (c) Analytic: The model is available in closed form, without the need for numerical simulations. (d) Low‐dimensional and parsimonious: Simple, global and universal rules specified with a few parameters, all global and apply at every node, are able to achieve sufficiently complex dynamics. The model is incredibly easy to evaluate, taking no more than a few seconds per brain and to infer model parameters directly from a subject's MEG data. The optimized model matches observed spectral and spatial patterns in MEG data quite well. No time‐consuming simulations of coupled neural masses or chaotic oscillators were needed; indeed, the latter greatly underperformed our model. We report several novel findings with potentially important implications, discussed below.

Recapitulating regional power spectra at all frequencies

Our main result is the robust demonstration of the model on 36 subjects' MEG data. The representative examples shown in Figures 3, 4, 5, 6 indicate that the graph model recapitulates the observed source localized MEG power spectra for the 68 parcellated brain regions, reproducing the prominent alpha and beta peaks. For each region, the model is also able to predict some characteristics of the full bandwidth power spectra, including what appears to be an inverse power law fall‐off over the entire frequency range of interest. However, this aspect will be quantitatively characterized in future work.

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 , driven by identically distributed random noise , 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ć et al., 2014; Mišić 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.

Relationship to other work

One can view the proposed generative model as a biophysical realization of a dynamic causal model (DCM) (Daunizeau, David, & Stephan, 2011; Friston, Preller, Mathys, et al., 2019; Pinotsis et al., 2017; Pinotsis, Hansen, Friston, & Jirsa, 2013; Razi, Kahan, Rees, & Friston, 2015) for whole brain electrophysiological activity but with very different goals, model dimensionality and inference procedures.

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 connectome—whether 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.

AUTHOR CONTRIBUTIONS

Ashish Raj conceived the study, wrote model algorithm, implemented software, contributed to data analyses, wrote and edited/improved the manuscript. Chang Cai contributed to experimental data acquisition and analysis, software implementation, and edited/improved manuscript. Xihe Xie contributed to experimental data acquisition and analysis, software implementation, wrote and edited/improved manuscript. Eva Palacios and Julia Owen contributed to experimental data acquisition and analysis and edited/improved manuscript. Pratik Mukherjee conceived the study and edited/improved manuscript. Srikantan Nagarajan conceived the study, contributed to data analyses, wrote and edited/improved manuscript.

Table S1 Parameters values and limits

Click here for additional data file.