Learning Dynamic Graph Representation of Brain Connectome with Spatio-Temporal Attention
Year: 2022

Abstract:
Functional connectivity (FC) between regions of the brain can be assessed by the degree of temporal correlation measured with functional neuroimaging modalities. Based on the fact that these connectivities build a network, graph-based approaches for analyzing the brain connectome have provided insights into the functions of the human brain. The development of graph neural networks (GNNs) capable of learning representation from graph structured data has led to increased interest in learning the graph representation of the brain connectome. Although recent attempts to apply GNN to the FC network have shown promising results, there is still a common limitation that they usually do not incorporate the dynamic characteristics of the FC network which fluctuates over time. In addition, a few studies that have attempted to use dynamic FC as an input for the GNN reported a reduction in performance compared to static FC methods, and did not provide temporal explainability. Here, we propose STAGIN, a method for learning dynamic graph representation of the brain connectome with spatio-temporal attention. Specifically, a temporal sequence of brain graphs is input to the STAGIN to obtain the dynamic graph representation, while novel READOUT functions and the Transformer encoder provide spatial and temporal explainability with attention, respectively. Experiments on the HCP-Rest and the HCP-Task datasets demonstrate exceptional performance of our proposed method. Analysis of the spatio-temporal attention also provide concurrent interpretation with the neuroscientific knowledge, which further validates our method. Code is available at https://github.com/egyptdj/stagin

heading: 1 Introduction
text: Neuroimaging modalities provide measurements of brain activity by capturing the signals of neural activity. Functional magnetic resonance imaging (fMRI) is a non-invasive imaging method that measures the blood-oxygen level dependence (BOLD) in order to estimate the neural activity of the whole brain over time [20]. Functional connectivity (FC) is defined as the degree of temporal
∗This work was conducted while the first author was at KAIST
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
correlation between regions of the brain. Based on the fact that these connectivities form networks that change over time, graph-based network analysis of brain connectome has been one of the key approaches to understanding how the brain works [7, 4, 43].
Graph neural networks (GNNs) are a type of deep neural networks that have recently been successful in learning the representation of graph-structured data [50]. The graph-structured nature of the brain has led to an increased interest in learning the reperesentation of the brain FC network with the GNNs. Learning the representation of the brain connectome can be linked to decoding trait or state from human brain signal measurements. Accordingly, the current trend in studies attempting to apply GNN to the brain connectome is to input the FC graph from either resting-state [25, 26, 2, 33, 23, 48, 49] or task fMRI data [29, 30, 32] and predict a particular phenotype of the subjects, such as gender [26, 2, 23, 22] or presence of a specific disease [26, 33, 29, 30, 32, 22]. While these studies have shown potential strengths and opportunities for learning the network representation of the brain, they also suggest limitations of current GNN-based methods.
One of the most common limitations with previous GNN-based FC network analysis methods is that most of them fail to take advantage of the dynamic properties of the FC network, which fluctuates over time. Incorporating the dynamic features of the FC network into the neuroimaging analysis has been an important direction in the field of functional neuroimaging [21, 37]. A work by [15] tried to address this issue by using the Spatial Temporal Graph Convolutional Network (ST-GCN) [54] model to incorporate dynamic features of the FC network. However, [15] reported lower accuracy than other non-dynamic GNN-based FC methods [23, 2] in the gender classification experiment, leaving a question about the effectiveness of the dynamic FC method. In addition, another limitation of the method is that no temporal explainability is provided from the model. This is a major drawback considering that the goal of applying GNNs to functional neuroimaging methods is not only to achieve high classification accuracy, but also to uncover the functional basis of the brain [23, 32]. Another recent work by [3], using GraphNets [5] and DiffPool [58] for the dynamic FC analysis, also suffers from the same limitations in terms of poor classification accuracy and lack of temporal explainability.
Here, we propose Spatio-Temporal Attention Graph Isomorphism Network (STAGIN) for learning the dynamic graph representation of the brain connectome with spatio-temporal attention. The proposed method exploits the temporal features of the dynamic FC network graphs to improve the classification accuracy of the model. In particular, we address the issue that the node features of the input dynamic graph should contain temporal information and concatenate encoded timestamp with the node features (Section 4.1). In addition, the proposed method includes novel attention-based READOUT modules (Section 4.2) and the Transformer encoder [46] (Section 4.3) in order to further improve the classification performance and provide spatial-temporal explainability at the same time. STAGIN achieves state-of-the-art performance with the Human Connectome Project (HCP) dataset [45] in gender classification for resting-state fMRI and task decoding for task fMRI. We inherit k-means clustering analysis of the resting-state dynamic FC [1] and general linear model (GLM) statistical mapping of task fMRI [14] for interpreting the spatio-temporal attention learned from STAGIN, which are widely accepted analysis methods for the fMRI data. The interpretation of the learned spatio-temporal attention replicates neuroscientific findings from previous large-scale fMRI studies in both resting-state and task fMRI, which further validates our proposed method.
Our work holds potential societal impact in that brain decoding methods can be linked to finding neural biomarkers of important phenotypes or diseases. However, potential negative impact related to privacy concerns that arise from abuse or misuse of accurate decoding methods should also be noted. Although our method is yet behind the decoding capability that can be abused or misused, our research cannot still be free from these ethical considerations.

