Heterogeneous Graph Convolutional Neural Network via Hodge-Laplacian for
  Brain Functional Data
This study proposes a novel heterogeneous graph convolutional neural network
(HGCNN) to handle complex brain fMRI data at regional and across-region levels.
We introduce a generic formulation of spectral filters on heterogeneous graphs
by introducing the $k-th$ Hodge-Laplacian (HL) operator. In particular, we
propose Laguerre polynomial approximations of HL spectral filters and prove
that their spatial localization on graphs is related to the polynomial order.
Furthermore, based on the bijection property of boundary operators on simplex
graphs, we introduce a generic topological graph pooling (TGPool) method that
can be used at any dimensional simplices. This study designs HL-node, HL-edge,
and HL-HGCNN neural networks to learn signal representation at a graph node,
edge levels, and both, respectively. Our experiments employ fMRI from the
Adolescent Brain Cognitive Development (ABCD; n=7693) to predict general
intelligence. Our results demonstrate the advantage of the HL-edge network over
the HL-node network when functional brain connectivity is considered as
features. The HL-HGCNN outperforms the state-of-the-art graph neural networks
(GNNs) approaches, such as GAT, BrainGNN, dGCN, BrainNetCNN, and Hypergraph NN.
The functional connectivity features learned from the HL-HGCNN are meaningful
in interpreting neural circuits related to general intelligence.
Jinghan Huang, Moo K. Chung, Anqi Qiu
# Heterogeneous Graph Convolutional Neural Network via Hodge-Laplacian for Brain Functional Data

###### Abstract

This study proposes a novel heterogeneous graph convolutional neural network (HGCNN) to handle complex brain fMRI data at regional and across-region levels. We introduce a generic formulation of spectral filters on heterogeneous graphs by introducing the \(k-th\) Hodge-Laplacian (HL) operator. In particular, we propose Laguerre polynomial approximations of HL spectral filters and prove that their spatial localization on graphs is related to the polynomial order. Furthermore, based on the bijection property of boundary operators on simplex graphs, we introduce a generic topological graph pooling (TGPool) method that can be used at any dimensional simplices. This study designs HL-node, HL-edge, and HL-HGCNN neural networks to learn signal representation at a graph node, edge levels, and both, respectively. Our experiments employ fMRI from the Adolescent Brain Cognitive Development (ABCD; n=7693) to predict general intelligence. Our results demonstrate the advantage of the HL-edge network over the HL-node network when functional brain connectivity is considered as features. The HL-HGCNN outperforms the state-of-the-art graph neural networks (GNNs) approaches, such as GAT, BrainGNN, dGCN, BrainNetCNN, and Hypergraph NN. The functional connectivity features learned from the HL-HGCNN are meaningful in interpreting neural circuits related to general intelligence.

## 1 Introduction

Functional magnetic resonance imaging (fMRI) is one of the non-invasive imaging techniques to measure blood oxygen level dependency (BOLD) signals [8]. The fluctuation of fMRI time series signals can characterize brain activity. The synchronization of fMRI time series describes the functional connectivity among brain regions for understanding brain functional organization.

There has been a growing interest in using graph neural network (GNN) to learn the features of fMRI time series and functional connectivity that are relevant to cognition or mental disorders [17, 23].
GNN often considers a brain functional network as a binary undirected graph, where nodes are brain regions, and edges denote which two brain regions are functionally connected. Functional time series, functional connectivity, or graph metrics (i.e., degree, strength, clustering coefficients, participation, etc.) are defined as a multi-dimensional signal at each node. A substantial body of research implements an convolutional operator over nodes of a graph in the spatial domain, where the convolutional operator computes the fMRI feature of each node via aggregating the features from its neighborhood nodes [23, 17]. Various forms of GNN with spatial graph convolution are implemented via 1) introducing an attention mechanism to graph convolution by specifying different weights to different nodes in a neighborhood (GAT, [9]); 2) introducing a clustering-based embedding method over all the nodes and pooling the graph based on the importance of nodes (BrainGNN, [17]); 3) designing an edge-weight-aware message passing mechanism [3]; 4) training dynamic brain functional networks based on updated nodes' features (dGCN, [23]). BrainGNN and dGCN achieve superior performance on Autism Spectrum Disorder (ASD) [17] and attention deficit hyperactivity disorder (ADHD) classification [23]. Graph convolution has also been solved in the spectral domain via the graph Laplacian [2]. For the sake of computational efficiency when graphs are large, the Chebyshev polynomials and other polynomials were introduced to approximate spectral filters for GNN [4, 10]. For large graphs, the spectral graph convolution with a polynomial approximation is computationally efficient and spatially localized [10].

