I. INTRODUCTION

H UMAN brains are at the center of complex neurobiolog- ical systems in which neurons, circuits, and subsystems interact to orchestrate behavior and cognition. Understanding the structures, functions, and mechanisms of human brains has been an intriguing pursuit for researchers with various goals, including neural system simulation, mental disorder therapy, as well as general artificial intelligence. Recent studies in neuroscience and brain imaging have reached the consensus that interactions between brain regions are key driving factors for neural development and disorder analysis [1,2]. Inspired by graph theory, brain networks composed of nodes and edges are developed to describe the interactions among brain regions.

The human brain can be scanned through various medical imaging techniques, including Magnetic-Resonance Imaging (MRI), Electrogastrography (EGG), Positron Emission Tomography (PET), and so on. Among all these acquisitions, MRI data are the most widely used for brain analysis research. There are also different modalities of MRI data such as functional MRI (fMRI) and Diffusion Tensor Imaging (DTI), from which functional and structural brain networks can be constructed respectively. Specifically, the connectivity in functional brain networks describes correlations between time-series signals of brain regions, while the connectivity in structural brain networks models the physical connectivity between gray matter regions [3]. Both functional and structural connections are widely acknowledged as valuable resources of information for brain investigation [4,5].

Previous work on brain network analysis has studied shallow models based on graph theory [5,6] and tensor factorization [7,8] extensively, which focuses on proposing neurobiologically insightful graph measures and approaches from the node, motif, and graph level to detect network communities or modules and identify central network elements. Methodological developments in graph research enable us to quantify more topological characteristics of complex systems, many of which have already been assessed in brain networks, such as modularity, hierarchy, centrality, and the distribution of network hubs. However, shallow modeling techniques can be inadequate for the sophisticated connectome structures of brain networks [9]. On the other hand, deep learning models have become extraordinarily popular in machine learning, achieving impressive performance on images [10,11], videos [12], and speech processing tasks [13]. These regular data are represented in 1D/2D/3D Euclidean spaces and can be suitably handled by traditional Recurrent (RNNs) or Convolutional Neural Networks (CNNs). In contrast, the irregular structural and functional brain connectivity networks constructed from neuroimaging data are more complex due to their non-Euclidean characteristics. In recent years, Graph Neural Networks (GNNs) have attracted broad interest due to their established power for analyzing graph-structured data [14][15][16]. Several pioneering deep models have been devised to predict brain diseases by learning graph structures of brain networks. For instance, Li et al. [1] propose BrainGNN to analyze fMRI data, where ROIaware graph convolutional layers and ROI-selection pooling layers are designed for neurological biomarker prediction. Kawahara et al. [17] design a CNN framework BrainNetCNN composed of edge-to-edge, edge-to-node, and node-to-graph convolutional filters that leverage the topological locality of structural brain networks. However, they mainly experiment w it h a t t e n t io n g n

Graph pooling

Output: gn

Multiple layers

Fig. 1: An overview of our BrainGB framework for brain network analysis with graph neural networks.

with their proposed models on specific private datasets. Due to the ethical issue of human-related research, the datasets used are usually not publicly available and the details of imaging preprocessing are not disclosed, rendering the experiments irreproducible for other researchers.

To address the aforementioned limitations, there is an urgent need for a public benchmark platform to evaluate deep graph models for brain network analysis. However, it is non-trivial to integrate different components within a unified benchmarking platform. Current brain network analyses are typically composed of two steps. The first step is to construct brain networks from neuroimaging data. Then, in the second stage, the resulting brain connectivity between all node pairs is used to classify individuals or predict clinical outcomes. The difficulties in the initial stage are mostly due to restricted data accessibility and sophisticated brain imaging preprocessing and network construction pipelines that differ across cohorts and modalities. The difficulty of the second stage is to establish a standard evaluation pipeline based on fair experimental settings, metrics, and modular-designed baselines that can be easily validated and extended for future research.

In this work, we propose Brain Graph Neural Network Benchmark (BrainGB)-a novel attempt to benchmark brain network analysis with GNNs to the best of our knowledge. The overview of BrainGB is demonstrated in Fig. 1 and the main contributions are four-fold: We conduct a variety of empirical studies and suggest a set of general recipes for effective GNN designs on brain networks, which could be a starting point for further studies. To foster future research, we release the source code of BrainGB at https://github.com/HennyJie/Brai nGB and provide an out-of-box package that can be installed directly, with detailed tutorials available on our hosted website at https://braingb.us. Preprocessing instructions and models are provided for standardized model evaluations. We enable the community to collaboratively contribute by submitting their own custom models, and we will maintain a leaderboard to ensure such efforts will be recorded.

II. PRELIMINARIES

A. Brain Network Analysis