heading: 2 Related works

heading: 2.1 Graph Neural Network on Dynamic Graphs
text: Many networks that arise around us are inherently dynamic, with changes in the existence of nodes and edges over time. Learning the representation of dynamic graphs has piqued the interest of researchers and has led to development of methods that can embed dynamic graphs using their time information [35]. Methods that incorporate attention for learning the representation of dynamic graphs have also been proposed [51, 40]. However, it is not easy to apply these techniques directly to the dynamic brain graphs because of the different inherent properties of the dynamic brain graphs that
do not include any addition or deletion of nodes and are sampled uniformly over time. Nonetheless, our work is inspired by these earlier studies, particularly for the encoding of temporal information and their concatenation to the node features, as proposed in Section 4.1 [51, 40].

heading: 2.2 Attention in Graph Neural Networks
text: Bringing attention to the GNNs is a topic that is being actively studied in the field of geometric deep learning [27]. One of the most successful uses of attention is to compute the attention at edges of the graph and scale the importance of the links when the features of the neighborhood node are aggregated [47, 6], often providing performance gain in learning the representation of input graphs. Another stream of applying attention to the GNNs comes with the motivation to define a pooling function on the graph domain. Since it is not straightforward to decide on what basis the coarsening should be carried out for graph structured data, works such as [16, 28, 38] have addressed this problem by selecting the nodes with top scores computed from projecting the node feature vectors into a learnable parameter vector, or from a GNN layer aggregated local graph features. Although the motivation may have been different, these graph pooling methods are closely related to the spatial attention modules that we propose in Section 4.2 in that they exploit learned relative scores across the vertices of the graph. While some works have already been aware that the appropriate use of node-wise attention can improve performance of downstream tasks [55, 11], we note that previous methods tend to score attention based on randomly initialized parameters or local graph structures which may be suboptimal for graph classification tasks that require taking the whole graph feature into account.

heading: 3 Theory

heading: 3.1 Problem definition
text: The goal of our study is to train a neural network
f : Gdyn → hGdyn ,
where Gdyn = (G(1), ..., G(T )) is the sequence of brain graphs with T timepoints and hGdyn ∈ RD is the vector representation of the dynamic graph G(t) with length D. The graph G(t) = (V (t), E(t)) at time t is a pair of vertex set V (t) = {x1(t), ...,xN (t)} of N nodes and edge set E(t) ={ {xi(t),xj(t)} | j ∈ N (i), i ∈ {1, ...N} } where N (i) denotes the neighborhood of the vertex
i. If f learns to extract a disentangled representation of the dynamic brain graph Gdyn, then the classification of a certain phenotypic characteristic (e.g. gender) from hGdyn can be performed with a linear mapping as a downstream task. Another important consideration in this work is to ensure the explainability of the model f , being able to inform us which part of the brain at which timepoint was considered important when extracting the meaningful representation hGdyn . Specifically, we formulate f = q ◦ g as a composition of the GNN g and the Transformer encoder q, where g outputs the set of graph representations hG(t) from each timepoint and q exploits self-attention to integrate hG(t) into the final representation hGdyn :
g : Gdyn → (hG(1), ...,hG(T )), (1) q : (hG(1), ...,hG(T ))→ hGdyn . (2)
We will omit timepoint notation (t) for brevity, whenever it is not of contextual importance.