Despite the success of the GNN techniques on cognitive prediction and disease classification [17, 23], the graph convolution aggregates brain functional features only over nodes and updates features for each node of the graph. Nevertheless, signal transfer from one brain region to another is through their connection, which can, to some extent, be characterized by their functional connectivity. The strength of the connectivity determines which edges signals pass through. Therefore, there is a need for heterogeneous graphs with different types of information attached to nodes, such as functional time series and node efficiency, and edges, such as functional connectivity and path length.

Lately, a few studies have focused on smoothing signals through the topological connection of edges [13, 12]. Kawahara et al. [15] proposed BrainNetCNN to aggregate brain functional connectivities among edges. However, brain functional connectivity matrices at each layer are no longer symmetric as the construction nature of the brain functional network. Jo et al. [13] employed a dual graph with the switch of nodes and edges of an original graph so that the GNN approaches described above can be applied (Hypergraph NN). But, the dual graph normally increases the dimensionality of a graph. To overcome this, Jo et al. [13] only considered important edges. Similarly, Jiang et al. [12] introduced convolution with edge-node switching that embeds both nodes and edges to a latent feature space. When graphs are not sparse, the computation of this approach can be intensive. The above-mentioned edge-node switching based model achieved great success on social network and molecular science [13, 12], suggesting that GNN approaches on graph edges have advantages when information is defined
on graph edges. Thus, it is crucial to consider heterogeneous graphs where multiple types of features are defined on nodes, edges, and etc. This is particularly suitable for brain functional data.

This study develops a novel heterogeneous graph convolutional neural network (HGCNN) simultaneously learning both nodes' and edges' functional features from fMRI data for predicting cognition or mental disorders. The HGCNN is designed to learn 1) nodes' features from their neighborhood nodes' features based on the topological connections of the nodes; 2) edges' features from their neighborhood edges' features based on the topological connections of the edges. To achieve these goals, the HGCNN considers a brain functional network as a simplex graph that allows characterizing node-node, node-edge, edge-edge, and higher-order topology. We develop a generic convolution framework by introducing the Hodge-Laplacian (HL) operator on the simplex graph and designing HL-spectral graph filters to aggregate features among nodes or edges based on their topological connections. In particular, this study takes advantage of spectral graph filters in [4, 10] and approximates HL-spectral graph filters using polynomials for spatial locations of these filters. We shall call our HGCNN as HL-HGCNN in the rest of the paper. Unlike the GNNs described above [12, 23], this study also introduces a simple graph pooling approach based on its topology such that the HL can be automatically updated for the convolution in successive layers, and the spatial dimension of the graph is reduced. Hence, the HL-HGCNN learns spectral filters along nodes, edges, or higher-dimensional simplex to extract brain functional features.

We illustrate the use of the HL-HGCNN on fMRI time series and functional connectivity to predict general intelligence based on a large-scale adolescent cohort study (Adolescent Brain Cognitive Development (ABCD), n=7693). We also compare the HL-HGCNN with the state-of-art GNN techniques described above and demonstrate the outstanding performance of the HL-HGCNN. Hence, this study proposes the following novel techniques:

1. a generic graph convolution framework to smooth signals across nodes, edges, or higher-dimensional simplex;
2. spectral filters on nodes, edges, or higher-dimensional simplex via the HL operator;
3. HL-spectral filters with a spatial localization property via polynomial approximations;
4. a spatial pooling operator based on graph topology.

## 2 Methods

This study designs a heterogeneous graph convolutional neural network via the Hodge-Laplacian operator (HL-HGCNN) that can learn the representation of brain functional features at a node-level and an edge-level based on the graph topology. In the following, we will first introduce a generic graph convolution framework to design spectral filters on nodes and edges to learn node-level and edge-level brain functional representation based on its topology achieved via the HL operator. We will introduce the polynomial approximation of the HL spectral
filters to overcome challenges on spatial localization. Finally, we will define an efficient pooling operation based on the graph topology for the graph reduction and update of the HL operator.

### Learning Node-Level and Edge-Level Representation via the Hodge-Laplacian Operator