Brain networks are complex graphs with anatomic Regions of Interest (ROIs) represented as nodes and connectivities between the ROIs as links [18]. In recent years, the analysis of brain networks has become increasingly important in neuroimaging studies to understand human brain organization across different groups of individuals [19][20][21][22][23]. Abundant findings in neuroscience research suggest that neural circuits are highly related to brain functions, with aberrations in these neural circuits being identified in diseased individuals [24][25][26].

Formally, in the task of brain network analysis, the input is a brain network dataset D = {G n , y n } N n=1 consisting of N subjects, where G n = {V n , E n } represents the brain network of subject n and y n is the subject's label of the prediction, such as neural diseases. In D, the brain network G n of every subject n involves the same set of M nodes defined by the ROIs on a specific brain parcellation, i.e., ∀n, V n = V = {v i } M i=1 . The difference across subjects lies in the edge connections E n among M brain regions, which are often represented by a weighted adjacency matrix W n ∈ R M ×M describing the connection strengths between ROIs. The edge weights in W are real-valued and the edges are potentially dense and noisy. The model outputs a prediction ŷn for each subject n, which can be further analyzed in terms of features and biomarkers.

Given brain networks constructed from different modalities such as Diffusion Tensor Imaging (DTI) and functional Magnetic Resonance Imaging (fMRI) [5,27,28], effective analysis of the neural connectivities of different label groups (e.g., disease, gender) plays a pivotal role in understanding the biological structures and functions of the complex neural system, which can be helpful in the early diagnosis of neurological disorders and facilitate neuroscience research [29][30][31][32][33][34][35]. Previous models on brain networks are mostly shallow, such as graph kernels [36] and tensor factorization [37,38], which are unable to model the complex graph structures of the brain networks [9].

B. Graph Neural Networks

Graph Neural Networks (GNNs) have revolutionized the field of graph modeling and analysis for real-world networked data such as social networks [14], knowledge graphs [39], protein or gene interaction networks [15], and recommendation systems [40]. The advantage of GNNs is that they can combine node features and graph structures in an end-to-end fashion as needed for specific prediction tasks. A generic framework of GNN could be represented in two phases. In the first phase, it computes the representation h i of each node v i ∈ V n by recursively aggregating messages from v i 's multi-hop neighborhood, where h 0 i is initialized with node features. After getting the last-layer node representation h (L) , an extra pooling strategy is adopted to obtain the graph representation. Thereafter, a Multi-Layer Perceptron (MLP) can be applied to make predictions on the downstream tasks.

It is worth noting that brain networks are different from other real-world graphs such as social networks or knowledge graphs, due to (1) the lack of useful initial node (ROI) features on brain networks represented by featureless graphs, (2) the real-valued connection weights that can be both positive or negative, and (3) the ROI identities and their orders are fixed across individual graph samples within the same dataset. The design of GNN models should be customized to fit the unique nature of brain network data. Recently, there have been emerging efforts on GNN-based brain network analysis [1,17,[41][42][43][44][45][46][47]. However, these models are only tested on specific local datasets, mainly due to the convention in neuroscience that researchers are more used to developing methods that are applicable to their specific datasets and the regulatory restrictions that most brain imaging datasets are usually restrictively public, meaning that qualified researchers need to request access to the raw imaging data and preprocess them to obtain brain network data, but they are not allowed to release the preprocessed data afterwards. These challenges largely prohibit the methodology development in computational neuroscience research.

III. BRAIN NETWORK DATASET CONSTRUCTION

A. Background: Diverse Modalities of Brain Imaging

Models of the human brain as a complex network have attracted increasing attention due to their potential for helping understand human cognition and neurological disorders. In practice, human brain data can be acquired through various scanning techniques [48], such as Magnetic-Resonance Imaging (MRI), Electroencephalography (EEG) and Magnetoencephalography (MEG), Positron Emission Tomography (PET), Single-Photon Emission Computed Tomography (SPECT), and X-ray Computed Tomography (CT). Among them, MRI is one of the most widely used techniques in brain research and clinical practice, due to its large range of available tissue contrast, detailed anatomical visualization, and high sensitivity to abnormalities [49].

1) MRI Data: In this paper, we focus on MRI-derived brain networks. Specifically, for different modalities of MRI data, we can reconstruct different types of brain networks. Functional MRI (fMRI) is one of the most popular modalities for investigating brain function and organization [31,32,50] by detecting changes in blood oxygenation and blood flow that occur in response to neural activity. Diffusion-weighted MRI (dMRI), on the other hand, can enable inference about the underlying connection structure in the brain's white matter by recording the diffusion trajectory of molecules (usually water). fMRI focuses on functional activity, while dMRI presents brain structural information from different perspectives. Specifically, two types of brain networks, functional and structural, can be constructed from the aforementioned modalities by following different connectivity generation paradigms [51].

