Gradients of connectivity as graph Fourier bases of brain activity.

The application of graph theory to model the complex structure and function of the brain has shed new light on its organization, prompting the emergence of network neuroscience. Despite the tremendous progress that has been achieved in this field, still relatively few methods exploit the topology of brain networks to analyze brain activity. Recent attempts in this direction have leveraged on the one hand graph spectral analysis (to decompose brain connectivity into eigenmodes or gradients) and the other graph signal processing (to decompose brain activity "coupled to" an underlying network in graph Fourier modes). These studies have used a variety of imaging techniques (e.g., fMRI, electroencephalography, diffusion-weighted and myelin-sensitive imaging) and connectivity estimators to model brain networks. Results are promising in terms of interpretability and functional relevance, but methodologies and terminology are variable. The goals of this paper are twofold. First, we summarize recent contributions related to connectivity gradients and graph signal processing, and attempt a clarification of the terminology and methods used in the field, while pointing out current methodological limitations. Second, we discuss the perspective that the functional relevance of connectivity gradients could be fruitfully exploited by considering them as graph Fourier bases of brain activity.

INTRODUCTION

Modern attempts at understanding brain function have leveraged the use of graph theory to grasp complex properties of neuronal networks, giving rise to the field of network neuroscience (Bassett & Sporns, 2017; Sporns, 2017). Modeling brain organization using graphs has led to fascinating results, such as the brain’s hypothetical rich-club organization (van den Heuvel & Sporns, 2011), the cortical organization in functionally relevant modules (Sporns & Betzel, 2016), as well as common wiring principles across species (Goulas, Majka, Rosa, & Hilgetag, 2019). Despite the tremendous progress that has been achieved in network neuroscience, surprisingly relatively few methods such as graph signal processing (GSP; Shuman, Narang, Frossard, Ortega, & Vandergheynst, 2013) exploit brain connectivity (i.e., take into account the topology of brain networks) to characterize brain activity (Ju & Bassett, 2020).

First steps in the direction of exploiting connectivity graphs in the analysis of brain signals have been made using spectral graph theory (Chung, 1997; Mhaskar, 2016). The underlying idea of this theory is to interpret the eigenvectors of graph Laplacians as harmonic components. Increasing evidence of the functional relevance of these spectral components of brain networks (i.e., , eigenmodes, or harmonics; Atasoy, Donnelly, & Pearson, 2016; Belkin & Niyogi, 2003; W. Huang et al., 2016; see Table 1), has recently been shown with a variety of approaches (Atasoy, Deco, Kringelbach, & Pearson, 2017; Margulies, Ghosh, Goulas, Falkiewicz, & Huntenburg, 2016; M. B. Wang, Owen, Mukherjee, & Raj, 2017).

 

Graph spectral analysis and GSP applied to neuroscience: Terminology. White: Based on spectral graph theory (spectral decomposition of brain networks). Light blue: Using GSP (graph Fourier transform: Spectral decomposition of a brain signal based on the underlying brain graph topology).

	Approach	Toolbox	Key references
Connectome Harmonics, Harmonic Brain Modes	Laplacian eigenvectors ui of large-scale brain networks estimated from DWI and anatomical MRI. They are interpreted as spectral components of spatiotemporal neural activity and compared with resting-state networks and oscillatory patterns (neural field model).		(Atasoy, Deco, et al., 2017; Atasoy et al., 2016; Atasoy et al., 2017)
Brain Gradients	Eigenvectors ui obtained applying diffusion map embedding on large-scale functional (Margulies et al., 2016), microstructural (Paquola et al., 2019), or spontaneous oscillation (Mahjoory et al., 2019) networks estimated respectively from resting-state fMRI data, myelin-sensitive MRI data, or MEG signals. They reveal macroscale axes of cortical organization with functional and neurodevelopment relevance.	BrainSpace	(Hong et al., 2019; Margulies et al., 2016; Mckeown et al., 2020; Paquola et al., 2020; Paquola et al., 2019)
Connectopies	Laplacian eigenvectors ui of the graph obtained computing the correlation between voxels within a selected ROI and the rest of the gray matter voxels. The approach reveals fine-grained topographic organization of a brain region’s connectivity (i.e., primary motor or visual cortex).	Congrads	(Haak, Marquand, & Beckmann, 2018; Schröder et al., 2015)
Brain Activity Eigenmodes	Excitatory or inhibitory neural activity expanded in terms of spatial eigenmodes of the cortex mesh obtained solving corticothalamic neural field theory equations. These brain eigenmodes show high similarity with spherical harmonics (cortical folding = 0) and DWI connectivity eigenvectors (graph Laplacian).		(Robinson et al., 2016)
Fourier/Harmonic Modes	Graph Fourier modes obtained applying graph Fourier transform to a signal (i.e., fMRI, W. Huang et al., 2016; Preti, Bolton, & Van De Ville, 2017; or EEG, Glomb et al., 2020) on a graph (i.e., structural connectivity graph estimated from DWI). This analysis reveals low (high) frequency modes that are aligned (liberal) with respect to the underlying graph structure.	PyGSP	(Glomb et al., 2020; B. W. Huang et al., 2018; Medaglia et al., 2017; Preti & Van De Ville, 2019)
Graph Neural Fields	Excitatory or inhibitory neural activity expressed as stochastic neural field equations on the human connectome graph. This approach combines Wilson-Cowan neural field equations and graph signal processing to model and analyze whole-brain activity.	Code available†	(Aqil et al., 2020)
Notes	†https://github.com/marcoaqil/Graph-Stochastic-Wilson-Cowan-Model