In this study, the brain functional network is characterized by a heterogeneous graph, \(G=\{V,E\}\) with brain regions as nodes, \(V=\{v_{i}\}_{i=1}^{n}\), and their connections as edges, \(E=\{e_{ij}\}_{i,j=1,2,\cdots,n}\), as well as functional time series defined on the nodes and functional connectivity defined on the edges. This study aims to design convolutional operations for learning the representation of functional time series at nodes and the representation of functional connectivity at edges based on node-node and edge-edge connections (or the topology of graph \(G\)).

Mathematically, nodes and edges are called \(0-\) and \(1-\)dimensional simplex. The topology of \(G\) can be characterized by _boundary operator_\(\boldsymbol{\partial}_{k}\). \(\boldsymbol{\partial}_{1}\) encodes how two \(0\)-dimensional simplices, or nodes, are connecting to form a \(1\)-dimensional simplex (an edge) [6]. In the graph theory [16], \(\boldsymbol{\partial}_{1}\) can be represented as a traditional incidence matrix with size \(n\times n(n-1)/2\), where nodes are indexed over rows and edges are indexed over columns. Similarly, the second order boundary operator \(\boldsymbol{\partial}_{2}\) encodes how \(1\)-dimensional simplex, or edges, are connected to form the connections among \(3\) nodes (\(2\)-dimensional simplex or triangle).

The goal of spectral filters is to learn the node-level representation of fMRI features from neighborhood nodes' fMRI features and the edge-level representation of fMRI features from neighborhood edges' fMRI features. The neighborhood information of nodes and edges can be well characterized by the _boundary operators_\(\boldsymbol{\partial}_{k}\) of graph \(G\). It is natural to incorporate the _boundary operators_ of graph \(G\) in the \(k\)-th Hodge-Laplacian (HL) operator defined as

\[\boldsymbol{\mathcal{L}}_{k}=\boldsymbol{\partial}_{k+1}\boldsymbol{\partial }_{k+1}^{\top}+\boldsymbol{\partial}_{k}^{\top}\boldsymbol{\partial}_{k}. \tag{1}\]

When \(k=0\), the \(0\)-th HL operator is

\[\boldsymbol{\mathcal{L}}_{0}=\boldsymbol{\partial}_{1}\boldsymbol{\partial}_{1} ^{\top} \tag{2}\]

over nodes. This special case is equivalent to the standard Graph Laplacian operator, \(\boldsymbol{\mathcal{L}}_{0}=\Delta\). When \(k=1\), the \(1\)-st HL operator is defined over edges as

\[\boldsymbol{\mathcal{L}}_{1}=\boldsymbol{\partial}_{2}\boldsymbol{\partial}_{2 }^{\top}+\boldsymbol{\partial}_{1}^{\top}\boldsymbol{\partial}_{1}. \tag{3}\]

We can obtain orthonormal bases \(\boldsymbol{\psi}_{k}^{0},\boldsymbol{\psi}_{k}^{1},\boldsymbol{\psi}_{k}^{2},\cdots\) by solving eigensystem \(\boldsymbol{\mathcal{L}}_{k}\boldsymbol{\psi}_{k}^{j}=\lambda_{k}^{j} \boldsymbol{\psi}_{k}^{j}\). We now consider an HL spectral filter \(h\) with spectrum \(h(\lambda_{k})\) as

\[h(\cdot,\cdot)=\sum_{j=0}^{\infty}h(\lambda_{k}^{j})\psi_{k}^{j}(\cdot)\psi_{ k}^{j}(\cdot). \tag{4}\]
A generic form of spectral filtering of a signal \(f\) on the heterogeneous graph \(G\) can be defined as

\[f^{\prime}(\cdot)=h*f(\cdot)=\sum_{j=0}^{\infty}h(\lambda_{k}^{j})c_{k}^{j}\psi_{ k}^{j}(\cdot)\, \tag{5}\]

where \(f(\cdot)=\sum_{j=0}^{\infty}c_{k}^{j}\psi_{k}^{j}(\cdot)\). When \(k=0\), \(f\) is defined on the nodes of graph \(G\). Eq. (5) indicates the convolution of a signal \(f\) defined on \(V\) with a filter \(h\).

Likewise, when \(k=1\), \(f\) is defined on the edges of graph \(G\). Eq. (5) then indicates the convolution of a signal \(f\) defined on \(E\) with a filter \(h\). Eq. (5) is generic that can be applied to smoothing signals defined on higher-dimensional simplices. Nevertheless, this study considers the heterogeneous graph only with signals defined on nodes and edges (0- and 1-dimensional simplices). In the following, we shall denote these two as "HL-node filtering" and "HL-edge filtering", respectively.