2) Challenges in MRI Preprocessings: The raw MRI data collected from scanners is not directly usable for brain network construction or imaging analysis. A complicated preprocessing pipeline is necessary to remove unwanted artifacts, transform the data into a standard format, and perform structure discovery. Although there are several widely-used neuroimaging data preprocessing tools, such as SPM 1 , AFNI 2 and FSL 3 , each of them still needs considerable training and learning efforts. Moreover, the functionality of these software varies, and for dMRI, no one software contains all the necessary preprocessing capabilities. In addition, many neuroimaging datasets cannot be made public due to privacy or ethical concerns. Due to the variety of preprocessing approaches and issues with making data publically available, there are difficulties in reproducibility in neuroimaging studies. Additionally, the preprocessing steps are distinctive across modalities. All these challenges make it difficult for deep learning researchers with little knowledge in medical imaging processing to get into the field.

B. Brain Network Construction from Raw Data

In this section, we provide a general overview of the standard preprocessing pipelines for the construction of brain networks of different modalities. Due to the regulation restrictions for direct sharing of the brain network data, we provide two complete pipelines, one for functional brain networks (ABCD 4 specifically) and one for structural brain networks (PPMI 5 specifically), with step-by-step commands and parameter settings on our hosted website for public access 6 . The idea

Construct Network

Calculate pairwise correlations between ROIs as edges

SPM 12 AFNI FSL Free Surfer CONN fMRI Prep ANTs Nilearn ✓

Eddy-current and Head Motion Correction Align all raw images to the b0 image to correct for head motion and eddy current distortions is, by following our instructions there, researchers with basic programming skills can easily construct similar brain networks after requesting the data access by themselves, so everyone in the community can compare their methods fairly on the same datasets.

1) Functional Brain Network Construction: The left side of Fig. 2 shows a standard preprocessing procedure for functional brain imaging, with the corresponding commonlyused toolboxes (i.e., SPM12 1 , AFNI 2 , FSL 3 , FreeSurferfoot_0 , CONN 8 , fMRI Prepfoot_2 , ANTsfoot_3 , Nilearnfoot_4 ) shown on the right side. Note that each step in the preprocessing and network construction pipeline needs quality control by the experts, and the specific order of preprocessing steps may change slightly based on the acquisition conditions of the dataset. Some representative functional neuroimaging datasets in literature to facilitate scientific research include ADHD 200 [52], ADNI (fMRI part) [53], HCP 900 [54], ABIDE [55], etc.

To measure functional connectivity, some preprocessing of the fMRI time series is often performed including detrending, demeaning, and whitening fMRI BOLD time series at each voxel [56]. To construct the brain networks, a brain atlas or a set of Regions of Interest (ROI) are selected to define the nodes. Then, the representative fMRI BOLD series from each node are obtained by either averaging or performing Singular Value Decomposition (SVD) on the time series from all the voxels within the node. Various measures have been proposed for assessing brain connectivity between pairs of nodes. One of the simplest and most frequently used methods in the neuroimaging community is via pairwise correlations between BOLD time courses from two ROIs. Other methods including partial correlations [56], mutual information, coherence, Granger causality, etc [57]. After selecting the Functional Connectivity (FC) measure, one can evaluate the strength of connectivity between each pair of nodes in the brain network. Often, some transformation, such as the Fisher's transformation, is performed to transform the original FC measures to improve their distribution properties. The transformed FC measures can then be utilized for the subsequent analysis of functional brain networks.

To facilitate public testing, we take Adolescent Brain Cognitive Development Study (ABCD) as an example and provide a step-by-step instruction for functional brain network construction on our hosted BrainGB website 12 . The ABCD-HCP BIDS 13 pipeline is used to preprocess the data. In brief, anatomical preprocessing included normalization, coregistration, segmentation, and brain extraction. Functional data preprocessing included slice-time correction, motion correction, distortion correction, co-registration, normalization, and spatial smoothing. Brain parcellation schemes were then applied to the functional data to obtain time courses for each ROI, and Pearson correlation was used to construct brain networks representing the connectivity between ROIs.

2) Structural Brain Network Construction: Structural brain networks provide a systematic perspective for studying the anatomical and physiological organization of human brains and help to understand how brain structure influences function. Some representative neuroimaging studies include diffusion MRI data are PPMI [58], ADNI [53], HCP [54], AIBL [59], OASIS [59], etc. The commonly-used toolboxes for dMRI include FSL 3 , AFNI 2 , FreeSurfer 7 , TrackVisfoot_7 , 3D Slicerfoot_8 , Tortoise 16 , MRtrix3 17 , DSI Studio 18 .

