Connectome harmonics predict resting state networks

To define the extension of the Laplace operator applied to the human connectome, we utilized its discrete counterpart, the graph Laplacian. We first created a graph representation for each of 10 human connectomes by combining the cortical surface extracted from MRI data of 10 subjects with the long-range (white matter) cortico-cortical and thalamo-cortical connections generated from cortical fibre tracts derived from diffusion tensor imaging (DTI) data of the same subjects (Fig. 1b). For each subject, we formed the graph representation of the modelled connectome, where the nodes −n being the total number of nodes—uniformly sample the curved anatomy of the cortical surface and the edges incorporate both local and long-range cortico-cortical and thalamo-cortical connections. By computing the graph Laplacian on the introduced representation , we define the discrete counterpart of the Laplace operator applied to the human connectome, the connectome Laplacian, for each individual and estimate its eigenvalue–eigenvector pairs , , connectome harmonics (see Methods). It is worth noting that by harmonics here we refer to spatial harmonics—as opposed to temporal harmonics. We will see later that these spatial harmonics can emerge from neuronally plausible dynamics but, at this stage, we are basing our Laplace eigenfunctions on static structural connectomes. It is also important to note that the introduced graph differs from previous graph representations of the human connectome used in population models10 in that it does not incorporate any parcellation of the cortical surface and cortico-cortical and thalamo-cortical connections. Thus, it provides a densely sampled connectome model with the minimum amount of discretization possible in the given resolution of the MRI and DTI data. Notably, when the number of uniformly sampled data points taken from the underlying manifold, such as the cortical surface, increases, the graph Laplacian converges to its continuous counterpart, the Laplace–Beltrami operator—the generalization of the Laplacian to non-euclidean geometries such as the curved anatomy of the human cortex25.

In the literature, Laplace eigenvalue–eigenvector pairs (eigenmodes) have received significant attention initially due to their relation to the excitation spectrum of a given geometry: the eigenvalues relate to the natural frequencies, the allowed frequencies of standing waves emerging on that particular geometry, whereas the eigenvectors yield the associated wave patterns18,23. Recent studies demonstrate also the relevance of Laplace eigenmodes for other physical phenomena including the phase extraction of an electron wave function18, the patterns emerging in electromagnetic interactions of ions15,16 and morphogenesis22,26 (Fig. 1a). Here, we utilized the eigenvectors of the connectome Laplacian to describe the spatio-temporal patterns of macroscale neural activity.

We found that the eigenvectors of the connectome Laplacian, the connectome harmonics, yield frequency-specific spatial patterns across distributed cortical regions (Fig. 1c; Supplementary Fig. 1). Here the red–blue patterns represent examples from the first 20 connectome harmonics in ascending order of frequency (wavenumber—shown in left) for one representative subject.

To quantitatively evaluate any similarity between the connectome harmonics and seven RSNs (Fig. 2a), commonly observed in the human brain; default mode, control, dorsal attention, ventral attention, visual, limbic and somato-motor network, we measured the mutual information between the reference RSNs27 and each of the connectome harmonic patterns (Fig. 2b: filled, coloured data points). As a control, we performed Monte Carlo simulations (2,000 simulations per subject randomly combined to 500,000 simulations for group average), in which we randomized the long-range white-matter cortico-cortical and thalamo-cortical connections while preserving the local anatomical structure of the subject’s cortical surface and computed the harmonics of each randomized network using the same methodology (Fig. 2b: black data points). Crucially, we found statistically significant similarity between the default mode network (DMN) and one particular connectome harmonic (in the range of the 9th connectome harmonic, ±2 due to individual differences) in the group average (Fig. 2b, top row, *P<0.0002, **P<0.0001; 500,000 Monte Carlo simulations for group average, corrected for multiple comparisons) as well as for all 10 subjects individually (Supplementary Fig. 2, P-values: *P<0.05, **P<0.02; 2,000 Monte Carlo simulations per subject, corrected for multiple comparisons). For all other resting state networks, we observed larger individual differences in the mutual information values (Supplementary Figs 2 and 3), while the Monte Carlo simulations still yielded significant similarity for different ranges of the connectome harmonics in the group analysis (Fig. 2b). In particular, we found that the visual, somato-motor and limbic networks showed significant similarity to only low-frequency connectome harmonics, while networks associated with higher cognitive functions, that is, control, dorsal attention and ventral attention, significantly matched a broad range of connectome harmonics distributed over the spatial frequency spectrum.

