# Benchmarking Graph Neural Networks
for FMRI analysis

Ahmed El Gazzar
Amsterdam University Medical Center
<EMAIL_ADDRESS>&Rajat Mani Thomas
Amsterdam University Medical Center
<EMAIL_ADDRESS>&Guido Van Wingen
Amsterdam University Medical Center
<EMAIL_ADDRESS>
###### Abstract

Graph Neural Networks (GNNs) have emerged as a powerful tool to learn from
graph-structured data. A paramount example of such data is the brain, which
operates as a network, from the micro-scale of neurons, to the macro-scale of
regions. This organization deemed GNNs a natural tool of choice to model brain
activity, and have consequently attracted a lot of attention in the
neuroimaging community. Yet, the advantage of adopting these models over
conventional methods has not yet been assessed in a systematic way to gauge if
GNNs are capable of leveraging the underlying structure of the data to improve
learning. In this work, we study and evaluate the performance of five popular
GNN architectures in diagnosing major depression disorder and autism spectrum
disorder in two multi-site clinical datasets, and sex classification on the
UKBioBank, from functional brain scans under a general uniform framework. Our
results show that GNNs fail to outperform kernel-based and structure-agnostic
deep learning models, in which 1D CNNs outperform the other methods in all
scenarios. We highlight that creating optimal graph structures for functional
brain data is a major bottleneck hindering the performance of GNNs, where
existing works use arbitrary measures to define the edges resulting in noisy
graphs. We therefore propose to integrate graph diffusion into existing
architectures and show that it can alleviate this problem and improve their
performance. Our results call for increased moderation and rigorous validation
when evaluating graph methods and advocate for more data-centeric approaches
in developing GNNs for functional neuroimaging applications.111Our code is
publicly available at https://github.com/elgazzarr/fMRI-GNNs

_Keywords_ graph neural networks $\cdot$ functional connectivity $\cdot$ major
depression disorder $\cdot$ autism spectrum disorder $\cdot$ functional
magnetic resonance imaging

##  1 Introduction