The left side of Fig. 3 summarizes the pipeline for reconstructing the structural brain network. Preprocessing steps for the dMRI data include removal of eddy current-induced distortions, brain extraction, and co-registration between diffusion and structural images. Next, some modeling strategies are applied to reconstruct the local diffusion patterns. Commonly adopted models include the DTI modeling, which fits a tensor model or muti-tensor model [60] to capture the local diffusion patterns, and the Ball and Sticks model [61]. After reconstructing the local diffusion patterns, a tractography algorithm is performed to computationally reconstruct fiber tract connections between brain regions. Commonly-used algorithms include the deterministic tractography [62] and the probabilistic tractography [63]. The deterministic tractography connects neighboring voxels from seed regions based on the major direction of the DTI tensor. The probabilistic tractography involves first estimating fiber orientation and its uncertainty at each voxel and building a diffusion path probability map based on the estimated orientation and uncertainty. While deterministic tractography is a more computationally efficient approach to reconstruct major fiber bundles in the brain, probabilistic tractography has become more popular because it is more robust to noise and allows tractography to progress beyond uncertain regions by taking into account uncertainty in fiber orientations at each voxel [64]. To construct the structural network, the structure connectivity for each node pair is calculated based on the empirical probability of fiber tracts connecting the two regions. Note that each step in the preprocessing and network construction needs quality control from experts.

Similarly to functional brain network construction, we also take Parkinson's Progression Markers Initiative (PPMI) as an example and provide an instruction pipeline for structural brain network construction on our hosted BrainGB website 19 . Specifically, the Diffusion Toolkit from TrackVis is used to reconstruct local diffusion patterns and tractography. The brain region parcellation is completed with both FSL and Freesurfer. Then local diffusion pattern reconstruction and the network computation are further performed by calculating the number of fibers within each ROI after removing the false positive ones.

C. Discussions

In addition to the mainstream methods of constructing connections in brain networks discussed above, there are also other ways to construct different types of edges. For example, directional connectivity that characterizes effective interactions for fMRI [65]; hybrid functional brain networks where different orders of relationships can be sensitive to different levels of signal changes [66]; and dynamic functional brain networks which include derivatives of windowed functional network connectivity in the identification of reoccurring states of connectivity [65,67]. Apart from fMRI and DTI, the most commonly used modalities to construct functional and structural brain networks, other neuroimaging modalities have also been explored in literature, such as metabolic brain network constructed from PET imaging [68], functional brain network constructed from EEG signals [69], etc. Recent studies have shown that the combination of both functional and structural neuroimaging modalities can be more effective than using only a single one, which can exploit complementary information across different modalities [4,70].

IV. GNN BASELINES FOR BRAIN NETWORK ANALYSIS

The process of applying GNNs to brain networks starts from initialization of the ROI features, followed by the forward pass which includes two phases, message passing, and pooling. The learned graph-level representation then can be utilized for brain disease analysis. In the machine learning domain, the rapid evolution of GNNs has led to a growing number of new architectures. Specifically for GNNs on brain network analysis, we decompose the design space of interest for basic message passing GNNs into four modules: node feature construction, message passing, attention enhanced message passing, and pooling strategies. An illustration of these modules is shown in the middle of Fig. 1.

A. Node Feature Construction

In neuroscience analysis, researchers mostly focus on brain connectivity represented by a featureless graph. To apply GNNs on non-attributed brain networks, researchers in the graph machine learning domain have studied several practical methods to initialize node features [71,72]. In this paper, we focus on the following node features that can be categorized as positional or structural:

• Identity: A unique one-hot feature vector is initialized for each node [73,74]. By giving each ROI in the brain network a unique high-dimensional vector, this identity node feature allows the GNN model to learn the relative positions of the nodes by memorizing their k-hop neighbors. They are essentially the same as random initialization considering the parameters in the first linear layer of the GNN are randomly initialized. • Eigen: Eigen decomposition is performed on the weighted matrix describing the connection strengths between ROIs and then the top k eigenvectors are used to generate a k-dimensional feature vector for each node [75][76][77].

The optimal value of k is decided by grid search. This feature is essentially dimension reduction and targets at grouping brain regions with respect to their positions, with global graph information condensed into a lowdimensional representation. • Degree: The degree value of each node is obtained as a one-dimensional vector as the node feature. This feature captures structural information of brain regions, meaning that neighborhood structural similarity of two regions will be partially recorded in the initialized node features. • Degree profile: This method takes advantages of existing local statistical measures on degree profiles [78], where each feature x i of node v i on graph G n is computed as

where

describes the degree values of node v i 's one-hop neighborhood and ∥ denotes concatenation.

• Connection profile: The corresponding row for each node in the edge weight matrix is utilized as the initial node feature, which contains connections with respect to all other nodes in the brain network. This feature aligns with the common practice of using pairwise connections to perform brain parcellation. Also, it reflects the whole picture of structural or functional information in the brain network.

B. Message Passing Mechanisms

The power of most GNNs to learn structures lies in their message passing schemes, where the node representation is updated iteratively by aggregating neighbor features through local connections. In each layer l, the node representation h l i is updated through two steps, namely message passing and update respectively. In the message passing step (Eq. 2), each node v i receives messages from all its neighbors, and then all the messages are aggregated with a sum function:

where N i denotes the neighbors of node v i in graph G, w ij represents the edge weights between node v i and v j , M l is the message function. In the update step (Eq. 3), the embedding of each node is updated based on the aggregated messages from Eq. 2 and optionally the previous embedding of node v i , where the update function can be arbitrary differentiable functions (e.g., concat the aggregated message with the previous node embedding and then pass them into another learnable linear layer).