Figure 2: Prediction of the RSNs by connectome harmonics. (a) Patterns of synchronous oscillations, i.e., the RSNs, of the human brain overlap with the established functional systems, i.e., groups of cortical regions, which coactivate during certain tasks6,7. For quantitative evaluation of any similarity between the RSNs and connectome harmonics, we use the seven RSNs (default mode, control, dorsal attention, ventral attention, visual, limbic and somato-motor networks) (shown in a) identified from 1,000 subjects’ intrinsic functional connectivity data27. (b) Similarity measured by mutual information and (c) predictive power measured by F-measure values between the connectome harmonics with 40 lowest frequencies and the reference RSNs in ref. 27 (shown in a) compared with those of randomized harmonics (*P<0.0002, **P<0.0001 estimated by Monte Carlo simulations with 2,000 simulations per subject and 500,000 simulations for group average, after multiple comparison correction by false discovery rate, error bars indicate standard error across 10 subjects). Full size image

To further assess the predictive power of connectome harmonics for the resting state networks, we utilized an information retrieval metric, F-measure, which simultaneously quantifies the recall and precision of connectome harmonics’ prediction of the RSNs. We computed F-measure values between the connectome harmonics (after applying a binary indicator function) and the binary patterns of the RSNs27 and compared it with the F-measure values of randomized harmonics derived from Monte Carlo simulations. The F-measure provides a stricter evaluation than the information theoretic mutual information, as the indicator function imposes a binary decision to each node. This evaluation (Fig. 2c; Supplementary Figs 4 and 5) confirmed our previous findings based on mutual information. Taken together, these results suggest that the visual, somato-motor, and limbic, as well as the default mode network, are well predicted by individual connectome harmonics of a narrow frequency range, whereas the higher cognitive networks rely on a broader frequency range of connectome harmonics distributed over the spatial spectrum.

Connectome harmonics as the Fourier basis on the connectome

Next, we investigated the use of connectome harmonics as a function basis to represent spatial patterns of cortical networks. The orthogonality of connectome harmonics means that a linear combination of these eigenfunctions can be used to recreate any spatial pattern of neural activity. The importance of the connectome harmonic basis lies in its close ties to the classical Fourier transform, which corresponds to the decomposition of a signal into a linear combination of the eigenfunctions of the Laplace operator applied to a circular domain, that is, sine and cosine functions with different frequencies23. Since connectome harmonics are defined as the eigenvectors of the connectome Laplacian—the discrete counterpart of the Laplacian applied to the human connectome—they extend the Fourier basis to the particular geometry of the human connectome. Hence, the spectral transform onto the connectome harmonic basis naturally extends the classical Fourier transform to the human connectome.

To analyse the spatial frequency content of the resting state networks27 (Fig. 2a), we performed a spectral transform to the connectome harmonic basis and reconstructed the spatial patterns of individual networks. Although the binary nature of the reference networks theoretically necessitates the use of the whole connectome harmonic spectrum for reconstruction—the same way that a square wave can only be reconstructed using the sine waves with infinite many frequencies—sharp decreases in the normalized reconstruction errors were observed by using just 0.1% of the connectome harmonic spectrum (low-frequency range; Fig. 3a). The steepest decrease of the DMN’s reconstruction error occurred for the frequency band that also showed significant similarity and predictive power in mutual information and F-measure values, respectively (highlighted by the red column in Fig. 3a, best matching connectome harmonic of each subject shown in Supplementary Fig. 6, reconstruction shown in Fig. 3c) while for the visual, somato-motor and limbic networks the decrease of the reconstruction error remained large but constant within 0.1% of the spectrum (Fig. 3a). Slower convergence was observed for the reconstruction errors of higher cognitive networks within the range of 1.2% of the spectrum (low-frequency range) suggesting the reliance of these networks on a broader range of frequencies (Fig. 3b). These results confirm our previous findings while providing a novel analytical language of cortical activity analogous to the classical Fourier transform that can be utilized to quantify any activity pattern including task-based event-related designs.

Figure 3: Reconstruction of the RSNs from connectome harmonic basis. Normalized reconstruction error of each resting state network using (a) 0.1% and (b) 1.2% of the connectome harmonics spectrum averaged across 10 subjects (error bars and shading indicate standard error across 10 subjects) compared to the reconstruction of a randomized binary pattern. Red band in a highlights the steepest decrease for the DMN. (c) Reconstruction of the DMN using (from top to bottom) 5, 0.5 and 0.05% of the spectrum and the best matching connectome harmonic of one subject’s data. Full size image

Biological mechanisms underlying connectome harmonics

We also investigated the biological mechanisms likely to underlie the self-organization of connectome harmonics on the cortex. Hitherto, we have assumed that the graph Laplacian based upon structural connectivity provides a plausible proxy for the effective connectivity (also known as the Jacobian-please see ref. 28) of an underlying neuronal dynamical system. In what follows, we explore the biological mechanisms likely to underlie these neuronal dynamics. In particular, we exploit the efficient description of these dynamics given by neural field equations in terms of connectome harmonics and show how they give rise to the emergence of connectome harmonics on the cortex.

