A plug-in graph neural network to boost temporal sensitivity in fMRI
  analysis
Learning-based methods have recently enabled performance leaps in analysis of
high-dimensional functional MRI (fMRI) time series. Deep learning models that
receive as input functional connectivity (FC) features among brain regions have
been commonly adopted in the literature. However, many models focus on
temporally static FC features across a scan, reducing sensitivity to dynamic
features of brain activity. Here, we describe a plug-in graph neural network
that can be flexibly integrated into a main learning-based fMRI model to boost
its temporal sensitivity. Receiving brain regions as nodes and
blood-oxygen-level-dependent (BOLD) signals as node inputs, the proposed
GraphCorr method leverages a node embedder module based on a transformer
encoder to capture temporally-windowed latent representations of BOLD signals.
GraphCorr also leverages a lag filter module to account for delayed
interactions across nodes by computing cross-correlation of windowed BOLD
signals across a range of time lags. Information captured by the two modules is
fused via a message passing algorithm executed on the graph, and enhanced node
features are then computed at the output. These enhanced features are used to
drive a subsequent learning-based model to analyze fMRI time series with
elevated sensitivity. Comprehensive demonstrations on two public datasets
indicate improved classification performance and interpretability for several
state-of-the-art graphical and convolutional methods that employ
GraphCorr-derived feature representations of fMRI time series as their input.
Irmak Sivgin, Hasan A. Bedel, Şaban Öztürk, Tolga Çukur
# A plug-in graph neural network to boost temporal sensitivity in fMRI analysis

###### Abstract

Learning-based methods have recently enabled performance leaps in analysis of high-dimensional functional MRI (fMRI) time series. Deep learning models that receive as input functional connectivity (FC) features among brain regions have been commonly adopted in the literature. However, many models focus on temporally static FC features across a scan, reducing sensitivity to dynamic features of brain activity. Here, we describe a plug-in graph neural network that can be flexibly integrated into a main learning-based fMRI model to boost its temporal sensitivity. Receiving brain regions as nodes and blood-oxygen-level-dependent (BOLD) signals as node inputs, the proposed GraphCorr method leverages a node embedder module based on a transformer encoder to capture temporally-windowed latent representations of BOLD signals. GraphCorr also leverages a lag filter module to account for delayed interactions across nodes by computing cross-correlation of windowed BOLD signals across a range of time lags. Information captured by the two modules is fused via a message passing algorithm executed on the graph, and enhanced node features are then computed at the output. These enhanced features are used to drive a subsequent learning-based model to analyze fMRI time series with elevated sensitivity. Comprehensive demonstrations on two public datasets indicate improved classification performance and interpretability for several state-of-the-art graphical and convolutional methods that employ GraphCorr-derived feature representations of fMRI time series as their input.

functional MRI, time series, neural network, graph, classification, connectivity

## 1 Introduction

The human brain comprises networks of regions that interactively process information during cognitive processing [1]. In turn, correlated activity within individual functional networks has been associated with unique mental states [2; 3]. Functional MRI (fMRI) is a powerful modality to examine functional networks as it can non-invasively measure whole-brain blood-oxygen-level-dependent (BOLD) signals consequent to neural activity at high spatio-temporal resolution [4; 5]. In fMRI studies, functional connectivity (FC) measures are used to assess similarity of BOLD signals among brain regions [6; 7; 8; 9; 10]. The traditional approach to map FC measures onto mental states is then based on conventional methods such as logistic regression and support vector machines (SVM) [11; 12; 13; 14]. Unfortunately, conventional methods are often insufficiently sensitive to the intricate information patterns in whole-brain fMRI time series [15].

In recent years, the success of deep learning (DL) models at exploring features in high-dimensional datasets has motivated their adoption for fMRI analysis as an alternative to conventional methods [16; 17; 18; 19; 20]. Earlier attempts in this domain have proposed shallow multi-layer perceptron (MLP) [21; 22] and Boltzmann machine (BM) models [23; 16]. Later studies have adopted deeper architectures based on convolutional neural network (CNN) [24; 17; 25], graph neural network (GNN) [26; 18; 27; 28; 29; 30], and transformer [31; 32; 33; 34; 35] models for improved performance. Typically, these models start by constructing construct a set of nodes corresponding to brain regions defined based on an anatomical or functional atlas [36; 37; 38], and receive input features at these nodes based on the FC strength among brain regions [39; 18]. A common approach has been to employ static FC features for derived from aggregate correlation measures across the entire duration of fMRI time series [40; 18]. Yet, this approach is insufficiently
sensitive to the dynamic inter-regional interactions in the human brain during resting-state or cognitive tasks [41]. While alternative strategies have recently been proposed to assess the temporal variability in FC features, these methods commonly consider instantaneous signal correlations across local time windows within the time series [42; 43; 44]. As such, they do not possess explicit mechanisms to capture delayed correlations between brain regions that can be present in fMRI times series due to hierarchical cognitive processing in the brain or hemodynamic lags in BOLD measurements [45].

In this study, we introduce a plug-in graphical neural network, GraphCorr, that provides enhanced input features to learning-based fMRI models so as to boost their sensitivity to dynamic, lagged inter-regional interactions. To capture dynamic changes in interactions, GraphCorr leverages a novel node embedder module based on a transformer encoder that computes hierarchical embeddings of windowed BOLD signals across the time series. To capture lagged interactions between brain regions, GraphCorr employs a novel lag filter module that computes nonlinear features of cross-correlation between pairs of nodes across a range of time delays. The graph model is initialized with node features taken as embeddings from the node embedder module, and with edge weights taken as lag features from the lag filter module. Afterwards, a message passing algorithm is used to compute enhanced node embeddings that account for dynamic, lagged inter-regional interactions.