### Laguerre Polynomial Approximation of the HL Spectral Filters

The shape of spectral filters \(h\) in Eq. (5) determines how many nodes or edges are aggregated in the filtering process. Our goal of the HL-HGCNN is to design \(h\) such as the representation at nodes and edges are learned through their neighborhood. This is challenging in the spectral domain since it requires \(h(\lambda)\) with a broad spectrum. In this study, we propose to approximate the filter spectrum \(h(\lambda_{k})\) in Eq. (5) as the expansion of Laguerre polynomials, \(T_{p}\), \(p=0,1,2,\ldots,P-1\), such that

\[h(\lambda_{k})=\sum_{p=0}^{P-1}\theta_{p}T_{p}(\lambda_{k})\, \tag{6}\]

where \(\theta_{p}\) is the \(p^{th}\) expansion coefficient associated with the \(p^{th}\) Laguerre polynomial. \(T_{p}\) can be computed from the recurrence relation of \(T_{p+1}(\lambda_{k})=\frac{(2p+1-\lambda_{k})T_{p}(\lambda_{k})-pT_{p-1}( \lambda_{k})}{p+1}\) with \(T_{0}(\lambda_{k})=1\) and \(T_{1}(\lambda_{k})=1-\lambda_{k}\).

We can rewrite the convolution in Eq. (5) as

\[f^{\prime}(\cdot)=h*f(\cdot)=\sum_{p=0}^{P-1}\theta_{p}T_{p}(\mathbf{ \mathcal{L}}_{k})f(\cdot). \tag{7}\]

Analog to the spatial localization property of the polynomial approximation of the graph Laplacian (the 0-th HL) spectral filters [10, 4, 21], the Laguerre polynomial approximation of the 1-st HL spectral filters can also achieve this localization property. Assume two edges, \(e_{ij}\) and \(e_{mn}\), on graph \(G\). The shortest distance between \(e_{ij}\) and \(e_{mn}\) is denoted by \(d_{G}(ij,mn)\) and computed as the minimum number of edges on the path connecting \(e_{ij}\) and \(e_{mn}\). Hence, \((\mathbf{\mathcal{L}}_{1}^{P})_{e_{ij},e_{mn}}=0\quad if\quad d_{G}( ij,mn)>P\), where \(\mathbf{\mathcal{L}}_{1}^{P}\) denotes the \(P\)-th power of the 1-st HL. Hence, the spectral filter represented by the \(P\)-th order Laguerre polynomials of the 1-st HL is localized within the \(P\)-hop edge neighborhood.
Therefore, spectral filters in Eq. (6) have the property of spatial localization. This proof can be extended to the \(k\)-th HL spectral filters. In Section of Results 3, we will demonstrate this property using simulation data.

### Topological Graph Pooling (TGPool)

The pooling operation has demonstrated its effectiveness on grid-like image data [22]. However, spatial graph pooling is not straightforward, especially for heterogeneous graphs. This study introduces a generic topological graph pooling (TGPool) approach that includes coarsening of the graph, pooling of signals, and an update of the Hodge-Laplacian operator. For this, we take an advantage of the one-to-one correspondence between the _boundary operators_ and graph \(G\) and define the three operations for pooling based on the _boundary operators_. As the _boundary operators_ encode the topology of the graph, our graph pooling is topologically based.

For graph coarsening, we generalize the Graclus multilevel clustering algorithm [5] to coarsen the \(k-\)dimensional simplices on graph \(G\). We first cluster similar \(k-\)dimensional simplices based on their associated features via local normalized cut. At each coarsening level, two neighboring \(k-\)dimensional simplices with maximum local normalized cut are matched until all \(k-\)dimensional simplices are explored [19]. A balanced binary tree is generated where each \(k-\)dimensional simplex has either one (i.e., singleton) or two child \(k-\)dimensional simplices. Fake \(k-\)dimensional simplices are added to pair with those singletons. The weights of \(k+1-\)dimensional simplices involving fake \(k-\)dimensional simplices are set as \(0\). The pooling on this binary tree can be efficiently implemented as a simple \(1\)-dimensional pooling of size \(2\). Then, two matched \(k-\)dimensional simplices are merged as a new \(k-\)dimensional simplex by removing the \(k-\)dimensional simplex with the lower degree and the \(k+1-\)dimensional simplices that are connected to this \(k-\)dimensional simplex. To coarsen the graph, we define a new _boundary operator_ by deleting the corresponding rows and columns in the boundary operator and computing the HL operators via Eq.2. Finally, the signal of the new \(k-\)dimensional simplex is defined as the average (or max) of the signals at the two \(k-\)dimensional simplices. Fig.1 illustrates the graph pooling of \(0\)-dimensional and \(1\)-dimensional simplices and the boundary operators of the updated graph after pooling.