heading: 3.2 Graph Isomorphism Network
text: The GNNs are generally composed of functions that (i) integrate the node features from its neighbors, and (ii) embed the integrated information with a nonlinear transformation to obtain the next layer node features. These functions are called AGGREGATE, and COMBINE functions, respectively, and the choice of these functions define many variants of the GNN,
a(k)v = AGGREGATE (k) ({ h(k−1)u : u ∈ N (v) }) , (3)
h(k)v = COMBINE (k) ( h(k−1)v ,a (k) v ) , (4)
where h(k)v denotes the feature vector of node v at layer k and h (0) v := xv .
The Graph Isomorphism Network (GIN) is a variant of the GNN suitable for graph classification tasks, which is known to be as powerful as the WL-test under certain assumptions of injectivity [52]. The GIN typically defines sum as the AGGREGATE and a multi-layer perceptron (MLP) with two layers as the COMBINE updating the node representation h(k)v at layer k [52] by :
h(k)v = MLP (k) ( (1 + (k)) · h(k−1)v + ∑
u∈N (v)
h(k−1)u
) , (5)
where is a learnable parameter initialized with zero. Equation (5) can be easily reformulated into the matrix form [23] by:
H(k) = σ ( ( (k) · I +A)H(k−1)W (k) ) , (6)
where H(k) = [ h
(k) 1 , · · · ,h (k) N
] ∈ RD×N
is the stack of node feature vectors, I is the identity matrix,A is the adjacency matrix between the node features,W is the network weights of the MLP, and σ is the nonlinearity function.
The READOUT function takes the updated node features h(k)v to compute the representation of the whole graph:
h (k) G = READOUT ( {h(k)v | v ∈ G} ) . (7)
In general, the READOUT function is defined simply as computing the sum or average of the input node features. This is equivalent to multiplication with the length N pooling vectors φ>sum = [1, ..., 1] or φ>mean = [1/N, ..., 1/N ] for the matrix form:
h (k) G = H (k)φmean. (8)

heading: 3.3 Encoder-decoder understanding of GNNs
text: Although formulating the GIN (5) as a combination of AGGREGATE and COMBINE function might not suggest its close relationship with convolutional neural networks (CNNs) at first glance, previous works by [23, 8] show that the matrix formulation of the GIN operation (6) can be thought of a CNN layer with shift operation of the convolution as the adjacency matrixA. We extend the understanding of encoder-decoder CNN as a framelet expansion [57, 56] to the GIN to formulate node feature vectorsH(k) at layer k with respect to the input node feature xi as:
Vec ( H(k) ) = Σ(k)E(k)> · · ·Σ(1)E(1)>x, where x := Vec ([x1, · · · ,xN ]) (9)
where Vec(·) denotes the vectorization operation, and the k-th layer encoder matrix E(k) is defined as E(k) = W (k) ⊗ ( (k) · I +AT ) where⊗ refers to the Kronecker product, and Σ(k) is the diagonal matrix with values 1 or 0 depending on the activation pattern of the nonlinearity. Now, φmean of equation (8) can be thought as the decoder at the k-th layer which yields the whole graph feature vector from the encoded node feature vectors. Proposition 1. The READOUT function φmean in (8) generates a decoder with fixed constant bases.
Proof. From the READOUT function (8), we have
h (k) G =Vec
( h
(k) G
) = Vec ( H(k)φmean ) = ( φTmean ⊗ I ) Vec ( H(k) ) = ( φTmean ⊗ I ) Σ(k)E(k)> · · ·Σ(1)E(1)>x
Now let bi and b̃i denote the i-th column of the encoder matrix E(1)Σ(1) · · ·E(k)Σ(k) and the decoder matrix ( φTmean ⊗ I ) , respectively. Then, it is straight to obtain the following representation:
h (k) G = ∑ i 〈bi,x〉b̃i
Therefore, we can see that although the encoder basis bi is a function of x, the decoder basis b̃i is a constant.
We address the issue that the decoder being a constant function can restrict the expressivity of the neural network, and explore adaptive READOUT functions with attention in Section 4.2.

heading: 4 STAGIN: Spatio-Temporal Attention Graph Isomorphism Network
text: In this section, we discuss the details of our main contribution. Specifically, we propose STAGIN with two novel attention-based READOUT modules for learning the dynamic graph representation of the brain connectome (Figure 1).