GSP, by modeling attributes of network nodes as signals onto a graph (Ortega, Frossard, Kovacevic, Moura, & Vandergheynst, 2018), takes a step further as it allows a joint analysis of brain activity and connectivity. The emergence of GSP is mostly due to the elegant and powerful analogy between graph Laplacian eigenvectors and classical Fourier analysis (Girault, 2015) and the possibility of decomposing a signal “living on a graph” as a combination of spatial harmonics. Recent works have exploited GSP to decompose brain activity in graph Fourier modes, with encouraging results (Glomb, Queralt, Pascucci, & Tourbier, 2020; B. W. Huang et al., 2018; Medaglia et al., 2017; Preti & Van De Ville, 2019).

The gradients/GSP framework is complementary to the classical approach of mapping cortical area functions (brain parcellations) and discrete networks associated with a particular condition or task. Exploiting the topology of brain networks, this framework allows for a decomposition of brain activity or structure as a continuum of spectral components, better describing subregional heterogeneity and multiplicity than parcellation approaches, which consider uniform brain regions (nodes or parcels; Haak & Beckmann, 2020). This novel approach has been successfully applied to the analysis of the healthy or pathological brain in an increasing number of studies (Glomb et al., 2020; Medaglia et al., 2017; Mortaheb et al., 2019; Preti & Van De Ville, 2019). However, readers that first approach the recent literature are often confronted with different notations, terminology, and methods (e.g., Laplacian embedding, diffusion maps, graph Fourier modes; see Table 2) that sometimes are not consistent between the GSP or “gradients” communities and may be difficult to unravel.

 

Graph spectral analysis and GSP applied to neuroscience: Methods. For each methodology, a brief description of the approach and main references are provided (in bold, the key references describing the algorithm in detail). White: Based on spectral graph theory (spectral decomposition of brain connectivity graphs). Light blue: Using GSP (graph Fourier transform: Spectral decomposition of a brain signal based on the underlying brain graph topology). Notations: .

	Notation	Method	Key references
Laplacian Embedding or Eigenmaps	Operator	Spectral decomposition of a graph G in eigenvectors of the graph Laplacian L. It is the discrete counterpart (on graph) of the Laplacian-Beltrami operator on continuous manifolds. The Laplacian eigenvectors associated with the lowest eigenvalues provide a dimensionality-reduction mapping that preserves locality.	(Haak et al., 2018; Schröder et al., 2015) (Belkin & Niyogi, 2003)
	L = D − W Laplacian
	Ln = D−1/2LD−1/2 Normalized Laplacian
	Lrw = D−1L Random Walk Laplacian
	Low-dimensional representation
	G = [u1, u2, u3, …]
Diffusion Maps	Operator	Diffusion map embedding treats the graph G as the basis of a diffusion process. Diffusion maps are a family of graph Laplacians that depend on a diffusion parameter α. They can be employed to embed the data into a Euclidean space where the probability of transition between nodes defines the Euclidean distance between the corresponding points in the embedding space.	(Hong et al., 2019; Huntenburg et al., 2018; Langs et al., 2016; Margulies et al., 2016; Mckeown et al., 2020) (Coifman & Lafon, 2006; Coifman et al., 2005)
	W(α) = D−1/αWD−1/α
	P(α) = D(α)−1W(α) Transition Probability
	Low-dimensional representation
	>G = [λ1Tu1, λ2Tu2, …]