where U l stands for the update function and the number of running steps L is defined by the number of GNN layers. The message passing mechanism can leverage both permutation equivariance and inductive bias towards learning local structures and achieve good generalization on new graphs. For brain networks, whether incorporating connections into the message function is beneficial for graph-level prediction tasks remains to be investigated. In this paper, we discuss the influence of different message function M l designs including:

• Edge weighted: The message m ij passed from node v j to node v i is calculated as the representation of node v j weighted by the corresponding edge weight w ij , that is

This is the standard message passing implementation in Graph Convolutional Network (GCN) [14] when w ij = 1 /Ni. With this message vector design, the update of each brain region representation is influenced by its neighbor regions weighted by the connection strength between them. • Bin concat: In this scheme, we map the edge w ij into one of the equally split T buckets based on its weight value. Each bucket corresponds to a learnable representation b t , t = {1 . . . T }. The total bucket number encompassing the entire value range of edge weights is determined by grid search and the representation dimension of each bin is set to the same as node features. Specifically, given the number of buckets is T , we first rank all the edge weights and then divide them into the equally divided T buckets from the lowest to the highest. All those edges falling into the same bucket will be mapped to the same learnable edge representation b t , so region connections with similar strength are binned together. In our experiment, we simply select from [5,10,15,20] as the possible number of buckets for grid search, which is a common practice in machine learning for hyperparameter tuning. The message m j passed from node v j to node v i is calculated as the concatenation of the representation of node v j and its corresponding bucket representation b t followed by a MLP,

The usage of bins helps to clusters region connections with similar strengths. By concatenating with the unique neighbor node representation, this message captures both common and peculiar characteristics of each neighbor. representation of node v j and the scaled edge weight d • w ij , followed by a MLP,

where d is a constant equals to the dimension number of node features. The motivation behind edge weight scaling is to increase the influence of edge features to the same scale as node features. Compared with bin concat design where edges with weight values falling into the same bin interval share the same initial edge representation, directly concatenating the scaled edge weights as the edge representations can retain the original edge information, therefore reserving more uniqueness on the pairwise connection when performing the aggregation from neighboring brain regions. • Node edge concat: To investigate the influence of preserving the brain region representation from the last time step while iterative updating the new representation, we design a message m j as the concatenation of both embeddings of node v i , v i and the edge weight w ij between them, followed by a MLP, that is

In this paradigm, every message passed from the local neighbors of each central node is reinforced with its representation from the last time step. This design may alleviate the over-smoothing problem of GNNs, where the feature distance between all nodes becomes too close and not distinguishable after layers of convolutions. • Node concat: Since the effect of involving connection weights into message passing is still unknown, we also include another message m ij design similar to node edge concat but without the concatenation of edge weights, where

C. Attention-Enhanced Message Passing

Attention is arguably one of the most important mechanisms in modern deep learning [79,80]. It is inspired by human cognitive systems that tend to selectively concentrate on the important parts as needed when processing large amounts of information. Various fields in deep learning communities such as natural language processing [81] and computer vision [82] have widely benefited from attention mechanisms in terms of model efficiency and accuracy. The attention mechanism can also be used to enhance the message passing scheme of GNNs, while also providing interpretations over the edge importance.

Specifically in brain network analysis, by utilizing the attention-enhanced version of message passing, the model updates the brain region representation in a data-driven way, where adjustable attention weights from each local neighbor perform as an additional influence factor besides the neural signals represented by edge weights. It is worth noting that the traditional designs of graph attention mechanisms on general graphs usually do not take the edge attributes (i.e., connection weights in the brain network scenario) into consideration.

However, for brain networks, the correlation between two regions contains meaningful biomedical information and might be helpful for graph-level tasks. In this paper, we design several attention-enhanced message passing mechanisms including:

• Attention weighted: This is the original GAT [16] on general graphs without involving edge attributes. The message from node v j to v i is calculated as the node representation v j weighted by the corresponding attention weights α ij ,

The α ij is calculated from a single-layer feed-forward neural network parameterized by a weight vector a, followed by the LeakyReLU nonlinearity σ,

, (10) where Θ represents a learnable linear transformation matrix.

• Edge weighted w/ attn: This is the attention-enhanced version of edge weighted message passing in Eq. 4. The message from v j to v i is obtained as the multiplication of node v j 's representation h j , the edge weight w ij and the attention score α ij in Eq. 10,

• Attention edge sum: This is another version of attentionenhanced edge weighted (Eq. 4) message passing. The edge weight w ij and the attention score α ij are first summed, then used as the impact factor on the node embedding h j ,

• Node edge concat w/ attn: This is the attention-enhanced version of node edge concat (Eq. 7) message passing, where the attention score α ij (Eq. 10) between node v i and v j is multiplied on the node representation h j before concatenation, followed by a MLP,