heading: 4.1 Dynamic graph definition
text: The sequence of input dynamic FC graphs is constructed from 4D fMRI data with 3D voxels across time. The ROI-timeseries matrix P ∈ RN×Tmax is extracted by taking the mean values within a predefined 3D atlas which consists of N ROIs at each timepoint. Values of each ROI are standardized across time. Constructing dynamic FC matrix follows the sliding-window approach, where the temporal window of length Γ is shifted across time with stride S to generate T = bTmax − Γ/Sc windowed matrices P̄ (t) ∈ RN×Γ (Figure 2 (a)). The FC at time t is defined as the correlation coefficient matrixR(t) of the windowed timeseries between p̄i(t) and p̄j(t):
Rij(t) = Cov(p̄i(t), p̄j(t))
σp̄i(t)σp̄j (t) ∈ RN×N ,
where the subscript i and j are the row and column indices of P̄ (t), Cov denotes the cross covariance, and σp denotes the standard deviation of p. The final binary adjacency matrix A(t) ∈ {0, 1}N×N is obtained from the FC matrixR(t) by thresholding the top 30-percentile values of the correlation matrix as connected, and otherwise unconnected following [23]. Other thresholds for binarizing the correlation matrix are also experimented and the results are provided in the Appendix Section C.2.
Unlike the adjacency matrix A(t), conventional definition of node feature vectors xv(t) at node index v as coordinates [29], mean-activation [29, 15], or one-hot encoding [23], do not change over t, disregarding any temporal variation. To address this issue, we concatenate encoded timestamp η(t) ∈ RD to the spatial one-hot encoding ev , followed by linear mapping with a learnable parameter matrixW ∈ RD×(N+D) to define the input node feature,
xv(t) = W [ev||η(t)]. (10) Here, the learnable timestamp encoder η is a Gated Recurrent Unit (GRU) [9] which takes ROItimeseries upto the endpoint of the sliding-window as the input. Both the vertex set V (t) and the
edge set E(t) of graph G(t) now incorporates temporal information at time t. See Figure 2 for an illustration of the dynamic graph definition.

heading: 4.2 Spatial attention with attention-based READOUT
text: As suggested from Proposition 1, conventional READOUT function of GNN can be thought of as a fixed decoder that decodes whole-graph feature from the node features with no learnable parameters. We address this issue by incorporating attention to the READOUT function, which the attention here refers to the scaling coefficient across the nodes learned by the model. Specifically, the spatial attention vector zspace(t) ∈ [0, 1]N is computed by taking theH as a prior:
zspace = s(H), (11)
h̃G = Hzspace, (12)
where s : RD×N → [0, 1]N is the attention function and h̃G denotes spatially attended graph representationhG. We propose two types of attention function s(·) for the attention-based READOUT, named Graph-Attention READOUT (GARO) and Squeeze-Excitation READOUT (SERO) inspired by the attention mechanisms of [46] and [19], respectively.

heading: 4.2.1 GARO: Graph-Attention READOUT
text: The GARO follows key-query embedding based attention of the Transformer [46]. However, the key embedding is computed from the matrix of node featuresH , while the query embedding is computed from the vector of unattended graph representationHφmean:
K = WkeyH,
q = WqueryHφmean, zspace = sigmoid (q>K√
D
) , (13)
whereWkey ∈ RD×D,Wquery ∈ RD×D are learnable key-query parameter matrices,K ∈ RD×N is the embedded key matrix, and q ∈ RD is the embedded query vector.

heading: 4.2.2 SERO: Squeeze-Excitation READOUT
text: The SERO follows MLP based attention of the Squeeze-and-Excitation Networks [19]. However, attention from the squeezed graph representation does not scale the channel dimension, but the node dimension in SERO:
zspace = sigmoid ( W2σ(W1Hφmean) ) , (14)
where σ is the nonlinearity function and W1 ∈ RD×D, W2 ∈ RN×D are learnable parameter matrices. This type of spatial dimension squeeze-excitation module has been shown to improve
performance of the CNN models [41], but was not easily applicable to general graphs which may vary in number of nodes for each graph. We exploit the fact that the brain graphs have fixed number of nodes N across participants based on the chosen atlas.

heading: 4.2.3 Orthogonal regularization
text: If we take a closer look at (8) and (12), computation of graph feature vector hG from the node feature matrix H can also be viewed as reconstructing signal hG from the basis frames H with vectors φmean and zspace, respectively. While zspace provides further expressivity of the model with adaptive coefficients when compared to φmean, we find it desirable to encourage the orthogonality of H as elaborated in the Appendix Section A. The orthogonal regularization Lortho is defined as:
Lortho = ∥∥1/m ·H>H − I∥∥
2 , (15)
where m = max(H>H). The scaling term 1/m ensures the columns of the matrix H become orthogonal to each other with the same length, while not restricting the specific length that the column vectors should follow.