Figure 1: Topological Graph Pooling (TGPool). Panels (a) and (b) illustrate the topological graph pooling of (a) \(0\)-dimensional (nodes) and (b) \(1\)-dimensional (edges) simplices. The color at each node or edge indicates features and their similarity across nodes or edges.
### Hodge-Laplacian Heterogeneous Graph Convolutional Neural Network (HL-HGCNN)

We design the HL-HGCNN with the temporal, node, and edge convolutional layers to learn temporal and spatial information of brain functional time series and functional connectivity. Each layer includes the convolution, leaky rectified linear unit (leaky ReLU), and pooling operations. Fig.2 illustrates the overall architecture of the HL-HGCNN model, the temporal, node, and edge convolutional layers.

**Filters.** Denote \(h_{t}\), \(h_{v}\), \(h_{e}\) to be temporal filters, HL-node filters, HL-edge filters, respectively. \(h_{t}\) is a simple 1-dimensional filter along the time domain with different kernel sizes to extract the information of brain functional time series at multiple temporal scales. \(h_{v}\) and \(h_{e}\) are defined in Eq. (6), where \(\theta_{p}\) are the parameters to be estimated in the HL-HGCNN. As mentioned earlier, \(P\) determines the kernel size of \(h_{v}\) and \(h_{e}\) and extracts the higher-order information of the brain functional time series and functional connectivity at multiple spatial scales.

**Leaky ReLU.** This study employs leaky rectified linear unit (ReLU) as an activation function, \(\sigma\), since negative functional time series and functional connectivity are considered biologically meaningful.

**Pooling.** In the temporal convolutional layer, traditional 1-dimensional max pooling operations are applied in the temporal dimension of the functional time series. In the edge and node convolutional layers, TGPool is applied to reduce the dimension of the graph and the dimension of the node and edge signals.

**Output Layer.** We use one more graph convolutional layer to translate the feature of each node or edge into a scalar. Then, we concatenate the vectorized node and edge representations as the input of the output layer. In this study, the output layers contain fully-connected layers.

Figure 2: HL-HGCNN architecture. Panel (A) illustrates the overall architecture of the HL-HGCNN model. Panels (B-D) respectively show the architectures of the HL-edge, temporal, and HL-node convolutional layers.
### Implementation

\(\mathbf{\mathcal{L}}_{0}\) **and \(\mathbf{\mathcal{L}}_{1}\).** Given a brain functional connectivity matrix, we first build a binary matrix while the element in the connectivity matrix with its absolute value greater than a threshold is assigned as one, otherwise zero. We compute the boundary operator \(\mathbf{\partial}_{1}\) with the size of the number of brain regions and the number of functional connectivities. The \(i\)-th row of \(\mathbf{\partial}_{1}\) encodes the functional connection of the \(i\)-th vertex and the \(j\)-th column of \(\mathbf{\partial}_{1}\) encodes how two vertices are connecting to form an edge [6, 7]. Hence, \(\mathbf{\mathcal{L}}_{0}=\mathbf{\partial}_{1}\mathbf{\partial}_{1}^{\top}\).

According to Eq. (1), the computation of \(\mathbf{\mathcal{L}}_{1}\) involves the computation of \(\mathbf{\partial}_{2}\) that characterizes the interaction of edges and triangles. The brain functional connectivity matrix does not form a triangle simplex so the second order boundary operator \(\mathbf{\partial}_{2}=0\). Hence, \(\mathbf{\mathcal{L}}_{1}=\mathbf{\partial}_{1}^{\top}\mathbf{\partial}_{1}\).