• Node concat w/ attn: This design corresponds to the attention-enhanced version of node concat (Eq. 8) message passing, where the attention score α ij (Eq. 10) between node v i and node v j is multiplied on the node representation h j before concatenation, followed by a MLP,

D. Pooling Strategies

In the second phase of GNNs, a feature vector for the whole graph g n is computed using the pooling strategy R, where

The pooling function R operates on the set of node vectors and is invariant to permutations of the node vectors. In this paper, we cover three basic global pooling operators [83,84]:

• Mean pooling: The graph-level representation is obtained by averaging node features. For each single graph G n , the graph-level representation is computed as

• Sum pooling: The graph-level representation is obtained by summing up all node features. For each single graph G n , the graph-level representation is computed as

• Concat pooling: The graph-level representation is obtained by concatenating node features of all nodes contained in the graph. For each single graph G n , the graph-level representation is computed as

Note that there are also other complex pooling strategies such as hierarchical pooling [85], learnable pooling [86] and clustering readout [87], which are usually viewed as independent GNN architecture designs that are not defined based on combinative modules. Here we include the representative method of DiffPool [85] to provide a view of the comparison between basic and more complex pooling methods.

V. EXPERIMENTAL ANALYSIS AND INSIGHTS

In this section, we show experimental results on brain networks generated from real-world neuroimaging studies with different GNN modular designs. Varying each design dimension under each module results in a total of 375 different architectures. Note that here we do not target at covering all the combinations, but with an ambition to quickly find a good design choice for a specific dataset or a downstream task. Furthermore, we emphasize that the design space can be expanded as new design dimensions emerge in state-of-the-art models.

A. Experimental Settings 1) Datasets:

To establish a benchmark for generic brain network analysis models, we include four datasets processed and constructed from different neuroimaging modalities, specifically fMRI (HIV [38], PNC 20 , ABCD 21 ) and dMRI (PPMI 22 ), based on different brain atlas. For the HIV and PPMI datasets, the task is to classify patients from healthy control (Patient, Normal Control); while for the PNC and ABCD datasets, the task is gender prediction (Male, Female). We intentionally cover such a diverse set of datasets from different modalities (and preprocessing procedures/parcellations/tasks), because our purpose is to establish a benchmark for generic brain network analysis models. Thus observations on a diverse set of datasets can be more instructive for methodology focused studies. All the datasets we used have been visually checked by imaging experts in our team for quality control. Among these four datasets, PNC, PPMI, and ABCD are restrictively publicly available ones that can be requested and downloaded from their official website. The dataset information is summarized in TABLE I. Since the datasets can be acquired from multiple sites, multisite issues need to be addressed when performing the analysis on the constructed networks. Over the past few years, ComBat techniques [88,89] from the microarray literature have started to be used more frequently to deal with multi-site batch effects. Since our benchmark focuses more on a comprehensive overview of brain network construction and effective GNN designs for brain networks, advanced methods for handling multi-site issues are out of the scope of this work. Interested readers can refer to [90][91][92][93][94] for more advanced multisite data handling methods. [98] and the brain network is reconstructed using the deterministic 2nd-order Runge-Kutta (RK2) whole-brain tractography algorithm [64].

This study recruits children aged 9-10 years across 21 sites in the U.S. Each child is followed into early adulthood, with repeated imaging scans, as well as extensive psychological and cognitive testing [99]. After selection, 7,901 children are included in the analysis, and 3,961 (50.1%) among them are female. We use rs-fMRI scans for the baseline visit processed with the standard and open-source ABCD-HCP BIDS fMRI Pipeline 23 . After processing, each sample contains a connectivity matrix whose size is 360 × 360 and BOLD time-series for each node. The region definition is based on the HCP 360 ROI atlas [100]. Structural connectivity and functional connectivity are different in their strength and sparsity, thus need to be handled differently. For structural connectivity, we normalize the edge weights by dividing each value by the maximum value in a sample. The processed edge weights are thus ranged from 0 to 1. For functional connectivity, we follow common practice to remove the negative values for GNNs that cannot handle negative values (like GCN), and keep them for GNNs that can handle negative values (like GAT).