Graph Signal Processing	Laplacian eigendecomposition	Expansion of a signal (or a stochastic function) X in terms of the eigenvectors of the underlying graph Laplacian L. Laplacian eigenvectors carry a notion of spatial frequency: i.e. eigenvectors corresponding to low eigenvalues vary smoothly across the graph; those corresponding to large eigenvalues have higher spatial frequencies (i.e. are more likely to have different values across adjacent nodes).	(Aqil et al., 2020; Glomb et al., 2020; B. W. Huang et al., 2018; Medaglia et al., 2017; Preti & Van De Ville, 2019) (Ortega et al., 2018; Shuman et al., 2013)
	L = UΛUT
	Graph Fourier transform of X
	X˜ = UTX

The goal of this paper is first to summarize recent contributions that have used connectivity gradients and GSP for neuroimaging. In particular, we aim to clarify terminology, compare different methodologies, and provide resources (see Box 1) and key references (Table 1 and 2). The second goal of this paper is to make a link between the two frameworks by discussing gradients as graph Fourier bases. We argue that using GSP for the analysis of multimodal neuroimaging data will pave the way to more interpretable analysis methods.

A variety of resources for graph analysis are available and extensively used in data science across multiple research domains. As graph structure and properties are encoded in matrices, in principle every toolbox manipulating arrays and matrices can be adapted to GSP and graph spectral theory; this includes for instance Python packages such as numpy, scikit, and pytorch. Here we list toolboxes and Python modules that were specifically designed or that can be easily applied for GSP and graph spectral analysis of brain networks and their visualization.

Nilearn: A widely used Python toolbox for machine learning in neuroimaging. It also includes useful functions for brain connectivity computation and visualization that can be easily adapted to plot gradients and signals on brain networks. https://nilearn.github.io/#

BrainSpace: A Matlab/Python software package that allows connectivity gradients computation and analysis specifically adapted to neuroimaging and connectome datasets (Vos de Wael et al., 2020). https://github.com/MICA-MNI/BrainSpace

PyGSP: A Python package specifically designed for graph signal processing that implements a variety of operations on graphs (computing graph Fourier transform, filtering or interpolating signals on graphs, plotting) that also scale to very large graphs. https://pygsp.readthedocs.io/en/stable/index.html

NetworkX: A Python package to analyze network structure, build network models, and visualize networks. https://networkx.github.io/

Congrads: A Python package to compute and map connectopies (connectivity topography) of a predefined region of interest. https://github.com/koenhaak/congrads

VB toolbox: A pair of packages (available in both MATLAB and Python) including connectivity gradient analysis pipelines and the computation of the Vogt-Bailey index (Bajada et al., 2020). https://github.com/VBIndex

GRADIENTS OF BRAIN CONNECTIVITY

Connectivity gradients are obtained from graph spectral decomposition of the connectivity matrix. As described in Box 2, they correspond to Laplacian eigenvectors u of connectivity graphs (the reader can also refer to Bajada et al. [2020] and Vos de Wael et al. [2020] for a detailed description and step-by-step tutorials of the “connectivity gradient analysis”). Connectivity gradients provide a representation of cortical organization as a continuum of spatial harmonics that can overlap in space (graph nodes) in the absence of hard boundaries. This approach is therefore complementary to a brain atlas that maps cortical areas onto a set of discontinuous functional or structural regions (brain parcellations; Eickhoff, Yeo, & Genon, 2018; Glasser et al., 2017; Thomas Yeo et al., 2011). A varied (and sometimes inconsistent) terminology has been used in relation to connectivity gradients, even though they are grounded in the same fundamental approach. While a comprehensive review of the growing gradients literature is outside the scope of this paper, Tables 1 and 2 aim to clarify terminology and describe corresponding methodologies, while providing some key references.

Graphs are tools that are ubiquitous in many fields of science, thanks to their generic nature and expressivity. They make it possible to efficiently represent relations between items, called nodes. These relations are modeled using edges, which are most often pairs of nodes. For example, if i and j are nodes, (i, j) is a potential edge in the graph. In many cases, relations are weighted with the convention that a weight 0 corresponds to the absence of an edge, and any nonzero value can be used otherwise. A concise way to represent a weighted graph consists of using its (weighted) adjacency matrix W, indexed by nodes. As such, W is 0 if and only if there is no edge between nodes i and j, and W represents the weight of the edge (i, j) otherwise. By summing a row of W, we obtain the strength (or weighted degree) of the corresponding node, which can be thought of as the importance it has, compared with the other nodes (Bollobas, 1998). These strengths can be arranged in a diagonal matrix, called the strength matrix D, which makes it possible to define the graph Laplacian: L = D − W. For an undirected graph with N nodes, the Laplacian is a real symmetric matrix and thus has an orthogonal basis of eigenvectors u associated with eigenvalues λ such as Lu = λu, with l = 0, 1, ⋯, N − 1. The graph Laplacian is key to many fundamental properties about graphs, which can be found in the literature about spectral graph theory (Spielman, 2012). For example, assuming W contains only nonnegative values, the spectrum of L is also nonnegative. It always contains the element 0. The magnitude of the second smallest eigenvalue of Lλ1 is an important indicator of the global connectivity of the graph; for example, λ1 = 0 if the graph is not connected (i.e., it is constituted of at least two separate subgraphs).