Here, we demonstrate GraphCorr for gender classification from fMRI scans in two public datasets: the Human Connectome Project (HCP) dataset [46] and the ID1000 dataset from Amsterdam Open MRI Collection (AOMIC) [47]. GraphCorr is coupled as a plug-in to state-of-the-art baseline models for fMRI analysis including SAGE [48], BrainGNN [18], BrainNetCNN [17] and GCN [29]. Significantly enhanced performance is obtained from each baseline model when coupled with GraphCorr. We devise an explanatory analysis approach for GraphCorr to interpret the time frames and brain regions that most significantly contribute to classification decisions. We show that GraphCorr improves explainability of baseline models, resulting in interpretations that are more closely aligned with prominent neuroscientific findings from the literature. We also demonstrate the benefits of GraphCorr-derived features against features extracted via plug-in recurrent neural networks (RNN) and dynamic FC features computed directly from BOLD signals.

## 2 Related Work

Cognitive processes elicit broadly distributed response patterns across the human brain [49]. In turn, process-related information such as stimulus or task variables can be decoded by analyzing resultant multi-variate BOLD signals [50; 51]. Initial studies in this domain employed relatively simpler, traditional machine learning (ML) methods for fMRI analysis [41; 52]. These traditional methods rely heavily on feature selection procedures to cope with the intrinsically high dimensionality of fMRI data [53]. Arguably, FC features among brain regions have been most commonly used to capture discriminative information about cognitive processes [21; 17; 39; 18]. Many studies have reported that external variables or disease states can be detected given FC features of individual subjects under resting state [54; 55], cognitive tasks [11; 18], or both [14].

Given their earlier success in fMRI analysis, FC features have also been pervasively adopted in recent DL methods that leverage more complex models to enhance performance [56; 57]. A pervasive approach in DL-based fMRI analysis relies on static FC features as model inputs, where FC between a pair of regions is taken as the aggregate correlation of their BOLD signals across the entire scan. To extract hierarchical latent representations of these features, earlier studies have proposed either relatively compact fully-connected architectures including BM and MLP models [23; 16; 21; 22], or computation-efficient deep architectures including CNN models [25; 17]. Later studies have considered GNN models given their natural fit to analyzing fMRI data that follows an intrinsic connectivity structure [58; 59; 60; 61; 18]. These DL methods have all enabled substantial performance improvements in fMRI analysis over traditional methods. Yet, analyses rooted in static FC features can still yield suboptimal sensitivity to fine-grained temporal information across the fMRI time series [12; 25].

To improve temporal sensitivity, several alternative strategies have been proposed that incorporate time-varying features to DL models for fMRI analysis. A first group of methods pre-compute FC features over moving windows across the time series based on standard correlation measures, and concatenate them across windows to form a higher-dimensional input [62; 63; 28; 43; 44]. While these dynamic FC features promise enhanced temporal sensitivity, they result in elevated complexity due to their intrinsic dimensionality that can degrade model performance. A second group of methods instead provide voxel-level BOLD signals spatially encoded via a CNN module. Spatially-encoded BOLD signals are then processed with RNN or transformer models to extract the time-varying information [32; 42]. Yet, CNN modules based on voxel-level inputs can be difficult to train from scratch under limited data regimes. A third group of methods retain static FC features as their input, albeit augment them with dynamic features captured by RNN modules that directly encode BOLD signals [42]. Besides elevated model complexity, these methods can suffer from intrinsic limitations of RNNs in terms of vanishing/exploding gradients over the extensive number of time steps in typical fMRI scans [64; 65]. Importantly, a common attribute of these previous approaches is that they primarily consider the temporal variation in instantaneous correlations among brain regions. However, they can elicit suboptimal sensitivity as they lack explicit mechanisms to capture delayed inter-regional interactions that can occur due to hierarchical processing or hemodynamic delays [45].

Here, we propose to improve the temporal sensitivity of downstream fMRI analysis models by integrating a novel plug-in GNN, GraphCorr. The proposed GraphCorr method uses a novel node embedder module to find contextual embeddings of dynamic FC features based on instantaneous correlations of windowed BOLD signals; and it uses a novel lag filter module to compute embeddings of cross-correlation features of windowed BOLD signals across various time delays. Following a message passing algorithm across the graph, GraphCorr provides
enhanced input features that preserve dynamic, delayed correlations among brain regions to a downstream analysis model so as to improve its performance. Unlike methods based on static FC features [17; 18], GraphCorr leverages dynamic FC features to capture the variability in connectivity among brain regions. Unlike methods that receive multiple sets of dynamic FC features across separate time windows [66; 42; 28], GraphCorr fuses its node features across time windows to lower model complexity without sacrificing performance. Unlike methods that employ recurrent architectures that involve sequential processing [42], GraphCorr leverages a transformer encoder on dynamic FC features that enables efficient parallel processing. Unlike methods that solely focus on instantaneous signal correlations [17], GraphCorr adopts an explicit lag filter mechanism to learn delayed cross-correlations among brain regions.

## 3 GraphCorr

Analysis procedures for fMRI time series typically start by defining a collection of \(R\) regions of interest (ROI) across the brain based on an anatomical atlas [18; 42]. Voxel-level BOLD signals within each ROI are then averaged to derive ROI-level signals, resulting in \(\mathbf{B}\in\mathbb{R}^{R\times T}\) as the matrix of BOLD signals where \(T\) denotes the number of time frames. Static FC features are conventionally computed based on Pearson's correlation coefficient of these BOLD signals across ROIs: \(\mathbf{sFC}_{i,j}=\text{Corr}(\mathbf{B}_{i,\cdot},\mathbf{B}_{j,\cdot})\), where \(\mathbf{sFC}\in\mathbb{R}^{R\times R}\) and \(i,j\) are ROI indices. Many previous traditional and learning-based methods use downstream classification models on static FC features, which result in suboptimal temporal sensitivity. Instead, here we propose to extract dynamic FC features of BOLD signals based on a novel GNN plug-in, and to use these enhanced features to improve the performance of downstream classification models. The proposed GraphCorr method forms a graph structure to represent brain connectivity features, leverages node embedder and lag filter modules to capture dynamic, lagged correlations in BOLD signals, and finally performs message passing on the graph to compute enhanced features (Fig. 1). The methodological components and procedures in GraphCorr are described below.