**Optimization.** We implement the framework in Python 3.9.13, Pytorch 1.12.1 and PyTorch Geometric 2.1.0 library. The HL-HGCNN is composed of two temporal, node, and edge convolution layers with \(\{8,8\}\), \(\{16,1\}\), and \(\{32,32\}\) filters, respectively. The order of Laguerre polynomials for the 0-th and 1-st HL approximation is set to 3 and 4, respectively. The output layer contains three fully connected layers with 256, 128 and 1 hidden nodes, respectively. Dropout with 0.5 rate is applied to every layer and Leaky ReLU with a leak rate of 0.33 are used in all layers. These model-relevant parameters are determined using greedy search. The HL-HGCNN model is trained using an NVIDIA Tesla V100SXM2 GPU with 32GB RAM by the ADAM optimizer with a mini-batch size of 32. The initial learning rate is set as 0.005 and decays by 0.95 after every epoch. The weight decay parameter was 0.005.

### ABCD Dataset

This study uses resting-state fMRI (rs-fMRI) images from the ABCD study that is an open-sourced and ongoing study on youth between 9-11 years old ([https://abcdstudy.org/](https://abcdstudy.org/)). This study uses the same dataset of 7693 subjects and fMRI preprocessing pipeline stated in Huang et al. [11]. A node represents one of 268 brain regions of interest (ROIs) [18] with its averaged time series as node features. Each edge represents the functional connection between any two ROIs with the functional connectivity computed via Pearson's correlation of their averaged time series as edge features. General intelligence is defined as the average of 5 NIH Toolbox cognition scores, including Dimensional Change Card Sort, Flanker, Picture Sequence Memory, List Sorting Working Memory, and Pattern Comparison Processing Speed [1]. General intelligence ranges from 64 to 123 with mean and standard deviation of \(95.3\pm 7.3\) among 7693 subjects.

## 3 Results

This section first demonstrates the spatial localization property of HL-edge filters in relation to the order of Laguerre polynomials via simulated data. We then demonstrate the use of HL-edge filtering and its use in GNN for predicting fluid intelligence using the ABCD dataset.
### Spatial Localization of the HL-Edge Filtering via Laguerre Polynomial Approximations

We illustrate the spatial location property of the HL-edge filtering by designing a pulse signal at one edge (Fig. 3 (a)) and smoothing it via the HL-edge filter. When applying the HL-edge filter approximated via the \(1^{st}\)-, \(2^{nd}\)-, \(3^{rd}\)-, \(4^{th}\)-order Laguerre polynomials, the filtered signals shown in Fig. 3 (b-e) suggest that the spatial localization of the HL-edge filters is determined by the order of Laguerre polynomials. This phenomenon can also be achieved using multi-layer HL-edge filters where each layer contains HL-edge filters approximated using the \(1^{st}\)-order Laguerre polynomial (see Fig. 3 (f)).

### HL-node vs. HL-edge filters

We aim to examine the advantage of the HL-edge filters over the HL-node filters when fMRI data by nature characterize edge information, such as the functional connectivity. When functional connectivities are defined at a node, they form a vector of the functional connectivities related to this node. In contrast, by nature, the functional connectivity represents the functional connection strength of two brain regions (i.e., edge). Hence, it is a scalar defined at an edge. We design the HL-node network with the two HL-node convolutional layers (see in Fig. 2D) and the output layer with three fully connected layers. Likewise, the HL-edge network with the two HL-edge convolutional layers (see Fig. 2B) and the output layer with three fully connected layers. We employ five-fold cross-validation six times to evaluate the prediction accuracy between predicted and actual general intelligence based on root mean square error (RMSE). Table 1 shows that the HL-edge network has smaller RMSE and performs better than the HL-node network (\(p=1.51\times 10^{-5}\) ). This suggests the advantage of the HL-edge filters when features by nature characterize the weights of edges.

Figure 3: Spatial localization of the HL-edge filtering. Panel (a) shows the simulated signal only occurring at one edge. Panels (b-e) show the signals filtered using the HL-edge filters with the \(1^{st}\)-, \(2^{nd}\)-, \(3^{rd}\)-, \(4^{th}\)-order Laguerre polynomial approximation, respectively. Panel (f) illustrates the signals generated from the HL-edge convolution networks with 4 layers. Each layer consists of the HL-filter approximated using the \(1^{st}\)-order Laguerre polynomial.
### Comparisons with existing GNN methods

We now compare our models with the existing state-of-art methods stated above in terms of the prediction accuracy of general intelligence using the ABCD dataset. The first experiment is designed to compare the performance of the HL-node network with that GAT [9], BrainGNN [17], and dGCN [23]. We adopt the architecture of BrainGNN and dGCN from Li et al. [17] and [23] as both methods were used for fMRI data. The GAT is designed with two graph convolution layers, each consisting of 32 filters and 2-head attention, which is determined via greedy search as implemented in our model. The functional connectivity vector of each region is used as input features. Table 1 suggested that the HL-node network performs better than the GAT (\(p=0.0468\)) and BrainGNN (\(p=0.0195\)), and performs equivalently with dGCN (\(p=0.0618\)).

The second experiment compares the HL-edge network with BrainNetCNN [15] and Hypergraph NN [13]. The Hypergraph NN comprises two graph convolution layers with 32 filters and one hypercluster layer after the first graph convolution layer. The BrainNetCNN architecture follows the design in [15]. Table 1 shows that the HL-edge network has smaller RMSE and performs better than the BrainNetCNN (\(p=4.49\times 10^{-5}\)) and Hypergraph NN (\(p=0.0269\)).

Finally, our HL-HGCNN integrates heterogeneous types of fMRI data at nodes and edges. Table 1 shows that the HL-HGCNN performs the best compared to all the above methods (all \(p<0.03\)).

\begin{table}
\begin{tabular}{c|c c c} \hline \hline  & **GNN model** & **RMSE** & \(p\)**-value** \\ \hline \multirow{4}{*}{GNN with node filtering} & **HL-Node network (ours)** & \(7.134\pm 0.011\) & \(4.01\times 10^{-6}\) \\  & **GAT**[9] & \(7.165\pm 0.020\) & \(1.91\times 10^{-5}\) \\  & **BrainGNN[17]** & \(7.144\pm 0.013\) & \(1.51\times 10^{-6}\) \\  & **dGCN[17]** & \(7.151\pm 0.012\) & \(9.83\times 10^{-6}\) \\ \hline \multirow{4}{*}{GNN with edge filtering} & **HL-Edge network (ours)** & \(7.009\pm 0.012\) & \(2.48\times 10^{-2}\) \\  & **BrainNetCNN[15]** & \(7.118\pm 0.016\) & \(5.34\times 10^{-6}\) \\ \cline{1-1}  & **Hypergraph NN[13]** & \(7.051\pm 0.022\) & \(3.74\times 10^{-5}\) \\ \hline GNN with node and edge filtering & **HL-HGCNN (ours)** & **6.972\(\pm\)0.015** & - \\ \hline \hline \end{tabular}
\end{table}
Table 1: General intelligence prediction accuracy based on root mean square error (RMSE). \(p\)-value is obtained from two-sample \(t\)-tests examining the performance of each method in reference to the proposed HL-HGCNN.

Figure 4: The saliency map of the brain functional connectivity. Red boxes highlight brain networks with higher weights, indicating greater contributions to the prediction of general intelligence.
### Interpretation

We use the graph representation of the final edge convolution layer of the HL-HGCNN to compute the saliency map at the connectivity level. The group-level saliency map is computed by averaging the saliency maps across all the subjects in the dataset. The red boxes in Fig.4 highlight the functional connectivities of the occipital regions with the prefrontal, parietal, salience, and temporal regions that most contribute to general intelligence. Moreover, our salience map also highlights the functional connectivities of the right prefrontal regions with bilateral parietal regions, which is largely consistent with existing findings on neural activities in the frontal and parietal regions [14, 20].

## 4 Conclusion

This study proposes a novel HL-HGCNN on fMRI time series and functional connectivity for predicting cognitive ability. Our experiments demonstrate the spatial localization property of HL spectral filters approximated via Laguerre polynomials. Moreover, our HL-node, HL-edge, and HL-HGCNN perform better than the existing state-of-art methods for predicting general intelligence, indicating the potential of our method for future prediction and diagnosis based on fMRI. Nevertheless, more experiments on different datasets are needed to further validate the robustness of the proposed model. Our method provides a generic framework that allows learning heterogeneous graph representation on any dimensional simplices, which can be extended to complex graph data. The HL-HGCNN model offers an opportunity to build high-order functional interaction among multiple brain regions, which is our future research direction.

**Acknowledgements.** This research/project is supported by the Singapore Ministry of Education (Academic research fund Tier 1) and A*STAR (H22P0M0007). Additional funding is provided by the National Science Foundation MDS-2010778, National Institute of Health R01 EB022856, EB02875. This research was also supported by the A*STAR Computational Resource Centre through the use of its high-performance computing facilities.