In the field of graph signal processing (GSP), authors are interested in manipulating graph matrices together with vectors indexed by the nodes of the graph, called graph signals. In the light of this formalism, it is possible to define ad hoc graph Fourier modes, tied to the specific considered graph structure. These Fourier modes simply consist of the eigenvectors of L. The corresponding eigenvalues λ exhibit behaviors that can be interpreted in terms of spatial frequencies over the graph structure (see Figure 1, panel D, for a schematic illustration of a signal on graph, graph spectrum, and examples of low-medium- and high-frequency graph Fourier modes, e.g., Laplacian eigenvectors). In the case of simple ring graph structures, graph Fourier modes and classical discrete Fourier modes become identical. For more complex and arbitrary graphs, the abundant GSP literature explains how to design filters and other operators adapted to the underlying structure (Shuman et al., 2013).

 

From graph theory to graph signal processing in brain imaging. (A) Different areas of the brain can be represented as nodes and structural and functional relationships between them as edges of a complex large-scale network, also known as the connectome (Sporns, Tononi, & Kötter, 2005). Various approaches exist to identify the nodes of the connectome (atlas and anatomic based, data-driven, etc.; Glasser et al., 2017; Thomas Yeo et al., 2011). Similarly, edges of a brain network can be assessed with a range of neuroimaging techniques (DWI, EEG, fMRI, MEG, PET) and methods (structural, functional, or effective connectivity; Friston, 2011). (B) Graph theory allows us to describe salient properties of network topology with matrices (i.e., adjacency, Laplacian, degree matrices, etc.) and graph indexes (i.e., efficiency, clustering, centrality; Bassett & Sporns, 2017; Bullmore & Sporns, 2009; Fornito, Zalesky, & Bullmore, 2016). (C) Graph spectral analysis (e.g., Laplacian eigenvectors) is used to extract low–dimensional representations of brain networks known as brain gradients (Margulies et al., 2016; see Tables 1 and 2). (D) Graph signal processing (GSP; Shuman et al., 2013) takes a step forward as it associates a signal with an underlying graph. It extends classical analysis methods from regular domains (discrete time signals) to nonregular graphs. GSP allows us to analyze brain activity taking into account the underlying topology of brain networks. GSP also allows for a spectral decomposition of brain activity based on the underlying graph Laplacian eigenvectors (graph Fourier transform; see Table 2). In the figure, a brain signal (whose amplitude is encoded in the height and color of the vertical bars) “lives” on a brain network (black) and can be decomposed in low (high) graph frequency harmonics corresponding to small (high) graph Laplacian eigenvalues. In this example the graph spectrum and corresponding Laplacian eigenvectors were obtained from the spectral analysis of an averaged structural graph from the Human Connectome Project.

In short, the graph Laplacian (and other versions of the Laplacian matrix, see Table 2) is ubiquitously used in both graph spectral analysis and GSP. In the first case, the properties of the Laplacian are typically exploited for dimensionality reduction of the graph (i.e., for the extraction of brain connectivity gradients from a brain network). In the case of GSP, Laplacian eigenvectors are used to decompose the signal in graph Fourier modes by defining a graph Fourier transform (see Table 2). It is worth clarifying that if these approaches are all grounded in a Laplacian-based spectral decomposition of a matrix, the adjacency matrix itself can be estimated with a variety of techniques (Figure 1, panel A). For instance, structural connectivity refers to anatomical connections between brain regions and is most commonly estimated using diffusion-weighted imaging (DWI). Functional connectivity is defined as the statistical dependence among measurements of neural activity (Friston, 2011) and it is usually inferred through correlations among neurophysiological time series. Effective connectivity estimates the influence that one neuronal system exerts on another and, because it refers to the notion of causality, is intrinsically directed. Functional and effective connectivity can be assessed by a range of neuroimaging techniques (electroencephalography; EEG, functional magnetic resonance imaging, fMRI; Magnetoencephalography, MEG).