2) Baselines: For comprehensiveness, we compare our modular design with competitors of both shallow and deep models. The shallow methods we consider include M2E [7], MPCA [101], and MK-SVM [102], where the output graphlevel embeddings are evaluated using logistic regression classifiers. Specifically, M2E is a partially-symmetric tensor 23 https://github.com/DCAN-Labs/abcd-hcp-pipeline factorization based method for brain network analysis, and it has been empirically compared with spectral embedding clustering methods such as SEC [103] or spectral learning frameworks such as AMGL [104]; MPCA is proposed for the feature extraction and analysis of tensor objects such as neuroimaging; multiple kernel SVM (MK-SVM) is essentially an extension of the conventional SVM algorithm and has been applied for the analysis of functional and structural connectivity in Alzheimer's disease. We also include two stateof-the-art deep models specifically designed for brain networks: BrainGNN [1] and BrainNetCNN [17]. The message passing in BrainGNN is Edge weighted and it further leverages additional regional information (such as coordinates or ROI ordering based one-hot embeddings) to assign a separate GCN kernel for each ROI where ROIs in the same community are embedded by the similar kernel and those in different communities are embedded in different ways, but this will introduce a lot of additional model parameters and make the model hard to train. On the other hand, BrainNetCNN models the adjacency matrix of a brain network as a 2D image and does not follow the message passing mechanism as we discussed in Section IV-B. Note that the purpose of our paper, and of most benchmark papers, is not to establish superior performance of a certain method, but rather to provide an effective and fair ground for comparing different methods.

3) Implementation Details: The proposed model is implemented using PyTorch 1.10.2 [105] and PyTorch Geometric 2.0.3 [106]. A Quadro RTX 8000 GPU with 48GB of memory is used for model training. The optimizer we used is Adam. We train all of our models through 20 epochs, and the learning rate is 1e-3. We use a weight decay of 1e-4 as a means of regularization. The loss function is cross entropy. Hyperparameters are selected automatically with an open-source AutoML toolkit NNI 24 . Please refer to our repository for comprehensive parameter configurations. When tuning the hyperparameters, we first split the dataset into a train set and a test set with the ratio of 8:2. The k-fold validation is performed on the train set, where we further divide the train set into 10 parts and take one in each run to use as the validation set. The selection of the best hyperparameter is based on the average performance of the model on the validation sets. The reported metrics in Table II, on the other hand, is the average performance on the test set, with each run trained on different train sets. The competing methods are also tuned in the same way. For BrainGNN, we used the author's open-source code 25 . For BrainNetCNN, we implemented it by ourselves with PyTorch, which is publicly available in our BrainGB package 26 . For the hyper-parameter tuning, we selected several important hyper-parameters and performed the grid search on them based on the provided best setting as claimed in their paper. To be specific, for BrainGNN, we searched for different learning rates in {0.01, 0.005, 0.001} with different feature dimensions in {100, 200} and the number of GNN layers in {2, 3}. For BrainNetCNN, we searched for different dropout rates in {0.3, 0.5, 0.7} with learning rates in {0.001, 0.0005, 0.0001} and the number of layers in MLP in {1, 2, 3}. The reported results of these two baselines in Table II are from the best performing groups, where for BrainGNN, the learning rate is 0.01, the feature dimension is 200 and the number of GNN layers is 2, and for BrainNetCNN, the dropout rate is 0.3, the learning rate is 0.0001 and the number of layers in MLP is 3. The metrics used to evaluate performance are Accuracy, F1 score, and Area Under the ROC Curve (AUC), which are widely used for disease identification. To indicate the robustness of each model, all the reported results are the average performance of ten-fold cross-validation conducted on different train/test splits.

B. Performance Report

1) Node Feature: On comparing node features, we set the other modules as the well-performed settings in individual tests. Specifically, we use node edge concat in Eq. 7 as the message passing scheme, and concat pooling in Eq. 18 as the pooling strategy. Our experimental results demonstrate that the connection profile which uses the corresponding row in the adjacency matrix as the node features achieves the best performance across all datasets, with up to 33.99% improvements over the second-best, degree, on ABCD. We believe this is because the connection profile captures the whole picture of structural information in the brain network, and preserves rich information on pairwise connections that can be used to perform brain parcellation. In general, the structure node features (e.g., degree, connection profile) perform better than the positional ones (e.g., identity, eigen), indicating that the overall structural information of graph and the structural role of each node are important in the task of brain network analysis. This conclusion is consistent with previous findings in the literature that structural artificial node features work well for graph-level tasks on general graphs [71].

2) Message Passing: To study the effectiveness of different message passing schemes, we initialize the node features with connection profile and apply the concat pooling to produce graph-level representations, which both perform best when examined separately in each module. Our results reveal that node concat (Eq. 8) message passing has the highest AUC performance across four datasets, followed by node edge concat (Eq. 7), which achieves a similar AUC performance with sometimes slightly better accuracy and F1 scores (ABCD). The performance superiority of the last two methods may arise from their advantage of reinforcing self-representation of the central node during each step of message passing. This helps to retain the original information from the last step and avoid overfitting towards a biased direction in the optimization process. Surprisingly, the edge involved node edge concat performs slightly worse than the pure node concat, though the gap gets closer on larger datasets. This indicates that encoding edge weights as a single value may not be useful when the global structure has already been used as the initial node features.