### Graph formation

As a learning substrate, GraphCorr first forms a graph \(G(N,E)\) with \(N\) and \(E\) denoting nodes and edges, respectively. The node set \(N=\{r_{i}\,|\,i=1,...,R\}\) includes ROIs defined according to the atlas, whereas the binary edge set is given as \(E=\{e_{i,j}=1\,|\,i=1,...,R;j\in\mathcal{N}(i)\}\) where \(\mathcal{N}(i)\) is the neighborhood of the \(i\)-th node. Edges are defined by thresholding to retain the strongest \(z\%\) of correlation coefficients in \(\mathbf{sFC}\) while excluding self connections, resulting in \(E_{n}\) number of edges. The node features \(\mathbf{F}=\{f_{i}\,|\,i=1,...,R\}\) are initialized as the time-windowed BOLD signals at each corresponding node to capture local dynamics in the fMRI times series. For this purpose, the scan containing \(T\) time frames is split into \(W\) windows of size \(T_{w}\) and stride value \(s\):

\[W=\lfloor\frac{T-T_{w}}{s}\rfloor, \tag{1}\]

resulting in a feature tensor of \(\mathbf{F}\in\mathbb{R}^{R\times T_{w}\times W}\).

### Architecture

To sensitively extract dynamic connectivity features, GraphCorr utilizes a novel node embedding module. To also capture delayed connectivity features, GraphCorr utilizes a novel lag filter module. The two modules are detailed below.

**Node embedder module:** Receiving as input time-windowed BOLD signals, this module computes latent representations of dynamic FC features (Fig. 2). First, dynamic FC features are extracted from the time-windowed BOLD signals as: \(\mathbf{FC}_{i,j,w}=\text{Corr}(\mathbf{F}_{i,\cdot,w},\mathbf{F}_{j,\cdot,w})\) where \(w\in\{1,...,W\}\) indicates window index, \(i\in\{1,...,R\}\), \(j\in\{1,...,R\}\) denote node indices. These FC features are then processed with a transformer encoder where windows across the fMRI time series correspond to the sequence of transformer tokens. Attention calculations are performed on window-specific keys \(K_{w}\in\mathbb{R}^{R\times d}\), queries \(Q_{w}\in\mathbb{R}^{R\times d}\) and values \(V_{w}\in\mathbb{R}^{R\times d}\) derived via learnable linear projections \(f_{q},f_{k}\) and \(f_{\cdot}\):

\[Q_{w} =U_{q}(|\mathbf{FC}_{1,\cdot,w},\mathbf{FC}_{2,\cdot,w},..., \mathbf{FC}_{\mathcal{K},w}|),\] \[K_{w} =U_{k}(|\mathbf{FC}_{1,\cdot,w},\mathbf{FC}_{2,\cdot,w},..., \mathbf{FC}_{\mathcal{K},w}|),\] \[V_{w} =U_{\tau}(|\mathbf{FC}_{1,\cdot,w},\mathbf{FC}_{2,\cdot,w},..., \mathbf{FC}_{\mathcal{K},w}|), \tag{2}\]

where \(d\) is the dimensionality of each attention head. The above computations can be performed separately for \(H\) attention heads. The window-specific attention matrix \(\mathbf{A}_{w}\in\mathbb{R}^{R\times R}\) is then derived as [67]:

\[\mathbf{A}_{w}=\text{Att}(Q_{w},K_{w},V_{w})=\text{Softmax}(\frac{Q_{w}K_{w}^{ \intercal}}{\sqrt{d}})V_{w}. \tag{3}\]

Attention matrices are concatenated across attention heads and propagated to an MLP block following layer normalization:

\[\mathbf{EMB}_{w}=\text{MLP}(\mathbf{A}_{w})=\text{GELU}(\mathbf{A}_{w}\mathbf{ M}_{1})\mathbf{M}_{2} \tag{4}\]

where \(\mathbf{M}_{1}\in\mathbb{R}^{(R\times D)}\) and \(\mathbf{M}_{2}\in\mathbb{R}^{(D\times D)}\) denote MLP model parameters, GELU is the Gaussian activation unit, and \(\mathbf{EMB}\in\mathbb{R}^{R\times D\times W}\) are window-specific node embeddings, and \(D\) is the embedding dimensionality with \(D<R\).

**Lag filter module:** Receiving as input time-windowed BOLD signals, this module computes cross-correlation features across a range of temporal delays (Fig. 3). For this purpose, initial node features from the graph formation stage are zero-padded across the time dimension:

\[\mathbf{X}_{i,w}=[\mathbf{0}_{(1\times m)},\mathbf{F}_{i,\cdot,w},\mathbf{0}_{ (1\times m)}]. \tag{5}\]

where \(\mathbf{X}\in\mathbb{R}^{R\times(T_{w}+2m)\times W}\), and \(m\) defines the range of delays \(\tau\in\{-m,-m+1,...,m-1,m\}\) that will be considered in the module. First, cross-correlations are computed between pairs of nodes connected by \(e_{ij}\) at each lag value separately:

\[\rho_{i,j,w,\tau}=\text{Corr}(\mathbf{X}_{i,\cdot,w},\mathbf{X}_{j,\cdot,w},\tau) \tag{6}\]

Afterwards, learnable lag filters \(\mathbf{P}_{LF}\in\mathbb{R}^{(2m+1)\times k}\) with \(k\) denoting the number of filters are used to map cross-correlations onto lag activations:

\[\mathbf{LAG}_{w}=\text{GELU}(\rho_{\cdot,\cdot,w},\mathbf{P}_{LF}) \tag{7}\]

where \(\mathbf{LAG}\in\mathbb{R}^{E_{w}\times W\times d}\) are window-specific lag activations.
### Graph learning

The node embedder produces time-windowed node embeddings, **EMB**, that reflect instantaneous inter-regional correlations. The lag filter produces time-windowed lag activations, **LAG**, that reflect reflect delayed inter-regional correlations. To consolidate these feature sets on the graph, node embeddings are taken as node features and lag activations are taken as edge weights (Fig. 1). A message passing algorithm is then run on the graph to compute enhanced FC features. To do this, a message tensor \(\textbf{MES}\in\mathbb{R}^{R_{\textit{c}}(D\times k)\times W}\) is computed between pairs of connected nodes \((r_{i},r_{j})\) as:

\[\textbf{MES}_{i,j,w}=\textbf{EMB}_{j,w}\textbf{LAG}_{i,j,w}^{\intercal} \tag{8}\]

Messages are first averaged across windows, and then propagated to a target node \(r_{i}\) that sums all messages from one-hop vicinity nodes \(j\in\mathcal{N}(i)\):

\[\textbf{AGG}_{i}=\sum_{j\in\mathcal{N}(i)}\frac{1}{W}\sum_{w=1}^{W}\textbf{MES }_{i,j,w}. \tag{9}\]

Figure 1: Overview of GraphCorr. **A.** GraphCorr utilizes two parallel modules to extract dynamic, lagged features of inter-regional correlations across the brain. The node embedder module receives as input time-windowed BOLD signals, and uses a transformer encoder to compute node embeddings of dynamic FC features \(\textbf{EMB}\in\mathbb{R}^{R\times D\times W}\). The lag filter module also receives as input time-windowed BOLD signals, and it computes lag activations due to cross-correlation across a range of lag values \(\textbf{LAG}\in\mathbb{R}^{R_{\textit{c}}\times W\times W}\). Cross-correlation is calculated only for connected node pairs \((e_{i,j}=1)\). **B.** To consolidate the extracted feature sets on a graph, node embeddings are taken as node features and lag activations are taken as edge weights. A message passing algorithm is then run on the graph to produce enhanced FC features in an output feature matrix, \(\textbf{OUT}\in\mathbb{R}^{R\times D\times(D\times k+1)}\).
This aggregate message is the concatenated with the window-averaged node embedding at \(r_{i}\):

\[\mathbf{OUT}_{i}=[\frac{1}{W}\sum_{v=1}^{W}\mathbf{EMB}_{i_{v},w},\mathbf{AGG}_{ i}], \tag{10}\]

where \(\mathbf{OUT}_{i}\in\mathbb{R}^{D\times(k+1)}\) denotes enhanced features for \(r_{i}\).

## 4 Methods

### Experimental procedures

Demonstrations were performed on fMRI data from the HCP S1200 release1[46] and ID1000 dataset from Amsterdam Open MRI Collection (AOMIC)2[47]. In the HCP dataset, preprocessed data from resting-state fMRI scans were analyzed. The first resting-state scan among four sessions was selected for each subject, excluding short scans with \(T<1200\). This resulted in a total of 1093 healthy subjects (594 female and 499 male). In the ID1000 dataset, preprocessed data from task-based fMRI scans recorded during movie watching were analyzed. All scans had a fixed duration of \(T=240\). A total of 881 healthy subjects were examined (458 female and 423 male). For both datasets, two alternative ROI definitions were considered, based on either the Schaefer atlas [68] or the AAL atlas [36]. The Schaefer atlas includes \(R=400\) ROIs within 7 intrinsic networks, whereas the AAL atlas defines \(R=116\) ROIs.

