{'fraction_non_alphanumeric': 0.029999904550096882, 'fraction_numerical': 0.016159668597936373, 'mean_word_length': 4.782855881216537, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Forecasting Brain Activity Based on Models of Spatio-Temporal Brain Dynamics: A Comparison of Graph Neural Network Architectures

Comprehending the interplay between spatial and temporal characteristics of neural dynamics can contribute to our understanding of information processing in the human brain. Graph neural networks (GNNs) provide a new possibility to interpret graph structured signals like those observed in complex brain networks. In our study we compare different spatio-temporal GNN architectures and study their ability to model neural activity distributions obtained in functional MRI (fMRI) studies. We evaluate the performance of the GNN models on a variety of scenarios in MRI studies and also compare it to a VAR model, which is currently often used for directed functional connectivity analysis. We show that by learning localized functional interactions on the anatomical substrate, GNN based approaches are able to robustly scale to large network studies, even when available data are scarce. By including anatomical connectivity as the physical substrate for information propagation, such GNNs also provide a multi-modal perspective on directed connectivity analysis, offering a novel possibility to investigate the spatio-temporal dynamics in brain networks.


Introduction
Distinct concepts of brain connectivity can provide different but complementary aspects of information processing in the brain (Lang et al., 2012). On the one hand, imaging modalities like functional magnetic resonance imaging (fMRI) allow us to temporally resolve dynamic neural activity patterns in distinct spatial locations in the human brain. Different statistical approaches that describe the coherency of activity profiles in brain networks have been proposed based on the notion of functional connectivity (FC). On the other hand, diffusion tensor imaging (DTI) can provide insights into the structural and relatively static aspects of the brain. By reconstructing white matter tracks from DTI data, the anatomical or structural connectivity (SC) between different brain areas can be estimated. Also directed and potentially causal relationships between regions become of interest in fMRI and are studied with respect to directed functional or effective connectivity. The latter is most often inferred from Granger causality or dynamic causal modeling (Friston et al., 2013;Bielczyk et al., 2019).
Based on these concepts, spatial and temporal relationships between brain areas can be represented by graphical models, which have recently received increasing attention in the field of machine learning (Wu et al., 2020;Bronstein et al., 2017). So-called graph neural networks (GNNs) allow us to effectively process signals in the non-Euclidean geometry of graphs, providing also novel possibilities for applications in brain connectivity research (Rosenthal et al., 2018;Wein et al., 2021b;Arslan et al., 2018;Kim and Ye, 2020). Given a decomposition of the brain into specified areas, their spatio-temporal neural activity patterns can be interpreted as graph structured signal distributions. Nodes in brain networks can be associated with variables like the temporal neuronal activity of neuron pools, while edges in such networks reflect the strength of interactions between neural populations (Bullmore and Bassett, 2011). As proposed in our recent study (Wein et al., 2021b), these complex signals exhibited in non-Euclidean geometries can be processed with a variant of GNN denoted as spatiotemporal graph neural network (STGNN). Such STGNNs can allow us to simultaneously models spatial and temporal dependencies in such graph structured signals and thereby provide a new possibility to combine DTI with fMRI data. Activity propagation based approaches made already various interesting contributions to brain connectivity research, and could for example explain how resting-state FC (Cole et al., 2016) or SC (Yan et al., 2021) are related to cognitive task activations observed in task-based fMRI. Moreover, they could give us insights in what way the hierarchical organization of the brain is related to the propagation of information along structural connections (Vézquez-Rodríguez et al., 2020). In our study we compare different approaches based on GNNs to study the spatio-temporal propagation of information in brain networks from a multi-modal and data-driven perspective.
Recently, several different GNN architectures have been proposed to model the flow of information across time and space in graphical signals (Wu et al., 2020). Convolution operations, often applied in deep learning, have recently been extended successfully to graphical models and allow us to capture inherent spatial dependencies on non-Euclidean network structures (Defferrard et al., 2016). They were later combined with recurrent neural networks (RNNs) (Rumelhart et al., 1986), which can detect sequential relations in signals. This combined spatio-temporal GNN framework was proposed in the notion of diffusion convolution recurrent neural network (DCRNN) (Li et al., 2018). However RNNs can have problems with long time series and, when combined with graph convolution operations, their gradients are more likely to explode or vanish (Seo et al., 2018;Li et al., 2018). This was the motivation for introducing approaches that combine spatial graph convolutions with standard one-dimensional temporal convolutions (Wu et al., 2020). But their receptive field size can only grow if many hidden layers are used (linear growth) or global pooling is applied. To alleviate such shortcomings so-called WaveNets (WNs) have been introduced that employ stacked dilated temporal convolutions, which provide a long-term temporal memory (van den Oord et al., 2016). They have recently been combined with graph convolution operations in an architecture denoted as graph WaveNet (GWN) (Wu et al., 2019). As an alternative method for the temporal processing, also attention mechanisms have been recently included in STGNN architectures (Zheng et al., 2020). Attention mechanisms select, from all inputs, information that is critical to the task at hand and modify edge connection strengths accordingly. They have been applied already to natural language processing, speech recognition and image processing (Xu et al., 2015;Vaswani et al., 2017;Liang et al., 2019), but applications to analyze the dynamics in neural signals are still missing. In this study we compare these different STGNN architectures with each other and evaluate their effectiveness in replicating functional activity distributions observed in brain networks. In addition to these distinct temporal models, we study different approaches to model the spatial information exchange between brain regions. At first we employ the structural connectivity as the substrate for information propagation between brain regions. Further we evaluate the effectiveness of employing connectome embeddings (CEs) of the anatomical network to characterize the node relations. In a recent study Rosenthal et al. (2018) have shown that embeddings of nodes in the anatomical network can inherently capture higher order topological relations between different structurally connected nodes in this network. Finally we compare it to the case when we incorporate no predefined spatial layout into the GNN models, trying to learn the spatial structure by gradient descent based optimization during model training. We demonstrate that by modeling the functional information exchange between regions in STGNNs based on structural connectivity, we can significantly increase the accuracy of predicting future neural signals. This points out that STGNN models are capable of learning informative functional interactions between areas in such brain networks. Based on these different comparisons, we at first try to identify the most effective STGNN architectures to investigate such spatial and temporal dynamics in brain networks.
In a subsequent step we then compare this STGNN based approach to a currently popular data-driven model for characterizing directed functional information exchange in brain networks. Until now methods for the inference of directed functional or effective connectivity are often based on Granger causality (Barnett and Seth, 2013) or dynamic causal modeling (Friston et al., 2013) and its recent extensions (Frässle et al., 2018;Prando et al., 2020). In addition various algorithmic and information theory based methods have been developed meanwhile in this field (Bielczyk et al., 2019;. In the following, we compare the performance of the STGNN prediction models to the forecasting model most often used in Granger causality analysis (Barnett and Seth, 2013). Granger causality is based on the idea that if one event A would cause another event B, then A should precede B and the occurrence of event A should contain information about the occurrence of event B (Friston et al., 2013). In the context of neuroimaging this is realized in a predictive framework, by testing if adding information on activity in a region A improves the prediction of activity in region B. For practical applications of this idea, the underlying predictive model in Granger causality is usually based on a vector auto regression (VAR) for multivariate timeseries inference (Friston et al., 2013;Barnett and Seth, 2013;Bielczyk et al., 2019). In a brain network with N regions of interest (ROIs) the parameters in a VAR model grow with N 2 , so for larger brain networks it can be challenging to accurately fit the model if only limited data are available. This can be problematic in fMRI, where the temporal sampling rate is relatively low, while its good spatial resolution would allow for a detailed, high-resolution network analysis. Therefore it would be desirable to have a predictive model that can learn interactions between all brain areas of interest, and in addition naturally scales to larger brain networks. In our study we compare the STGNN approaches to a classical VAR model and test their accuracy on a variety of network sizes and data set sizes. We show that by learning localized functional interactions based on the anatomical network, the prediction accuracies of STGNN models remain significantly more accurate, even when brain networks become very complex and only limited data are available to fit the models. This points out that the STGNN approaches are robust among a large variety of MRI study scenarios, and are therefore also suitable for the analysis of smaller subject cohorts, like in studies of patients with rare neurological diseases.
Finally we analyze the spatial interactions between brain regions, which are learned by the STGNN models. By integrating prior knowledge on the brain anatomy in form of structural connectivity or based on connectome embeddings, such models can provide multi-modal perspective on directed relations between brain areas. So far, a variety of approaches were proposed to study the structure-function relation in the human brain (Suarez et al., 2020), which are based on computational modeling (Honey et al., 2009;Deco et al., 2012Deco et al., , 2013Chen and Wang, 2018;Messé et al., 2014), graph theory (Liang and Wang, 2017;Abdelnour et al., 2018;Becker et al., 2018;Lim et al., 2019) and machine learning (Sarwar et al., 2021;Deligianni et al., 2016;Rosenthal et al., 2018). Biophysically inspired models for example could describe how functional connectivity patterns emerge from neural dynamics with static couplings characterized by anatomical connections (Deco et al., 2013;Messé et al., 2014;Honey et al., 2009), and were recently used to also study spatial heterogeneity of local circuit properties across the cortex (Demirtaş et al., 2019;Wang et al., 2019). Methods from graph the-ory can supplement such computational frameworks by specifically pointing out how indirect structural connections contribute to the inference of FC strength (Liang and Wang, 2017;Becker et al., 2018). Also hybrid methods could demonstrate that the typology of structural brain networks supports in neuromorphic networks their memory capacity (Suarez et al., 2021). Such insights into structural and functional connectivity can then provide a basis to better understand the cognitive information processing in the human brain (Ito et al., 2020). The vast majority studies on structure-function relations focuses currently on inferring overall FC patterns from their SC, although static coherency based measures of FC might have limitations in their ability to capture the rich nature of dynamic brain activity (Wein et al., 2021a). To the contrary, STGNNs are able to directly predict the measured BOLD dynamics, and their interactions between brain regions, without relying on the indirect representation of functional dynamics based on coherency. This characteristic of STGNNs allows us to additionally investigate the structure-function coupling in brain networks from a novel perspective. To study this relation further on the individual brain region level, we demonstrate how a perturbation based approach can be utilized to reconstruct how dynamic functional interactions are mediated by their structural substrate in STGNN models. By inferring how information is propagated between individual regions in STGNNs, these predictive models have the potential to reveal directed relationships between individual areas in brain networks from a multimodal perspective. In general due to the low temporal sampling rate and physiological artifacts, fMRI can have several limitations in detecting directed relations in brain networks (Smith et al., 2011;Webb et al., 2013). Still some recent fMRI studies and computational simulations could demonstrate that also lag-based methods like Granger causality can be useful for detecting such directed dependencies in fMRI data (Seth et al., 2012;Wang et al., 2014;Mill et al., 2017;Duggento et al., 2018).
The possibility to combine structural and functional imaging data in STGNNs can make these models as well interesting for several practical applications in brain connectivity research. For instance they can be used to investigate differences in the structure-function relationship between resting-state and task-based fMRI. Furthermore, in clinical applications these models could be employed to study how lesions in the structural connectome affect the functional organization of the brain network. In our current study we compare therefore such mechanisms for spatial and temporal modeling in STGNNs, with the objective to establish their methological foundation for brain connectivity research, and thereby providing a basis for future applications of STGNNs in multi-modal neuroimaging studies.

Graph neural network models
In our context of MRI the goal of the spatio-temporal GNNs will be to forecast the observed BOLD signal as accurately as possible in order to precisely capture the spatio-temporal dynamics of the underlying mechanisms in the brain. The learning objective can be formalized by introducing a graph signal x (t) ∈ R N , representing the BOLD signal measured at timestep t in N different brain regions. The goal of the models is to predict from an input sequence of T p past neural activity states t = 1, . . . , T p a sequences of future states t = 1, . . . , T f . In addition to the temporal information, also spatial dependencies are included in the GNN architectures. The spatial relations between the N brain regions can be represented in the notion of a graph G = (V, E, A w ), containing vertices (nodes) V, with |V| = N , and edges E. The structure of the graph is characterized by a weighted adjacency matrix A w ∈ R N ×N , where an entry w nn describes the connection strength between brain region n and n . An illustration of such a graphical representation of a dynamic brain state is provided in figure 1. Based on this concept, the task of the GNN models is to derive a function h(·) that best predicts T f future activity states from an input sequence of T p past states: (1) Figure 1: The spatio-temporal representation of a signal in a brain network is illustrated. The temporal component of the graph-like signal is given by the BOLD signal x (t) ∈ R N in N brain regions sampled at different timesteps t, as shown on the left side. The strength of the edges in the brain network are defined by a weighted adjacency matrix A w ∈ R N ×N , as illustrated on the right side. One entry w nn of the matrix A w characterizes thereby the spatial relation between brain region n and n in the network.
Until now various spatio-temporal GNN architectures have been proposed to account for spatial and temporal dependencies of such graph structured signals (Wu et al., 2020). An overview of different possibilities for the temporal modeling is given in figure 2 (A). In timeseries analysis, recurrent neural networks (RNNs) (Rumelhart et al., 1986) provide one efficient way to detect patterns in sequential data structures, like in our context the BOLD signal, subsequently sampled at different timesteps t ( Figure  2 A1). This approach can be extended to a RNN based sequence-to-sequence architecture, whereby an encoder recursively processes an input sequence of T p past neural activity states x (t) and encodes the temporal information into a hidden state H (Tp) (Sutskever et al., 2014). Next a decoder network uses the information in H (Tp) to generate a prediction for T f future activity states. To account for vanishing gradients during training, the encoder and decoder consist of gated recurrent unit (GRU) cells (Chung et al., 2014). An alternative for detecting repetitive patterns in sequential data is provided by convolutional neural networks (CNNs) (Figure 2 A2). By employing one-dimensional convolutions in the time domain they are used in our context to process temporal dynamics of neural activity. To more efficiently capture long-term dependencies in temporal data the WaveNet (WN) architecture has been proposed (van den Oord et al., 2016). This model introduces dilated causal convolution operations to generate a large receptive field when using only relatively few network layers, which alleviates the processing of long temporal input horizons. The growth of the receptive field of a neuron (marked in green) in a WN layer is illustrated in figure 2 (A2). More recently, also attention mechanisms have been proposed to detect underlying hidden correlations in sequential data structures (Vaswani et al., 2017). In timeseries analysis, the idea of a temporal attention (TAtt) architecture is thereby to adaptively focus on the most important temporal features in a sequence (Figure 2 A3).
These different fundamental approaches for temporal dependency modeling have been recently combined with techniques to additionally capture spatial relationships in graph structured signals (Li et al., 2018;Wu et al., 2019;Zheng et al., 2020). Graph convolutional neural networks can be incorporated to model the propagation of information between adjacent nodes in the graphical representation of the sig- Figure 2: Overview on the different techniques used by spatio-temporal GNNs to learn the spatial and temporal dynamics in brain networks. The temporal component of the STGNN tries to infer dependencies between the activity at different timesteps t = 1, 2, 3, . . . in a certain brain region n, as illustrated in (A). A recurrent neural network (RNN) architecture processes in a recursive manner the activity states at the different timesteps t (A1). In our study this principle will be implemented in a RNN based sequence-to-sequence architecture for timeseries forecasting. As an alternative, convolutional neural network (CNN) architectures employ one-dimensional convolutions in the time-domain (A2). This idea is picked up in the WaveNet (WN) architecture, which introduces dilated convolutions to obtain exponentially growing receptive fields, as highlighted here in green. Following the idea of an attention mechanism, temporal relevance scores can be dynamically computed to weight the importance of a temporal feature observed at timestep t (A3). A temporal attention (TAtt) based architecture is thereby composed of a multitude of stacked attention layers. Based on these temporal modeling approaches, STGNNs additionally detect spatial dependencies between a node n and other interconnected nodes in the brain network, as illustrated in (B). The substrate for the spatial interactions can be based on the structural connectivity matrix A SC , as reconstructed from DTI data (B1). Alternatively, the node2vec algorithm can be used to detect higher order relationships between brain regions in the structural connectome, for characterizing the spatial similarity based on connectome embeddings A CE (B3). Finally, spatial connections can be tried to be learned by the model itself, by freely adapting the edges in an adjacency matrix A Adap at different iterations i, i + 1, i + 1 . . . during the model training (B2).
nal (Defferrard et al., 2016). The neighborhoods of the vertices/nodes V in the network are characterized by the adjacency matrix A w . In our study we additionally investigate different possibilities for defining the spatial layout for the information propagation between brain regions, as illustrated in figure 2 (B). As a first choice for the adjacency matrix we will employ the structural connectivity A SC between the N brain areas, as it could be reconstructed from DTI data (Figure 2 B1). This choice is motivated by the idea that white matter connections obtained from this modality would establish the anatomical backbone for information exchange between brain areas. In a recent study Rosenthal et al. (2018) demonstrated that connectome embeddings (CE) can be utilized for projecting the structural connectome into a continuous vector space, which captures meaningful correspondences between different brain areas. This technique can thereby allow us to additionally account for long range and inter-hemispheric homotopic connections, which are usually only weakly expressed in DTI based anatomical connectivity (Thomas et al., 2014). In our study we utilized this technique to represent the edge weight w nn in the adjacency matrix as the similarity between the vector representations of two nodes n and n , which will be denoted as A CE (Figure 2 B3). The information is accordingly propagated between brain regions which possess high similarity based on their neighborhood role within the anatomical network. Finally we compare these techniques to the case when the model is given the freedom to learn spatial dependencies between the N regions itself. In this setup the adjacency matrix is represented by a self-adaptive matrix A Adap ∈ R N ×N , which is learned during the model optimization (Figure 2 B2). A detailed formal description of the model architectures and the training involved is outlined in the Materials and methods section. In the following we will assess the effectiveness of the different spatial and temporal modeling approaches by comparing their predictive performance on a MRI dataset from the Human Connectome Project (HCP) (Van Essen et al., 2013).

Data description
For the different evaluations in this study, resting-state fMRI data provided by the HCP S1200 release was incorporated (Smith et al., 2013b). To define the nodes of the brain network, the multi-modal parcellation proposed by Glasser et al. (2016) was applied, which is composed of 180 segregated regions within each hemisphere. The average of the BOLD signal was computed within each brain region, so for each restingstate session, N = 360 time courses were obtained (180 per hemisphere). During one session T = 1200 fMRI images were collected, so that the ROI timeseries can be represented by a data matrix X ∈ R N ×T . We filtered the resting-state fMRI timeseries data with a 0.04 − 0.07Hz Butterworth bandpass filter, because this frequency band has shown to be most reliable and functionally relevant for gray matter activity (Glerean et al., 2012;Buckner et al., 2009;Deco et al., 2017;Biswal et al., 1995;Achard et al., 2006).
For learning the predictions of the BOLD signal, samples of input and output sequences were generated from the timeseries data in X (Wein et al., 2021b). This was achieved by selecting windows of length T p to obtain input sequences of neural activity states [x (1) , . . . , x (Tp) ], and respective target sequences of length T f denoted as [x (Tp+1) , . . . , x (Tp+T f ) ]. The time index t was propagated through each fMRI session, where in total T − T p − T f + 1 input-output pairs were generated per session. For each fMRI session the first 80% of those time window samples were used as the training data for the models, the subsequent 10% as a validation set, and the last 10% were employed for testing. For the following comparisons, the length of the input and output sequences were selected to be T p = T f = 60, which corresponds to a time span of roughly 43 s, based on a sampling interval of T R = 0.72 s . This time window has been shown to be long enough to be sufficiently challenging for the models and to make clear the differences in their performance. Likewise, the time window of 60 timepoints is short enough for them to make reasonable non-random forecasts of the BOLD signal.
In addition to the functional dynamics in the different brain regions derived from fMRI, the structural connectivity between those regions was reconstructed from DTI data. For this purpose the DTI dataset in the HCP S1200 release was processed using the multi-shell, multi-tissue constrained spherical deconvolution model (Jeurissen et al., 2014), made available in the MRtrix3 software package (Tournier et al., 2019). White matter tractography was performed to estimate the anatomical connection strength between the regions defined by the multi-modal parcellation atlas (Glasser et al., 2016). The number of the streamlines which connect two atlas regions was used to determine the structural connectivity values between the N brain regions, which were then collected in a structural connectivity matrix A SC ∈ R N ×N . A detailed description of the MRI datasets and their preprocessing is provided in the Dataset section. In addition the embeddings of the nodes within the structural network A SC ∈ R N ×N were generated using the node2vec algorithm (Grover and Leskovec, 2016). The parameters for this algorithm are outlined in detail in the Spatial dependencies section and further Pearson correlation was used to quantify the degree of similarity of structural nodes in their connectome embedding space. The pairwise similarities between the N nodes were then collected in the matrix A CE ∈ R N ×N .

Comparison of GNN architectures
Before evaluating the performance of different models on a larger variety of MRI study scenarios, we will first focus on the effects of different temporal and spatial modeling techniques. For this purpose a dataset with a sample size of a medium sized fMRI study including 25 subjects will be incorporated. Each resting-state fMRI session was decomposed into pairs of input and output samples, as discussed in the Data description section, and the generated training, validation and test samples were then aggregated across the 25 fMRI sessions. The neural signal in regions within the right hemisphere (Glasser et al., 2013), consisting of N = 180 ROIs, will be included in the following comparison. At first we evaluate the prediction accuracy of the different temporal modeling strategies. For this purpose we compare the recurrent neural network (RNN) model, with the WaveNet (WN) model and the temporal attention (TAtt) model. The influence of the model hyperparameters, which are used in the following comparisons are described in the Model training section and discussed in detail in Supplement I. The BOLD signal data was scaled to zero mean and unit variance for the evaluations, to obtain values of a magnitude that is easier to interpret. Figure 3 (A) shows the test mean absolute error (MAE) between the predicted and the true neural activity. By generating windowed input-output pairs of activity values from the fMRI data, the last 10% of samples from each session correspond to 108 of such input-output pairs per session for testing, each containing 60 output time points (corresponding to roughly 43s of activity). The overall test errors were computed as the average across all these test samples from the 25 subjects and across all 180 brain regions. The comparison shows that RNN and WN have very similar capabilities in predicting the BOLD signal, while the TAtt model exhibits a worse performance. To test the significance of this difference between the models, we further computed the test MAE of each individual subject as an average across predicted timepoints, brain regions and test samples per subject. By applying a paired t-test, the differences between the WN and RNN model to the TAtt model were shown to be both highly significant with p ≤ 0.0001 across subjects (Cohen's d = 10.68 and d = 10.66 respectively). In addition to the MAE, we evaluated these models in Supplement III using scale-free measures like R-squared (R 2 ) and the similarity of the predicted FC states. Despite their conceptual differences the results show that the RNN and WN based approach can both recover a comparable and consistent amount of temporal information from the fMRI data. In comparison to these, the TAtt architecture appears to be less suitable to accurately predict the BOLD signal with this limited amount of data. For this reason in the following we will focus on RNN and WN based approaches for identifying suitable models to model functional dynamics in brain networks. The overall test error was computed as an average across samples, brain regions and subject sessions. The errorbars represent the standard deviation of the test MAE across subjects. Due to their high accuracy in the temporal domain, we focus on RNN and WN based approaches for forecasting the spatio-temporal dynamics in the following. Spatial relations are added to the temporal models in form of graph convolutions, and the spatio-temporal extension of the RNN and WN models are respectively denoted as diffusion convolution recurrent neural network (DCRNN) and graph WaveNet (GWN) (Li et al., 2018;Wu et al., 2020). Spatial transitions are based on the relations of network nodes captured in a weighted adjacency matrix, which is either based on structural connectivity (A SC ), connectome embedding similiarity (A CE ), or adapted during model training (A Adap ). In (B) the adjacency matrix A SC based on structural connectivity within the 180 regions of the right hemisphere is illustrated, together with the adjacency matrix A CE derived from structural connectome embedding similarities. The regions in this illustration are ordered according to the atlas proposed by Glasser et al. (2016). Figure (C), (D) and (E) show the prediction accuracies of the DCRNN and GWN model in dependence of the walk order K. In figure (C) the overall test MAE is shown when incorporating the SC as an adjacency matrix A SC , figure (D) illustrates the test MAE when employing CEs in an adjacency matrix A CE to define spatial relationships, and (E) displays the case when using a self-adaptive weight matrix A Adap .
In the next step we will study the impact of adding information on spatial relations between the different regions in the brain network. This will be implemented by invoking graph convolution operations in the predictive models, as outlined in detail in the Spatial dependencies section. The definition of an adjacency matrix determines how information is propagated between the different nodes in our brain network, and in our evaluations we investigate three conceptually different possibilities. In the first approach we use the structural connectivity as derived from DTI as the substrate for information exchange between different ROIs. The SC based adjacency matrix A SC is illustrated in figure 3 (B). The information can propagate along direct connections in the graph (K = 1), but also higher orders (K = 2, 3, . . .) expressing the influence of indirect connections can considerably contribute to interactions between different areas in the brain (Becker et al., 2018;Liang and Wang, 2017;Bettinardi et al., 2018). A walk order of K = 0 denotes the case when including no spatial information exchange between network areas, exclusively incorporating temporal information for the predictions. Figure 3 (C) depicts the test MAE in dependence of the walk order K when using the SC derived from DTI as a basis for information propagation in space. The RNN based model in combination with graph convolution operations is referred to as DCRNN (Li et al., 2018) and the MAE of its predictions, averaged across test samples, brain regions and predicted timepoints is depicted in blue. Figure 3 (C) shows that it has the lowest test MAE when incorporating walks on the structural graph up to a order of K = 2. The WN incorporating graph convolution operations is denoted as GWN (Wu et al., 2019) and its average test MAE is shown in red in figure 3 (C). The influence of the walk order K on the GWN accuracy suggests that its performance can be successively improved by including first-order connections, followed by the second-and third-order connections. As an alternative thereto, the structural similarity between ROIs can be based on their CE similiarity A CE , as illustrated in figure 3 (B). The comparison between A CE and the structural connectivity matrix A SC highlights that in the adjacency relation defined by the structural embeddings, long range connections between brain regions are considerably more pronounced. Figure 3 (D) shows the test MAE of the models when incorporating A CE in the graph convolution operations. In this case we can observe for both models a sharp drop in the error at walk order K = 1. This suggests that the node embeddings already inherently capture higher order relations between nodes in the brain network. Finally in figure 3 (E) the test MAE is shown when treating the connections between nodes as learnable weights. In this case we do not observe an improvement in the test error. This observation indicates that it is rather challenging to learn all N 2 connections between brain regions without any prior knowledge. In general both STGNN models could profit the most when using CEs to characterize the spatial layout for functional interactions between brain regions. For the DCRNN the test error was M AE = 0.1388 when incorporating no information from other brain regions in the network, and could be reduced to M AE = 0.1158 (for K = 1) when using CEs to model the information exchange within the brain network. To test whether incorporating information about structural connections significantly increases the prediction accuracy of our models, we at first recomputed the overall test MAE for each subject again. Then by using a paired t-test, we find that, for both STGNN models (DCRNN and GWN) and both adjacency types (A SC and A CE ), the impact of structural modeling is positive (Cohen's d > 1 for all comparisons) and significant (p ≤ 0.0001 for all comparisons), compared to the case in which it is not considered. Although the performance differences between the GWN and DCRNN are quite small in general, the DCRNN slightly outperformed with a test error of M AE = 0.1158 the GWN with a test error of M AE = 0.1211 at K = 1 (significant with p ≤ 0.0001, Cohen's d = 0.49). In addition the distribution of the test error across subjects and ROIs, with and without the structural modeling in STGNNs is illustrated in Supplement III in figure S6. This demonstrates that around 17% more information on functional dynamics can be directly retrieved from nodes with similar context within the structural network. Using the SC to model transitions could only reduce the MAE of the DCRNN by 5% at K = 1. This observation supports the idea that structural node embeddings can strengthen the relationship between structural data derived from DTI with functional data observed in fMRI (Rosenthal et al., 2018). When applying a paired t-test, the improvement of the prediction accuracy when using the CE similarity in comparison to the SC became for both, the DCRNN and GWN model, significant with p ≤ 0.0001 at K = 1 (Cohen's d = 1.45 for the DCRNN and d = 0.95 for the GWN). By inherently capturing higher order transitions in A CE , only a low walk order K is required to capture information from structurally connected ROIs. In this manner, this technique can help to efficiently reduce the number of necessary parameters to account for spatial dependencies in STGNN models.

Model accuracy and network scaling
In this section we study the prediction accuracy of the above introduced STGNN based approaches and compare it to the VAR model, which is currently most often used for directed functional connectivity analysis (Friston et al., 2013;Barnett and Seth, 2013). In practicable applications the amount of available fMRI data may vary depending on the project size and on the recruited subject cohort. Also the size of the brain network of interest can range from a few specific areas in a single functional network to a large-scale whole brain analysis. For this purpose we consider different scenarios in our following evaluations, by analyzing the models accuracies in dependence of the brain network size and the fMRI dataset size. We consider one larger subject dataset consisting of resting-state fMRI sessions from 50 different subjects, one medium sized dataset of 25 subjects and one small dataset including data from 10 subjects. In addition we vary the size of the analyzed brain network. The first network consists of 22 ROIs per hemisphere involved in visual processing as defined by the Glasser parcellation (Glasser et al., 2016) (a complete list of selected ROIs is provided in the Supplement II). The second network includes the regions within one hemisphere, and for that purpose the 180 regions within the right hemisphere included in the Glasser atlas were selected (Glasser et al., 2016). Finally the whole brain network of in total 360 regions was incorporated. As discussed in the Data description section, windowed input and output time sequence pairs were created from the data and the goal of the different models is accordingly to predict T f = 60 TRs of neural activity from the past T p = 60 activity values. We fitted the VAR model using the ordinary least squares (OLS) method as implemented in the multivariate Granger causality (MVGC) toolbox (Barnett and Seth, 2013), and for each dataset we selected the VAR model with order p that achieved the best MAE on the test set, as outlined in more detail in the Vector auto regressive model section. The hyperparmeters used for the STGNNs are described in the methods part in the Model training section. Further for this comparison the CE similarity A CE with transition order of K = 1 was used to define the structural relations in the STGNN models, which has shown to improve the GNNs forecasting accuracy with low computational cost, as discussed in the section Comparison of GNN architectures. Figure 4 shows the test accuracy of the VAR, DCRNN and GWN model in dependence on the dataset size and brain network size. It can be observed in figure 4 (A) that if a large dataset of 50 subjects is available, all models are able to accurately predict the BOLD signal with a low test MAE, and a notable increase in the test error only appears for the VAR model, when it is fitted to the whole brain network. Figure 4 (B) shows the test MAE when data from 25 subjects is incorporated. In this case the test error of the VAR model starts to increase noticeably when modeling activity distributions within one hemisphere and becomes quite large when including the whole brain network. In contrast to these, the prediction accuracies of the DCRNN and GWN models remain stable in all cases. Finally when only 10 subject datasets are available, the test MAE of VAR model is highly dependent on the analyzed network size, as illustrated in figure 4 (C). The DCRNN and GWN model can still achieve a high accuracy, also when a limited amount of data are available and the network size is relatively large. In addition, this comparison of the models is replicated in Supplement IV using additional measures like R 2 and the similarity of predicted FC states. After applying a paired t-test, the differences between the DCRNN and GWN to the VAR were shown to be in all cases highly significant with p ≤ 0.0001 (Cohen's d 1), except when the VAR is only fitted to the single visual network, where it still could make reliable forecasts. The test MAE of the VAR is here depicted in orange, the MAE of the DCRNN in blue and the error of the GWN in red. The overall error was computed as an average across brain regions, timesteps and test samples. The error bars represent the standard deviation of the test error across subjects. In (A) the test MAE using a dataset of 50 subjects is shown for the visual network, the network within the right hemisphere and the whole brain network (Glasser et al., 2016). Figure (B) and (C) show the test performances in dependence of the network size using the 25 and 10 subject dataset, respectively.
To illustrate the prediction accuracies of the different models in more detail, an example of the predictions using the dataset including 25 subjects, and modeling the activity within one hemisphere is shown in figure 5. Figure 5 (A) shows the MAE of the models computed as an average across test samples and ROIs in dependence of the forecasting horizon. In this case within the first 15 predicted timesteps all three models can generate very accurate predictions, but after that period the error of the VAR model starts to accumulate, while the GNN based approaches remain considerably more stable and precise. The predicted BOLD signals of the different models in a few representative samples are shown in figures 5 (B), (C) and (D). Additionally, the predicted FC states A F C ∈ R N ×N were computed as the Pearson correlation between predicted BOLD signals of all N brain regions, and a comparison of representative predictions with the true FC state is illustrated in figure 5 (E). The average correlation of the predicted FC state to the true FC state was for the VAR model r F C = 0.635 on this dataset, while the GWN could achieve a correlation value of r F C = 0.948, and the DCRNN a value of r F C = 0.950. The overall FC similarity for all different datasets of all prediction models is given in Supplement IV figure S8. Furthermore, in Supplement V in figure S9 we performed the analysis on the same dataset using a more liberal frequency filtering within the 0.01 − 0.1Hz frequency range. In this range the signal dynamic becomes more complex and we can observe an increase in the prediction error of the different models accordingly.
In addition we evaluated in more detail how the prediction errors are distributed across different subjects and different ROIs. Figure 6 shows the distribution of the test MAE of the DCRNN, GWN and VAR model across subjects, and in dependence of the brain region within the right hemisphere. For all three models we observe a consistently greater prediction error in the posterior cingulate cortex and medial orbitofrontal cortex, which could point towards a more complex BOLD dynamic in those regions. Alternatively, the prediction accuracy might be also affected by a lower signal-to-noise ratio observed in medial brain regions (Olman et al., 2009). Figure 5: The prediction accuracy of the different models is presented in more detail for the 25 subject dataset and the brain network including the ROIs within the right hemisphere (Glasser et al., 2016). In (A) the test MAE in dependence of the forecasting horizon is shown, computed as an average across test samples and brain regions.

Multi-modal directed connectivity
In the Comparison of GNN architectures section different approaches have been investigated to model functional interactions between segregated regions in the brain network. The results showed that in- Figure 6: The distribution of the test error across subjects and brain regions is shown. In (A) the MAE across subjects and brain regions of the DCRNN is first visualized in a boxplot on the left side. Additionally on the right side of the figure, the MAE values are projected onto the cortical surface within the right hemisphere, where the colormap was linearly scaled between 0 and 0.18. In (B) the distribution of the test MAE of the GWN is shown and in (C) the MAE distribution of the VAR model. For the VAR the colormap was adjusted to account for larger error values by scaling it between 0 and 0.6. corporating information on the spatial relation between ROIs in the form of the structural connectivity or connectome embedding similarity could considerably improve the prediction accuracy of the GNNs models. This points out that the GNNs are able to learn relevant and functional informative transitions of neural activity on their structural spatial layout. Based on the idea of Granger causality (Granger, 1969) that the observation of one event A carries information about the occurrence of a future event B, this might represent initial evidence for a potentially causal relation between A and B. Due to the relatively low temporal sampling rate and physiological artifacts in fMRI (Smith et al., 2011;Webb et al., 2013), it still is a matter of discussion to what extent we can observe a causal relationship between brain regions in this imaging modality (Wang et al., 2014;Mill et al., 2017;Bielczyk et al., 2019). But the observation that activity in one region carries additional information among activity in another region in the brain can also go beyond simple undirected FC or SC, and we therefore refer to this kind of relation as directed connectivity in the following. Propagating the information between ROIs related to their SC or structural CE similarity has the potential to give us in this manner a multi-modal perspective of such a directed relationship between different brain areas. For this purpose we choose a perturbation base approach to reconstruct the amount of information individual ROI carries about other ROIs (Zeiler and Fergus, 2013). By learning a function h(·), the GNN models try to infer from an input sequence of neural activity states [x (1) , . . . , x (Tp) ] a sequence of future activity states [x (Tp+1) , . . . ,x (Tp+T f ) ], whereby x (t) ∈ R N denotes the activity at timestep t in all regions n = 1, . . . , N . To induce an artificial perturbation into the system of neural dynamics, we remove all activity in a certain ROI n by setting its activity values to the sample mean x n = 0. By using the perturbed timeseries as an input for our trained model h(·) the model generates then a prediction [x (Tp+1) , . . . ,x (Tp+T f ) ]. To reconstruct the directed influence of ROI n on ROI n we evaluate the overall difference between the original prediction and the prediction with perturbation in the input: where I n (n ) denotes the impact of ROI n on n. Furtherx n (s) denote the predictions in ROI n with and without the perturbation in n of one test sample s at time step t. Note that this artificial perturbation approach has only the goal of making spatial relations between ROIs in STGNN models explainable, and should not be equated with the effect of an experimental perturbation applied to the human brain, like for example induced by transcranial magnetic stimulation (TMS).
In the following we compare this proposed measure of directed influence I(n ) to the classical undirected types of brain connectivity. First we compare it to structural connectivity as derived from DTI, characterized by the number of fiber tracks connecting two brain regions. Then we incorporate functional connectivity, defined as the Pearson correlation of functional activity timecourses between two areas. We employ the above introduced GWN model to obtain a multi-modal measure of directed connectivity I(n ), first using the SC as substrate for information propagation, captured in A SC , and then also employing the similarity of CEs, represented by A CE . In the following example we study the connectivity of V1 within the right hemisphere by incorporating data of 25 subjects. For this purpose a perturbation was induced into the target region V1, and the impact of this perturbation on all other 179 regions within the right hemisphere was then characterized by computing the measure of directed influence I(n ), as defined in equation 2. These values for directed connectivity strength can then be visualized by projecting them onto the 179 other areas of the cortical surface. For the following comparison, all connectivity values were rescaled by normalizing them between 0 and 100. At first in figure 7 (A) the structural connectivity between V1 and all other 179 regions within the right hemisphere is depicted. The target region V1 is marked here in light blue, and the strength of connectivity to all other regions is encoded in red. Figure 7 (A) shows that we can mainly observe a pronounced structural connectivity between V1 and V2 and some structural connections leading to V3. Figure 7 (B) shows the undirected functional connectivity in resting-state. In this type of connectivity we can observe predominantly correlations to the functional activity in V2 and V3, but also a considerable connection strength to V3, V4 and V6. In figure 7 (C) the directed connectivity strength I(n ) is depicted, when using the SC as spatial backbone for the information exchange between brain regions in the GNN. In comparison to the SC, in this variant of brain connectivity we can observe in addition to V2 also a more pronounced relationship to areas V3 and V4, and to some anatomically more distant areas like V6 and the ventromedial visual area VMV1. This multi-modal type of connectivity also reflects the role of indirect structural connections by modeling higher order transitions on the structural scaffold captured by the STGNN model. As an alternative to the SC, in figure 7 (D) the directed connectivity patterns when using CE similarity as the spatial layout in the GWN are displayed. Here we can see an even stronger integrity of V1 within the visual network, which is in agreement with the observation that CEs capture higher order topological information of anatomical connectivity (Rosenthal et al., 2018). Figure S10 in Supplement VI shows Figure 7: Different types of connectivity are illustrated between V1 and all other regions within the right brain hemisphere. In (A) the structural connectivity is shown, whereby the target region V1 is marked in light blue and the connectivity strength is encoded in red. In (B) the correlation based functional connectivity is illustrated, which was computed as an average across subjects. Further (C) shows the measures of influence I(n ), derived from the GWN model using the SC for information propagation and figure (D) depicts the influence when incorporating CEs for the information exchange between ROIs. The values of the connectivity measure were linearly mapped between 0 and 100 (and between −100 and 100 for FC). The default scaling of the color values provided by the connectome workbench (version 1.4.2) was used, adjusting the colormap between the 2th and 98th percentile of the values respectively. additionally the spatial relations learned by the DCRNN model. Here we can observe a pronounced similarity to the directed connectivity pattern learned by the GWN architecture, showing additionally strong relations to areas like V3 and V4. Based on this observation, such a GNN based connectivity approach can serve as a link between structural and functional connectivity, and as such they can provide a multi-modal perspective on directed influences between individual areas in brain networks. So far, we have only studied the effect of a single artificial perturbation in V1 directed to all other regions within the right hemisphere. This approach can further be extended to sample a full connectivity network, by systematically inducing perturbations into all regions of interest for the analysis, and then systematically observing the effect on the other network regions. In figure S11 in Supplement VI we studied the effects of perturbations induced in some different additional areas of the visual network based on this approach. In this manner, this technique can allow us to reconstruct the directed spatial relations between brain areas captured in STGNN models, and could be applied in practical applications by for example comparing these connectivity patterns between different conditions in task-based fMRI, or studying the difference between healthy and diseased brain states.
A study of Shinn et al. (2021) demonstrated that FC topology in resting-state fMRI is shaped by and can be predicted from spatial and temporal autocorrelations, as typically observed in fMRI data. To also investigate to what extend STGNN based connectivity patterns are related to such correlations, we computed the temporal autocorrelation as the Pearson correlation between the BOLD signal values and its lagged values with different lag orders j, depicted in figure S12 in Supplement VI. We could observe a relatively high temporal correlation of 0.89 around lag order j = 14, and of 0.62 around lag order j = 27, what shows that within these first 30 TRs of the signal, such temporal autocorrelations can still be detected. In section Model accuracy and network scaling we could show that STGNNs were also able to make reliable long-term predictions of the BOLD signal up to a horizon of 60 TRs, what demonstrates that STGNNs capture properties of neural activity dynamics that go clearly beyond the range that is shaped by temporal autocorrelations. In addition, spatial autocorrelations can also play a distinctive role in shaping FC network properties (Shinn et al., 2021). The analysis in section Comparison of GNN architectures showed that by modeling the information exchange between structurally connected brain regions in the spatial domain, the prediction accuracy of the STGNNs could be significantly improved, in comparison to the null-models, which incorporated no spatial dynamics. The observation that spatially connected regions contain some additional relevant information on functional dynamics points out that spatial interactions captured in STGNNs go beyond simple correlation based spatial network relations. A more detailed comparison of the individual differences and similarities between correlation based FC and the STGNN based connectivity pattern is additionally provided in a barplot in Supplement VI (figure S13).

Conclusion
In this study we have compared different STGNN architectures for learning the spatio-temporal dynamics in brain networks. First in the section Comparison of GNN architectures we studied different mechanisms for learning the temporal dynamics in the BOLD signal. We could show that a RNN based model and a WN based model exhibit very similar capabilities in learning the temporal characteristics in neural activity timeseries. Despite their conceptual differences in their architectures, they demonstrated almost the exact same prediction accuracy, which indicates that they are both very consistent in capturing the temporal information in the data. As an alternative, we also studied TAtt mechanisms to learn temporal characteristics of neural signals. The TAtt model showed to be less suitable to model the dynamics in the BOLD signal with a limited amount of fMRI data. Despite incorporating techniques into the TAtt model that in general stabilize the learning, like residual connections and batch normalization (He et al., 2016;Ioffe and Szegedy, 2015), its prediction error was considerably higher in comparison to the RNN and WN based approach. This indicates that the geometric assumptions which are realized by the temporally structured inference in the RNNs and WNs based on either recurrent computations or causal convolutions can contribute to the learning of the temporal characteristics of the BOLD signal. We then studied the impact of adding spatial dependencies to the temporal models, realized by invoking graph convolution operations. We have compared different spatial layouts for information propagation between ROIs, and therefore included either the structural connectivity (A SC ), the CE similarity (A CE ), or a self-adaptive adjacency matrix (A Adap ) into the STGNN models. While the model performance of the GWN and DCRNN steadily improved with higher walk orders K on the anatomical substrate, we could observe a more pronounced improvement already when using CEs with a walk order of only K = 1. This embedding strategy turns out to be therewith also interesting in applications of STGNNs, because it helps to effectively incorporate indirect structural connections with low computational cost. In addition, the observed characteristics of CEs in our application support the ideas of Rosenthal et al. (2018), which showed in their study that embeddings of the structural network can naturally capture higher order topological relations between ROIs within the structural layout. In our context of modeling spatio-temporal dynamics this method also proved to strengthen the relationship between brain structure and functional dynamics. Learning all N 2 connections of the underlying structural graph during the model training has been shown to be challenging for the STGNN models, in case no prior knowledge is provided to them in form of the anatomical brain connectivity. While such highly parameterized artificial neural network models can be in theory quite flexible in learning complex relations (Hornik et al., 1989;Brüel Gabrielsson, 2020), often the decisive limitation is the successful optimization of the parameters during model training (Dauphin et al., 2014). In the discussed applications of STGNNs in fMRI, where the amount of training data is most often quite limited, prior knowledge in form of the anatomical graph structure has been shown to considerably support the learning of spatial relations between brain areas captured in STGNN models.
So far methods based on biophysical modeling (Honey et al., 2009;Deco et al., 2012Deco et al., , 2013Messé et al., 2014), graph theory (Liang and Wang, 2017;Abdelnour et al., 2018;Becker et al., 2018;Lim et al., 2019) or machine learning (Deligianni et al., 2016;Rosenthal et al., 2018;Sarwar et al., 2021) have contributed already numerous valuable insights into the structure-function relation in brain networks, and could highlight the role and importance of the structural connectome in shaping functional connectivity patterns. While the majority of approaches studying the structure-function relationship infer brain dynamics by fitting the models to empirically observed FC patterns, STGNNs provide us with a possibility to directly predict the observed neural activity states. Similar to some other recently proposed predictive models (Singh et al., 2020;Suarez et al., 2021), this principle can allow us to investigate additional interesting aspects of dynamic brain functions. As discussed in section Comparison of GNN architectures this could enable us for example to study directly the amount of information on the activity of one ROI that is contained in the activity of other structurally connected ROIs. For a comparison with other currently used approaches investigating SC-FC mappings, the predicted BOLD signal states of STGNN can be used then again to reconstruct predictions for FC states, as shown in figure 5 (E). The relatively high accuracy in predicting empirical FC states already point out the potential of STGNNs in this field. Moreover, in comparison to other currently popular approaches used in this area , by learning localized graph filters in STGNNs, their forecasting accuracy is also robust with regards to the brain network size. While such highly parameterized artificial neural network models appear to be promising for achieving high prediction accuracies of FC states (Sarwar et al., 2021), they cannot provide us with the same mechanistic insights into physiological processes as biophysically inspired models. Still they can be used to supplement current biophysically inspired models, for studying different aspects of the structure-function relationship from a novel data-driven perspective. A more comprehensive comparison of these different new approaches, evaluating in detail their interrelations like in the study of , could be thereby interesting for future studies in this area.
In the Model accuracy and network scaling section we have compared the STGNN models to a VAR model, which is currently most often used in Granger causality analysis for inference of directed relationships between brain regions (Barnett and Seth, 2013). We evaluated the accuracy of the different approaches on a variety of brain network sizes and data set sizes to account for different possible scenarios in their application in fMRI studies. The results showed that if a sufficiently large cohort of 50 subjects is available, also a VAR model is able to make very reliable long-term predictions, and only for a large network consisting of N = 360 there is a notable increase in the prediction error. But the dependency of the accuracy on the network size N becomes more apparent when data from only 25 subjects are used to fit the VAR model, and when only 10 subjects are available, the error grows strongly with N . This demonstrates that a VAR is a very reliable and fast model for fMRI studies with a sufficiently large test subject size and for connectivity studies including a limited amount of pre-defined regions. However it can be desirable in some cases to include a larger amount of brain areas into the connectivity analysis, in order to avoid omitting relevant areas in the network of interest. Also in MRI studies it can be very costly and time-consuming to collect a large amount of data, which is, for example, especially challenging in studies on rare neurological disorders. A classical VAR model fits a parameter for every possible connection between the N regions in a network, so that the number of parameters in a VAR based approach grow strongly with an order of N 2 . In contrast thereto, STGNNs utilize prior information in the form of the anatomical connectivity, and then model the functional information exchange based on this underlying structural substrate. By incorporating graph convolution operations in STGNNs, the amount of parameters only linearly scale with walk order K, which can even be chosen to be K = 1, if higher order structural relations are already expressed in an adjacency matrix derived from connectome embeddings (A CE ). This property allowed STGNNs to make very robust inferences also on large networks and when only limited data are available, thereby providing a flexible method for various connectivity analysis scenarios.
Finally in the Multi-modal directed connectivity section we studied the individual spatial interactions within the brain network which were learned by the STGNN models. By integrating information on the anatomical connectivity into the GNN based models, we could derive a multi-modal connectivity measure for directed relationships between brain regions. When comparing this measure of influence to the original structural connectivity, we can observe that STGNNs have learned to include transitions along higher order structural connections in the network. The models could infer links between V 1 and V 2, but additionally strong connections to V 3 and V 4. Especially when incorporating the CE based similarity A CE to define spatial node relations in the STGNN models, we can observe a high integration of V 1 within the visual system. However due to the relatively low temporal sampling rate in fMRI (Friston et al., 2013), and the indirect measurement of neural signals based on their hemodynamic response (Webb et al., 2013), one should also be aware of these limitations in the inference of directed and potentially causal connections in fMRI studies (Smith et al., 2011). Our lag-based predictive approach based on STGNN models might therefore also be affected by the same limitations as classical Granger causality in fMRI. On the other hand, a combined fMRI-MEG study by Mill et al. (2017) and different computational simulations of fMRI data (Wang et al., 2014;Seth et al., 2012;Wen et al., 2013;Duggento et al., 2018) could establish evidence that Granger causality is still able to identify meaningful directed relationships between brain areas in fMRI, despite the indirect measurements based on the hemodynamic response. As an alternative, deconvolution based approaches can have the potential to infer from the measured BOLD signals the underlying neural timeseries (Bush et al., 2015;Mill et al., 2017) for assessing effective brain connectivity, rather than only estimating directed functional connectivity. But the estimation of the underlying hemodynamic response from the data might come with the cost of introducing additional assumptions and uncertainties (Roebroeck et al., 2009;Bielczyk et al., 2019). A more detailed discussion on considerations concerning Granger causality, and in general causal inference in fMRI, is provided in the comprehensive review of Bielczyk et al. (2019), as well as in the perspective on FC and its variants by Reid et al. (2019). Despite these current limitations in fMRI, a multi-modal GNN based approach can allow us to join structural and functional imaging data in a new manner, and reveals thereby potential for supplementing current analysis methods in brain connectivity research by studying such directed relations under a novel perspective (Reid et al., 2019).
In conclusion, in our study we found that the DCRNN and GWN architecture are both suitable for the task of functional dynamics inference. Using CEs to characterize the structural similarities between brain regions could further improve their prediction accuracy. Their robust scaling properties and the possibility to combine the information in structural and functional MRI data reveal the potential of STGNNs in the field of brain connectivity analysis. Besides their applications in fMRI, other functional neuroimaging techniques like electroencephalography (EEG) or magnetoencephalography (MEG) might be interesting for analyzing temporal dynamics with STGNNs in the high frequency range. While in this presented approach, we only incorporated a single temporal feature (the BOLD signal) into the STGNNs, in general such a flexible data-driven approach could be expanded to account for different types of data and annotations. For example the activity measured in a combined EEG-fMRI experiment (Abreu et al., 2018;Mele et al., 2019) could be also simultaneously integrated in STGNNs as different temporal features, or adding the temporal response of a subject could be helpful to better predict activity patterns in task-based fMRI. Also alternative structural imaging techniques like neurite orientation dispersion and density imaging (NODDI) (Zhang et al., 2012) might capture additional aspects of the brain structure, which could be included as structural information in STGNN based models. In clinical applications multi-modal STGNNs could be interesting for studying how the relationship between structure and function is affected in the diseased brain (Panda et al., 2021), or which impact a structural lesion might have on the functional organization of the brain network (Alstott et al., 2009). Still the research on GNNs is a relatively new field in machine learning and recent developments in this field can make interesting contributions to our understanding of information processing in brain networks (de Haan et al., 2020;Schnake et al., 2021).

Dataset
The MRI data set used in our study is provided by the HCP data repository (Hodge et al., 2015;Van Essen et al., 2013). As part of the HCP protocol, the study participants gave written informed consent to the HCP consortium. The MRI scanning protocols were approved by the Institutional Review Board at Washington University in St. Louis. We incorporated data of the S1200 release, which provides data from resting-state fMRI sessions, each with a duration of 14.4 minutes, whereby 1200 volumes were sampled per session. The data was acquired with customized Siemens Connectome Skyra magnetic resonance imaging scanners with a field strength of B 0 = 3T , using multi-band (factor 8) acceleration Feinberg et al., 2010;Setsompop et al., 2012;Xu et al., 2012). A gradient-echo echo-planar imaging (EPI) sequences with a repetition time T R = 720 ms and an echo time T E = 31.1 ms was used. The field of view of the fMRI sequence was F OV = 208 mm × 180 mm and in total N s = 72 slices with a slice thickness of d s = 2 mm were collected, containing voxels with an isotropic size of 2 mm. The preprocessing of the HCP fMRI data includes corrections of gradient-nonlinearity-induced distortions, registration to a single-band reference image to account for subject motion, and registration to the structural T1w image (Glasser et al., 2013;Jenkinson et al., 2002Jenkinson et al., , 2012Fischl, 2012). Further ICA-FIX was applied to automatically classify and remove artifactual components in the resting-state fMRI data (Smith et al., 2013b;. Finally the volumetric fMRI images are mapped into the CIFTI grayordinate space and Gaussian surface smoothing with a FWHM of 2 mm is performed. A detailed description of the standard minimal preprocessing pipelines of the HCP can be found in Glasser et al. (2013). In a next step to define our brain network, the multi-modal parcellation proposed by Glasser et al. (2016) was applied, which divides the cortical surface into 180 segregated areas per hemisphere. The BOLD signal within each area was averaged, to obtain the temporal activity evolution for each node in our brain network. For this study we considered it useful to apply global signal regression in our preprocessing (Power et al., 2016). Firstly, in a systematic comparison of different preprocessing methods to address motion artifacts, Ciric et al. (2017) could show that an ICA-based de-nosing in combination with global signal regression is among the most effective methods to reduce movement artifacts. This result is in line with the study of Burgess et al. (2016), investigating the effect of ICA-FIX in combination with global/grayordinate signal regression on resting-state fMRI data provided by the HCP. Furthermore in our study of functional interactions between specific brain regions, the objective was to extract the additional information, which certain regions contain about the activity in other regions. Therefore local interactions rather than global modulations in the signal were of main interest for our analysis (Power et al., 2016). The time courses were then bandpass filtered in the 0.04 − 0.07Hz frequency range. In a summary of several different studies that account for different artifacts in the BOLD signal related to MRI scanner drift in the frequency range below 0.015Hz (Smith et al., 1999), respiratory and cardiac frequencies around 0.3Hz and 1 − 2Hz respectively (Biswal et al., 1996), and fluctuations in arterial carbon dioxide level around 0.0−0.05Hz (Wise et al., 2004), the study of Glerean et al. (2012) identified the 0.04 − 0.07Hz frequency band to be most reliable and relevant for gray matter activity in resting-state fMRI (Achard et al., 2006;Zou et al., 2008;Buckner et al., 2009). To additionally ensure that the low frequency signals are not mainly related to respiratory artifacts, we studied the respiratory signals recorded with a Siemens respiratory belt during resting-state fMRI, as provided by the HCP (Smith et al., 2013b). The average respiratory frequency spectrum is depicted in Supplement VII in figure S14, and we can observe that respiratory fluctuations are mainly present in the higher frequency range around 0.28Hz in this resting-state fMRI dataset. In addition, the different models were tested on data incorporating a more liberal frequency filtering within the 0.01 − 0.1Hz range, as presented in Supplement V.
In the S1200 release, diffusion MRI data was collected in 6 runs, whereby approximately 90 directions were sampled during each run, using three shells of b = 1000, 2000, and 3000 s/mm 2 , with additionally 6 b = 0 images (Sotiropoulos et al., 2013a). A spin-echo EPI sequence was incorporated with repetition time T R = 5520 ms, echo time T E = 89.5 ms, using a multi-band factor of 3. In total N s = 111 slices were collected, with field of view F OV = 210 mm × 180 mm and an isotropic voxel size of 1.25 mm. The minimal preprocessing pipeline of the HCP includes intensity normalization across runs, EPI distortion correction using the FSL5 "topup" tool, correction of eddy-current induced field inhomogeneities and head motion artifacts using the FSL5 "eddy" tool, and finally includes gradient non-linearity corrections and registration to the structural T1w image (Glasser et al., 2013;Sotiropoulos et al., 2013b;Andersson et al., 2003;Andersson and Sotiropoulos, 2015a,b). More details on the minimal preprocessing of the HCP diffusion MRI are described in Glasser et al. (2013). To reconstruct the anatomical connection strengths between regions within the multi-modal parecellation (Glasser et al., 2016), the MRtrix3 software package was incorporated (Tournier et al., 2019). Multi-shell multi-tissue constrained spherical deconvolution (Jeurissen et al., 2014) was applied to obtain response functions for fiber orientation distribution estimation (Tournier et al., 2004(Tournier et al., , 2007. Then 10 million streamlines were created using anatomical constrained tractography . Finally spherical-deconvolution informed filtering was used (Smith et al., 2013a), reducing the number of streamlines to 1 million. The strength of SC was defined as the number of streamlines connecting two brain regions, normalized by the region volumes. The group structural connectivity matrix A SC was obtained as the average SC across the first 10 subjects, because the variance in the SC strength was relatively low across subjects (Zimmermann et al., 2019), while probabilistic tractography methods are computationally demanding. For the HCP dataset, including only young healthy subjects, the similarity of the SC across subjects was quite high, and the Pearson correlation coefficient between SC values of the 10 subjects was on average 0.91. But when comparing very different subject cohorts, like healthy and diseased subjects, the anatomical connectivity can differ considerably between those cohorts, and the SC matrix should then be computed for every studied group individually.

Graph neural networks
Different brain areas communicate via bioelectrical signals transmitted along neuronal axons and collected by neuronal dendrites. Spatio-temporal GNNs provide a novel possibility to incorporate such a structural scaffold into a graph-based prediction model (Wein et al., 2021b). Due to cognitive information processing in the brain, the spatial interactions of the activity distribution changes dynamically. Spatiotemporal GNNs thus encompasses both the information about the layout of the physical scaffold encoded by the graph structure and the dynamical information about temporal activity correlations. Recently we used a DCRNN architecture to model the spatio-temporal brain dynamics in resting-state fMRI (Wein et al., 2021b). In this study, spatial dependencies of brain activities were modeled via diffusion convolution operations based on the anatomical connectivity and the temporal dynamics of the graph signal were captured in an RNN based model architecture (Li et al., 2018). In our current study we evaluate some alternative spatial and temporal approaches to model dynamics in brain networks. In addition to RNNs, a CNN based architecture for temporal modeling has been introduced by Wu et al. (2019). These authors built upon the WaveNets (van den Oord et al., 2016) and stack dilated causal convolution layers to capture long-range temporal dependencies. Dilated convolutions support exponentially growing receptive fields in deeper layers of the network and allow us to handle long-range temporal sequences efficiently (van den Oord et al., 2016). In addition to the temporal processing based on RNNs and CNNs, we also follow ideas expressed in attention networks and incorporated a relevance score that was computed in temporal attention layers (Vaswani et al., 2017;Zheng et al., 2020).
Based on these temporal approaches, we further study different concepts for representing the spatial dependency between brain regions. First we integrated the SC reconstructed from DTI to represent the anatomical substrate for information propagation in graph convolution operations. Then we additionally incorporated CEs of the structural graph to inherently capture higher order relations between ROIs. Finally we used no predefined spatial layout and treated the spatial connection strengths between ROIs as free parameters. These different spatio-temporal GNN architectures have to the best of our knowledge not been applied yet to analyze the dynamics of brain networks, and in our study we investigate their effectiveness in spatio-temporal modeling of functional MRI.

Preliminaries
Let us represent the brain network as a graph. Every specific brain area or region of interest (ROI), then forms a node in the graph. Let these N ROIs form an graph G = (V, E, A w ) encompassing N vertices, i.e. the meta-voxels or ROIs, and a set E of edges connecting the vertices v n , v n . The graph structure can then be captured in a weighted adjacency matrix A w ∈ R N ×N , whose entries w nn provide the connection strengths between vertices v n and v n , and implicitly define the spatial structure of the graph. As introduced above, in our study we compared three different variants to define the spatial relationship between ROIs. Once we incorporated the SC derived from DTI data as an adjacency matrix A SC , we next employed CE to additionally capture higher order topological features in SC represented by A CE , and finally we treated the spatial relations as adaptive learnable parameters A Adap in the GNN models. The dynamics of the graph signal is then represented by the time-varying neural activity obtained from functional imaging data. Let us first assume that each node of the graph is associated with a single feature represented by the BOLD activity. By considering voxel time series of brain activity maps, then all data can be collected into a data matrix X = (x (1) , . . . , x (T ) ) ∈ R N ×T with x (t) ∈ R N . Given N ROIs, taken from a brain atlas and each represented by a meta-voxel, and considering T time points for each meta-voxel time series, which represents the activation time course of one of the ROIs, then we have, for the BOLD feature represented at node n a related graph signal matrix or BOLD feature matrix: Note that the columns x (m) t ∈ R N of the data matrix describe the activation of all ROIs at any given time point 1 ≤ t ≤ T , while its rowsx (m) n represent the meta-voxel time course of every single ROI 1 ≤ n ≤ N . More generally, if nodes not only represent a single feature m, like the input BOLD signal, but an M -dim feature vector X nt: ∈ R M , then we obtain a feature tensor X ∈ R N ×T ×M , whose frontal, lateral (vertical) and horizontal slices, respectively, read X ::m ∈ R N ×T , X :t: ∈ R N ×M and X n:: ∈ R T ×M .
In addition to above frontal slices X ::m ≡ X (m) of the data tensor X , we thus have the lateral tensor slices: and the horizontal tensor slices: Note that the column fibers of the data tensor X :tm , denoted as x (m) t represent, at every time point t, the distribution of the activity of feature m across all nodes n of the graph. Correspondingly, the row fibers of the tensor X n:m , denoted as x (m) n , represent the time course of every feature m at node n. Finally the tube fibers of the tensor X nt: , denoted as x (n) t , represent the distributions of features at every node n and time point t. This notation will in the following provide the framework to introduce the different techniques to model either dependencies between nodes n, time t or features m.

Spatial dependencies.
Diffusion convolution: In the following we provide a short introduction on a variant of graph convolution denoted as diffusion convolution in the context STGNNs (Li et al., 2018;Wu et al., 2020). The information flow in the underlying graph G = (V, E, A w ) is considered as a stochastic random walk process modeled by a state transition matrix T = D −1 A w = (ŵ 1 . . .ŵ N ) where A w represents a weighted adjacency matrix. The diagonal node degree matrix is given by: whereŵ n = (ŵ 1n . . .ŵ N n ) T ∈ R N ∀ n = 1, . . . , N withŵ nn = w nn / n w nn denoted normalized edge strengths. State transitions were modeled as a diffusion process on an unstructured graph. The former was represented by a random walk Laplacian: where the transition operator T was replaced by its eigen-decomposition with U the matrix of eigenvectors and Λ the diagonal matrix of eigenvalues. Hence, the set of eigenvectors provided an orthogonal basis system for the spatial representation of the brain graph. With the help of these eigenvectors u n the spatial structure of the graph could be implemented. A spectral representation in combination with the convolution theorem then provided a definition of the graph convolution operator G C (Shuman et al., 2013), which served to compute the spatial convolution of the input signal and a spatial filter kernel to yield the output of the -th convolution layer as: Here the approximation resulted from a power series expansion of the convolution kernel with respect to the eigenvalue matrix Λ of the transition operator T (Wein et al., 2021b;Defferrard et al., 2016). Finally considering a CNN architecture and applying the graph convolution operator G C , the filtered input signal y (q) t was transformed with an activation function σ(·) to yield the output h (q) t of each of the q ∈ {1, . . . , Q} graph convolution layers as follows: Hereby x (m) t ∈ R N denotes the m-th input feature component at time t, h (q) t ∈ R N the corresponding output component of the q-th convolution channel, Θ (q) k ∈ R N parameterizes the q-th convolutional kernel of order k and σ(·) denotes any suitable activation function. Note that for deeper convolution layers > 1, = 1, . . . , L, the input to the convolution layer is given by the output component of the convolution layer − 1 instead of the input signal. In summary, these graph convolution layers can learn to represent graph structured data and could be trained with gradient descent based optimization techniques.
Structural connectivity: One possibility to define the spatial layout of the brain network characterized by the weighted adjacency matrix A w is to directly incorporated the structural connection strength as reconstructed from DTI data. The weights w nn in our adjacency matrix would accordingly reflect the number of fibers connecting two brain regions n and n , derived from probabilistic fiber tracking (Tournier et al., 2004). This type of structural adjacency relation is denoted as A SC ∈ R N ×N . The acquisition parameters of the DTI data and the structural connectome generation are outlined in detail in the Dataset section.

Connectome embeddings:
As an alternative to the original SC, connectome embeddings (CEs) can generate node embeddings that capture also higher order topological features of the structural layout (Rosenthal et al., 2018). The idea of such a graph embedding is to represent each node in the graph by a M -dimensional feature vector. This technique is originally inspired by the word2vec algorithm introduced by Mikolov et al. (2013) who proposed a technique to learn vector-valued representations for words in a text which preserve linguistic regularities in their embedding space. Similarly the node2vec algorithm can be used to embed vertices of a graph into a subspace where similar embeddings capture the k-step (k = 1, 2, . . . , K) relation between the vertices and their k-step neighbors (Rosenthal et al., 2018;Grover and Leskovec, 2016). We used this technique to embed each brain region n in the SC graph into a 64-dimensional vector representation. We therefore employed the gensim python package (Řehůřek and Sojka, 2010) using the skip-gram model to learn the node representations (Mikolov et al., 2013). Briefly, in this context the idea of the skip-gram model is to predict from a target node in a network its neighboring nodes, whereby a sequence of neighboring nodes is created by performing a biased random walk on the structural graph (Grover and Leskovec, 2016). To generate the node sequences in total 100 random walks were performed for each node with walk a length of 80 nodes. The return parameter of the random walk was set to p = 2 and the in-out parameter to q = 1. The similarity between the N brain regions in their embedding space was computed using the Pearson correlation coefficient, yielding a connectivity matrix denoted with A CE ∈ R N ×N . As illustrated in figure 3 (B), the embeddings could yield meaningful representations that revealed long-range connections between regions which were not present in the original SC (Rosenthal et al., 2018).
Adaptive adjacency matrix: So far the spatial layout of the brain graph has been represented with the help of the orthogonal eigenbasis system U of the transition operator proportional to the random walk Laplacian. This presupposed a thorough knowledge about the spatial structure of the underlying brain network that entered the related adjacency matrix. Remember that the weights of the adjacency matrix were deduced from DTI measurements based on SC or their CE similarity. However there may exist hidden activity correlations that are not represented in the original adjacency matrix used to construct the random walk Laplacian. Hence, one may wish to introduce an additional self-adaptive, normalized adjacency matrix A Adap ∈ R N ×N (Wu et al., 2019). The latter has been constructed as a matrix of trainable weights V Adap ∈ R N ×N , which were at first initialized as zero and then again optimized via gradient descent (Kingma and Ba, 2014). Inspired by the study of Wu et al. (2019), a normalized self-adaptive adjacency matrix was computed as: The transformation function σ(·) ≡ tanh(·) confined the adaptive weights to the range [−1, 1], which then were normalized by the number of nodes N in the network. This self-adaptive adjacency matrix can help to uncover any hidden, still unknown dependencies between ROIs of a given graph structure. Thus it may extend any graph diffusion convolution layer to yield its output activity as: Note that the normalized self-adaptive adjacency matrix A Adap may be considered as an additional transition operator here. In an attempt to decouple the temporal processing from any underlying spatial layout of the graph connectivity, the first term within parentheses may be skipped and the self-adaptive adjacency matrix may possibly identify the underlying graph structure from the data alone. This may be applicable to situations, where no predefined graph structure is known or involved. The output of the WaveNets: Rather than incorporating diffusion convolution layers into RNNs, dilated causal convolution (DCC) layers (van den Oord et al., 2016) have been instead employed in the GWN architecture (Wu et al., 2019). The full GWN model is illustrated in figure 9. The DCC was defined through a dilated causal convolution operator D C : whereby d denoted the dilation factor and Θ (q) t represented the filter kernel. DCC could be implemented by sliding over the input time series x (m) t while skipping input values while, from layer to layer, increasing step size d · r. This procedure leads to an exponential growth of the receptive field with increasing layer depth as is schematically illustrated in figure 9 (C). The information flow was controlled by a gated temporal convolution network (TCN) as shown in figure 9 (B), which is obtained as: Here tanh(·) denotes the output activation function, and Θ (q) 1 , Θ 1 , and b 1 , b 2 represent the convolution filers and bias terms respectively. Further D C represents the causal convolution operator, the Hadamard product and σ(·) denotes the logistic function, which controls the information passed to the next layer. To achieve large receptive fields, the layers in a WN architecture are organized in blocks, whereby in each block the dilation factor d is doubled with d = 1, 2, 4, . . . up to a certain limit and then repeated in the same manner in the next block (van den Oord et al., 2016). After each such dilated convolution layer a diffusion convolution layer G C (equation 9) is subsequently applied to account for the spatial dependencies, as illustrated in figure 9 (A). Figure 9: An overview of the complete GWN model is provided in (A). The GWN model consists of L layers. For the temporal modeling the GWN applies first the gated TCN mechanism and then for the spatial aspects utilizes graph convolution operations (G C ) in each layer. Each layer additionally incorporates residual connections to stabilize the gradient during learning (He et al., 2016). The information in each layer is combined by using skip connections, and the final predictions are generated by passing the output of the skip connections through two fully connected layers. The gated temporal convolution network (TCN) mechanism (B) applies a dilated causal convolution (D C ) in combination with a tanh(·) and a σ(·) activation function to control the information flow. In (C) the dilated causal convolutions are illustrated. In each layer a temporal convolution is applied whereby the dilation factor can be increased in subsequent layers. These dilations lead to exponentially growing receptive fields for neurons in higher layers. The receptive field of a neuron in layer is highlighted in blue.
Temporal relevance: Yet another approach to solve spatio-temporal time series prediction problems considers attention mechanisms in spatial and temporal domains to capture dynamic correlations (Vaswani et al., 2017;Zheng et al., 2020). In this study we therefore additionally explore non-linear temporal correlations via a temporal relevance mechanism for modeling temporal fluctuations in the BOLD signal. Let the temporal state of the brain network be represented by the multivariate signal tensor X ∈ R N ×T ×M such that the temporal states of any node n be collected in the signal matrix X n:: ≡ X (n) ∈ R T ×M . The activity at any node n and at any time t was then represented by the tube fibersx (n) t ∈ R M , where M denoted the number of features characterizing the node activity. Temporal correlations between different node states could be estimated by filtering the multivariate signals in a cascade of temporal relevance blocks, as illustrated in figure 10. The queries and keys are computed from the input in the -th block at time point t with a simple non-linear transformation g r (x (n) For any node n and any time point t i the relevance of its statesx (n) t i at time points t j < t i with respect to the considered statex (n) t j could then be assessed by computing the inner product between the queries and keys: A normalized temporal relevance scoreδ (n) (t i ,t j ) could then be computed according to: Finally, t j < t i , j ∈ {1, . . . , T p } denoted a set of time steps before time point t i . After computing the temporal relevance scoreδ (n) t i ,t j , the hidden state of node n at time t i could be derived as: Whereby g r (·) denotes a non-linear projection again. Note that all parameters W r and b r to be learned were shared across all nodes and time steps. In total L layers of temporal attention mechanisms were stacked to generate a final prediction for the BOLD signal. After each layer batch normalization was applied and additionally residual connections were incorporated to stabilize the gradient (He et al., 2016).
Figure 10: An overview over the temporal relevance or attention model. The single feature input X (in) , representing the BOLD signal, is first projected by a fully connected layer onto M output features. The temporal relevance scores are computed in each of the L attention layers, and to further account for vanishing gradients, additionally residual connects are incorporated (He et al., 2016). The output of the final layer L is then projected back onto a single feature, representing the predicted neural signal.

Model training
In this section we outline the training procedures which were used for the different neural network models to learn the temporal and spatial dynamics in the BOLD signal. Before training, the fMRI data of each session was linearly scaled between 0 and 1, to get gradients of a small magnitude during the backpropagation learning, what facilitates the fine-tuning of the learning rate. For all models, the mean absolute error (MAE) was used as an objective function to quantify the overall difference between the true BOLD signal x (t) and predicted signalx (t) in all N brain regions: DCRNN The DCRNN model, based on a RNN architecture, was trained with backpropagation through time (Werbos, 1990), with the objective to maximize the likelihood of generating the target timeseries.
To additionally account for a mismatch between training and testing distributions of stimuli, a scheduled sampling strategy was used . The probability of using a true label as a decoder input decayed according to: with τ > 0 the decay parameter and i ∈ N counting the iterations. During supervised learning, instances to be predicted were, of course, known. For this optimization problem, the Adam algorithm (Kingma and Ba, 2014) was employed, and the model was trained for 70 epochs on mini-batches of 16 training samples. To further improve convergence, an annealing learning rate was used, initialized as η = 0.1, and decreased by a factor of 0.1 at epochs 20, 40 and 60, or if the validation error did not improve for more than 10 epochs. Before lowering the learning rate, the weights with lowest validation error were restored, in order to avoid getting stuck in local optima. For the training data set including only 10 subjects used in the Model accuracy and network scaling section, the number of training epochs was increased to 140 and the learning rate decay applied at epochs 40, 80 and 120. The influence of the DCRNN model hyperparameters are discussed in Supplement I (figure S1) and were chosen to yield a reasonable trade-off between accuracy and computational requirements. The encoder and decoder of the sequence-to-sequence architecture consist to two diffusion convolution GRU layers each, and the hidden state size was set to 64. TAtt The TAtt model was trained using the Adam optimizer (Kingma and Ba, 2014) for in total 40 epochs, minimizing the MAE defined in equation 19, using a batch size of 16. The learning rate was initialized with η = 0.1 and decreased by a factor of 0.1 at epochs 10, 20 and 30. The influence of the hyperparmeters is evaluated in Supplement I ( figure S3). The number of neurons in the temporal attentions were set to 32 thereby using 4 attention heads in the 4 TAtt layers. With this setup of hyperparameters one epoch of the TAtt model took around 7.7 minutes.

Vector auto regressive model
Granger causality (Granger, 1969) is currently most often based on linear vector autoregressive (VAR) models for stochastic time series data. Therefore we compare our GNN based approach with a VAR model, as implemented in the multivariate Granger causality (MVGC) toolbox (Barnett and Seth, 2013). An autoregressive process (AR) is based on the idea that a signal x (t) can be described as a linear superposition of the first T p of its lagged values (Luetkepohl, 2005): with coefficients or weights α 1 , . . . , α p , an intercept β and the error term u (t) . This univariate formulation can be then extended to a multivariate VAR model including N time series N ] (Luetkepohl, 2005): whereby the coefficients are now collected in matrices A ∈ R N ×N , and intercepts and errors are characterized by vectors b ∈ R N and u (t) ∈ R N respectively. In our study the multivariate time series x (t) reflect the BOLD signal strength in the N brain regions, sampled at different timesteps t.
To estimate parameters of the VAR model we used the ordinary least squares (OLS) fit provided in the MVGC toolbox (Barnett and Seth, 2013). As outlined in the Data description section we used the first 80% of the data from each fMRI session to fit the model. Then for the comparison to the GNN approaches in the Model accuracy and network scaling section we tested the model order p in steps of 5 with p = 5, 10, . . . , T p and chose the model with highest accuracy individually for each dataset. To check for stationarity of the signals, an augmented Dickey-Fuller test for unit roots was applied to the BOLD timecourses (Hamilton, 1994;Mackinnon, 1994), using a p-value of p < 0.01. For the 25 subject dataset, around 10.0% of the BOLD time courses do not fulfill the stationarity criteria of the augmented Dickey-Fuller test (p > 0.01) when using such a high lag order of T p = 60. But as the objective criterion of the evaluation in the Model accuracy and network scaling section was to assess the capabilities of the models to predict empirically observed neural activity patterns, we chose the VAR model with best prediction accuracy for comparisons with the GNNs.

Data availability
A demo version for MRI data preparation and training the DCRNN model is provided under: https://github.com/simonvino/DCRNN_brain_connectivity. In addition, a demo version for the GWN model is provided under: https://github.com/simonvino/GraphWaveNet_brain_connectivity. Preprocessed HCP data is publicly available under: https://db.humanconnectome.org.

Supplementary Information Supplement I
In this supplement the influence of the model hyperparameters for the different neural network architectures is discussed. The hyperparameters are chosen as described in the Model training section and held constant, while only the hyperparameter of interest is varied in the following evaluations. Figure S1 and S2 show that the performance of he DCRNN and GWN in general still could slightly improve with a larger number of model parameters. Because the computation time and memory requirements linearly grow with the number of parameters, we chose the model hyperparameters as described in the Model training section to yield a good trade-off between model performance and computational requirements. Also the TAtt model in figure S3 shows some improvement with a larger number of parameters, however the MAE is still considerably higher compared to the RNN and WN based GNN architectures.    Table 1: List of ROIs involved in visual processing based on the multi-modal parcellation proposed by Glasser et al. (2016). The table shows the index of the region in the atlas for the right/left hemisphere including the name of the region.

Supplement III
In this supplement the different spatial and temporal modeling approaches are evaluated using the scalefree R-squared (R 2 ) measure and the similarity of predicted connectivity states. The R 2 measure can be obtained with: where x (t) n represents the true BOLD signal in brain region n at timestep t. Correspondinglyx  figure S4. To test the significance also based on this measure, the R 2 values were computed for the individual subjects, and by applying a paired t-test the RNN and WN model both outperformed the TAtt model with a p-value of p ≤ 0.0001 (Cohen's d 1) based on this measure. Also the impact of structural modeling showed to be significant with p ≤ 0.0001 (Cohen's d > 1) for both STGNN models.
In addition, we evaluated the similarity of the predicted FC states to the ground truth FC states. For this purpose we computed a FC matrix A F C ∈ R N ×N from the true ROI timecourses x ∈ R N ×T f , and respectively the predicted FC matrixÂ F C from predicted timecoursesx based on Pearson correlation. Then we sorted the elements below the diagonal of the symmetric FC matrix A F C into a vector v ∈ R 1 2 N (N −1) , and defined a correlation based measure of FC similarity as: withv c denoting the average and σ v the variance of FC values. Based on this metric we replicated the comparison of spatial and temporal GNN approaches in figure S5. Also using this similarity measure, the RNN and WN model outperformed the TAtt model with p ≤ 0.0001 (Cohen's d 1) across subjects. The impact of structural modeling showed to be significant with p ≤ 0.0001 (Cohen's d > 1) across subjects for both STGNN models. In addition the distributions of the test MAE across different subjects, timepoints and brain ROIs are visualized for the STGNNs with and without incorporating spatial modeling in figure S6.  CEs in an adjacency matrix A CE to define spatial relationships, and (D) displays the case when using a self-adaptive weight matrix A Adap . Figure S6: The distribution of the test MAE error across timepoins, subjects and ROIs is illustrated for the STGNNs, with and without incorporating spatial modeling. In (A) the prediction error in dependence of the forecasting horizon is depicted. Further (B) shows the MAE across subjects and brain regions of the RNN, visualized in boxplots on the left side. On the right side the MAEs in dependence of the brain regions were projected onto the cortical surface. The same analysis was carried out for the DCRNN, WN and GWN model in figure (C), (D) and (E) respectively. Figure S7: The figure shows a comparison of the model performances, based on the R 2 measure as defined in equation 23. The errorbars represent the standard deviation of the R 2 values across subjects. In (A) the model accuracies when using a dataset of 50 subjects are shown for the visual network, the network within the right hemisphere and the whole brain network. Figure (B) and (C) show the test performances in dependence of the network size using the 25 and 10 subject dataset, respectively. Except for the single network dataset, the improvements of accuracy of the DCRNN and GWN in comparison to the VAR model became in all cases highly significant with p ≤ 0.0001 (Cohen's d 1). Figure S8: The figure shows a comparison of the model performances, based on the similarity of predicted FC states, as defined in equation 24. The errorbars represent the standard deviation of the values across subjects. In (A) the model accuracies when using a dataset of 50 subjects are shown for the visual network, the network within the right hemisphere and the whole brain network. Figure (B) and (C) show the test performances in dependence of the network size using the 25 and 10 subject dataset, respectively. Except for the single network dataset, the improvements of accuracy of the DCRNN and GWN in comparison to the VAR model became in all cases highly significant with p ≤ 0.0001 (Cohen's d 1). Figure S9: In this figure the prediction accuracy of the different models is presented when using a bandpass filter within the 0.01 − 0.1Hz frequency range for the 25 subject dataset and the brain network including the ROIs within the right hemisphere (Glasser et al., 2016). In (A) the test MAE in dependence of the forecasting horizon is shown, computed as an average across test samples and brain regions.   , and then averaged across all subjects and brain regions. The errorbars represent the standard deviations across subjects. Figure S13: The FC connectivity of V1 to all other regions within the right hemisphere is shown in (A). In (B) and (C) the GWN based directed connectivity is illustrated when using the SC (A SC ) and CE similarity (A CE ) for the information propagation respectively. Figure S14: The distribution of respiratory frequencies is shown in this figure. For the 25 subjects resting-state fMRI dataset, physiological data was available from 22 subjects, and the depicted spectrum was computed as the average across the individual spectra of the 22 subjects.