The dynamics of the oscillatory cortical networks is thought to emerge from the interplay of excitation; for instance mediated by glutamatergic principal cells, and inhibition; for instance mediated γ-aminobutyric acid GABAergic interneurons29. To describe macroscale cortical dynamics, we extend a variant of neural field models based on Wilson–Cowan equations, the most commonly utilized mathematical description of the excitatory–inhibitory neural dynamics30. Neural field models based on Wilson–Cowan equations are a variant of reaction-diffusion systems31 originally introduced by Turing as a mathematical model for morphogenesis19. Based on the principle that mutual interactions between a diffusing activator and inhibitor can result in self-organizing pattern formation, reaction-diffusion models have provided valuable insights into the mechanisms underlying the emergence of non-linear waves in several biological processes including morphogenesis19,20, formation of ocular dominance patterns in the visual cortex32, visual hallucinations33,34 and interactions of excitatory and inhibitory activity in neural populations of the cortical and thalamic tissue30,35. Crucially, the pattern formation in reaction-diffusion systems is caused by the exponential growth of certain eigenfunctions of the Laplacian applied to the patterning domain (see refs 21, 22 and Supplementary Notes 1–4). The selection of which harmonic patterns are ‘activated’ (grow) is determined by the diffusion parameters of excitation and inhibition (Supplementary Fig. 7). Hence, the Laplace eigenfunctions provide the building blocks of complex patterns emerging in reaction-diffusion systems.

On a two-dimensional continuous idealization of the cortex, the Wilson–Cowan equations lead to self-organization of neural oscillatory patterns when short-range excitation is coupled with broader lateral inhibition36. This type of functional circuitry, known as ‘Mexican hat’ organization37 or centre-on and surround-off connectivity38 is well-observed experimentally in the early visual cortex, known as cortical surround suppression39,40 and is likely to extend throughout the neocortex41,42.

Recent experimental evidence showed that a plausible mechanism underlying the cortical surround suppression in V1 is the activity of the somatostatin-expressing inhibitory neurons in the superficial layers of the mouse visual cortex and a similar neural circuit is also likely to underlie surround suppression in other cortical areas42. These findings are further supported by the report of broader spatial extent of inhibition compared with excitation in the primary visual cortex of awake mice43. Furthermore, layer-specific suppression and facilitation also generates the necessary circuits for lateral inhibitory interactions29,44: coordinated modulation of superficial (L2/3) and deep cortical layers (L5) gives rise to competition between neighbouring domains and lateral inhibition, although the spatial extent of excitation and inhibition across cortical domains show overlapping distributions vertically (across cortical layers) and horizontally (within layers)44. Notably, it has been shown that the Mexican hat type of functional circuitry can also be generated in anatomical circuits with short-range inhibition when a fraction of the total excitatory conductance is slower than the inhibition37. A potential biological source of slower excitatory conduction is the slow synaptic transmission caused by the N-methyl-D-aspartate receptors contributing mostly to excitatory currents37.

Taken together, this converging evidence suggests that various biological mechanisms can give rise to the functional circuitry equivalent to short-range excitation coupled with broad inhibition, which indeed is the necessary condition for the self-organization of oscillatory patterns in Wilson–Cowan-type neural field models30,33,35 (Supplementary Notes 1–5). Next, we extend the Wilson–Cowan equations to the full structural connectivity of the thalamo-cortical system by incorporating the connectome Laplacian and demonstrate the relation between the emerging oscillatory patterns and the connectome harmonics.

Neural field model of connectome-wide neural dynamics

We extended a variant of the neural field model35 based on the Wilson–Cowan equations30 to the three-dimensional connectome model by incorporating the connectome Laplacian into the spatial propagation (diffusion) term (Fig. 4a; Supplementary Notes 2 and 3). Numerical simulations were then performed by combining the network diffusion on the human connectome—modelled by iterative application of symmetric graph Laplacian23,45 —with the excitatory–inhibitory reaction dynamics (Fig. 4a). This allows us to extend continuous form neural field models by incorporating the connectivity of the human connectome. This neural field approach differs from previous macroscale simulations10,46 in that spatial propagation is modelled by network diffusion, as opposed to discrete coupling, yielding a (spatially) near-continuous model of cortical dynamics.