Footnote 1: [https://db.humanconnectome.org](https://db.humanconnectome.org)

Experiments were conducted on a single NVIDIA Titan Xp GPU using the PyTorch framework. A nested cross-validation procedure was performed with 5 outer and 1 inner folds Domain adaptation procedures can be employ to improve reliability [69]. Data were three-way split into a training set (70%), a validation set (10%) and a test set (20%) with no subject overlap between the sets. For fair comparison, all models were trained, validated and tested on identical data splits. All models were trained based on cross-entropy loss. For each model, hyperparameters were selected to maximize the average performance across the validation sets. A common set of hyperparameters that were observed to yield near-optimal performance were used across datasets and atlases [70]. Details regarding model implementations are discussed in Section 4.2.

### Comparative analysis

GraphCorr was demonstrated on several learning-based methods taken as downstream classification models including SAGE [48], GCN [29], BrainGNN [18], and BrainNetCNN [17]. For each method, a vanilla downstream model was trained by providing static FC features as model input, and an augmented downstream model was separately trained where GraphCorr was employed as a plug in to provide model input. Vanilla and augmented models were obtained with identical training procedures. In all graph models, ROIs in a given brain atlas were taken as nodes, and edge selection was then performed based on correlations of BOLD signals. Edges whose correlation coefficients were in the top \(z=2\%\) were retained, while remaining edges were discarded. The implementation details of the downstream models and GraphCorr are discussed below.

**SAGE:** A GNN model was built based on a module with graph convolution pooling, and fully-connected layers [48]. SAGE comprised a cascade of two graphical modules with hidden dimension of 250 and dropout rate of 0.5. Cross-validated hyperparameters were a learning rate of \(3\times 10^{-3}\), 20 epochs, a batch size of 12.

**GEN:** GCN is a GNN model based on graph convolution, pooling and fully-connected layers [29]. GCN comprised a cascade of two graphical modules with hidden dimension of 100 and dropout rate of 0.5. Cross-validated hyperparameters were a learning rate of \(5\times 10^{-3}\), 30 epochs, a batch size of 12.

**BrainGNN:** BrainGNN is a GNN model based on ROI-aware graph convolution, pooling and fully-connected layers [18]. A single graphical module with hidden dimension of 100 and dropout rate of 0.5 was used. Cross-validated hyperparameters were a learning rate of \(8\times 10^{-4}\), 80 epochs, a batch size of 16.

**BrainNetCNN:** BrainNetCNN is a CNN model based on convolutional layers with edge-to-edge and edge-to-node filters [17]. The convolutional layers had hidden dimension of 32 and dropout rate of 0.1. Vanilla BrainNetCNN expects a 2D input of size \(R\times R\) taken as the static FC matrix. When it was augmented with GraphCorr, its input dimensionality was modified as \(R\times D(k+1)\) for compatibility. Cross-validated hyperparameters were a learning rate of \(2\times 10^{-4}\), 20 epochs, a batch size of 16.

Figure 2: The node embedder module. The time-windowed FC features \(\mathbf{FC}_{\cdot,w}\in\mathbb{R}^{R\times R}\) at window \(w\) are processed with a transformer encoder with multi-head self-attention (MHSA), layer normalization, and multi-layer perceptron (MLP) layers. The output is a node embedding matrix \(\mathbf{EMB}_{\cdot,w}\in\mathbb{R}^{R\times D}\) where \(D<R\) denotes the embedding dimensionality.
**GraphCorr:** The node embedder module was built with a single-layer transformer encoder. Because the scan durations differed across HCP and ID1000, dataset-specific \(T_{w}\) (window size) and \(s\) (stride) were selected while common \(m\) (maximum lag) and \(k\) (filter count) were used. Accordingly, cross-validated parameters were (\(T_{w}\)=50, \(s\)=30, \(m\)=5, \(k\)=3) for HCP, (\(T_{w}\)=40, \(s\)=15, \(m\)=5, \(k\)=3) for ID1000.

### Explanatory analysis

To assess the influence of GraphCorr on interpretability, the vanilla and augmented versions of trained downstream models were examined. An explanation procedure was devised to identify the brain regions within the fMRI times series that most saliently contribute to the model decisions. First, a gradient-based approach was used to compute a saliency tensor summarizing inter-regional interactions [71; 39]. For vanilla models, gradients were computed with respect to static FC features **sFC**:

\[\mathbf{SAL}_{i,j}^{\text{sum}}=|\nabla_{\mathbf{sFC}_{i,j}}y_{\text{sum}}| \tag{11}\]

where \(\mathbf{SAL}^{\text{sum}}\in\mathbb{R}^{Re\times R}\) and \(y_{\text{sum}}\) denotes the model prediction. For augmented models, gradients were computed with respect to time-windowed FC features \(\mathbf{FC}_{i,j\omega}\):

\[\mathbf{SAL}_{i,j\omega}^{\text{ang}}=|\nabla_{\mathbf{FC}_{i,j\omega}}y_{ \text{sum}}| \tag{12}\]

where \(\mathbf{SAL}^{\text{ang}}\in\mathbb{R}^{Re\times R\times W}\) and \(y_{\text{ang}}\) denotes the model prediction. Afterwards, an ROI-specific saliency score was computed by aggregating values across windows and interacting ROI dimensions of the saliency tensor:

\[\mathbf{rSAL}^{\text{sum}}=\sum_{j=1}^{R}\mathbf{SAL}_{\cdot,j}^{\text{sum}} \tag{13}\]

\[\mathbf{rSAL}^{\text{ang}}=\sum_{j=1}^{R}(\frac{1}{W}\sum_{w=1}^{W}\mathbf{SAL }_{\cdot,j\omega}^{\text{ang}}) \tag{14}\]

where \(\mathbf{rSAL}^{\text{sum},\text{ang}}\in\mathbb{R}^{R}\). For saliency assessment at the level of functional brain networks, the seven intrinsic brain networks defined within the Schaefer atlas were used. For each network, ROI-specific saliency scores were averaged across the regions within the network to obtain a network saliency score per hemisphere.

While unsigned ROI-specific saliency scores reflect the relative importance of each region on the model decision, they do not indicate whether the model output is driven by an increase or decrease in BOLD signals within the ROI. To address this question, a post-hoc logistic regression analysis was conducted. First, important windows in the fMRI time series were determined by aggregating values in the saliency tensor across ROI dimensions:

\[\mathbf{wSAL}^{\text{ang}}=\sum_{i=1}^{R}\sum_{j=1}^{R}\mathbf{SAL }_{i,j,\cdot}^{\text{ang}} \tag{15}\] \[w^{*}=\underset{w}{\arg\max}\ \mathbf{wSAL}^{\text{ang}} \tag{16}\]

Here, \(\mathbf{wSAL}^{\text{ang}}\in\mathbb{R}^{W}\) denotes the window-specific saliency score used for important window selection. BOLD signals within the most important window were extracted, and thresholded according to intensity to select the top 5 time frames [72; 73; 35]. A logistic regression model was then fit to map the BOLD signal vector across ROIs onto the output class, i.e., performing the same task as the downstream model. The logistic model returns a weight for each ROI: a positive weight indicates that an increase whereas a negative weight indicates that a decrease in the ROI's BOLD signal elicits the downstream model's decision.

\begin{table}
\begin{tabular}{l l l l l l} \hline \hline \multirow{2}{*}{Model} & \multicolumn{3}{c}{HCP} & \multicolumn{2}{c}{ID1000} \\ \cline{3-6}  & & Acc (\%) & ROC (\%) & Acc (\%) & ROC (\%) \\ \hline \multirow{2}{*}{SAGE} & Vanilla & 75.2 \(\pm\) 2.84 & 85.29 \(\pm\) 1.67 & 62.39 \(\pm\) 2.17 & 68.59 \(\pm\) 3.76 \\  & Augmented & **89.57**\(\pm\) 0.68 & **94.27**\(\pm\) 1.98 & **81.7**\(\pm\) 2.67 & **87.02**\(\pm\) 2.10 \\ \hline \multirow{2}{*}{GCN} & Vanilla & 79.14 \(\pm\) 2.93 & 86.00 \(\pm\) 1.47 & 67.84 \(\pm\) 2.95 & 71.78 \(\pm\) 3.88 \\  & Augmented & **89.94**\(\pm\) 2.18 & **94.52**\(\pm\) 1.61 & **80.80**\(\pm\) 0.98 & **87.90**\(\pm\) 1.75 \\ \hline \multirow{2}{*}{BrainGNN} & Vanilla & 72.83 \(\pm\) 1.98 & 78.85 \(\pm\) 2.45 & 62.50 \(\pm\) 1.80 & 65.63 \(\pm\) 2.81 \\  & Augmented & **84.72**\(\pm\) 1.33 & **92.97**\(\pm\) 0.98 & **79.32**\(\pm\) 1.96 & **88.10**\(\pm\) 2.19 \\ \hline \multirow{2}{*}{BrainNetCNN} & Vanilla & 82.52 \(\pm\) 2.80 & 91.23 \(\pm\) 1.21 & 75.45 \(\pm\) 2.01 & 83.65 \(\pm\) 2.05 \\  & Augmented & **88.47**\(\pm\) 2.63 & **94.71**\(\pm\) 1.86 & **82.73**\(\pm\) 1.63 & **89.85**\(\pm\) 2.24 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance of downstream models on the HCP and ID1000 datasets with the Schaefer atlas. Results are listed as mean+std across test folds for vanilla and GraphCorr-augmented versions. Boldface indicates the better performing version of each model.

Figure 3: The lag filter module. Cross-correlation of time-windowed BOLD signals at window \(w\) are computed for delays \(\tau\in\{-m,-m+1,...,m-1,m\}\), where \(m\) defines the range. This computation is only performed for pairs of connected nodes (\(e_{i,j}=1\)). Afterwards, cross-correlation values \(\rho_{\cdot,\cdot,\cdot,\cdot}\in\mathbb{R}^{Re\times(2m+1)}\) are linearly transformed onto with a learnable filter \(M_{LF}\in\mathbb{R}^{(2m+1)\times 4}\) onto window-specific lag activations \(\mathbf{LAG}_{\cdot,\cdot,w}\in\mathbb{R}^{Re\times R\times R}\).
## 5 Results

### Comparative analysis

GraphCorr was demonstrated on downstream classification models based on SAGE [48], GCN [29], BrainGNN [18], and BrainNetCNN [17]. A gender detection task was performed given resting-state fMRI scans in individual subjects. Performances of vanilla and augmented versions of downstream models on HCP and ID1000 datasets are listed in Table 1 for the Schaefer atlas, and in Table 2 for the AAL atlas. In all examined cases, augmentation with GraphCorr significantly enhances the performance of downstream models (p\(<\)0.05, Wilcoxon signed-rank test). When ROIs are defined via the Schaefer atlas, GraphCorr enables (accuracy, ROC)% improvements of (14.37, 8.98)% for SAGE, (10.80, 8.52)% for GCN, (11.89, 14.12)% for BrainGNN, and (5.95, 3.48)% for BrainNetCNN on HCP; and it enables improvements of (19.31, 18.43)% for SAGE, (13.32, 16.12)% for GCN, (16.82, 22.47)% for BrainGNN, and (7.28, 6.20)% for BrainNetCNN on ID 1000. When ROIs are defined via the AAL atlas, GraphCorr enables improvements of (17.19, 15.54)% for SAGE, (14.46, 13.14)% for GCN, (14.91, 17.53)% for BrainGNN, and (15.83, 16.49)% for BrainNetCNN on HCP; and it enable improvements of (14.66, 17.69)% for SAGE, (13.98, 14.47)% for GCN, (12.84, 16.72)% for BrainGNN, and (3.98, 4.72)% for BrainNetCNN on ID1000. We observe that for vanilla versions of the relatively simpler GNN models perform poorly against the more complex BrainNetCNN model. However, GraphCorr-augmented versions of these GNN models start outperforming the augmented BrainNetCNN. Thus, our results suggest that the feature extraction capabilities of vanilla GNN models might be suboptimal in comparison to CNN-based architectures, albeit a powerful feature extractor on the input side can mitigate this deficit in favor of GNN models.

### Explanatory analysis

To assess the influence of GraphCorr on interpretability, an explanatory analysis was conducted separately on the trained vanilla and augmented downstream models. For this analysis, the HCP dataset and the Schaefer atlas were selected that have been broadly studied in the literature for intrinsic brain networks during resting state [74]. First, network saliency scores obtained in each hemisphere were compared between vanilla and augmented versions of SAGE, which generally maintains the highest performance after GraphCorr augmentation. Literature reports that BOLD signals across the sensorimotor network (SMN), the default mode network (DMN) and the visual network bilaterally across the two hemispheres carry discriminative information on subject gender [75; 39]. Accordingly, we reasoned that a successful downstream classification model for gender detection should focus on these networks. Vanilla SAGE shows somewhat heterogeneous results with significant salience in the DMN in the left hemipshere (LH); the attention network, the SMN, and the visual network in the right hemisphere (RH); but unexpectedly it also yields strong salience in the limbic network in both hemispheres (\(p<0.05\), Wilcoxon signed-rank test) although it is not considered to carry information on gender. In contrast, augmented SAGE shows significant salience across the SMN, the DMN and the attention network in the RH; the visual network in both hemispheres (\(p<0.05\)), without any salience in the limbic network. These results imply that GraphCorr helps improve interpretability of the downstream classification model by allowing it to focus on brain regions that carry task-relevant information.

Significant saliency in a network indicates that BOLD signals across that network carry information about subject gender. Yet, it does not explain whether an increase or decrease in BOLD signals is evoked for individual genders. To address this question, a logistic regression analysis was conducted on BOLD signals extracted from the most important time window determined according to ROI saliency scores. Specifically, a logistic regression model was fit to detect subject gender given important BOLD signals. Fig. 4 illustrates the ROI weights in the logistic model, where a positive weight indicates that elevated BOLD signals in the ROI are associated with female subjects, and a negative weight indicates that elevated BOLD signals in the ROI are associated with male subjects. Accordingly, ROI weights were inspected for the logistic regression analyses based on the augmented SAGE model. In females, elevated BOLD signals are identified in RH parietal DMN areas including posterior cingulate cortex (PCC), RH prefrontal DMN areas, LH prefrontal control areas and LH-RH SMN areas. In males, elevated BOLD signals are identified in LH prefrontal DMN areas, LH-RH dorsal attention areas (DAN), LH parietal and prefrontal control areas, and RH parietal and extrastriate visual areas. These findings are consistent with evidence that females have relatively higher activations across DMN areas including PCC, and that males have relatively higher activations in visual and attentional areas [76; 77]. Our results are also consistent with recent studies suggesting that SMN and prefrontal regions show discriminate activation patterns across the two genders [39].

### Ablation studies

Ablation studies were performed to assess the contribution of the individual design elements in GraphCorr to model performance. These analyses were conducted based on the SAGE model using the HCP dataset and the Schaefer atlas, i.e., the

\begin{table}
\begin{tabular}{l l l l l l} \hline \hline \multirow{2}{*}{Model} & \multicolumn{3}{c}{HCP} & \multicolumn{3}{c}{ID1000} \\ \cline{3-5}  & \multicolumn{1}{c}{Acc (\%)} & \multicolumn{1}{c}{ROC (\%)} & \multicolumn{1}{c}{Acc (\%)} & \multicolumn{1}{c}{ROC (\%)} \\ \hline \multirow{2}{*}{SAGE} & Vanilla & 68.26 \(\pm\) 3.31 & 75.65 \(\pm\) 1.29 & 62.84 \(\pm\) 1.81 & 67.23 \(\pm\) 2.71 \\  & Augmented & **85.45**\(\pm\) 3.57 & **91.19**\(\pm\) 2.63 & **77.50**\(\pm\) 3.68 & **84.92**\(\pm\) 1.92 \\ \hline \multirow{2}{*}{GCN} & Vanilla & 69.90 \(\pm\) 1.35 & 75.85 \(\pm\) 1.06 & 65.45 \(\pm\) 1.36 & 70.96 \(\pm\) 1.23 \\  & Augmented & **84.36**\(\pm\) 3.17 & **88.99**\(\pm\) 2.83 & **79.43**\(\pm\) 3.45 & **85.43**\(\pm\) 2.49 \\ \hline \multirow{2}{*}{BrainGNN} & Vanilla & 65.69 \(\pm\) 3.00 & 71.79 \(\pm\) 2.97 & 62.27 \(\pm\) 3.44 & 66.42 \(\pm\) 4.24 \\  & Augmented & **80.60**\(\pm\) 2.67 & **89.32**\(\pm\) 2.67 & **75.11**\(\pm\) 0.56 & **83.14**\(\pm\) 1.41 \\ \hline \multirow{2}{*}{BrainNetCNN} & Vanilla & 68.16 \(\pm\) 3.53 & 74.76 \(\pm\) 2.08 & 75.00 \(\pm\) 2.19 & 81.44 \(\pm\) 3.27 \\  & Augmented & **83.99**\(\pm\) 2.92 & **91.25**\(\pm\) 2.65 & **78.98**\(\pm\) 1.83 & **86.16**\(\pm\) 0.84 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Performance of downstream models on the HCP and ID1000 datasets with the AAL atlas. Results are listed as mean\(\pm\)std across test folds for vanilla and GraphCorr-augmented versions. Boldface indicates the better performing version of each model.
setting that yields the highest overall performance for gender detection. First, we assessed contributions of the node embedder module, lag filter module, and time windowing in GraphCorr. To ablate the node embedder module, node embeddings prior to message passing were initialized with the unlearned time-windowed FC matrix derived via conventional correlation measures on BOLD signals. To ablated the lag filter module, a single filter at zero lag was used within the module to consider only instantaneous correlations. To ablated time windowing, the entire fMRI time series was provided to GraphCorr with a single window of size equal to scan duration. Table 3 lists performance metrics for ablated variants of GraphCorr. We find that the node embedder module, the lag filter module and time windowing enable (accuracy, ROC)% improvements of \((5.31,3.2)\)%, \((0.65,0.17)\)%, and \((8.41,5.61)\)%, respectively.

Next, we assessed the benefits of GraphCorr over alternative plug-in approaches to improve the temporal sensitivity of downstream model. In particular, we considered providing the downstream model pre-computed dynamic FC features across time windows via conventional correlation measures [66], an RNN model based on LSTM layers [78], and an RNN model based on GRU layers [79]. The feature dimensionality at the output of all plug-in models were identical to that for GraphCorr. Table 4 lists performance metrics for different plug-in methods. GraphCorr outperforms all other plug-in methods, with \((5.77,2.15)\)% higher performance than the top-contending GRU method.

## 6 Discussion

Here we reported a novel plug-in GNN method, GraphCorr, to improve the performance of downstream classification models in fMRI analysis by capturing dynamic, lagged FC features of BOLD signals. Demonstrations were provided on two large-scale resting-state fMRI datasets, where substantially improved performance was achieved following model augmentation with GraphCorr. The proposed method can be trivially combined with classification models to detect other categorical variables related to cognitive task or disease [18; 17]. Alternatively, it can be employed as a plug-in to downstream regression models to boost sensitivity in predicting continuous variables related to stimulus or task features [1].

A mainstream approach in neuroimaging studies rests on prediction of experimental variables typically related to stimulus or task from BOLD signals [41; 51]. Here we adopted this approach to build decoding models that predict subject gender from resting-state fMRI scans. An alternative procedure to examine cortical function rests on encoding models that instead predict BOLD signals from experimental variables [45; 80; 81]. It may be possible to adopt GraphCorr to improve sensitivity of such downstream encoding models. In this case, GraphCorr would receive as input the time course of experimental variables during an fMRI scan. In turn, it would learn dynamic, lagged correlations among experimental variables to better account for their distribution. Learned correlations might help improve performance of downstream regression models that aim to predict measured BOLD signals. Future work is warranted to investigate the potential of GraphCorr in building encoding models for fMRI.

In conjunction with downstream models, GraphCorr was directly trained end-to-end on the HCP or ID1000 datasets that contained data from several hundred subjects. While the lag filter module has low complexity, the node embedder module uses a transformer encoder with a relatively large number of parameters. To improve learning on limited datasets, transfer learning can be performed where the encoder is initialized with pre-trained weights [82]. Data augmentation procedures that can produce a large variety of realistic samples from a learned distribution might further facilitate learning [83; 84]. GraphCorr forms an initial graph where edges are retained in a single-hop neighborhood based on static FC values between corresponding

\begin{table}
\begin{tabular}{c c c} \hline \hline Plug-in & Accuracy (\%) & ROC (\%) \\ \hline Dynamic FC & \(81.06\pm 2.94\) & \(89.04\pm 1.53\) \\ LSTM & \(83.26\pm 2.05\) & \(90.29\pm 1.40\) \\ GRU & \(83.80\pm 2.29\) & \(92.12\pm 2.01\) \\ GraphCorr & \(\textbf{89.57}\pm 0.68\) & \(\textbf{94.27}\pm 1.98\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Performance of competing plug-in methods on the HCP dataset with the Schaefer atlas. Results are listed as mean\(\pm\)std across test folds for the downstream SAGE model. Boldface indicates top-performing plug-in.

Figure 4: Salient ROIs for gender detection assessed via the logistic regression analysis. Results are shown for the GraphCorr-augmented SAGE model on the HCP dataset with the Schaefer atlas. ROIs with the top 2% saliency scores are marked. Red color indicates ROIs whose BOLD signals are elevated in female subjects, whereas blue color indicates ROIs whose BOLD signals are elevated in male subjects.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Node & Lag Filter & Windowing & Accuracy (\%) & ROC (\%) \\ Embedder & & & & \\ \hline ✗ & ✗ & ✗ & \(75.20\pm 2.83\) & \(85.29\pm 1.67\) \\ ✓ & ✗ & ✗ & \(80.51\pm 1.92\) & \(88.49\pm 1.54\) \\ ✓ & ✓ & ✗ & \(81.16\pm 1.69\) & \(88.66\pm 2.07\) \\ ✓ & ✓ & ✓ & \(\textbf{89.57}\pm 0.68\) & \(\textbf{94.27}\pm 1.98\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Performance for ablated variants of GraphCorr on the HCP dataset with the Schaefer atlas. Results are listed as mean\(\pm\)std across test folds for the downstream SAGE model. Boldface indicates top-performing variant.
nodes. This structure is kept fixed during subsequent training procedures. To improve performance, an adaptive structure can be used instead where the edge weights are taken as learnable parameters.

Here each individual subject's fMRI scans were aligned to an anatomical template, and brain regions were then defined with guidance from a brain atlas. The mean BOLD signals in each ROI were then processed in downstream models. Benefits of this approach include computational efficiency due to relatively lower model complexity, and consistency in region definitions across subjects [85]. Meanwhile, information losses naturally occur during registration of individual-subject fMRI data onto a standardized template. To alleviate these losses, ROI definitions in the template space could instead be backprojected onto the brain spaces of individual subjects. This way ROI definitions can be performed while leaving fMRI data in its original space [51].

## 7 Conclusion

In this study, we introduced a novel plug-in graph neural network to improve the performance of downstream models for fMRI classification. The proposed GraphCorr method employs node embedder and lag filter modules to sensitively extract dynamic and lagged functional connectivity features from whole-brain fMRI time series. As such, it transforms raw BOLD signals into a graph representation where neighboring nodes are taken as brain regions with correlated signals and node features are extracted via message passing on connectivity features from the two modules. This procedure restores the fine-grained temporal information that can otherwise be diminished in conventional functional connectivity features. As augmenting downstream classification models with GraphCorr significantly improves their performance and interpretability, GraphCorr holds great promise for analysis of fMRI time series.

## Acknowledgments

This study was supported in part by a TUBITAK BIDEB scholarship awarded to H.A. Bedel, by TUBA GEBIP 2015 fellowship, and TUBITAK 121N029 grant awarded to T. Cukur.