See also Table 1 and 2 for more details on the terminology and methods of different approaches using Laplacian decomposition of brain signals.

Connectivity gradients estimated from functional (Margulies et al., 2016), structural (Robineau et al., 2017), or microstructural (Paquola et al., 2019) brain data have been applied to study the hierarchical organization of brain structure and function. The pioneering work of Margulies and colleagues (Huntenburg, Bazin, & Margulies, 2018; Margulies et al., 2016) introduced the concept of gradients to indicate eigenvectors obtained applying diffusion map embedding to a resting-state fMRI connectivity matrix. Interestingly, the first nontrivial eigenvector of the spectral decomposition (corresponding to eigenvalue λ1) revealed a macroscale cortical organization spanning from unimodal sensory areas to transmodal association areas. Atasoy and colleagues (Atasoy et al., 2016) obtained connectome harmonics through Laplacian decomposition of diffusion-weighted imaging (DWI) connectivity and showed that they can predict resting-state functional network activity. Results indicated that visual, sensory-motor, and limbic networks more closely matched low-frequency harmonics, while higher cognitive function networks spanned over a broader range of brain modes. Another recent work on brain structure/function relation (Park et al., 2021) compared structural connectome gradients and functional states estimated from fMRI using hidden Markov-autoregressive models. Results from Park and colleagues revealed a strong structure-function coupling for functional states anchored in sensorimotor areas, whereas functional states associated with transmodal areas more largely relied on rich-club nodes and long-range connectivity. In Paquola et al. (2019), functional gradients were shown to only partially align with microstructural gradients obtained from histology and high-field MRI: If both axes of variance originate in primary sensory areas, the functional gradient identified in Margulies et al. (2016) arches towards default mode and frontoparietal networks, while its microstructural counterpart extends to limbic cortices.