Figure 4: Neural field model. (a) Left: dynamics of excitatory (E) and inhibitory (I) activity. Right: time evolution of the excitatory E(v i , t) and inhibitory I(v i , t) activity at the cortical location at time t where , , and describe the diffusion processes of E and I activity acting on excitatory (EE, IE) and inhibitory (EI, II) neural populations and τ s refers to the units of system time, that is, characteristic time scale. (b) Linear stability analysis of the neural field model in terms of connectome harmonics. The red regions correspond to the diffusion parameters in the phase space that algebraically satisfy the necessary condition for oscillations, that is, the critical Hopf regime, plotted as a function of the analysed diffusion parameter vertical axis and the eigenvalue of the connectome harmonic horizontal axis. (c) Power spectrum of the temporal oscillations in (a total of 267) numerical simulations averaged over all nodes. (d) Spatial pattern for an arbitrary time slice and the temporal profile of four seed locations shown in e. (f) Seed-based correlation analysis of the neural field patterns demonstrates the decoupling between the posterior and anterior midline nodes of the DMN for the same set of parameters leading to slow cortical oscillations. Full size image

We found that structured oscillatory patterns naturally self-organize on the human connectome for a wide-range of diffusion parameters in the model (Fig. 4a–d; Supplementary Figs 8 and 9; Supplementary Movies 1–4). Linear stability analysis of the neural field model revealed that a wide range of connectome harmonics could be activated for different diffusion speeds of excitation and inhibition (Fig. 4b; Supplementary Notes 3–5), rendering the neural field model a plausible neural mechanism for the self-organization of connectome harmonics. In particular, we observed a decrease in the frequency of coherent oscillations when excitatory activity is decreased, modelled by slower diffusion of excitation, or when inhibitory activity is increased, modelled by faster diffusion of inhibition in the neural field model (Fig. 4c,d). This relationship between the frequency of temporal oscillations and the excitation–inhibition balance shows remarkable overlap with the neurophysiological changes observed during the loss and recovery of consciousness (for the analysed parameter range). Neurophysiological evidence suggests that drug- or sleep-induced loss of consciousness is associated with increasing inhibitory or decreasing excitatory activity, which is accompanied by a transition from the low amplitude, high-frequency patterns to low-frequency coherent oscillations in cortical activity47. Recent work also shows gradual decoupling between the posterior and anterior midline nodes of the DMN during loss of consciousness48,49. We observed this decoupling in seed-based correlation analysis of the neural field patterns for the exact parameters, which resulted in slower cortical oscillations (Fig. 4e,f; Supplementary Fig. 9).

Finally, we tested the stability of the emerging oscillatory patterns to external perturbations such as noise using Lyapunov stability analysis50. This method involves perturbing the system at some time t* and observing whether the perturbed system converges to the original system. In dynamical systems, such as described by neural field models, the eigenvalues of the linearized system around a known fixed point immediately reveals the local behaviour of the system; that is, it would allow one to identify the unstable ‘growing’ eigenmodes that dominate observed fluctuations (Supplementary Notes 3–5). However, due to the nonlinearities inherent in the Wilson–Cowan equations as well as the high dimensionality of the modelled system, the continuum of the fixed points, that is, the trajectory, was not readily obtainable analytically. Therefore, we numerically integrated the system to determine the stable solution, in our case a periodic solution, for certain parameters. We first perturbed the system separately 10 times by white noise at time t* and analysed the largest difference between the unperturbed trajectory and each of the perturbed trajectories over time by computing the L-infinity norm L(t) (see Methods). Fig. 5a shows that L(t) is bounded and converges to ɛ≈0, whereas Fig. 5b,c illustrate the convergence of the limit cycle and the temporal oscillations of the perturbed system to that of the unperturbed system for two example vertices v 1 and v 2 . These results demonstrate that for the analysed parameter range the extended Wilson–Cowan equations (Fig. 4a) are robust to external perturbations such as noise.

Figure 5: Stability analysis of the neural field model. Stability of the emerging oscillations to external perturbation is tested by perturbing the neural field model for a sample oscillatory parameter set ( , , and ). (a) The distance measure L(t) versus the time from perturbation t−t* where t* denotes the time of perturbation by white noise. Before the perturbation L(t) is identical to 0. L(0) shows how far the perturbed system altered from the underlying oscillatory stable state. After perturbation, that is, when t−t*>0, L(t) approaches 0 showing the system is Lyapunov stable. (b) The convergence of the limit cycle; that is, excitatory activity of two non-adjacent nodes plotted over time, of the perturbed (black solid line) to unperturbed (orange dashed line) system. (c) The original (unperturbed) trajectory (dashed line) and the perturbed trajectory (solid line) for two sample vertices (green and blue lines) after perturbation t−t*∈[0, 5] unit time. The perturbed trajectory converges back to the original trajectory (with a small phase shift) demonstrating that the state corresponding to the original trajectory is Lyapunov stable to small perturbations. Full size image

Taken together, the extension of the Wilson–Cowan-type neural field models to the particular structural connectivity of the human connectome provides a biologically plausible, robust mechanism likely to underlie the self-organization of connectome harmonics in the thalamo-cortical system.