heading: 4.3 Temporal attention with Transformer encoder
text: For attention across time, we employ a single-headed Transformer encoder [46] upon the sequence of graph features (h̃G(1), ..., h̃G(T )). The temporal attention can be measured by the self-attention weights Ztime ∈ [0, 1]T×T after the softmax function of the Transformer encoder. Per-layer dynamic graph representation h(k)Gdyn is computed by summing the temporally attended feature output from the Transformer encoder across time at each layers, where the final representation:
hGdyn = concatenate({h (k) Gdyn | k ∈ {1, ...,K}}) (16)
is the concatenation of dynamic graph representation of all K layers following [53].

heading: 5 Experiment

heading: 5.1 Dataset
text: Publicly available2 fMRI data from the HCP S1200 release [45] was used for our experiments. The data was collected from voluntary participants with informed consent and was fully anonymized. We constructed two datasets, the HCP-Rest and the HCP-Task, depending on whether the subject was resting or performing specific tasks during the acquisition of the image. The HCP-Rest dataset consisted of pre-processed and ICA denoised resting-state fMRI data [17], which the subjects were instructed to rest for 15 minutes during the data acquisition. We used first run data of the four sessions, and excluded data with short acquisition time with Tmax < 1200. There were 1093 images finally included in the dataset, which consisted of 594 female and 499 male subjects. The gender of each subject served as the labels of the HCP-Rest dataset letting the number of classes C = 2. The HCP-Task consisted of pre-processed task fMRI data [17], which the subjects were instructed to perform specific tasks during data acquisition. For example in the "Motor" task fMRI, participants were told to perform one of the subtasks during the acquisition to make motor movements on one’s left hand, left foot, right hand, right foot, or tongue. There were seven types of tasks including working memory, social, relational, motor, language, gambling, and emotion. After excluding the fMRI data with short acquisition time, there were 7450 images included in the dataset. The task type during the data acquisition served as the labels of the HCP-Task dataset, letting C = 7. A more detailed description of the experiment datasets with a note on the twin subjects of HCP can be found in the Appendix Section B.

heading: 5.2 Experimental settings
text: Experiments were performed on a workstation with two NVIDIA GeForce GTX 1080 Ti GPUs. The STAGIN model f is trained end-to-end in a supervised manner with the loss L = Lxent + λ · Lortho
2https://db.humanconnectome.org
where Lxent is the cross entropy loss and λ is the scaling coefficient of the orthogonal regularization. We set the number of layersK = 4, embedding dimensionD = 128, window length Γ = 50, window stride S = 3, and regularization coefficient λ = 1.0×10−5. The window length and stride correspond to capturing the FC within 36 seconds every 2.16 seconds, which follows the standard setting of the sliding-window dFC analyses [59, 37]. Dropout rate 0.5 is applied to the final dynamic graph representation hGdyn , and rate 0.1 is applied to the attention vectors zspace and ztime during training. For nonlinearity σ in (6) and (14), GELU [18] is used instead of ReLU with batch normalization before each σ. One-cycle learning rate policy is employed, which the learning rate is gradually increased from 0.0005 to 0.001 during the early 20% of the training, and gradually decreased to 5.0× 10−7 afterwise. Thirty training epochs were run for the HCP-Rest dataset with minibatch size 3, while ten epochs were run with minibatch size 16 for the HCP-Task dataset. We performed 5-fold stratified cross-validation of the dynamic graphs from the dataset, and report mean and standard deviation across the folds. To extract the ROI-timeseries, the Schaefer atlas [42] with 400 regions (N = 400) labelled with 7 intrinsic connectivity networks (ICNs) was used. The time dimension of ROI-timeseries matrix P was randomly sliced with a fixed length (600 for HCP-Rest, 150 for HCP-Task) at each steps during training for (i) relieving computational overload, (ii) stochastic augmentation of the training dataset, (iii) mitigating unwanted memorization of the specific timing of subtask onset, and (iv) matching the number of timepoints T across different task labels for the HCP-Task dataset. Unsliced full matrix P was used for inference at test time. The end-to-end inference from the construction of the dynamic graph to the acquisition of the final prediction required 1.68 seconds per sample with given experimental settings.