Earlier work applied Laplacian decomposition to the cortical surface mesh (instead of graphs of connectivity) to extract anatomically relevant features. In Germanaud et al. (2012), Laplacian-based spectral analysis was applied to the cortical mesh to investigate gyrification complexity. Using a segmentation of the cortex based on Laplacian spectral decomposition, the authors were able to identify developmentally relevant features such as primary, secondary, and tertiary folds. Spatial eigenmodes of the curved cortex were also used in Robinson et al. (2016) to solve neural field corticothalamic equations and estimate so-called activity eigenmodes of the brain (Deco, Jirsa, Robinson, Breakspear, & Friston, 2008). The authors also show that excitatory or inhibitory states can be reconstructed through a finite number of eigenfunctions, and that these largely overlap with the Laplacian eigenvectors of a connectivity matrix estimated from DWI. Interestingly enough, and in line with our attempt at giving a unified view of the gradients and GSP methodologies, this approach linking neural field equations and Laplacian decomposition has been revisited in a GSP perspective in a recent publication (Aqil, Atasoy, Kringelbach, & Hindriks, 2020 introducing the graph neural fields framework (see Table 1).

Different from the whole-brain approaches described hitherto, Haak and colleagues (Haak, Marquand, & Beckmann, 2018; Schröder, Haak, Jimenez, Beckmann, & Doeller, 2015) proposed a framework to map regional connection topography (i.e., connectopies) using Laplacian eigenmaps (Belkin & Niyogi, 2003). For every voxel in the selected region, a connectivity fingerprint was obtained computing correlation with the fMRI time series from all the other gray matter voxels; nonlinear manifold learning based on Laplacian decomposition of the connectivity matrix was then applied to extract corresponding eigenvectors named connectopies (i.e., connections topography). The authors underline the utility of using nonlinear dimensionality reduction based on connectivity (instead of linear approaches such as PCA and ICA) to identify more biologically plausible maps of functional organization. Following this perspective, recent work has investigated the relationship between connectivity gradients and cognitive processes, showing that these are altered depending on the ongoing cognitive experience or psychological state (Atasoy, Deco, et al., 2017; Ito, Hearne, & Cole, 2020; Lanzoni et al., 2020; Mckeown et al., 2020; X. Wang, Margulies, Smallwood, & Jefferies, 2020). In Mckeown et al. (2020), for example, macro-scale gradients of functional connectivity at rest are related to ongoing thoughts, and different gradient profiles are associated with different cognitive tasks (e.g., problem-solving vs. past-related thoughts).

Taken together, these studies underline the potential of data-driven dimensionality reduction based on brain networks to reveal principles of large-scale cortical organization and identify gradual changes in functional, white matter, and cytoarchitectonic architecture in different conditions. They also suggest that connectivity gradients can yield meaningful functional relevance, and thus it might be particularly sensible to use them as a basis to analyze brain activity using GSP.

GSP FOR NEUROIMAGING

The fundamental difference between GSP and graph theory is that while the latter provides tools to analyze and manipulate graphs, GSP focuses on analyzing signals on graphs. In other words, GSP leverages concepts developed in spectral graph theory to translate Fourier analysis to signals on graphs. A general overview of the GSP framework illustrating how to perform operations on graphs (spectral analysis, convolution, filtering) is presented in Ortega et al. (2018) and Shuman et al. (2013), and a first review on GSP for neuroimaging was recently proposed (W. Huang et al., 2017). Using a graph estimated from the white matter tracts of the brain, the authors detail how to apply graph Fourier transform to analyze brain activity, and exploit the corresponding graph frequency bands to interpret the data. Precisely, low frequencies correspond to brain activity that follows closely (i.e., is “aligned with”) the underlying white matter connectivity, while high frequencies are characterized by brain activity patterns that can be seen as “liberal” with respect to the network structure.

This approach provides a simple and interpretable framework, which has successfully been applied to study the decoupling between brain structure and function in several recent studies (Glomb et al., 2020; Medaglia et al., 2017; Preti & Van De Ville, 2019; Sihag et al., 2020). In Medaglia et al. (2017), the application of GSP to model fMRI activity onto a DWI connectivity graph allowed a novel analysis of brain functional-structural coupling. The authors investigated to what extent fMRI time series are constrained by the underlying structure. Results indicate that aligned signals were concentrated within default mode and frontoparietal systems, while the subcortical system included both aligned and liberal modes. Interestingly, these findings are interpreted in terms of brain dynamics flexibility and linked to cognitive performance, showing that GSP can discriminate behaviorally relevant signals. Preti and Van De Ville (2019) also decomposed resting-state fMRI time series using Laplacian eigenvectors of structural connectivity, revealing a gradient of large-scale organization of the cortex spanning from sensory-motor areas (with high functional-structural alignment) to higher cognitive function areas (whose activity is more decoupled from underlying structure). In seminal work by Glomb and colleagues (Glomb et al., 2020), the GSP framework allowed for a sparser representation of source EEG data than the conventional individual regions analysis. Few structural connectivity harmonics were shown to capture EEG task dynamics and, more importantly, revealed significant patterns of activation involving the entire cortex, which were disregarded in the classical region-by-region analysis. Together with a high-density EEG study in patients with consciousness disorders (Mortaheb et al., 2019), this is the only work applying GSP to the analysis of EEG signals, and it indicates that network harmonics also have functional significance, as they can be considered as an orthogonal basis of large-scale EEG dynamics.

Going beyond metrics and inference-based approaches, other studies have combined GSP and machine learning to derive features from graph Fourier transform (Brahim & Farrugia, 2020; Ménoret, Farrugia, Pasdeloup, & Gripon, 2017; Pilavci & Farrugia, 2019; Xu, Li, Zhu, & Wu, 2019). For instance, in Brahim and Farrugia (2020), a combination of GSP and machine learning was proposed for autism spectrum disorder classification. More specifically, the authors revealed that the analysis of fMRI data could be enriched by projecting resting-state fMRI (rs-fMRI) time series on a structural brain graph, as shown by substantial classification performance gains. In Pilavci and Farrugia (2019), the authors presented a feature extractor approach based on machine learning and spectral wavelets on brain graphs. In Zhang, Tetrel, Thirion, and Bellec (2020), a functional graph Laplacian embedding of deep neural networks (graph convolutional networks) is used to classify task fMRI time series, in a joint GSP-deep learning framework. Finally, other approaches have taken advantage of graphs to denoise brain signals, such as in Kim et al. (2019), where the authors simultaneously clean brain signals and learn the associated graph.

DISCUSSION AND PERSPECTIVES

In this paper we provided an overview of work that applied Laplacian spectral decomposition and GSP to analyze brain signals. We attempted to elucidate terminology and related approaches used in the “brain gradients” and GSP communities, systematically describing recent promising results. In this section we discuss the potential of an integrated gradients/GSP framework to reveal a spectral basis of brain activity grounded in brain connectivity topology. We will also bring up open questions and methodological challenges in this novel approach.

Connectivity Gradients as a Fourier Basis of Brain Activity

Recent work applying Laplacian-based spectral decomposition of brain networks has revealed functionally, developmentally, and anatomically relevant organizational hierarchies of the brain that have also been related to cognitive performance. Connectivity gradients obtained by such graph spectral decomposition not only have been shown to represent relevant axes of brain organization, but within the GSP framework can be seen as a Fourier basis to decompose brain activity. In this sense they provide a new spatial-frequency language to characterize patterns of neural activity and a novel perspective for probing brain dynamics.

A few studies have investigated the relationship between the principal (Ito et al., 2020; Lanzoni et al., 2020; X. Wang et al., 2020) or second (Mckeown et al., 2020) gradient of rs-fMRI connectivity during various tasks, suggesting that the cortical organization encoded in connectivity gradients supports specific cognitive or semantic functions. The GSP framework takes the analysis further by using the whole spectrum of brain gradients as a Fourier basis to decompose brain activity. For instance, the large-scale fMRI connectivity gradients identified in Margulies et al. (2016) could be exploited to disentangle brain activity measured during a complex cognitive task in unimodal (e.g., related to sensory processing) and transmodal (e.g., related to ongoing thoughts) patterns. This analysis would extend the work of Ito et al. (2020) relating local and distributed processes to the organization of the cortex in unimodal and transmodal areas. Moreover, by decomposing brain signals as a function of eigenvectors (or gradients) of the underlying connectivity, the GSP framework uses the information encoded in higher order connectivity gradients whose functional relevance has scarcely been explored in literature. For instance, in Preti and Van De Ville (2019), the spectrum of structural connectivity eigenvectors is split into low and high frequency components to define a binary decoupling index. Low frequencies correspond to signals coupled to structural connectivity, while high-frequency components are considered to be decoupled. The potential of the GSP framework could be further developed by considering the whole set of connectivity gradients, instead of a partition in low versus high frequencies. This is similar to the way classical Fourier analysis is used to decompose signals in the time domain. Continuing the analogy with the classical frequency domain, fundamental operations such as filtering and denoising can be generalized to brain signals on graphs by taking into account the full graph spectrum. For instance, artifactual components of brain activity (i.e., balistocardiogram artifact for EEG-fMRI simultaneous acquisition; Allen, Polizzi, Krakow, Fish, & Lemieux, 1998; Lioi et al., 2020) could be reduced by filtering out the graph frequency component or band that best represents the artifact (i.e., that maximally correlates with the electrocardiogram signal).

Using a continuous set of dimensions (graph Fourier modes) for the analysis of brain dynamics is an approach that complements (rather than excluding) the more classical hard boundaries parcellation. Some processes may be best characterized in terms of nonoverlapping fixed regions, others in terms of delocalized, overlapping eigenmodes. In this sense, this graph modal approach may be more appropriate than modular analysis in describing complex cognitive states depending on multiple overlapping phenomena. Neural patterns of ongoing activity could be seen as location in a multidimensional state-space constructed out of large-scale brain gradients (Margulies et al., 2016). The “biological” validity of this approach can be also found in the intrinsic organization of brain tissues. Brain structure (and function) are organized into overlapping hierarchical components (Betzel & Bassett, 2017). It is well known, for instance, that the visual and auditory cortices are organized into topographic maps that reflect how sensory information is processed (i.e., retinotopic or tonotopic mapping; Silver, Ress, & Heeger, 2005). As a result, in the same cortical area multiple and heterogeneous modes coexist. Moreover, while the cytoarchitecture of a region can be considered uniform, the same region can be heterogeneous in terms of function, gene profile, or axonal projections (i.e., connectivity topography; Haak et al., 2018). This intrinsic complexity may be better represented by a continuum of functions rather than a mosaic of brain areas (Haak & Beckmann, 2020).

The GSP approach also allows the efficient integration of different neuroimaging techniques (EEG, fMRI, DWI, MEG), thus exploiting complementary measurements of brain properties. Only a few studies have analyzed fMRI and EEG signals at rest using graph Fourier modes of the underlying structural graph (i.e., estimated from DWI). In the future, this framework could be applied to jointly decompose electrophysiological and functional time series using underlying structural topology, thus integrating different temporal and spatial scales in a multimodal analysis, with potential to shed new light on the complex interplay between brain function and structure. In addition, GSP could be extended to the analysis of brain signals during different tasks and (pathological) conditions, holding promise for developing more sensitive markers of disease.

Methodological Challenges and Future Directions

As described in Table 2 and discussed in Bajada et al. (2020), even if grounded in the same fundamental approach, various algorithms and similarity measures have been applied for graph inference and spectral analysis. These methodologies have different properties and it is not clear which one should be used and for which specific applications. For instance, different Laplacian matrices can be used for spectral decomposition. The normalized Laplacian has the useful property that its spectrum is limited to the [0, 2] interval. However its first eigenvector (associated with λ0 = 0) is not constant as for the combinatorial Laplacian, which is less intuitive for an interpretation of eigenvalues in terms of spatial frequency (Shuman et al., 2013). Similarly, some influential studies have used diffusion map embedding (Margulies et al., 2016; Mckeown et al., 2020; Paquola et al., 2019) with specific parameter choices (diffusion time and anisotropic diffusion parameter α) because of a series of advantageous properties for brain connectivity analysis. Future work should explore how different choices could affect results and which metrics are more adapted to the analysis of specific graph structures or brain signals.

Another point to clarify is the impact of the choice of parcellation on the connectivity gradient topography. In the pioneering work of Margulies et al. (2016), high-resolution rs-fMRI connectivity was computed between fMRI voxels, while in other work (Mckeown et al., 2020; Preti & Van De Ville, 2019; D. Wang, Yu, & Zou, 2020) different atlases were used to identify brain ROIs and then compute functional connectivity between them. It would be interesting to assess how the use of different parcellations affects the computation of connectivity gradients and whether it introduces bias in the characterization of connectivity topography, which is supposed to describe a continuous pattern of organization.

One important limitation of the studies we have reviewed is the lack of accounting for temporal dependencies in the graph model. Indeed, in most of the literature “constant” spatial dependencies between brain regions are considered (i.e., structural graphs built using white matter tracts, or functional graphs estimated using statistics of brain activity). Several theoretical breakthroughs have been made to address dynamic properties (i.e., time-varying graphs), and their application to neuroimaging data could prompt the understanding of cognitive process dynamics. Some promising frameworks to model time-varying aspects in graphs include graph slepians (Van De Ville, Demesmaeker, & Preti, 2017), sparseness of temporal variation (Yamada, Tanaka, & Ortega, 2019), or lapped Fourier transform (Lu & Ortega, 2019). Another avenue to model spatiotemporal dynamics could be to use deep learning models adapted to sequence modeling, combined with graph convolutional networks (J. Chen, Xu, Wu, & Zheng, 2018).

Finally, GSP is an active research field and there are a few recent theoretical proposals in GSP that have not yet been applied to neuroimaging data, but which could potentially bring interesting breakthroughs. Considering generalized signal processing operations on graphs such as graph filters (Segarra, Mateos, Marques, & Ribeiro, 2017), graph wavelets (Hammond, Vandergheynst, & Gribonval, 2011; Pilavci & Farrugia, 2019), multiscale graphs (Tremblay & Borgnat, 2014), graph slepians (W. Huang et al., 2017), graph sampling (Marques, Segarra, Leus, & Ribeiro, 2016; Puy, Tremblay, Gribonval, & Vandergheynst, 2016), or locating and decomposing signals on graphs (S. Chen, Yang, Moura, & Kovačević, 2016) could allow richer interpretations, and potentially a unified perspective on graph signals. GSP also includes novel methods for graph inference that are centered on signal representations on graphs (Dong, Thanou, Rabbat, & Frossard, 2019; Giannakis, Shen, & Karanikolas, 2018; Mateos, Segarra, Marques, & Ribeiro, 2019; Pasdeloup, Gripon, Mercier, Pastor, & Rabbat, 2018). GSP-based methods aim to infer graphs from measurements while enforcing specific properties of these signals on the learned graph (e.g., spectral distribution, smoothness), with potential to uncover specific structural (Hu et al., 2015) and functional (Shen, Giannakis, & Baingana, 2019) interactions between brain areas.

In conclusion, this work reviewed recent studies applying connectivity gradients and GSP for the analysis of brain signals, clarified terminology and methods, and related these two approaches grounded in the eigenvector decomposition of connectivity matrices. We point out that, given the increasing relevance connectivity gradients are taking in the understanding of brain macroscale organization, the application of GSP to neuroimaging is an exciting avenue towards a deeper understanding of brain organization. We also identify methodological challenges and suggest that future work should address multimodal and time-varying modeling and further explore the use of different metrics.

AUTHOR CONTRIBUTIONS

Giulia Lioi: Conceptualization; Resources; Writing – original draft. Vincent Gripon: Supervision; Writing – review & editing. Abdelbasset Brahim: Conceptualization; Writing – review & editing. François Rousseau: Supervision; Writing – review & editing. Nicolas Farrugia: Conceptualization; Supervision; Writing – original draft; Writing – review & editing.

FUNDING INFORMATION

Nicolas Farrugia, Région Bretagne (http://dx.doi.org/10.13039/501100011697), Award ID: SAD-2019.