Mapping the underlying functional brain activity to observed phenotypes and
neurological disorders has been a major goal of neuroimaging research. The
emergence of non-invasive imaging techniques such as functional magnetic
resonance imaging (fMRI) that estimates the neural activity of the brain by
measuring the blood-oxygen level dependence (BOLD) has been paramount in
advancing our understanding of the functional organization of the
brain(Huettel et al., 2004). Yet, modeling the fMRI signal is a challenging
feat due to the high dimensionality of the signal coupled with the low sample
size, especially in clinical datasets, and the low signal-to-noise ratio of
the scans. To build learning models from such data, researchers must opt for
tools that incorporate priors from the data structure to guide modeling. This
is also called finding the correct inductive bias. Given our knowledge that
the brain functions a network Sporns et al. (2005), graph-based methods have
been instrumental in modeling the brain with the field of network neuroscience
utilizing graph theory approaches to find interesting insights about the brain
function and dysfunction (Sporns, 2018). With the rise of deep learning, GNNs
have been developed to extend the merits of neural networks to graph
structured data (Bronstein et al., 2017). Naturally, the neuroimaging
community has picked up GNNs and applied it to wide range of tasks such as
phenotype prediction (Ktena et al., (2017; Parisot et al., (2018), disease
diagnosis (Gallo et al., (2021; Li et al., 2020) and task inference (Li et
al., 2019). This adoption is still in its infancy, where most of the work is
focused on adapting a specific architecture proposed on other domains for a
specific neuroimaging application. It therefore remains unclear which graph
architecture should be selected. Should we use graph convolution (Kipf and
Welling, 2016) as (Gallo et al., (2021) or Graph isomorphic networks (Xu et
al., 2018) as in (Kim and Ye, 2020). Should we use graph attention (Velic´ et
al., 2017) as in (Hu et al., 2021). Should the input graph be static or
dynamic (El-Gazzar et al., 2021)? How to create the edges? What is the sample
size suitable for training these models? And most importantly, are these
techniques capable of exploiting the underlying structure of the data to
improve learning?
Despite the theoretical appeal of the application of GNNs to learn from
functional brain activity, we find that similar to others domains, the
experimental settings for GNNs in neuroimaging applications are in many cases
ambiguous or not reproducible (Errica et al., 2019). Some of the most common
reproducibility problems concern hyperparameters selection, where
hyperparameters optimization is conducted for the proposed model and not the
baselines, or they are optimized based on their performance on the final test
set. Another main issue is the incorrect use of data splits for model
selection versus model assessment. These limitations could lead to over-
optimistic and biased estimates of the true performance of the model, making
it difficult to assess the advantage of adopting GNNs for this application or
to identify the key components in developing more powerful architectures.
With this premise, our primary aim is to provide the neuroimaging community
with a fair performance comparison among some of the popular GNN approaches
recently adopted for classification of phenotypes or psychiatric disorders
from resting-state fMRI scans. We conducted a large number of experiments
within a standardized and rigorous model selection and assessment framework
for three different large-scale dataset encompassing three different
classification tasks. Namely, the Rest-meta-mdd (Yan et al., 2019) dataset for
the classification of major depression disorder, the Autism Brain Imaging Data
Exchange (ABIDE I+II) (Martino et al., 2014) dataset for the classification of
autism spectrum disorder and the UkbioBank (Sudlow et al., 2015) dataset for
sex classification. Secondly, we aim to identify if and to what extent current
approaches can effectively exploit graph structures. To this end, we implement
three structure-agnostic baselines, whose purpose is to disentangle the
advantage of graph structure information from the nodes features. To ensure a
fair comparison against GNN methods, we conduct our baseline experiments under
a similar framework for hyperparameter optimization and model selection and
assessment. Our third aim is to study the effect of scaling the sample size on
the GNNs and the baselines. We evaluate the performance of the models when
trained on different sample sizes on the UkbioBank for sex classification and
tested on an independent subset. Finally and most importantly moving forward,
our objective is to apply insights from our implementations and empirical
results to highlight strengths and shortcomings of current methods and propose
a set of recommendations to help utilize the potential of this powerful set of
tools.

##  2 Background

Graph neural networks are a class of neural networks that extend deep learning
to graph structured data, and have been utilized in several domains for tasks
such as node classification, edge prediction and graph classification
(Bronstein et al., 2017; Zhou et al., 2020; Wu et al., 2019). While graph
neural networks have been applied for node classification tasks in functional
neuroimaging where the graph represented a population graph (Parisot et al.,
(2018), a more befitting application is graph classification, where a brain is
represented as a mathematical graph $G=(V,E,X)$ whose nodes $v_{i}\in V$
represent brain regions and edges $E$ constitute their connections. Input
graphs for GNN models are often attributed graphs with the attributes
representing the node features $x_{i}\in X$. A graph can be stationary or
dynamic depending on the evolution of the features and/or connections across
time. This work addresses the task of graph classification from both static
and dynamic functional graphs and using five popular GNN architectures. Our
criteria for selection are i) popularity i.e. the architecture has been
extensively validated in multiple domains for applications and is commonly
used as a benchmark in graph neural networks surveys; ii) strong architectural
differences; iii) original code publicly available; iv) previous adaptation on
fMRI data. Specifically, we adopt Graph Convolutional Networks (GCNs), Graph
Attention Networks (GATs), Graph Isomorphic Networks (GINs) as architectures
to learn from static graphs. While for dynamic graphs, we adopt Spatio-
Temporal Graph Convolutional Networks (ST-GCNs) and Adaptive Spatio-Temporal
Graph Convolutional Networks (AST-GCNs). Several architectural choices for
each model are represented as parameters in a search space and are selected
using an exhaustive hyperparameter search. Thus each of the models can be
viewed as a general class that encompasses several variants of previously
proposed GNN models. We briefly describe the main idea behind each class and
some of their previous applications for learning from fMRI data. We discuss
the specific details of each method under a uniform framework in Section 3.2.

##### Graph convolutional networks

A GCN (Kipf and Welling, 2016) can be viewed as the counterpart of a
convolutional neural network (CNN) model for graph-structured data that uses a
graph spectral approach to perform convolution. A graph convolutional layer
extracts features through a first-order spectral approximation of the graph by
restricting the filters (limiting the order of the Chebyshev polynomial) to
operate in the neighborhood at one step away from each node. GCNs are the most
popular method applied for graph classification from fMRI data and have been
utilized for the diagnosis of psychiatric disorders such as ASD (Arya et al.,
2020; Li et al., 2020; Ktena et al., (2017), MDD (Gallo et al., (2021; Qin et
al., 2022) and Schizophrenia (Oh et al., 2022; Lei et al., 2022).

##### Graph attention networks

GAT (Velic´ et al., 2017) is a graph neural network architecture that uses the
attention mechanism to learn weights between connected nodes. In contrast to
GCN, which uses predetermined weights for the neighbors of a node
corresponding to the normalization coefficients, GAT modifies the aggregation
process of GCN by learning the strength of the connection between neighboring
nodes through the attention mechanism on the features of the connected nodes.
(Li et al., 2019, 2020) have adopted GATs for the diagnosis of ASD.

##### Graph isomorphic network

GINs were proposed by (Xu et al., 2018) as a special case of spatial GNN
suitable for graph classification tasks. The network uses the sum operation to
aggregate features from the neighboring nodes and employs an extra fully
connected layer in a non-linear multi-layer perceptron (MLP) to transform the
output at each layer. The authors argue the GINs can learn an invective
function which deems the model to be possibly as powerful as the Weisfeiler-
Lehman test (Weisfeiler and Leman, 1968) for graph classification tasks. (Kim
and Ye, 2020)) have applied GINs for sex classification from rs-fMRI data.

##### Spatio-temporal graph convolution network

ST-GCNs (Yu et al., 2017) has been first proposed to model dynamic data for
the application of traffic forecasting. They were developed to extract
features from graphs with dynamic attributes that evolve with time and have
been popular in several applications such as weather forecasting and skeleton-
based activity recognition. A typical ST-GCN block consists of a temporal
layer (e.g. 1D Conv layer, LSTM, etc.) to extract temporal features followed
by a graph convolution layer for spatial feature extraction. Stacking ST-GCN
blocks enables spatio-temporal feature extraction.(Gadgil et al., 2020) and
(Azevedo et al., 2020) have applied ST-GCNs on fMRI data the tasks of sex and
age classification, and (Kong et al., 2021) for the task of MDD diagnosis.

##### Adaptive spatio-temporal graph convolution network

AST-CGNs (Bai et al., 2020) fundamentally share the same concept as ST-GCNs in
feature extraction however they differ in the input graph structure. While ST-
GCN relies on a predefined adjacency matrix, AST-GCN starts with a randomly
initialized graph structure and adaptively learns the optimal graph structure
through gradient descent in an end-to-end fashion. While computationally more
expensive, this setup enables learning task-specific graph structure in
domains where there is no known graph structure or the predefined graph
structure is not sufficient. (El-Gazzar et al., 2021) have applied AST-GCN for
the tasks of sex and age classification.

##  3 Methods

###  3.1 Graph creation

This work studies both static and dynamic representation of graphs. Static
graphs assume no temporal evolution of the graph structure and nodes features.
They are typically created solely using static measures of similarity (e.g.
correlations) for both the edges and features. Currently, the most popular
approach for static graph creation is using the Pearson correlation of the
ROIs time-courses (commonly refereed to as functional connectivity) to define
the nodes features as the connectivity profile of each ROI (the respective row
of the correlation matrix) and the graph edges are represented as a sparse
version of the correlation matrix. Note that while some works use other
approaches such as clustering, structural connectivity and anatomical location
for representing graph edges, the use of functional connectivity (FC) to
define graph structure remains the most popular approach in the large majority
of GNN brain studies. In this work we follow the same popular strategy to
create static graphs and use GCN, GAT, and GIN models for learning from static
graphs while treating the sparsfication threshold as a tunable hyperparamter.
The threshold is set based on a proportional thresholding (i.e. selecting
strongest PT % connections) rather than absolute threshold as it has shown in
graph to have improved test-retest reliability and is less sensitive to
variation in graph density (Braun et al., 2012; Van Wijk et al., 2010). Static
graphs are relatively easier to compute, have more training architectures to
select from, have a significantly lower number of dimensions and a more robust
to physiological noise and inter-scanner variability. However, they fail to
represent the dynamic nature of brain activity which is hypothesized to
contain the biomarkers for mental disorders (Hutchison et al., (2013).

Dynamic graphs enable modelling such dynamic signal by representing the nodes
features as ROIs time-courses. Similar to static graphs, dynamic GNNs employ
thresholded FC to define the graph adjacency matrices. Another strategy
adopted by (El-Gazzar et al., 2021), is to randomly initialize the adjacency
matrix and adaptively learn its weights in an end-to-end manner simultaneous
with the GNN weights via gradient descent. We adopt the former strategy for
our ST-GCN model and the latter for the AST-GCN model. Figure 1 summarizes the
graph creation strategy employed for all 5 models in this study.

Figure 1:  Graph creation strategy for the 5 models employed in this study.
Input graphs to GCNs, GATs and GINs solely rely on static functional
connectivity to create the nodes features and graph edges. Input graphs to ST-
GCNs and AST-GCNs utilize ROIs timecourses to define the nodes features. ST-
GCNs graph use thresholded FC to define the graph structure, where AST-GCNs
start with a randomly initialized graph and learn the edge weights via
gradient descent.

###  3.2 Graph Neural Networks

For graph classification, the goal of a GNN is to learn a function $f(.)$
which maps a graph to a feature vector that summarizes the graph and contains
information necessary for classification. The forward function thus can be
viewed as :

$\centering y=g(f(G))\@add@centering$ (1)

where $y\in R$ is the label of the graph and $g(.)$ is a function
parameterized by a neural network to map the output feature vector of the GNN
to the label. In this work, we choose $g(.)$ as a single-layer fully connected
neural network for all the different GNN architectures. Thus, the empirical
results would only reflect how powerful are different parameterizations of
$f(.)$ in representing the graph signal and structure. In this section, we
explain how each GNN architecture models this function.

####  3.2.1 GCNs

A GCN learns graph representations via an approximation of spectral
convolutions on the graphs. Spectral convolutions on graphs is defined as the
multiplication of the the graph signal $x$ with a filter $s_{\theta}$ in the
Fourier domain.

$s_{\theta}\ast x=Us_{\theta}U^{T}x$ (2)

where $U$ is the matrix of eigenvectors of the normalized graph Laplacian
$L=I-D^{-1/2}AD^{-1/2}=U\Lambda U^{T}$ with a diagonal matrix of its
eigenvalues $\Lambda$ and $U^{T}x$ being the graph Fourier transform of $x$.
Evaluating this spectral convolution in its current form is computationally
expensive because of the multiplication by $U$ and the eigendecomposition of
$L$. (Hammond et al., 2011) suggested that $s_{\theta}$ can be well-
approximated by a truncated expansion in terms of Chebyshev polynomials up to
K-th order. Follow up work by (Kipf and Welling, 2016) proposed to restrict
the polynomials to their first order for each graph convolutional layer and
stacking $K$ graph convolutional layer to recover a rich class of
convolutional filters. The layer-wise propagation rule of the GCN can thus be
defined in the spatial domain as :

$H^{k+1}=\sigma(\hat{D}^{-1/2}\hat{A}\hat{D}^{-1/2}H^{k}W^{k})$ (3)

where $H^{k}$ is the graph signal representation at layer $k$ ; $H(0)=X$,
$\hat{A}$ is the adjacency matrix with self-connections, $\hat{D}$ is the
degree matrix, $W(k)$ is a trainable weight matrix and $\sigma$ denotes a non-
linear activation function. (Gilmer et al., 2017) later formalized this
definition as message passing and update operations on the nodes. The node-
wise update rule for a GCN can thus be defined as

$h^{k+1}_{i}=\phi(h^{k}_{i},\sum_{j\in\mathcal{N}_{i}}c_{ij}\psi(h^{k}_{j}))$
(4)

where $\psi$ and $\phi$ are the message and update functions respectively,
$c_{ij}$ is the edge weight between nodes $i$ and $j$ (set to 1 in unweighted
graphs) and $\mathcal{N}_{i}$ is the neighborhood of the node $i$. In this
work, we set $K=2$ as we observed no additional gain from convolving higher-
order neighborhoods of the graph nodes. $\psi$ is a linear fully connected
layer $W$ and $\phi$ is a is a single-layer feed forward network with a hidden
dimension $d_{h}$ set as a hyperparameter with a non-linear ReLU(·) activation
function.

####  3.2.2 GATs

The main difference between a GCN and GAT layer lies in the permutation
invariant update function $\phi(h_{i},h_{\mathcal{N}_{i}})$ applied to each
node i to combine the node features hi with its aggregated neighborhood
features $h_{\mathcal{N}_{i}}$ (transformed by some function $\psi$) and
update its representation. In GCNs, importance of node j to node i’s
representation is defined using $c_{ij}$, the entry $\hat{A}_{ij}$on the
normalized adjacency matrix. For GATs, the importance is learned implicitly
using a self-attention mechanism. Formally, the update rule for each node in a
GAT layer can be defined as:

$\displaystyle
h^{k+1}_{i}=\phi(h^{k}_{i},\sum_{j\in\mathcal{N}_{i}}a(h^{k}_{i},h^{k}_{j})\psi(h^{k}_{j}))$
(5) $\displaystyle
a(h^{k}_{i},h^{k}_{j})=\textit{Softmax}(\sigma(m^{T}[Wh^{k}_{i}||Wh^{k}_{j}]))$
(6)

$||$ denotes a concatenation operation, $\sigma$ is a non-linear activation
function and $m$ is a trainable single-layer feed-forward neural network. This
attention mechanism can be extended to a multi-head attention via learning M
different copies of $a(h^{k}_{i},h^{k}_{j})$ and calculating their average. In
this work, we treat $M$ as a hyperparameter with an option of setting $M=1$
given that the sample size is limited and increasing the model trainable
parameters by a large factor can lead to over-fitting. We follow the author’s
recommendation of choosing $\sigma$ as LeakyReLU( ) and parameterizing the
linear transformation of the feature vectors via a linear layer. For the rest
of the parameters of the GAT model (Number of layers, size of the hidden
dimension and graph readout layer), we follow the same setting as described in
Section 4.1.1

####  3.2.3 GINs

The Graph Isomorphism Network (GIN) is a variant of the GNN suitable for graph
classification tasks, which is known to be as powerful as the WL-test under
certain assumptions of injectivity. Unlike GCNs and GATs which can use any
permutation invariant function to aggregate neighborhood features, GINs
typically only use the sum operation and do not include a transformation
function $\psi$ before combining the neighborhood features with the original
node features. Rather, a multi-layer perceptron (MLP) is applied directly on
the combined features. The update rule for each node in a GIN layer be defined
as:

$h^{k}_{i+1}=\textit{MLP}((1+\epsilon)h^{k}_{i}+\sum_{j\in\mathcal{N}_{i}}h^{k}_{j})$
(7)

where $\epsilon$ is a learnable parameter initialized to zero. The purpose of
this parameter is to learn the amount of original information in the node that
it should retain in the next representation. We use a two layer $MLP$ with
ReLU() as non-linearity and as recommended by the paper. For the READOUT
layer, we experiment with the three permutation invariant functions as in
Section 4.1.1.

####  3.2.4 ST-GCNs

ST-GCNs model graphs with dynamic features by utilizing spatio-temporal blocks
that alternate between sequence-to-sequence layers for learning temporal
features, and graph convolutions for spatial feature extraction. While there
has been several variants of ST-GCNs that differ either in the number or order
of layers in the ST-GCN block, or the the temporal learning layer (e.g., gated
temporal convolutions, dilated convolutions, LSTMs), we follow the work by (Wu
et al., 2019) to implement our variant of ST-GCN for fMRI as we have observed
it provided the most stable results across tasks. Namely, we utilize a cascade
of $M$ St-GCN blocks for spatio-temporal feature extraction, where each ST-GCN
block consists of a one layer of gated temporal convolution (GATED-TCN) for
temporal feature extraction followed by a spatial graph convolution as
described in Section for spatial feature extraction. For a node input
representation $x_{i}\in R^{T\times C}$ where $C$ is the number of channels, a
Gated-TCN extracts temporal feature representation $h_{i}$ as follows:

$h_{i}=\textit{tanh}(W_{a}\ast x_{i}+b)\odot\textit{Sigmoid}(W_{b}\ast
x_{i}+c)$ (8)

where $W_{a}$ and $W_{b}$ are the weights of two parallel 1D convolutions for
controlling the flow of information into the next layer, $b$ and $c$ are the
biases and $\odot$ represent the Hardmard product operation. The extracted
feature representations are then treated as the nodes input representation in
a GCN and is processed similar to Equation 4. Another key difference between
the ST-GCN employed this work and prior implementations of ST-GCNs for
modelling fMRI data is that we use the entire timecourses for the duration of
the scan as the node features as opposed to previous works by (Gadgil et al.,
2020) (Azevedo et al., 2020) who use crops of $t=30-60s$. We found that
cropping the timecourses for training and aggregating the predictions during
inference drastically deteriorates the performance for the clinical tasks.

####  3.2.5 AST-GCNs

In multiple spatio-temporal datasets, we know that the underlying interaction
between the temporal features can be modelled as a graph, but we don’t have
access to the ground truth structure of the graph. Adaptive learning of graph
structure is a growing area of research in the field of GNNs and several
variants have been proposed by the research community in recent years (Bai et
al., 2020), (Wu et al., 2019), (Wu et al., 2020). AST-GCNs can be generalized
as ST-GCNs equipped with a graph structure learning module. In this work, we
follow the implementation of (Bai et al., 2020) which defines the normalized
adjacency matrix as

$\hat{A}=I_{N}+\textit{Softmax}(\textit{ReLU}(E\cdot E^{T}))$ (9)

where $E\in R^{N\times d_{e}}$ is a learnable node embedding dictionary
initialized randomly with dimension $d_{e}$ as a hyperparameter. Multiplying
the dictionary by its transpose results in a dense symmetric $N\times N$
matrix, which is then passed through a $ReLU$ activation function to eliminate
weak connections and then row-normalized using the Softmax operation. We also
experimented with the Sparsemax operation (Martins and Astudillo, 2016) and
found that it improved stability during learning. The rest of our AST-GCN
implementation is kept similar to our ST-GCN to disentangle the effect of
adaptive learning of the graph structure.

###  3.3 Graph Readout

Following feature extraction using any of the proposed GNN architectures, the
readout phase computes a feature vector for the whole graph using some
permutation-invariant readout function to be passed as an input to the final
classification layer. In this work we experiment with 3 different permutation-
invariant functions to summarize the graph features after the last graph
layer. Namely Mean(), Mean()$||$Max() and Sum(). The selection of which
function to use for each architecture/dataset combination is set as a
hyperparameter. We empirically observe that Mean() and Mean()$||$Max() provide
the most stable results to all GNN architectures with exception of the GIN
architecture where the Sum() operation marginally improve the performance of
the model.

##  4 Data & Experimental Setup

###  4.1 Datasets

Throughout this study we evaluate the performance of the models on three
different datasets, each with a different classification target and a
different prepossessing pipeline. This setup enables us to make a generalized
conclusion from the empirical results that are robust to variation in
demographics, classification target, pre-processing pipeline and data
acquisition protocols.

ABIDE I+II contains a collection of rs-fMRI brain images aggregated across 29
institutions. It includes data from participants with a clinical diagnosis of
autism spectrum disorder and typically developing participants (TD). In this
study, we select a subset of the dataset with N = 1207 (558 TD/649 ASD) using
the largest nine sites in terms of number of participants (N$\geq$50) to
conduct our analysis. The 4D volumes were pre-processed according to the CPAC
pipeline (Craddock et al., (2013b) with the same parameters as in (Craddock et
al., (2013a). Further, the craddock-200 (Craddock et al., 2012) atlas was used
to parcellate the volume into 200 regions of interest (ROIs) and the
timecourses of each region were calculated as the mean signal between all the
voxels of the ROI.

Rest-meta-MDD is currently the largest open-source rs-fMRI database for
studying major depression disorder including clinically diagnosed patients and
healthy controls from 25 cohorts in China. In this study we select a subset of
the dataset with N = 1453( 628 HC/825 MDD) using the largest 9 sites (N $\geq$
40). Standard pre-processing of the data was done at each site using the Data
Processing Assistant for Resting-State fMRI (DPARSF) (Yan and Zang, 2010) and
112 timecourses of cortical and sub-cortical regions as defined by the
Harvard-Oxford atlas (Makris et al., (2006) were extracted.

UkBioBank is a large-scale biomedical database and research resource,
containing in-depth genetic and health information from half a million UK
participants. In this work we use a subset with (N=5500 (2750 M/ 2750 F))
randomly sampled from the first release. The dataset provides the pre-
processed 4D volumes in MNI format. These volumes were then parcellated to 116
anatomical regions using the automated anatomical labeling (AAL) atlas
(Tzourio-Mazoyer et al., (2002) and the mean timecourse of each region was
extracted.

###  4.2 Model selection & assessment

For a fair evaluation of the methods, we constructed a systematic pipeline to
find the optimal performance of each of the methods. This constitutes using a
configuration of hyperparameters for each model that are optimal for the given
task. A common pitfall in hyperparameter tuning is conducting the search based
on the performance of the final test set. This results in an over-optimistic
estimate of the model generalization capability. In this work, we conducted
hyperparameter search for each GNN model on a held out validation set for each
dataset. We conducted an exhaustive hyperparameter search on the key
parameters that could influence the convergence of the model and its
performance. While some of the parameters are exclusive to some architectures
(e.g. the number of attention heads in the GAT), the majority of
hyperparameters and their search spaces were shared across methods (See Table
S1 for details on the hyperparameters and their search spaces). This enables a
relatively fair evaluation of each architecture. After finding the best
hyperparameter configuration for each model and the respective tasks, we
followed a 5-fold cross validation scheme with model selection using the
lowest loss value on an inner validation split (85/15 training/validation
split) for the Rest-meta-MDD and ABIDE I+II datasets. For the UkBioBank
dataset, the models are trained on subsets of the full sample (N = [500, 1000,
2000, 5000]) with a similar inner validation split for model selection and the
test results are reported on a fixed independent set of N = 500 to evaluate
the effect of training sample size on the performance. (See Supplementary
information for details about training procedure)

###  4.3 Non-graph baselines

To fairly judge whether GNNs are fully capable of leveraging the graph
structure, their performance need to be evaluated against baselines that are
not inductively aware of the structure of the data, and a similar framework
should be followed for hyperparameter search and model selection. In this work
we employ three different structure-agnostic baselines. Specifically, we adopt
a multi layer perceptron (MLP), support vector machines with radial basis
function kernel (SVM-rbf) as baselines for learning from static
representations. The input to these models are the flattened lower-triangular
matrix of the correlation matrix. For a baseline that learns from dynamic
representations, we used 1D CNNs similar to (El Gazzar et al., 2019). The
inputs to the 1D CNN are the prepossessed ROI time-courses stacked as
channels. (See Supplementary information for details about training the
baselines and their hyperparameters)

Table 1: 5-fold cross validation results (balanced accuracy, sensitivity, specificity) for GNNs and baselines on the two clinical tasks of MDD and ASD diagnosis. |  | Rest-meta-MDD | ABIDE I+II
---|---|---|---
| Models | Acc. % | Sens. % | Spec. % | Acc. % | Sens. % | Spec. %
| GCN | 62.6 ± 3 | 57.6 ± 8 | 67.8 ± 8 | 67.6 ± 2 | 69.5 ± 4 | 65.7 ± 3
| GIN | 62.7 ± 3 | 67.3 ± 5 | 58.1 ± 3 | 68.5 ± 1 | 81.7 ± 3 | 55.4 ± 4
Static GNNs | GAT | 62.9 ± 2 | 69.5 ± 3 | 56.4 ± 5 | 68.2 ± 2 | 63.5 ± 5 | 72.9 ± 3
| ST-GCN | 58.2 ± 4 | 48.6 ± 8 | 67.8 ± 9 | 65.3 ± 2 | 67.2 ± 4 | 63.4 ± 4
Dynamic GNNs | AST-GCN | 60.7 ± 2 | 44.6 ± 2 | 75.8 ± 8 | 67.8 ± 2 | 70.8 ± 5 | 64.9 ± 3
| 1D CNN | 63.2 ± 3 | 68.5 ± 5 | 59.6 ± 6 | 70.6 ± 2 | 72.7 ± 3 | 68.5 ± 3
| SVM-rbf | 61.9 ± 2 | 79.0 ± 5 | 44.9 ± 5 | 69.5 ± 2 | 75.9 ± 3 | 63.0 ± 3
Baselines | MLP | 58.4 ± 4 | 60.1 ± 2 | 56.8 ± 3 | 66.6 ± 2 | 72.1 ± 1 | 61.2 ± 2

##  5 Results

###  5.1 Disorder diagnosis

The results in Table 1 show the test performance of the GNNs and baselines on
the clinical datasets. On the clinical tasks (N=1207 and N=1453), all the GNNs
perform competitively where no one architecture significantly superior to the
others. A notable observation is that static GNNs marginally outperformed
dynamic GNNs. This contradicts the main hypothesis that dynamic
representations are more faithful to natural dynamics of brain activity over
summarized static correlation values and thus would result in better results.
This is however true for the dynamic baseline, where the 1D CNN outperformed
all the methods on both tasks. In general, graph-agnostic baselines performed
on par with or marginally outperformed GNN methods despite having no awareness
of the structural organization of the data. These empirical results raise even
more questions; do GNNs simply fail to leverage the graph structure in
learning better representations and hence obtain superior results, or do the
heterogeneity of the data and possible noise in the clinical labels create an
upper bound on the performance of the models, rendering quantitative
evaluation of the merits of GNNs challenging? How does the sample size affect
the performance of these models?

###  5.2 Scaling performance on the UKBioBank

To shed light on some of these questions, we can look at the scaling
performance of these models on the UkBioBank at different training samples on
the task of sex classification. Figure 2 shows that, similar to the clinical
datasets, 1D CNN outperform the rest of the methods at all samples instances.
This is especially clear on the smaller sample sizes (N=500 and N=1000). This
gap is then reduced at the larger sample sizes, where the dynamic AST-GCN
model show improved scaling performance and eventually catch up to the 1DCNN.
Again, while these results emphasize the benefits of incorporating the dynamic
signal, it questions the sample-efficiency of GNN models in clinical
applications, especially since the motivation behind adopting GNNs is that
incorporating inductive bias about the structure of the data should alleviate
the need for larger sample sizes to learn this solely from data. An advantage
we unfortunately do not observe in our empirical results.

A key property of AST-GCN is that it does not require specifying the graph
structure in advance, and is jointly learned using gradient descent along with
the rest of the model at the cost of extra training parameters. The scaling
behaviour of AST-GCN is thus empirically justified and similar scaling
performance was previously reported for adaptive spatio-temporal GNNs (El-
Gazzar et al., 2021; Wu et al., 2020). The performance of AST-GCN however,
attest against the correctness of the graph structure of the data used as
input to the rest of the GNN models. As we highlighted in Section 3.1, unlike
several graph-structured data, there is no ground truth network structure for
the brain and thus defining adjacency matrix of fMRI graphs is rather an
arbitrary process with most works opting for FC for this role. In this work,
we followed a similar strategy in effort to evaluate the performance of GNN
models. The results of our experiments invites us to rethink about this aspect
and to find more principled approaches to accommodate this problem. The
results of AST-GCN guide toward end-to-end learning of the graph structure as
a potential solution to this problem. This however, comes at the expense of
large training samples, a luxury not available in most of clinical datasets.

Figure 2: Scaling performance of GNNs and Baselines on the task of sex
classification at different available training samples of the UkBioBank
dataset. Figure 3: Average test accuracy of the GCN model on the two clinical
datasets at different sparsfication thresholds for creating the graph with and
without graph diffusion.

###  5.3 Graph diffusion improves learning on noisy graphs

In network neuroscience, data-driven analysis on thresholded connectivity
graphs are acknowledged to be potentially problematic, where different
thresholding techniques and values and can lead to the inclusion of spurious
connections confounding the analysis (van den Heuvel et al., 2017; Van Wijk et
al., 2010; Fornito et al., 2012). Yet we argue that this remains an overlooked
problem when using GNNs. During the course of this work, we treated the
sparsfication proportional threshold as a hyperparameter that is selected
based on the best performance on a validation set. While conducting the
experiments, we observed a high sensitivity of the results to the sampled
sparsfication threshold. Further, we observed different optimal thresholds for
each task and GNN class. In Figure 3, we present the results of the GCN model
on the two clinical datasets with all parameters fixed and using sparsfication
thresholds in the range of [50, 95] % with a step size of 5 %. The empirical
results show that selecting an optimal threshold plays an important role in
the final evaluation performance and attest against arbitrary threshold
selection (e.g. 50 % as is the current practice) in creating input graphs for
GNNs. Further, improved performance at higher thresholds (70-85%) highlights
that proportional thresholding of FC results in noisy graphs. In graph
theoretical analysis, proposals to mitigate this problem are to include graph
density as a covariate or to adjust the null distribution to account for
random edges resulting from thresholding low density graphs (van den Heuvel et
al., 2017). An open question is how to deal with such noisy FC graphs while
training GNNs on fMRI data.

A possible solution to dealing with noisy graphs is graph diffusion
convolution (GDC) (Klicpera et al., 2019). Where instead of aggregating
information only from the first-hop neighbors, GDC models message passing
between nodes as a generalized continuous diffusion process defined by the
diffusion matrix

$S=\sum_{k=0}^{\infty}\theta_{k}T^{k}$ (10)

with the weighting coefficients $\theta_{k}$, and the generalized transition
matrix $T$. In principle, one can choose different methods to define the
coefficients and the transition matrix provided that Eq. 10 converges.
Diffusion coefficients can be defined using well-studied functions such as the
heat kernel (Kondor and Lafferty, 2002), personalized page rank (PPR) (Brin
and Page, 1998) or even learned using gradient descent. Popular examples of
defining the transition matrix are the random walk transition matrix
$T_{rw}=AD^{-1}$ and the symmetric transition matrix
$T_{sym}=D^{1/2}AD^{-1/2}$. Essentially by restricting the sum in Eq. 10 to a
finite number of steps $K$, and choosing suitable diffusion parameters, GDC
transforms the original potentially noisy adjacency matrix $A$ into a new
graph $S$ where any GNN model can be applied. This transformation has been
shown to consistently improve the performance of GNN models on several real
world graph datasets. In Fig 3, we show the effect of applying graph diffusion
using the heat kernel coefficient and the symmetric transition matrix with
$K=2$ on the GCN model on the two clinical datasets at different
sparsification thresholds compared to using the original defined graph
structure without any transformation (i.e. setting $K=1$, $\theta_{k}=1$ and
$T=T_{sym}$). The results show that by simply incorporating this
transformation as a prepossessing step, the performance of the model is less
sensitive to the sparsfication threshold at lower thresholds and can even
improve the empirical results without the need for expensive hyperparameter
search for the optimal sparsfication threshold.

##  6 Discussion

The rapid pace of adoption of GNNs in neuroimaging applications urges us to
pause and ponder about the rationale and success of this adoption. While we
agree that the motivation is well-grounded in the field and in fact necessary
for sample-efficient factual learning models, we question whether the current
methods are living up to the promises of the motivation. We argue that it is
difficult to answer this question without a systematic evaluation of the
methods and their current alternatives under a uniform fair setup. In this
study, our aim was to build such a framework by benchmarking popular classes
of GNN architectures through the applications of brain disorder diagnosis and
phenotype prediction. The choice of the application is motivated by the trend
and interest of the researchers to use GNNs for such tasks and the open-source
nature of the datasets. In building our framework, we strive for a rigorous
fair and setup through i) uniform model selection and assessment strategies
across methods, ii) fair allocation of computational resources to all the
methods by providing an exhaustive range of parameter search spaces for each
architecture/task combination to obtain the best possible performance for each
method iii) uniform training setup (e.g. loss functions, regularization,
initialization schemes) for all deep learning methods to disentangle only
architectural effects on the performance. Our empirical results show that, on
average, GNNs fail to outperform current learning alternatives, and that and
no one GNN architecture was consistently superior over the others. A simple 2
layer 1D CNN where all the ROIs stacked as channels with no structural
awareness of the graph consistently outperformed sophisticated structurally
aware GNN methods. The results of AST-GCN – a dynamic model that learns the
graph structure at the cost of more training parameters – showed an improved
performance at the largest sample size over FC-driven graphs. This opens the
question whether graph structures defined using thresholded functional
connectivity, as is done throughout this work and as is the current practice
in the field, is the most optimal strategy for constructing the graphs. We
highlighted that thresholded FC is a major issue when computing graph theory
measures in network neuroscience, yet, remains an unacknowledged problem in
GNNs research. We proposed to integrate graph diffusion as a prepossessing
step and show that it can partially alleviate the problem and consistently
improve the results. While not the only solution, this calls for a more data-
centeric approach when developing learning models from fMRI graphs over a
model-centeric approach where we observed that architectural differences did
not play an important role in the final outcome in the studied applications.
On a more general note, the empirical results open the debate of the utility
of deep learning in neuroimaging applications (He et al., 2020; Schulz et al.,
2020). On the clinical datasets, the SVM model performed competitively with
the DL models without the need for engineering complex neural network
architectures. We argue that the competitiveness of DL methods at this scale
is a positive rather than a criticism. On the UkBioBank, the dynamic DL models
(1D CNN, AST-GCN) do show improved performance over the SVM and other static
baselines. This guides towards developing DL models that utilize minimally
pre-processed dynamic signal over static feature engineered summaries such as
correlation values. Further, to objectively identify the merits of applying
GNNs for single-subject prediction we have to look further than the test
metrics and consider the advantage of enabling the application of deep
learning methods in this domain. DL is a versatile class of algorithms with
applications that extend beyond supervised classification and would
potentially be of remarkable value when applied on fMRI data. For example,
graph explainability techniques (Ying et al., 2019; Huang et al., 2022) could
be utilized to find potential functional connectivity biomarkers for
phenotypes and psychiatric disorders. The application of graph self-supervised
learning techniques Fedorov et al. (2021); Wang et al. (2022); Peng et al.
(2022) and graph normative modelling (Gazzar et al., 2022) on non-clinical
datasets would assist in learning representations that could further be
transferred to a downstream clinical tasks. Graph augmentation methods Zhao et
al. (2021) could boost the statistical power of the clinical datasets and
assist in solving the problem of labeled data scarcity and regularize
overfitting. These are some of the promising directions that support the
utility of GNNs in neuroimaging. That said, the empirical results presented in
this work advocate for more moderation and more rigorous validation when
reporting GNN results that do not clearly outperform graph-agnostic baselines.
Finally, we hope that framework developed throughout this work could serve as
a baseline for future novel GNN architectures or applications on fMRI data and
assist in identifying key components that could improve existing methods.

## Acknowledgments

This work was supported by the Netherlands Organization for Scientific
Research (NWO; 628.011.023); Philips Research; ZonMW (Vidi; 016.156.318). The
access of the UKbioBank data was granted under the application number 30091.

## References

  * Huettel et al. [2004] S. A. Huettel, A. W. Song, and G. McCarthy.  Functional magnetic resonance imaging.  _vol._ , 1, 2004.
  * Sporns et al. [2005] O. Sporns, G. Tononi, and R. Kotter.  The human connectome: a structural description of the human brain.  _PLoS Computational Biology_ , 1, 2005.
  * Sporns [2018] O. Sporns.  Graph theory methods: applications in brain networks.  _Dialogues in clinical neuroscience_ , 20:111, 2018.
  * Bronstein et al. [2017] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst.  Geometric deep learning: going beyond euclidean data.  _IEEE Signal Processing Magazine_ , 34:18–42, 2017.
  * Ktena et al. [(2017] S. I. Ktena, S. Parisot, E. Ferrante, M. Rajchl, M. Lee, B. Glocker, et al.  Distance metric learning using graph convolutional networks: Application to functional brain networks.  pages 469–477. Conference on Medical Image Computing and Computer-Assisted Intervention, (2017).
  * Parisot et al. [(2018] S. Parisot, S. I. Ktena, E. Ferrante, M. Lee, R. Guerrero, B. Glocker, et al.  Disease prediction using graph convolutional networks: application to autism spectrum disorder and alzheimer’s disease.  _Medical image analysis_ , 48:117–130, (2018).
  * Gallo et al. [(2021] S. Gallo, A. ElGazzar, P. Zhutovsky, R. M. Thomas, N. Javaheripour, L. Meng, et al.  _Thalamic hyperconnectivity as neurophysiological signature of major depressive disorder in two multicenter studies_.  (2021).
  * Li et al. [2020] X. Li, N. C. Dvornek, J. Zhuang, P. Ventola, and J. Duncan.  Graph embedding using infomax for asd classification and brain functional difference detection.  _In Medical Imaging_ , 1317:11317, 2020.
  * Li et al. [2019] X. Li, N. C. Dvornek, Y. Zhou, J. Zhuang, P. Ventola, and J. S. Duncan.  Graph neural network for interpreting task-fmri biomarkers.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 485–493, 2019.
  * Kipf and Welling [2016] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.arxiv.  preprint, 2016.
  * Xu et al. [2018] K. Xu, W. Hu, J. Leskovec, and S. Jegelka.  How powerful are graph neural networks? arxiv.  preprint, 2018.
  * Kim and Ye [2020] Byung-Hoon Kim and Jong Chul Ye.  Understanding graph isomorphism network for rs-fmri functional connectivity analysis.  _Frontiers in neuroscience_ , page 630, 2020.
  * Velic´ et al. [2017] P. Velic´, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio.  Graph attention networks.  preprint, 2017.
  * Hu et al. [2021] J. Hu, L. Cao, T. Li, S. Dong, and P. Li.  Gat-li: a graph attention network based learning and interpreting method for functional brain network classification.  _BMC bioinformatics_ , 22:1–20, 2021.
  * El-Gazzar et al. [2021] A. El-Gazzar, R. M. Thomas, and G. van Wingen.  Dynamic adaptive spatio-temporal graph convolution for fmri modeling.  _In International Workshop on Machine Learning in Clinical Neuroimaging._ , pages 125–134, 2021.
  * Errica et al. [2019] F. Errica, M. Podda, D. Bacciu, and A. Micheli.  A fair comparison of graph neural networks for graph classification. arxiv.  preprint, 2019.
  * Yan et al. [2019] C.-G. Yan, X. Chen, L. Li, F. X. Castellanos, T.-J. Bai, Q.-J. Bo, et al.  Reduced default mode network functional connectivity in patients with recurrent major depressive disorder.  _Proceedings of the National Academy of Sciences_ , 116:9078–9083, 2019.
  * Martino et al. [2014] A. D. Martino, C.-G. Yan, Q. Li, E. Denio, F. X. Castellanos, K. Alaerts, et al.  _Molecular psychiatry_ , 19:659–667, 2014.
  * Sudlow et al. [2015] C. Sudlow, J. Gallacher, N. Allen, V. Beral, P. Burton, J. Danesh, et al.  Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  _PLoS medicine e1001779_ , 12, 2015.
  * Zhou et al. [2020] Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  _AI Open_ , 1:57–81, 2020.
  * Wu et al. [2019] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chengqi Zhang.  Graph wavenet for deep spatial-temporal graph modeling.  _arXiv preprint arXiv:1906.00121_ , 2019.
  * Arya et al. [2020] Devanshu Arya, Richard Olij, Deepak K Gupta, Ahmed El Gazzar, Guido Wingen, Marcel Worring, and Rajat Mani Thomas.  Fusing structural and functional mris using graph convolutional networks for autism classification.  In _Medical imaging with deep learning_ , pages 44–61. PMLR, 2020\.
  * Qin et al. [2022] Kun Qin, Du Lei, Walter HL Pinaya, Nanfang Pan, Wenbin Li, Ziyu Zhu, John A Sweeney, Andrea Mechelli, and Qiyong Gong.  Using graph convolutional network to characterize individuals with major depressive disorder across multiple imaging sites.  _EBioMedicine_ , 78:103977, 2022.
  * Oh et al. [2022] Kang-Han Oh, Il-Seok Oh, Uyanga Tsogt, Jie Shen, Woo-Sung Kim, Congcong Liu, Nam-In Kang, Keon-Hak Lee, Jing Sui, Sung-Wan Kim, et al.  Diagnosis of schizophrenia with functional connectome data: a graph-based convolutional neural network approach.  _BMC neuroscience_ , 23(1):1–11, 2022.
  * Lei et al. [2022] Du Lei, Kun Qin, Walter HL Pinaya, Jonathan Young, Therese van Amelsvoort, Machteld Marcelis, Gary Donohoe, David O Mothersill, Aiden Corvin, Sandra Vieira, et al.  Graph convolutional networks reveal network-level functional dysconnectivity in schizophrenia.  _Schizophrenia Bulletin_ , 2022.
  * Weisfeiler and Leman [1968] B. Weisfeiler and A. Leman.  The reduction of a graph to canonical form and the algebra which appears therein.  _NTI, Series_ , 2:12–16, 1968.
  * Yu et al. [2017] B. Yu, H. Yin, and Z. Zhu.  Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. arxiv.  preprint, 2017.
  * Gadgil et al. [2020] S. Gadgil, Q. Zhao, A. Pfefferbaum, E. V. Sullivan, E. Adeli, and K. M. Pohl.  Spatio-temporal graph convolution for resting-state fmri analysis.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 528–538, 2020.
  * Azevedo et al. [2020] T. Azevedo, L. Passamonti, P. Lio, and N. Toschi.  Towards a predictive spatio-temporal representation of brain data. arxiv.  preprint, 2020.
  * Kong et al. [2021] Youyong Kong, Shuwen Gao, Yingying Yue, Zhenhua Hou, Huazhong Shu, Chunming Xie, Zhijun Zhang, and Yonggui Yuan.  Spatio-temporal graph convolutional network for diagnosis and treatment response prediction of major depressive disorder from functional connectivity.  _Human brain mapping_ , 42(12):3922–3933, 2021\.
  * Bai et al. [2020] L. Bai, L. Yao, C. Li, X. Wang, and C. Wang.  Adaptive graph convolutional recurrent network for traffic forecasting. arxiv.  preprint, 2020.
  * Braun et al. [2012] Urs Braun, Michael M Plichta, Christine Esslinger, Carina Sauer, Leila Haddad, Oliver Grimm, Daniela Mier, Sebastian Mohnke, Andreas Heinz, Susanne Erk, et al.  Test–retest reliability of resting-state connectivity network characteristics using fmri and graph theoretical measures.  _Neuroimage_ , 59(2):1404–1412, 2012.
  * Van Wijk et al. [2010] Bernadette CM Van Wijk, Cornelis J Stam, and Andreas Daffertshofer.  Comparing brain networks of different size and connectivity density using graph theory.  _PloS one_ , 5(10):e13701, 2010.
  * Hutchison et al. [(2013] R. M. Hutchison, T. Womelsdorf, E. A. Allen, P. A. Bandettini, V. D. Calhoun, M. Corbetta, et al.  Dynamic functional connectivity: promise, issues, and interpretations.  _Neuroimage_ , 14(630):360–378, (2013).
  * Hammond et al. [2011] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval.  Wavelets on graphs via spectral graph theory.  _Applied and Computational Harmonic Analysis_ , 30(2):129–150, 2011.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In _International conference on machine learning_ , pages 1263–1272. PMLR, 2017.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang.  Connecting the dots: Multivariate time series forecasting with graph neural networks.  In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_, pages 753–763, 2020.
  * Martins and Astudillo [2016] Andre Martins and Ramon Astudillo.  From softmax to sparsemax: A sparse model of attention and multi-label classification.  In _International conference on machine learning_ , pages 1614–1623. PMLR, 2016.
  * Craddock et al. [(2013b] C. Craddock, S. Sikka, B. Cheung, R. Khanuja, S. S. Ghosh, C. Yan, et al.  _Frontiers in Neuroinformatics_ , 42:10–3389, (2013b).
  * Craddock et al. [(2013a] C. Craddock, Y. Benhajali, C. Chu, F. Chouinard, A. Evans, A. Jakab, et al.  _Frontiers in Neuroinformatics_ , 7, (2013a).
  * Craddock et al. [2012] R. C. Craddock, G. A. James, P. E. H. Iii, X. P. Hu, and H. S. Mayberg.  A whole brain fmri atlas generated via spatially constrained spectral clustering.  _Human brain mapping_ , 33:1914–1928, 2012.
  * Yan and Zang [2010] C. Yan and Y. Zang.  Dparsf: a matlab toolbox for" pipeline" data analysis of resting-state fmri.  _Frontiers in systems neuroscience_ , 4:13, 2010.
  * Makris et al. [(2006] N. Makris, J. M. Goldstein, D. Kennedy, S. M. Hodge, V. S. Caviness, S. V. Faraone, et al.  _Decreased_ , 83:155–171, (2006).
  * Tzourio-Mazoyer et al. [(2002] N. Tzourio-Mazoyer, B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, et al.  Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain.  _Neuroimage_ , 15:273–289, (2002).
  * El Gazzar et al. [2019] Ahmed El Gazzar, Leonardo Cerliani, Guido van Wingen, and Rajat Mani Thomas.  Simple 1-d convolutional networks for resting-state fmri based classification in autism.  In _2019 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–6. IEEE, 2019.
  * van den Heuvel et al. [2017] Martijn P van den Heuvel, Siemon C de Lange, Andrew Zalesky, Caio Seguin, BT Thomas Yeo, and Ruben Schmidt.  Proportional thresholding in resting-state fmri functional connectivity networks and consequences for patient-control connectome studies: Issues and recommendations.  _Neuroimage_ , 152:437–449, 2017.
  * Fornito et al. [2012] Alex Fornito, Andrew Zalesky, Christos Pantelis, and Edward T Bullmore.  Schizophrenia, neuroimaging and connectomics.  _Neuroimage_ , 62(4):2296–2314, 2012.
  * Klicpera et al. [2019] Johannes Klicpera, Stefan Weißenberger, and Stephan Günnemann.  Diffusion improves graph learning.  _arXiv preprint arXiv:1911.05485_ , 2019.
  * Kondor and Lafferty [2002] Risi Imre Kondor and John Lafferty.  Diffusion kernels on graphs and other discrete structures.  In _Proceedings of the 19th international conference on machine learning_ , volume 2002, pages 315–322, 2002.
  * Brin and Page [1998] Sergey Brin and Lawrence Page.  The anatomy of a large-scale hypertextual web search engine.  _Computer networks and ISDN systems_ , 30(1-7):107–117, 1998.
  * He et al. [2020] Tong He, Ru Kong, Avram J Holmes, Minh Nguyen, Mert R Sabuncu, Simon B Eickhoff, Danilo Bzdok, Jiashi Feng, and BT Thomas Yeo.  Deep neural networks and kernel regression achieve comparable accuracies for functional connectivity prediction of behavior and demographics.  _NeuroImage_ , 206:116276, 2020.
  * Schulz et al. [2020] Marc-Andre Schulz, BT Yeo, Joshua T Vogelstein, Janaina Mourao-Miranada, Jakob N Kather, Konrad Kording, Blake Richards, and Danilo Bzdok.  Different scaling of linear models and deep learning in ukbiobank brain images versus machine-learning datasets.  _Nature communications_ , 11(1):1–15, 2020.
  * Ying et al. [2019] Rex Ying, Dylan Bourgeois, Jiaxuan You, Marinka Zitnik, and Jure Leskovec.  Gnn explainer: A tool for post-hoc explanation of graph neural networks.  _Neural Information Processing Systems_ , 2019.
  * Huang et al. [2022] Qiang Huang, Makoto Yamada, Yuan Tian, Dinesh Singh, and Yi Chang.  Graphlime: Local interpretable model explanations for graph neural networks.  _IEEE Transactions on Knowledge and Data Engineering_ , 2022.
  * Fedorov et al. [2021] Alex Fedorov, Lei Wu, Tristan Sylvain, Margaux Luck, Thomas P DeRamus, Dmitry Bleklov, Sergey M Plis, and Vince D Calhoun.  On self-supervised multimodal representation learning: an application to alzheimer’s disease.  In _2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI)_ , pages 1548–1552. IEEE, 2021.
  * Wang et al. [2022] Xuesong Wang, Lina Yao, Islem Rekik, and Yu Zhang.  Contrastive functional connectivity graph learning for population-based fmri classification.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 221–230. Springer, 2022.
  * Peng et al. [2022] Liang Peng, Nan Wang, Jie Xu, Xiaofeng Zhu, and Xiaoxiao Li.  Gate: Graph cca for temporal self-supervised learning for label-efficient fmri analysis.  _arXiv preprint arXiv:2203.09034_ , 2022.
  * Gazzar et al. [2022] Ahmed El Gazzar, Rajat Mani Thomas, and Guido Van Wingen.  Improving the diagnosis of psychiatric disorders with self-supervised graph state space models.  _arXiv preprint arXiv:2206.03331_ , 2022.
  * Zhao et al. [2021] Tong Zhao, Yozen Liu, Leonardo Neves, Oliver Woodford, Meng Jiang, and Neil Shah.  Data augmentation for graph neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 11015–11023, 2021.


# A nonequilibrium system on a restricted scale-free network

R. A. Dumer<EMAIL_ADDRESS>Instituto de Física - Universidade
Federal de Mato Grosso, 78060-900, Cuiabá, Mato Grosso, Brazil.  M. Godoy
<EMAIL_ADDRESS>Instituto de Física - Universidade Federal de Mato
Grosso, 78060-900, Cuiabá, Mato Grosso, Brazil.

###### Abstract

The nonequilibrium Ising model on a restricted scale-free network has been
studied with one- and two-spin flip competing dynamics employing Monte Carlo
simulations. The dynamics present in the system can be defined by the
probability $q$ in which the one-spin flip process simulate the contact with a
heat bath at a given temperature $T$, and with a probability $(1-q)$ the two-
spin flip process mimics the system subjected to an external flux of energy
into it. The system network is described by a power-law degree distribution in
the form $P(k)\sim k^{-\alpha}$, and the restriction is made by fixing the
maximum, $k_{m}$, and minimum, $k_{0}$, degree on distribution for the whole
network size. This restriction keeps finite the second and fourth moment of
degree distribution, allowing us to obtain a finite critical point for any
value of $\alpha$. For these critical points, we have calculated the
thermodynamic quantities of the system, such as, the total
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ and staggered
$\textrm{m}_{\textrm{N}}^{\textrm{AF}}$ magnetizations per spin,
susceptibility $\chi_{\textrm{N}}$, and reduced fourth-order Binder cumulant
$\textrm{U}_{\textrm{N}}$, for several values of lattice size $N$ and exponent
$1\leq\alpha\leq 5$. Therefore, the phase diagram was built and a self-
organization phenomena is observed from the transitions between
antiferromagnetic $AF$ to paramagnetic $P$, and $P$ to ferromagnetic $F$
phases. Using the finite-size scaling theory, we also obtained the critical
exponents for the system, and a mean-field critical behavior is observed,
exhibiting the same universality class of the system on the equilibrium and
out of it.

##  I Introduction

The dynamic evolution of equilibrium systems is related to the fact that the
transition rates of its states obey the principle of microscopic
reversibility. Otherwise, without the advanced tooling as proposed by Gibbs in
the equilibrium scene [1], nonequilibrium systems have aroused the interest of
researchers in finding out phase transitions with the particularities of
continuous phase transitions of reversible systems. One kind of the
nonequilibrium system is those subjected to two dynamics in competition [2,
3]. These systems are described by a master equation that involves the sum of
the operators on each present process and generally each of these processes
separately obeys the principle of microscopic reversibility. However, the
combination of these processes may not satisfy the detailed balance and the
system will be forced out of equilibrium.

In the last decades, the computerization of data acquisition on large
networks, make raised the possibility of understanding the dynamical and
topological stability of its networks. From that databases, the result is that
large networks that span fields as diverse as the World Wide Web (WWW) or
actors that have acted in a movie together, self-organize into a scale-free
state [4, 5]. This means that independent of the system and its constituents,
the probability $P(k)$ that a vertex interacts with $k$ other vertices in the
network, decay as a power law, i.e., $P(k)\sim k^{-\alpha}$. Barabási and
Albert [5] incorporating growth and preferential attachment on its network
model, were able to obtain this scale invariance, not present in the previous
random [6] and small-world networks [7]. These models and their interesting
ability to describe real networks instigated the curiosity of researchers to
know what would be the behavior of physical systems in complex networks [8, 9,
10, 11, 12]. Among these, we can highlight the simple but powerful Ising
model, comprising both exact [13, 14] and computational [15, 16, 17, 18] or
approximate [19, 20, 21] results for the critical behavior on arbitrary
networks.

In the same way, the study of nonequilibrium physical systems has been
spreading and continuous phase transitions, characteristic of equilibrium
systems is observed [22, 23]. Moreover, the same critical exponents have been
obtained in reversible and irreversible systems, that is, they belong to the
same universality class, acting as proof of what was conjectured by Grinstein
_et al_. [24], in which says that any nonequilibrium __ stochastic spin system
with spin-flip dynamics and up-down symmetry belongs to the same universality
class. The Ising model with complex networks is already being studied with
competing dynamics, analytically in 1D [25], by Monte Carlo simulations in 2D
[26], and by Gaussian model in 3D [27]. However, these studies were made only
for small-world networks, and by Monte Carlo simulations a mean-field critical
behavior is obtained, characteristic of equilibrium systems with random
interactions and convergent fourth moment of its network degree distribution
[13, 14, 15, 17]. Another interesting feature of that nonequilibrium systems
is the self-organization phenomena between antiferromagnetic $AF$ to
paramagnetic $P$, and $P$ to ferromagnetic $F$ phase transitions, as a
function of competition parameter [2, 3, 22, 23].

With this in mind, in the present work, we have investigated the Ising model
on a restricted scale-free network, where each site of the network is occupied
by a spin variable that can assume values $\pm 1$. Divided into two
sublattices, the connections between them in the network are made by the site
interactions, and the degree distribution of the network obey a power-law
distribution, with fixed values of minimum and maximum degree. The system is
in a nonequilibrium regime by competing between two reagent dynamic processes
that do not conserve the order parameter: with competition probability $q$,
the one-spin flip process simulates the system in contact with a heat bath at
temperature $T$, and with probability $1-q$, the two-spin flip process mimics
the system subjected to an external flux into it. Thus, here we have
investigated the phase transitions of the system and verified if the phase
diagrams present the same topology of systems with these same dynamics [23,
26], and in addition, the critical exponents carrying the universality class
of the system, is compared with previous works at equilibrium system [16].

This article is organized as follows: In Section II, we describe the network
used and the Hamiltonian model of the system. In Section III, we present the
Monte Carlo simulation method, some details concerning the simulation
procedures, and the thermodynamic quantities of the system, also necessary for
the application of FSS analysis. The behavior of thermodynamic quantities,
phase diagrams, and critical exponents are described in Section IV. Finally,
in Section V, we present our conclusions.

##  II Model

Figure 1: Schematic representation of the restricted scale-free network. Red
circles indicate the sites on one of the sublattices, blue circles are the
sites on the other sublattice, and the black solid lines are the connections
between the two sublattices. The size of the circles is proportional to the
degree of sites, varying from $k_{0}=2$ to $k_{m}=8$ in the distribution with
$\alpha=3$, and $N=10^{2}$.

The Ising model studied in this work has $N$ spins $\sigma_{i}=\pm 1$ on a
restricted scale-free network and ferromagnetic interaction of strength
$J_{ij}$. The degree distribution on the network follows the power-law
$P(k)\sim k^{-\alpha}$ and to distribute the connections between the sites, we
have used the same procedures shown in the paper [16]. In order to construct a
scale-free network with always convergent second and fourth moments on its
degree distribution and arbitrary value of $\alpha$. For that, we first define
minimum $k_{0}$ and maximum $k_{m}$ degree, and the exponent $\alpha$ of the
distribution. The next procedure is to calculate the normalization constant of
the distribution, $A=\sum_{k=k_{0}}^{k_{m}}k^{\alpha}$, and found the smaller
network size that we can use and guarantee the degree distribution,
$N_{0}=k_{m}^{\alpha}/A$. With these values, we create a set of site numbers,
$\\{N_{k}\\}$, and that will have the respective degrees $k$, where
$N_{k}=AN/k^{\alpha}$. On that distribution of connections, we have divided
the network into two sublattices, where one sublattice plays the role of
central spins, while the other sublattice contains the spins in which the
central spins can connect. Thus, starting with the lowest degree $k_{0}$,
connections of each $N_{k_{0}}$ sites are randomly created connecting the two
sublattices, and it was made until reach degree $k_{m}$ and the whole set
$\\{N_{k}\\}$ will be visited. An example of that construction can be seen in
Fig. 1 which was chosen $\alpha=3$, $k_{0}=2$, $k_{m}=8$ and $N=10^{2}$. In
Fig. 1, the sites in the middle of the figure are the more connected, while
the peripheral sites are the less connected, and sites from the blue
sublattice are only connected with sites from the red sublattice.

Based on this construction, in the course of this work, we have selected the
integer values of $1\leq\alpha\leq 5$, $k_{0}=4$, $k_{m}=10$, and network size
$(32)^{2}\leq N\leq(256)^{2}$ to study the nonequilibrium Ising model. The
ferromagnetic Ising spin energy is described by the Hamiltonian on the form

$\mathcal{H}=-\sum_{\left\langle i,j\right\rangle}J_{ij}\sigma_{i}\sigma_{j}$
(1)

where the sum is over all pair of spins, and $J_{ij}$ is the ferromagnetic
interaction, assuming the value of unity if sites $i$ and $j$ interact between
the sublattices.

In the nonequilibrium system presented here, let $p(\\{\sigma\\},t)$ be the
probability of finding the system in the state
$\\{\sigma\\}=\\{\sigma_{1},...,\sigma_{i},...,\sigma_{j},...\sigma_{N}\\}$ at
time $t$, the motion equation for the probability of states evolves in time
according to the master equation

$\frac{d}{dt}p(\\{\sigma\\},t)=qG+(1-q)D,$ (2)

where $qG$ represents the one-spin flip process, relaxing the spins in contact
with a heat bath at temperature $T$, favoring the lowest energy state of the
system, and has probability $q$ to occur. On the other hand, the $(1-q)D$
denotes the two-spin flip process, in which the energy of the system increases
by one external flow of energy into it, and has a probability $(1-q)$ to
occur. $G$ and $D$ are described as follows:

$\begin{array}[]{ccc}G=&\sum_{i,\\{\sigma^{\prime}\\}}\left[W(\sigma_{i}\to\sigma_{i}^{\prime})p(\\{\sigma\\},t)+\right.\\\
&\left.-W(\sigma_{i}^{\prime}\to\sigma_{i})p(\\{\sigma^{\prime}\\},t)\right]&,\end{array}$
(3)

$\begin{array}[]{ccc}D=&\sum_{i,j,\\{\sigma^{\prime}\\}}\left[W(\sigma_{i}\sigma_{j}\to\sigma_{i}^{\prime}\sigma_{j}^{\prime})p(\\{\sigma\\},t)+\right.\\\
&\left.-W(\sigma_{i}^{\prime}\sigma_{j}^{\prime}\to\sigma_{i}\sigma_{j})p(\\{\sigma^{\prime}\\},t)\right]&,\end{array}$
(4)

where $\\{\sigma^{\prime}\\}$ is the spin configuration after spin flipping,
$W(\sigma_{i}\to\sigma_{i}^{\prime})$ is the transition rate between the
states in the one-spin flip process, and
$W(\sigma_{i}\sigma_{j}\to\sigma_{i}^{\prime}\sigma_{j}^{\prime})$ the
transition rate between the states in the two-spin flip process.

##  III Monte Carlo simulations

Figure 2: Behavior of the thermodynamic quantities $\textrm{m}_{\textrm{L}}$,
$\textrm{U}_{\textrm{L}}$, and $\textrm{$\chi$}_{\textrm{L}}$as a function of
$q$ for different network sizes as presented in the figures. In the $AF-P$
transition: (a) $\textrm{m}_{\textrm{L}}^{\textrm{AF}}$, (b)
$\textrm{U}_{\textrm{L}}^{\textrm{AF}}$, and (c)
$\textrm{$\chi$}_{\textrm{L}}^{\textrm{AF}}$curves. In the $F-P$ transition:
(d) $\textrm{m}_{\textrm{L}}^{\textrm{F}}$ , (e)
$\textrm{U}_{\textrm{L}}^{\textrm{F}}$, and (f)
$\textrm{$\chi$}_{\textrm{L}}^{\textrm{F}}$curves. Here, we have fixed values
of $\alpha=1$, $k_{0}=4$, $k_{m}=10$, $T=1$.

In the simulation of the system specified by the Hamiltonian in Eq. (1), we
always have chosen the initial state of the system with all spin states at
random, and a new configuration is generated by the following Markov process:
for a given temperature $T$, competition probability $q$, distribution
exponent $\alpha$, network size $N$, and minimum $k_{0}$ and maximum $k_{m}$
degree, we choose at random a spin $\sigma_{i}$ in network, and generate a
random number $\xi$ between zero and one. If $\xi\leq q$, we choose the one-
spin flip process, in which the flipping probability is dependent of
$W(\sigma_{i}\to\sigma_{i}^{\prime})$ and given by the Metropolis
prescription:

$W(\sigma_{i}\to\sigma_{i}^{\prime})=\left\\{\begin{array}[]{cccc}e^{\left(-\Delta
E_{i}/k_{B}T\right)}&\textrm{if}&\Delta E_{i}>0\\\ 1&\textrm{if}&\Delta
E_{i}\leq 0&,\end{array}\right.$ (5)

where $\Delta E_{i}$ is the change in energy after flipping the spin,
$\sigma_{i}\to\sigma_{i}^{\prime}$, $k_{B}$ is the Boltzmann constant, and $T$
the temperature of the system. Thus, the acceptance of a new state is
guaranteed if $\Delta E_{i}\leq 0$, but, in the case where $\Delta E>0$ the
acceptance is pondered by the probability $\exp\left(-\Delta
E_{i}/k_{B}T\right)$ and just guaranteed if by choosing a random number,
$0<\xi_{1}<1$, it is $\xi_{1}\leq\exp\left(-\Delta E_{i}/k_{B}T\right)$. On
the other hand, if none of these conditions are satisfied, we do not change
the state of the system. Now, if $\xi>q$ the two-spin flip process is chosen,
and in addition to the spin $\sigma_{i}$ we also randomly choose one of its
neighbors $\sigma_{j}$, and these two spins are flipping simultaneously
according to transition rate
$W(\sigma_{i}\sigma_{j}\to\sigma_{i}^{\prime}\sigma_{j}^{\prime})$ given by

$W(\sigma_{i}\sigma_{j}\to\sigma_{i}^{\prime}\sigma_{j}^{\prime})=\left\\{\begin{array}[]{c}0\\\
1\end{array}\begin{array}[]{c}\textrm{if}\\\
\textrm{if}\end{array}\begin{array}[]{cc}\Delta E_{ij}\leq 0\\\ \Delta
E_{ij}>0&,\end{array}\right.$ (6)

where $\Delta E_{ij}$ is the change in the energy after flipping the spins
$\sigma_{i}$ and $\sigma_{j}$, and consequently, in this process, the new
state is only accepted if $\Delta E_{ij}>0$.

Repeating the Markov process $N$ times, we have one Monte Carlo Step (MCS). In
our simulations, we have waited for $10^{4}$ MCS to the system reach the
stationary state, in the whole the network sizes and adjustable parameters. In
order to calculate the thermal averages of the interest quantities, we used
more $4\times 10^{4}$ MCS, and the average over samples was done using $10$
independent samples for any configuration.

The measured thermodynamic quantities in our simulations are: magnetization
per spin $\textrm{m}_{\textrm{N}}^{\textrm{F}}$, staggered magnetization per
spin $\textrm{m}_{\textrm{N}}^{\textrm{AF}}$, magnetic susceptibility
$\textrm{$\chi$}_{\textrm{N}}$ and reduced fourth-order Binder cumulant
$\textrm{U}_{\textrm{N}}$:

$\textrm{m}_{\textrm{N}}^{\textrm{F}}=\frac{1}{N}\left[\left\langle\sum_{i=1}^{N}\sigma_{i}\right\rangle\right],$
(7)

$\textrm{m}_{\textrm{N}}^{\textrm{AF}}=\frac{1}{N}\left[\left\langle\sum_{i=1}^{N}(-1)^{(r+c)}\sigma_{i}\right\rangle\right],$
(8)

$\textrm{{\hbox{\chi}}}_{\textrm{N}}=\frac{N}{k_{B}T}\left[\left\langle
m^{2}\right\rangle-\left\langle m\right\rangle^{2}\right],$ (9)

$\textrm{U}_{\textrm{N}}=1-\frac{\left[\left\langle
m^{4}\right\rangle\right]}{3\left[\left\langle
m^{2}\right\rangle^{2}\right]},$ (10)

where $\left[\ldots\right]$ representing the average over the samples, and
$\left\langle\ldots\right\rangle$ the thermal average over the MCS in the
stationary state. To facilitate the calculation of
$\textrm{m}_{\textrm{N}}^{\textrm{AF}}$, the sites on the network are labeled
as if we had a square lattice, $N=L^{2}$, in this way, $r$ and $c$ are the row
and column of the site $i$, respectively. In Eqs. (9) and (10), $m$ can be
used to represent $\textrm{m}_{\textrm{N}}^{\textrm{F}}$ or
$\textrm{m}_{\textrm{N}}^{\textrm{AF}}$.

In the vicinity of the stationary critical point $\lambda_{c}$, the Eqs. (7),
(8), (9) and (10) obey the following finite-size scaling relations [28]:

$\textrm{m}_{\textrm{N}}=N^{-\beta/\nu}m_{0}(N^{1/\nu}\epsilon),$ (11)

$\textrm{{\hbox{\chi}}}_{\textrm{N}}=N^{\gamma/\nu}\chi_{0}(N^{1/\nu}\epsilon),$
(12)

$\textrm{U}_{\textrm{N}}^{\prime}=N^{1/\nu}\frac{U_{0}^{\prime}(N^{1/\nu}\epsilon)}{\lambda_{c}},$
(13)

where $\epsilon=(\lambda-\lambda_{c})/\lambda_{c}$ ($\lambda$ and
$\lambda_{c}$ can be used $T$ or $q$), and $\beta$, $\gamma$ and $\nu$ are the
critical exponents related the magnetization, susceptibility and length
correlation, respectively. The functions $m_{0}(N^{1/\nu}\epsilon)$,
$\chi_{0}(N^{1/\nu}\epsilon)$ and $U_{0}(N^{1/\nu}\epsilon)$ are the scaling
functions.

Using the data from simulations for the network sizes $(32)^{2}\leq
N\leq(256)^{2}$ in the Eqs. (11), (12) and (13), we have obtained the critical
exponents ratio $\beta/\nu$, $\gamma/\nu$, and $\nu^{-1}$ from the slope of
the straight lines in the log-log plot of
$\textrm{m}_{\textrm{N}}(\lambda_{c})$,
$\textrm{$\chi$}_{\textrm{N}}(\lambda_{c})$, and
$\textrm{U}_{\textrm{N}}^{\prime}(\lambda_{c})$ (derivative of
$\textrm{U}_{\textrm{N}}$) as a function of $N$. Besides that, we also used
data collapse from scaling functions to estimate the critical exponent values.

##  IV Results

Figure 3: Phase diagrams in the plane temperature $T$ versus $q$ for some
values of the exponent $\alpha$ and fixed values of $k_{0}=4$ and $k_{m}=10$.
On the top, we have the color bar of $m_{L}$ where the left side represent the
staggered magnetization $\textrm{m}_{\textrm{N}}^{\textrm{AF}}$ illustrating
the $AF-P$ transition and the right side the magnetization per spin
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ illustrating the $F-P$ transition. The
magenta circles are the critical points estimated by the crossing of
$\textrm{U}_{\textrm{L}}$ curves and the black solid lines are just a guide
for the eyes indicating the phase transition lines.

In this section, we present and discuss the results of the nonequilibrium
Ising model on a restricted scale-free network. For the two dynamic processes,
we have an adjustable parameter $q$ that controls the dynamic competition in
the system. If $0<q<1$, the two dynamic processes have a non-null probability
to be chosen and acting in the system, making it irreversible with respect to
the temporal evolution of its states. As these processes favor the states of
higher and lower energy of the system, with the competition is possible to
find stationary states in the $AF$, $F$, and $P$ phases, based on the
Hamiltonian of the system, Eq. (1). With this, it is worth noting that to
obtain a self-organization phenomenon passing from a $F$ to $P$ and from $P$
to $AF$ phases, the division of the network into two sublattices is essential,
once that for nonfrustrated antiparallelism we must to have well-defined who
the central spins are, and to whom they can connect in the network.

Figure 4: Data collapse of $\textrm{m}_{\textrm{L}}^{\textrm{AF}}$ (a),
$\textrm{$\chi$}_{\textrm{L}}^{\textrm{AF}}$ (b),
$\textrm{m}_{\textrm{L}}^{\textrm{F}}$ (c) and
$\textrm{$\chi$}_{\textrm{L}}^{\textrm{F}}$ (d) for different network sizes as
presented in the figures. In (a) and (c) from the left to the right side we
have respectively $\alpha=1,2,3,4$ and $5$, and
$\textrm{m}_{\textrm{L}}^{\textrm{AF}}$ and
$\textrm{m}_{\textrm{L}}^{\textrm{F}}$color bars as shown in the figures. In
(b) and (d), with $\epsilon=(q-q_{c})/q_{c}$ color bar and from the top to
bottom we have the collapses with $\alpha=1,2,3,4$ and $5$, respectively. In
these figures, we have changed the positions curves for $1<\alpha\leq 5$ in
order to compare all the collapses obtained. The critical exponents and
critical points used here can be seen in Table 1 and Table 2 respectively.
Here, we have fixed $T=1$, $k_{0}=4$, $k_{m}=10$.

Therefore, the first results can be seen in Fig. 2, where we have displayed
the thermodynamic quantities obtained with Eqs. (7), (8), (9) and (10). These
quantities were calculated as a function of the competition parameter $q$, in
which is verified that for lower values of $q$ we found an $AF$ phase, and for
higher values of $q$, an $F$ phase is observed. These phases are easily
explained when we look at the dynamics, once that for lower values of $q$, the
two-spin flip mechanism prevails and this favors the state of high energy in
the system, which based on the ferromagnetic Ising model Hamiltonian is the
one where the spin states are antiparallel, i. e., $AF$ phase. This $AF$ phase
is made explicit in Fig. 2(a) with the $\textrm{m}_{\textrm{N}}^{\textrm{AF}}$
curves, and with this magnetization is calculated
$\textrm{U}_{\textrm{N}}^{\textrm{AF}}$ present in Fig. 2(b), and its
susceptibility $\textrm{$\chi$}_{\textrm{N}}^{\textrm{AF}}$ in Fig. 2(c). On
the other hand, for higher values of $q$, the one-spin flip mechanism
prevails, and as it favors the states of lower energy in the system, i.e., all
spins in the same state, a ordered phase is also observed, $F$ phase. The
quantities related to this phase is specifically the magnetization
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ curves in Fig. 2(d), and the
$\textrm{U}_{\textrm{N}}^{\textrm{F}}$ and
$\textrm{$\chi$}_{\textrm{N}}^{\textrm{F}}$ curves in Figs. 2(e) and 2(f),
respectively.

Figure 5: In log-log plots is presented linear fits (black dotted lines) of
thermodynamic quantities in the critical point,
$\textrm{m}_{\textrm{L}}^{\textrm{AF}}$ $\left(\square\right)$ and
$\textrm{m}_{\textrm{L}}^{\textrm{F}}$ $\left(\Circle\right)$ in (a),
$\chi_{\textrm{L}}^{\textrm{AF}}$ $\left(\square\right)$ and
$\textrm{$\chi$}_{\textrm{L}}^{\textrm{F}}$ $\left(\Circle\right)$ in (b), and
$\textrm{U}_{\textrm{L}}^{\textrm{AF}}$ $\left(\square\right)$ and
$\textrm{U}_{\textrm{L}}^{\textrm{F}}$ $\left(\Circle\right)$ in (c), both has
a function of the network size $N$, and different values of $\alpha$ as shown
in the figures. Here we have fixed $T=1$, $k_{0}=4$ and $k_{m}=10$. The
critical points used for these fits can be seen in Table 2.

We have used the curves of the fourth-order Binder cumulants for different
network sizes to identify the critical points and order phase transition [29,
30, 31, 32]. The intersection point of the $\textrm{U}_{\textrm{N}}$ curves
indicates the phase transition point on a second-order phase transition. With
the critical point in hand for several values of adjustable parameters, a
phase diagram was built, which can be seen in Fig. 3. Therefore, for these
diagrams and later results, we will limit the values $k_{0}=4$ and $k_{m}=10$,
once we can build all networks with sizes ($32)^{2}\leq N\leq(256)^{2}$,
integer exponent $1\leq\alpha\leq 5$, and compare with others equilibrium [13,
14, 16] and nonequilibrium [23, 26] Ising model results. Fig. 3 presents the
phase diagrams of temperature $T$ as a function of competition parameter $q$
for some values of $\alpha$, in which we can see the $AF$, $F$, and $P$
phases.

Figure 6: Average static critical exponents $\beta$, $\gamma$, and $\nu$,
obtained form the slope of scaling functions and the data collapse of
magnetization and susceptibility curves as a function of $\alpha$. (a) For
$AF-P$ transition and (b) for the $F-p$ transition. These exponents were
obtained with fixed values of $T=1$, $k_{0}=4$ and $k_{m}=10$.

In these diagrams (Fig. 3), we have illustrated the self-organization
phenomena with the transitions between $AF$ to $P$ phases, and $P$ to $F$
phases. Since the scale is fixed in all the figures, we can also see that when
we decrease the value of $\alpha$, the region of ordered phases, $F$ and $AF$,
increases. This change in the topology of the diagram is related to the degree
distribution, once the lower values of exponent $\alpha$ mean a high
probability of having more connected sites on network, i.e., more sites with a
degree $k_{m}$. Consequently, knowing that more connected sites on the
stationary ordered state require more energy to override its interactions,
larger are the regions of the ordered phases. Another interesting observation
that we can do, is regarding the shape of the regions in the ordered phases.
The ferromagnetic phases are driven by the one-spin flip mechanism described
by Metropolis prescription, which is very dependent on $T$, and for high $T$
we observe the disordered phase $P$. On the other hand, $AF$ phases is driven
by the two-spin flip mechanism, in which is a simpler process and little
influenced by temperature.

All systems belonging to a given universality class share the same set of
critical exponents. The critical points can be used to describe the critical
behavior in the sense of universality class with the set of critical
exponents. Here, we have computed the exponents $\beta$, $\gamma$ and $\nu$,
by two methods. The first one is based on the data collapse, in which we use
the scaling relations, Eqs. (11) and (12), to obtain the scaling functions of
magnetization and susceptibility with its collapsed curves. This is possible
because in the proximity of the critical points the scaling relations are
independent of network size with the correct critical exponents and critical
point of the system [29, 30, 31]. To obtain the critical exponents by this
method and using the already estimated critical points, we have plotted the
scaling functions $m_{0}(N^{1/\nu}\epsilon)$ and $\chi_{0}(N^{1/\nu}\epsilon)$
as a function of $|\epsilon|N^{1/\nu}$ for different network sizes and in the
proximity of critical points. Therefore, for $\epsilon\to 0$ and adjusting the
involved critical exponents, when the curves of different network sizes
collapse better into a single curve, these exponents used are considered the
critical exponents of the system.

$\alpha$ | $\beta_{F}$ | $\gamma_{F}$ | $\nu_{F}\left(\textrm{m}_{\textrm{L}}\right)$ | $\nu_{F}\left(\textrm{$\chi$}_{\textrm{L}}\right)$ | $\beta_{AF}$ | $\gamma_{AF}$ | $\nu_{AF}\left(\textrm{m}_{\textrm{L}}\right)$ | $\nu_{AF}\left(\textrm{$\chi$}_{\textrm{L}}\right)$
---|---|---|---|---|---|---|---|---
$1$ | $0.51\pm 0.04$ | $0.98\pm 0.03$ | $1.98\pm 0.03$ | $2.02\pm 0.04$ | $0.50\pm 0.03$ | $1.00\pm 0.02$ | $2.00\pm 0.04$ | $2.00\pm 0.03$
$2$ | $0.54\pm 0.04$ | $0.92\pm 0.03$ | $2.00\pm 0.05$ | $2.00\pm 0.04$ | $0.51\pm 0.04$ | $0.95\pm 0.03$ | $2.00\pm 0.05$ | $2.00\pm 0.03$
$3$ | $0.52\pm 0.04$ | $0.93\pm 0.05$ | $2.00\pm 0.02$ | $1.96\pm 0.02$ | $0.50\pm 0.04$ | $1.00\pm 0.04$ | $2.00\pm 0.03$ | $2.02\pm 0.04$
$4$ | $0.56\pm 0.03$ | $0.90\pm 0.03$ | $1.96\pm 0.03$ | $2.06\pm 0.04$ | $0.56\pm 0.03$ | $0.96\pm 0.04$ | $2.02\pm 0.04$ | $2.04\pm 0.03$
$5$ | $0.54\pm 0.04$ | $0.95\pm 0.04$ | $1.96\pm 0.04$ | $2.06\pm 0.03$ | $0.55\pm 0.04$ | $0.95\pm 0.03$ | $2.01\pm 0.02$ | $1.92\pm 0.03$
Table 1: Critical exponents obtained by the data collapse method for several
values $\alpha$. The $F-P$ transitions are denoted by $F$ subscript and the
$AF-P$ transitions are denoted by $AF$ subscript. In these transitions, the
data collapse of magnetization and susceptibility curves returns us
respectively $\nu(\textrm{m}_{\textrm{L}})$ and
$\nu(\textrm{$\chi$}_{\textrm{L}})$ estimates using the $\nu$ exponent. For
these exponents, we have fixed $T=1$ and $k_{0}=4$ and $k_{m}=10$. The
collapsed curves can be seen in Fig. 4.

Fig. 4 display the scaling functions $m_{0}(N^{1/\nu}\epsilon)$ and
$\chi_{0}(N^{1/\nu}\epsilon)$ collapsed in the log-log plot to obtain its
asymptotic behavior. In these figures, we have fixed $T=1$ to obtain the
critical exponents of the system both in the $F-P$ transition and in the
$AF-P$ transition, for all values of $\alpha$. In Fig. 4(a), we can see the
function $m_{0}(N^{1/\nu}\epsilon)$ in the log-log plot, produced with the
collapse of $\textrm{m}_{\textrm{N}}^{\textrm{AF}}$ curves, and with this was
obtained the exponents $\beta$ and $\nu$. In the same way, in Fig. 4(b) the
scaling function $\chi_{0}(N^{1/\nu}\epsilon)$ is presented in the log-log
plot with $\textrm{$\chi$}_{\textrm{N}}$ curves based on staggered
magnetization, in which with the best data collapse we have obtained the
exponent $\gamma$, and another estimated value for the $\nu$ exponent. On the
other hand, in the $F-P$ transition, Figs. 4(c) and 4(d), respectively,
contain the log-log plot of the scaling functions based on
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ and its susceptibility,
$\textrm{$\chi$}_{\textrm{N}}^{\textrm{F}}$. The asymptotic behavior, away
from the critical point of these functions, is predicted to a slope $\Theta$
related to the obtained critical exponents, once that for the magnetization
curves starting from the ordered phase, below from the critical point
$\Theta=\beta$, and above it $\Theta=\nu/2-\beta$, and for the susceptibility
curves we only have $\Theta=-\gamma$. The critical exponents obtained by this
first method are presented in Table 1 and the critical point used for them can
be seen in Table 2.

$\alpha$ | $q_{c}^{F}$ | $q_{c}^{AF}$ | $\left(-\beta/\nu\right)_{F}$ | $\left(\gamma/\nu\right)_{F}$ | $\left(1/\nu\right)_{F}$ | $\left(-\beta/\nu\right)_{AF}$ | $\left(\gamma/\nu\right)_{AF}$ | $\left(1/\nu\right)_{AF}$
---|---|---|---|---|---|---|---|---
$1$ | $0.8267\pm 0.0002$ | $0.1320\pm 0.0002$ | $0.24\pm 0.04$ | $0.51\pm 0.05$ | $0.47\pm 0.06$ | $0.24\pm 0.04$ | $0.51\pm 0.04$ | $0.47\pm 0.06$
$2$ | $0.8396\pm 0.0003$ | $0.1156\pm 0.0002$ | $0.28\pm 0.04$ | $0.48\pm 0.05$ | $0.47\pm 0.07$ | $0.29\pm 0.05$ | $0.48\pm 0.05$ | $0.47\pm 0.05$
$3$ | $0.8534\pm 0.0004$ | $0.0995\pm 0.0005$ | $0.27\pm 0.06$ | $0.49\pm 0.05$ | $0.47\pm 0.06$ | $0.27\pm 0.05$ | $0.46\pm 0.05$ | $0.49\pm 0.06$
$4$ | $0.8664\pm 0.0004$ | $0.0842\pm 0.0004$ | $0.28\pm 0.05$ | $0.49\pm 0.06$ | $0.47\pm 0.05$ | $0.26\pm 0.04$ | $0.48\pm 0.05$ | $0.48\pm 0.05$
$5$ | $0.8788\pm 0.0004$ | $0.0714\pm 0.0004$ | $0.27\pm 0.06$ | $0.46\pm 0.05$ | $0.50\pm 0.06$ | $0.30\pm 0.05$ | $0.47\pm 0.03$ | $0.48\pm 0.04$
Table 2: Critical points obtained by the intersection of
$\textrm{U}_{\textrm{L}}$ curves, and the critical exponents by the second
method for several values $\alpha$. The $F-P$ transitions are denoted by $F$
subscript and the $AF-P$ transitions are denoted by $AF$ subscript. The slope
of the straight lines are presented in Fig, 5.

Now, let use a second method to calculate the critical exponents and also
using the scaling relations, but, with the log-log plot of
$\textrm{m}_{\textrm{N}}^{\textrm{AF}}$ and
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ at its respective
$\textrm{$\chi$}_{\textrm{N}}$ and $\textrm{U}_{\textrm{N}}$ in the proximity
of the critical point as a function of $N$. The slope on this set of point
returns us specific ratios between the critical exponents. Fig. 5(a) shown the
points of $\textrm{m}_{\textrm{N}}^{\textrm{AF}}$ and
$\textrm{m}_{\textrm{N}}^{\textrm{F}}$ in the vicinity of the critical point
as a function of network sizes $N$. With the best fit of these points and its
slope based on the scaling relation of Eq. (11), gives us the estimate of the
ratio $-\beta/\nu$. In the same way, but for the susceptibility of these
magnetizations, on the vicinity of the critical point as a function of $N$ in
the log-log plot, is presented in Fig. 5(b). The best fit with the points in
this figure gives us the slope related to the ratio $\gamma/\nu$ presented in
the scaling relation of Eq. (12). The ratio between these critical exponents
is interesting but does not reveal the correct value of the exponents
separately. Thus, to solve this, we used the scaling relation in Eq. (13), in
which the derivative of $\textrm{U}_{\textrm{N}}$ in the vicinity of the
critical point and different network sizes gives us the ratio $1/\nu$. This
ratio is illustrated in Fig. 5(c) by its log-log plot. All the ratio between
the critical exponents obtained on this method can be found in Table 2.

From these two used methods are obtained equivalent exponents. But, we have to
pay attention that as we are dealing with random interactions on the network,
we do not have a well-defined dimension, and consequently, it was necessary to
use scaling relations dependent only of the system size, Eqs. (11), (12) and
(13). Therefore, the expected mean-field finite-size scaling exponents due to
these equations are $\beta=1/2$, $\gamma=1$, and $\nu=2$ [28]. If compared to
the usual Ising model mean-field exponents $\beta=\widetilde{\nu}=1/2$, and
$\gamma=1$, the only exponent affected by the dimension of the system is the
related to the correlation length, $\nu$, but, can be derived by the relation
$\nu=d_{u}\widetilde{\nu}$, where $d_{u}$ is the upper critical dimension,
that in the Ising model is $d_{u}=4$. With these information, we have computed
$\widetilde{\nu}$ based on the exponents $\nu$ obtained here. The critical
exponents of the system, obtained by the two methods are compiled in Fig. 6(a)
for $AF-P$ transitions, and in Fig. 6(b) for $F-P$ transitions, both as a
function of $\alpha$. Comparing these obtained critical exponents with the
mean-field critical exponents, we can see that for $\alpha=1$ is obtained the
more accurate mean-field critical exponents, however, as $\alpha$ increase,
the critical exponents are still of mean-field but with a little deviation.
This deviation was explained in the work with the equilibrium Ising model on a
restricted scale-free network [16], and is due to the increase of degree-
degree correlations [33] with the decreasing of more connected sites.

For the sake of curiosity, our network was labeled as a square lattice, which
we always use $N=L^{2}$ sites. Thus, changing $N$ in Eqs. (11), (12) and (13)
by $L$, we have the scaling relations depending on the dimension of the
system, that in our case is two dimensions. With these new scaling relations,
we have computed again the critical exponents of the system and we have
obtained the same critical exponents of systems upper the Ising model critical
dimension, by adding long-range interactions on a regular square lattice [15,
26]. It indicates that with our selected network sizes $N=L^{2}$, we also
could use the scaling relations depending on the dimension of the system to
calculate the critical exponents. However, when dealing with complex networks,
this dimensioning possibility is not always available, once that the objective
is to model real networks [12, 34, 35]. In this case, Hong _et al._ [28]
proposed scaling relations for complex networks independent of system
dimension, and from them, obtained the set of mean-field finite-size-scaling
exponents.

##  V Conclusions

Here, we have employed Monte Carlo simulations to study the thermodynamic
quantities and the critical behavior of the nonequilibrium Ising model on a
restricted scale-free network. By using one- and two-spin flip competing
dynamics, we reach the stationary state of the system at the nonequilibrium
regime. Fixing the maximum and minimum degree values for the whole network
size and by using FSS analysis, we are able to always find a finite critical
point even being in a network with power-law degree distribution, since we
always have second and fourth convergent moments based on its distribution
$P(k)$. As a result, we have obtained the critical points from the second-
order phase transitions based on the intersection of $\textrm{U}_{\textrm{N}}$
curves and built a phase diagram of temperature $T$ as a function of the
competition parameter $q$. In these diagrams, we have verified the self-
organization phenomena in the transitions from $AF$ to $P$ phases in lower
values of $q$, and from $P$ to $F$ phases in higher values of $q$ and lower
$T$. Because we are dealing with a power-law degree distribution on the
network, $P(k)\sim k^{-\alpha}$, decreasing the value of $\alpha$, increase
the number of more connected sites, and as consequence, also increase the
region of the ordered phases in the diagram. Topologies equivalent to these
diagrams were also obtained in previous works with the same dynamics, but in
different networks and models [23, 26]. Through FSS arguments, we calculated
the critical exponents $\beta$, $\gamma$, and $\nu$ for the system, and as a
function of $\alpha$, because we have a restricted scale-free network in which
its second and fourth moments of degree distribution are convergent. In this
case, we have always found the mean-field critical exponents and a slight
deviation from them with the increasing degree-degree correlations. This mean-
field behavior follows the predicted and observed critical behavior in other
complex networks [13, 14, 15, 16, 17], in addition to being another agreement
of what was conjectured by Grinstein _et al_. [24], i.e., we obtained the same
universality class of the Ising model on a restricted scale-free network both
in the equilibrium regime [16] as out of it.

## References

  * [1] J. W. Gibbs. Elementary Principles in Statistic Mechanics. (Yale Universality Press, New Haven, 1902);
  * [2] T. Tomé and M. J. de Oliveira. Phys. Rev. A, 40, 6643 (1989);
  * [3] T. Tomé, M. J. de Oliveira, and M. A. Santos. J. Phys. A: Math. Gen., 24, 3677 (1991);
  * [4] R. Albert and A.-L. Barabási. Rev, Mod. Phys., 74, 47 (2002);
  * [5] A.-L. Barabási and R. Albert. Science, 286, 509 (1999);
  * [6] P. Erdös and A. Rényi. Mathematical Institute of the Hungrian Academy of Sciences, 5, 17 (1960);
  * [7] D. J. Watts and S. H. Strogatz. Nature, 393, 440 (1998);
  * [8] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes. Eur. Phys. J. B, 38, 177 (2004);
  * [9] G. Bianconi and A.-L. Barabási. Phys. Rev. Lett., 86, 5632 (2001);
  * [10] S. Aparicio, J. Villazón-Terrazas, and G. Álvarez. Entropy, 17, 5848 (2015);
  * [11] A. L. M. Vilela, B. J. Zubillaga, C. Wang, M. Wang, R. Du, and H. E. Stanley. Sci. Rep., 10, 8255 (2020);
  * [12] T. Gradowski and A. Krawiecki. A. Phys. Pol. A, 127, A-55 (2015);
  * [13] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes. Phys. Rev. E, 66, 016104 (2002);
  * [14] A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes. Phys. Rev. E, 67, 026123 (2003);
  * [15] R. A. Dumer and M. Godoy. Eur. Phys. J. B 95, 159 (2022);
  * [16] R. A. Dumer and M. Godoy. Physica A, 612, 128795 (2023);
  * [17] C. P. Herrero. Phys. Rev. E, 65, 066110 (2002);
  * [18] C. P. Herrero. Phys. Rev. E, 69, 067109 (2004);
  * [19] A. Aleksiejuk, J. A. Holyst, and D. Stauffer. Physica A, 310, 260 (2002);
  * [20] G. Bianconi. Phys. Lett. A, 303, 166 (2002);
  * [21] M. Leone, A. Vazquez, A. Vespignani, and R. Zecchina. Eur. Phys. J. B, 28, 191 (2002);
  * [22] W. Figueredo and B. C. S. Grandi. Braz. J. Phys., 30, 58 (2000);
  * [23] M. Godoy and W. Figueredo. Phys. Rev. E, 65, 026111 (2002);
  * [24] G. Grinstein, C. Jayaprakash, and Yu He. Phys. Rev. Lett., 55, 2527 (1985);
  * [25] W. Liu, W.-Y. Xiong, and J.-Y. Zhu. Phys. Rev. E, 71, 056123 (2005);
  * [26] R. A. Dumer and M. Godoy. Phys. Rev. E, 107, 044115 (2023);
  * [27] J.-Y. Zhu, W. Liu and H. Zhu. Eur. Phys. J. B, 33, 545 (2003);
  * [28] H. Hong, M. Ha and H. Park. Phys. Rev. Lett., 98, 258701 (2007).
  * [29] K. Binder and D. W. Heermann. Monte Carlo Simulation in Statistical Physics. An Introduction, 6rd ed. (Springer, Cham, Switzerland, 2019);
  * [30] K. Binder and D. P. Landau. A Guide to Monte Carlo Simulations in Statistical Physics, 4rd ed. (TJ International Ltd, Padstow, UK, 2015);
  * [31] L. Böttcher and H. J. Herrmann. Computational Statistical Physics, 1rd ed. (Cambridge University Press, NewYork, EUA, 2021);
  * [32] S.-H. Tsai and S. R. Salinas. Braz. J. Phys., 28, 1, (1998);
  * [33] M. E. J. Newman. Phys. Rev. Lett., 89, 208701 (2002);
  * [34] F. W. S. Lima. Entropy, 18, 81 (2016);
  * [35] A. L. M. Vilela, B. J. Zubillaga, C. Wang, M. Wang, R. Du, and H. E. Stanley. Sci. Rep., 10, 8255 (2020);