heading: 5.3 HCP-Rest: Gender classification
text: We first validate our proposed method by gender classification on the HCP-Rest dataset. The two proposed methods, named STAGIN-GARO and STAGIN-SERO based on the type of the spatial attention module, resulted in 87.01% and 88.20% mean accuracy on the 5-fold cross validation, respectively (Table 1). The mean area under receiver operator characteristic curve (AUROC) were 0.9151 and 0.9296. Classification performance of STAGIN is compared with other GNN methods for reprensentation learning of dynamic/static FC network, including ST-GCN [15], MS-G3D [10], GIN [23], GCN [24], GraphSAGE [31], and ChebGCN [2]. We used the code by the authors of [15]3 and [10]4 but modified the cross validation scheme to avoid early stopping based on the test dataset for fair comparison. It can be seen from Table 1 that our proposed method outperforms other GNN based methods. The results of the ablation study are shown in Table 3 in the Appendix.
We use STAGIN-SERO, which showed the best accuracy, for analyzing temporal and spatial attention of the dynamic FC networks. We define the temporal attention vector z(k)time ∈ [0, 1]T at layer k as the average of row elements in the self-attention weight matrix z(k)time[j] = 1 T ∑T i=1 Zij where
3https://github.com/sgadgil6/cnslab_fmri 4https://github.com/metrics-lab/ST-fMRI
z (k) time[j] and Zij are j-th element of z (k) time and (i, j)-th element of Z (k) time for the resting-state data, respectively. To employ k-means clustering to the resting-state dynamic FC analysis [1], we first define a set of attended timepoints T̃ = {t | ztime[t] > α · σztime} where α is the cutoff coefficient, and σ
z (k) time denotes the standard deviation of z(k)time. Defining the threshold based on standard deviation inherits the practice of the point-process analysis for dynamic FC, so we set α = 1.0 following [44]. Pattern of the FC matrices at attended timepoints AT̃ = {A(t) | t ∈ T̃} for each subject can now be analyzed with the k-means clustering. Specifically, we fit 7 template cluster centroids from the dynamic FC matricesA(t) over all subjects, and assign elements of AT̃ into one of the 7 template clusters. The ratio of each clusters from AT̃ with respect toA can then be analyzed with the subset of AT̃ including only the female or male subjects.
Evidences from large scale studies suggest that female subjects show hyperconnectivity of the DMN [34, 39] and hypoconnectivity of the SMN when compared to male subjects [39, 13]. We accordingly hypothesized that the FC at attended timepoints will show higher values for the DMN and lower values for the SMN in female participants. Figure 3 demonstrates that the clusters mainly attended by female participants show a trend of hyperconnectivity of the DMN and hypoconnectivity of the SMN. This can be interpreted to mean that the STAGIN is properly trained to take the dynamic state of the FC networks into account for predicting the phenotype of the subject.
The spatial attention across regions of the brain is analyzed with the z(k)space averaged across time z̃ (k) space := 1 T ∑T t=1 z (k) space(t). The regions with top 5 percentile attention values of z̃ (k) space are plotted with respect to the seven ICNs in Figure 9 in the Appendix. It can be seen that the majority of the top attended regions are from the SMN, which further suggests gender difference of resting-state FC within the SMN. A notable limitation here is that the threshold for determining the top attended region is heuristically set. Statistically determining the spatially attended regions from the resting-state data would further provide validity of the method, which is left as a future work.