3) Attention Enhanced Message Passing: When evaluating the effectiveness of different attention-enhanced message passing schemes, we set the node features as connection profile and apply the concat pooling strategy, just as for the evaluation of message passing without attention mechanisms. It is shown that the node concat w/ attn (Eq. 14) and node edge concat w/ attn (Eq. 13) yield very close results across four datasets and they alternately perform the best. Furthermore, the attentionenhanced version achieves better outcomes most of the time (up to 5.23% relative improvements) vs. the corresponding message passing architecture without an attention mechanism. This demonstrates the effectiveness of utilizing learnable attention weights in the GNN aggregation and update process in addition to the fixed edge weights. Also, the node edge concat w/ attn surpasses node concat w/ attn on the larger dataset (e.g., ABCD), which may imply potential advantages of involving edge weights into message design when there are enough training samples.

4) Pooling Strategies: For studying pooling strategies, we employ the node edge concat (Eq. 7) as the message passing scheme and connection profile as the initial node features. Our findings reveal that the concat pooling strategy (Eq. 18) consistently outperforms the other two methods across all four datasets. This is likely because when concat is used, the final node representations of all the brain regions are kept in the graph-level representation for classifiers. The other two paradigms, on the other hand, obtain a graph-level embedding with the same dimension of node features. Thus they lose some information that could be helpful for graph-level prediction tasks. Though concat does not ensure permutation invariance, it is actually not needed for brain network analysis since the node order given a parcellation is fixed. The compared hierarchical pooling method DiffPool demonstrates some advantages on the small HIV dataset but fails to surpass the simple concat pooling on three other larger datasets.

5) Other Baselines: In general, we expect deep models like GNNs to perform better on larger datasets. For example, the performance of GNN models on the ABCD dataset clearly surpasses all shallow models by about 2 percent. However, this trend should not prohibit one from experimenting with GNN models on smaller datasets. GNNs do perform well on some small datasets, such as the HIV dataset. Despite running on a small dataset, GNN models in BrainGB have an over 5 percent advantage over all shallow models. As for the deep baselines, BrainGNN can be out-of-memory (OOM) on large datasets. The best combination based on our modular design outperforms BrainGNN on all four datasets (HIV, PNC, PPMI and ABCD) and achieves comparable results with BrainNetCNN in most cases especially on smaller datasets (HIV, PPMI). These findings prove the need to carefully experiment with our modular designs of GNNs before further developing more complicated architectures, which might just overfit certain datasets.

6) Insights on Density Levels: Functional connectivity and structural connectivity have distinctive differences in sparsity levels. Functional networks like ABCD are fully connected. Structural networks like PPMI contain approximately 22.64% edges on average. Through our experiments, we found sparsity levels do have an impact on the choices of hyperparameters. For example, GNNs on the sparser structural networks of PPMI reach the maximum performance with a hidden dimension of 64, whereas on the functional network of ABCD, they have an optimal hidden dimension of 256, which indicates that GNN models should more complicated with more learnable parameters when the input networks are denser. This observation can be instructive for designing GNN architectures on brain networks constructed from different modalities.

VI. OPEN SOURCE BENCHMARK PLATFORM

To foster future research, we provide an out-of-box package that can be directly installed through pip, with installation and tutorials on our hosted BrainGB website https://brai ngb.us. The BrainGB package is also open-sourced at ht tps://github.com/HennyJie/BrainGB. We provide examples of GNN-based brain network analysis, trained models, and instructions on imaging preprocessing and functional and structural brain networks construction from raw fMRI and dMRI respectively. It is noted that due to the modular designs, BrainGB can also be extended to other tasks, by adding taskspecific functions in each module.

VII. DISCUSSION AND EXTENSIONS

In this paper, we first present BrainGB, a unified, modular, scalable, and reproducible framework for brain network analysis with GNNs. While the dataset generation, baselines, and evaluations we provide in BrainGB are thorough, we consider several limitations in the current paradigm:

• The aggregation mechanism in GNN is known to be effective for node-level tasks with the effect of node feature smoothing, and for graph-level tasks due to its capability in structure differentiation. However, for brain networks, what kinds of graph structures (e.g., communities, subgraphs) are effective beyond the pairwise connections are still unknown. • The small size of neuroimaging datasets may limit the effectiveness and generalization ability of complex deep learning models. Towards these two limitations, we envision several future directions that can be potentially helpful to fully unleash the power of GNNs for brain network analysis:

• Neurology-driven GNN designs: to design the GNN architectures based on neurological understandings of predictive brain signals, especially disease-specific ones. • Pre-training and transfer learning of GNNs: to design techniques that can train complex GNN models across studies and cohorts [107]. Besides, information sharing across different diseases could lead to a better understanding of cross-disorder commonalities.

https://abcdstudy.org/consortium_members/. ABCD consortium investigators designed and implemented the study and/or provided data but did not necessarily participate in the analysis or writing of this report. This manuscript reflects the views of the authors and may not reflect the opinions or views of the NIH or ABCD consortium investigators. The ABCD data repository grows and changes over time. The ABCD data used in this report came from NIMH Data Archive Release 4.0 (DOI 10.15154/1523041). DOIs can be found at https://nda.nih.gov/abcd.