heading: 5.4 HCP-Task: Task decoding
text: Task decoding refers to classifying which of the seven tasks the subject was performing during the acquisition of the brain fMRI. The STAGIN-GARO and STAGIN-SERO showed 99.02% and 99.19% mean accuracy for the task decoding experiment, respectively (Table 1). It can be seen that the proposed methods outperform the previous state-of-the-art model BAND and BAnD++ [36], which applied self-attention of the Transformer encoder directly to 3D ResNet extracted representation vectors of the fMRI without considering the network property of the brain. To account for the possible statistical disadvantage of voxel-based feature extraction, we further implemented a new region-based BAnD (r-BAnD) by using GIN without attention-based READOUT instaed of the 3D ResNet. Accuracy of r-BAnD resulted in an accuracy of 98.90%, suggesting that our method shows superior performance even when the statistical disadvantages are matched. Experiment on
other models including ST-GCN [15], GIN [23], GCN [24], GraphSAGE [31], and ChebGCN [2] demonstrate exceptional performance of our proposed method for HCP-Task (Table 1). The fact that subtask timing information is completely lost may reflect the reason behind poor performance of static FC methods, which can be a critical disadvantage in task classification.
We interpret the result from the working memory task for spatio-temporal attention analysis, where the subtask consists of either performing an n-back memory task or rest. Our key expectation of the temporal attention analysis was that if STAGIN learns to accurately attend to temporal features of the dynamic FC graphs, then Ztime should represent which subtask the subject was upto. Surprisingly, it can be clearly seen that the Transformer encoder of STAGIN learns to attend to the timing of subtasks from Figure 4 (a), which demonstrates mean temporal attention Ztime across all subjects. Notice that no supervision is provided to the STAGIN model regarding the subtask timing during training.
To analyze the spatially attended regions zspace of STAGIN, we construct a GLM [14] to statistically evaluate how much each region is responsible for performing the subtasks. The parameter vectors βtask ∈ RN and βrest ∈ RN are estimated with the sequence of spatial attention vectors and the subtask timing design matrix M ∈ {0, 1}T×2 by solving the following with least-squares estimation:
[zspace(0), · · · , zspace(T )]> = M [βtask,βrest]> + ,
where denotes residual error. The contrast of the estimated parameters β̂task and β̂rest was set to c = [1,−1] so the rejection of null hypothesis indicates β̂task[i] > β̂rest[i] at the i-th ROI. Multiple comparisons of the N ROIs are family-wise error (FWE) corrected.
Figure 4 (b) shows the proportion of statistically significant regions within the 7 ICNs for each layers. Interestingly, the layer 1 and 2 share a similar trend that the regions from SMN, visual network (VN), and salience/ventral attention network (SVN) are dominant. In contrast, layer 3 and 4 suggest a dominance of the regions from DMN and cognitive control network (CCN). We denote the layer 1 and 2 as the low-order layers (LoL) and the layer 3 and 4 as the high-order layers (HoL). The dominance of SMN and VN at LoL can be understood as the low-level sensorimotor function for perceiving the task is being processed within the short-range 1- or 2-hop connection of the networks. On the other hand, the dominance of DMN and CCN at HoL reflects the high-level cognitive integration for executing and controlling the given task being processed within the long-range 3- or 4-hop connection of the networks. Considering that the SVN is a network for integrating the low-level sensorimotor networks and the high-level executive networks to provide dynamic balancing between the two functions, the significant regions of SVN being present at both LoL and HoL is not surprising. Temporal and spatial attention plot of other six tasks are further provided in the Appendix Section D.2.
Acknowledgments and Disclosure of Funding
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2021M3E5D9025019, NRF-2020R1A2B5B03001980). This work was also supported by Institute of Information & communications Technology Planning & Evaluation (IITP) grant funded by the Korea government(MSIT) (No.2019-0-00075, Artificial Intelligence Graduate School Program(KAIST)) and the KAIST Key Research Institute (Interdisciplinary Research Group) Project.

Review summaries:
The paper proposes an attention based GNN method, which it names STAGIN. The new method handles dynamic FC matrices for gender and task classification using fMRI data and learns to focus on important time points and brain regions based on attention. The method achieves better accuracy and AUC for classification compared to the relevant competing methods that STAGIN is compared against. The temporal attention matches with the task-relevant time points in the input data. 

This paper proposes a method, STAGIN, to model dynamic connectivity. The architecture combines graph isomorphic network (GIN) with spatial and temporal attention. STAGIN is applied to HCP data to classify gender as well as tasks with higher accuracy than state-of-the-art methods shown.

This paper proposes a GNN model (STAGIN) for learning from dynamic, spatio-temporal fMRI data. The model, based on GIN, takes a temporal sequence of brain graphs as input and includes a time encoding in the node features. The standard readout function is modified to incorporate spatial attention in two possible ways (scaled dot product attention and squeeze-excitation model). The series of graph features is then passed to a transformer-style self-attention encoder to incorporate temporal attention. Orthogonality of node feature vectors is also encouraged in the loss function. The method was tested on the public HPC datasets (~1100 subjects) to classify gender from resting state fMRI and decode task from task-based fMRI, and analysis of the temporal and spatial attention results was performed.