Do graph neural network states contain graph properties?

## Abstract

Abstract Graph learning models achieve state-of-the-art performance on many tasks, but this often requires increasingly large model sizes. Accordingly, the complexity of their representations increase. Explainability techniques (XAI) have made remarkable progress in the interpretability of ML models. However, the non-relational nature of Graph Neural Networks (GNNs) make it difficult to reuse already existing XAI methods. While other works have focused on instance-based explanation methods for GNNs, very few have investigated model-based methods and, to our knowledge, none have tried to probe the embedding of the GNNs for well-known structural graph properties. In this paper we present a model agnostic explainability pipeline for Graph Neural Networks (GNNs) employing diagnostic classifiers. This pipeline aims to probe and interpret the learned representations in GNNs across various architectures and datasets, refining our understanding and trust in these models.

## 1 Introduction

In the last decade, significant progress has been made towards modelling non-Euclidean, graph-structured data (Kipf & Welling, ) on the one hand, and on interpreting the predictions of deep neural networks (DNN) on the other hand.
We often qualify DNNs as black box as their predictions are not inherently interpretable.
Occlusion, gradient, perturbation, layer-wise relevance propagation, and attention mechanisms have been proposed to solve this problem (Zeiler & Fergus, ; Denil et al., ; Li et al., ; Sundararajan et al., ).
These methods focus on highlighting
the importance of different input features.
They can, however, not be directly applied on GNNs due to the lack of a regular structure (e.g. vertices can have different degrees).
In this case, explaining a prediction means identifying important parts of the relational structure, or input features of nodes.
An issue is that finding the explanation is itself a combinatorial problem, making XAI (explainable AI) methods for GNN intractable (Longa et al., ; Ying et al., ; Lucic et al., ).

Previous surveys (Agarwal et al., ; Dai et al., ) highlighted the lack of comprehensive, robust and model-agnostic explainability methods. We also identified that there are very few model-level explainability methods. As an alternative to these more traditional XAI methods, we propose to apply probing techniques for graph properties (as developed for Natural Language Processing Giulianelli et al. (), Belinkov ()) to GNN embeddings.
In our pipeline (see fig. 1), we investigate both local properties like betweennes centrality, as well as global properties like average path length.
To our knowledge, this is the first work to explore this direction.

Figure: Figure 1: An example of the probing pipeline. First, a GNN is trained on a specific task, for example detecting whether a graph contains a grid or house shaped pattern. Then, we extract embeddings from the internal layers of the network. We use these embeddings to train the probing model; in this example a binary classifier which can detect whether the embedding contains predictive information for the diameter of the graph. If a linear probe has good performance ($R^{2}$ score) then there exists a hyperplane in the representation space that separates the inputs based on the property
Refer to caption: x1.png

Findings (^1^11All results and experiments accessible on [github](https://github.com/TomPelletreauDuris/Probing-GNN-representations)):

- •
We demonstrate the ability of diagnostic classifiers to effectively highlight known graph-theoretic and domain-specific properties in GNN learned latent representations (fig. 5).
- •
We explore how different regularization techniques (none, $L_{2}$ weight decay, dropout) affect the representation of graph properties within the same GNN architecture (fig. 10).
- •
We compare how various GNN architectures (GCN, R-GCN, GIN, GAT) differ in their ability to represent graph properties, analyzing whether these differences align with their mathematical frameworks (table 6).
- •
We apply this pipeline to a toxicity dataset showing that probed graph properties align with chemical knowledge (table 8) before exploring the pipeline’s inferential power on fMRI datasets, uncovering structural properties that might not yet have been extensively studied (table 23).

## 2 Background

### 2.1 Graph Neural Networks

Nowadays we have some theoretical understanding of the representational restrictions and capabilities of Graph Neural Networks (GNNs) with regard to the Weisfeiler-Lehman test (Akhondzadeh et al., ). We know that this cannot capture certain graph properties, such as connectivity or triangle-freeness (Franks et al., ; Kiefer, ; Kriege et al., ), due to its reliance on local structure. This constraint is also present in (message passing) GNNs.

Graph Convolutional Network (GCN) (Kipf & Welling, ) are GNNs where for a single layer, the node representation is computed as:
$\bm{X}^{\prime}=\sigma\left(\tilde{\bm{D}}^{-1/2}\cdot\tilde{\bm{A}}\cdot
\tilde{\bm{D}}^{-1/2}\cdot\bm{X}\cdot\bm{W}\right)$.
We know that GNNs which rely solely on local information, like the GCN and its relational variant (R-GCN) (Schlichtkrull et al., ), cannot compute important graph properties, such as girth and diameter or eigenvector centrality Garg et al. ().
We are therefore also investigating more globally aware networks
like GAT (Graph Attention Network) (Veličković et al., ) and GIN (Graph Isomorphism Network) (Xu et al., ).

GAT makes use of self-attention and is thereby more expressive than the GCN. However, its reliance on feature-dependent weights and structure-free normalisation limits its ability to capture specific structural properties that do not directly depend on edges.
This is particularly true for tasks where node features alone are not enough, and global graph structures are crucial (e.g., tasks requiring knowledge of subgraphs or non-local patterns).
GIN aggregates node features in a way that mimics the Weisfeiler-Lehman test for graph isomorphism, and with its strong inductive learning capabilities, it is likely to excel at encoding complex graph properties and solving classification tasks.

### 2.2 Graph properties

Graph theory is a branch of mathematics that studies the properties and relationships of graphs. Graphs can be undirected or directed and analysed through both local and global properties. Local properties like node degree which count the number of connections a node has, identifying highly connected hub nodes, or the clustering coefficient which measures how well a node’s neighbours are interconnected, capturing the local density of connections and giving the node a score, are based on a node with regard to its neighbour. In contrast, global properties such as diameter and characteristic path length assess the overall structure. They indicate how far nodes are from one another and how efficiently information can spread through the network. Global graph properties can be associated with higher level complex systems’ characteristics like the presence of some repeated motifs in the sub structures of the graphs or information-flow properties.

We can distinguish different global properties, basic ones like the number of nodes a graph has, clustering and centrality ones, graph motifs and substructures, spectral and small-world properties.

As an higher-order analysis, the recurrence of specific motifs within network substructures—such as triangles, cliques, or feed-forward loops can be seen as the fundamental building blocks that dictate the system’s functionality and resilience.
Small-worldness (^2^22We are using the Small-World Index, $SWI=\left(\frac{L-L_{l}}{L_{r}-L_{l}}\right)\times\left(\frac{C-C_{r}}{C_{l}-C _{r}}\right)$ in our experiment because it provides a more balanced and robust measure of small-world properties. Unlike the Small-World Quotient: $Q=\frac{C/C_{r}}{L/L_{r}}$, which can be sensitive to network size and degree, $SWI$ normalises both the clustering coefficient and average path length with respect to both random and lattice reference graphs. This dual normalisation approach ensures that $SWI$ is less prone to false positives or negatives, making it a more reliable metric for our analysis Neal ().), as characterised by Barabási Albert & Barabási (), reveal how networks can maintain short path lengths despite their expansive size and sparse connectivity. This kind of higher order properties are very interesting in order to understand how the macroscopic behaviour of complex systems emerges from the intricate interplay of their microscopic components Barabási et al. (). For example how diseases spread in social networks, how neurons interact in the brain, or how information propagates through the Internet.
GNNs synthesise local topological features into global structures, abstract these representations into higher-order graph attributes. Probing their learnt representations should act as a scalable proxy to investigate how global arrangement and connectivity patterns influence a system’s function. In other terms, by dissecting these learned embeddings, we can possibly delve into the intricate relationships between a network’s macroscopic arrangement and its emergent behaviours.

A Graph $G=(V,E)$, $V$ the set of vertices, $E$ the set of edges, can be analysed through both local and global properties. Local properties (like node degree or clustering coefficient) are based on the neighbors of a node.

In contrast, global properties (such as diameter and characteristic path length) assess the overall graph structure. Global graph properties can be associated with higher level complex systems’ characteristics like the presence of repeated motifs in the graphs or information-flow properties. See the appendix B for a list of local and global properties used in our experiments.

GNNs synthesise local topological features into global structures and then abstract these representations into higher-order graph attributes. Probing their learnt representations should act as a scalable proxy to investigate how global arrangement and connectivity patterns influence a system’s function. Based on the message passing paradigm in GNNs, as layers progress, one would expect an increased abstraction in the selection of graph properties. Initially, local features like node degree dominate, but deeper layers progressively capture more global properties, such as connectivity patterns and centrality.

Through hierarchical pooling or readout mechanisms, GNNs can aggregate node embeddings into a single, global graph-level embedding.
Graphs that share structural similarities or patterns of interaction among nodes are organised closely in the embedding space, allowing the model to differentiate between classes of graphs, such as those with and without long paths.

### 2.3 Probing classifiers

In prior work (Hupkes et al., ) probing classifiers have been used for linguistic properties. Here, we adapt them for graph features. Unlike unsupervised techniques such as Principal Component Analysis (PCA) or T-SNE, which are useful to visualise input data with regard to the embedding latent space, we adopt a supervised framework to quantitatively assess how specific properties are encoded within the embedding space of DNNs.
Let $g:f_{l}(x)\mapsto\hat{z}$ represent a probing classifier, used to map the learned intermediate representations from the original model $f$ to a specific property $\hat{z}$. The choice of a linear classifier for $g$ is motivated primarily by its simplicity. If a linear probe performs well, it suggests the existence of a hyperplane in the representation space that separates the inputs based on their properties, indicating linear separability.

Another advantage of a simple linear probe is avoiding the risk that a more complex classifier might infer features that are not actually used by the network itself Hupkes et al. (). While other non-linear probes have been explored in the literature Belinkov (), even studies showing improved performance with complex probes maintain the same logic: $\operatorname{Perf}(g,f_{1},\mathcal{D}_{O},\mathcal{D}_{P})>\operatorname{
Perf}(g,f_{2},\mathcal{D}_{O},\mathcal{D}_{P})$ holds across representations $f_{1}(x)$ and $f_{2}(x)$ when evaluated by a consistent probe $g$. This consistency ensures valid comparison, underscoring that if a property can be predicted well by a simple probe, it is likely relevant to the primary classification task.

From an information-theoretic perspective, training the probing classifier $g$ can be viewed as estimating the mutual information between the learned representations $f_{l}(x)$ and the property $z$. This mutual information is denoted as $\mathrm{I}(\mathbf{z};\mathbf{h})$, where $\mathbf{z}$ refers to the property and $\mathbf{h}$ represents the intermediate representations Belinkov ().

This supervised approach allows us to define hyperplanes or higher-dimensional decision boundaries that partition the embedding space according to the chosen graph property. The $R^{2}$ score serves as this information-theoretic measure indicating how well the hyperplane divides the inputs in the embedding space. A $R^{2}$ near 1 indicates that the embeddings are highly informative about $\hat{z}$, suggesting that the neural model has internalized this property in a linearly accessible manner.

By defining specific properties that could divide the embedding space
and assessing how well the corresponding hyperplanes make the embedding space linearly separable, we gain quantitative insights into the abstract features aggregated within the embeddings.
This method moves beyond mere hypothesis generation based on clustering patterns observed through techniques like PCA, providing a rigorous framework for understanding how well the embedding space represents complex graph properties. It can also be thought as complementary from the T-SNE and PCA visualisation techniques, as it provides a quantitative measure of the separability of the embeddings based on hypothesised properties of interest.

The best illustration of this comes with fig. 3. We illustrate the evolution of the separability of graphs in the embeddings also in fig. 5 using a T-SNE visualisation and the corresponding separability with the properties thanks to the probing.
This highlight the most interesting results of the paper, showing that the separability of house-only and grid-only graphs in the negative class (purple) match with the presence of the property number of triangles in the 5th layer of the GIN architecture.

## 3 Related Work

Existing post-hoc GNN explanations methods can be classified into two main categories: instance-level and model-level methods Barredo Arrieta et al. (). See Agarwal et al. (); Dai et al. () for nice reviews on the subject.
In the realm of instance based methods, gradient-based methods use the gradients of the output with respect to the input or intermediate features to measure the importance of each component of the graph. Decomposition-based methods try to decompose the input graph into smaller subgraphs or paths that can account for the output. Surrogate-based methods use a simpler, more interpretable model to approximate the behaviour of the original GNN and provide explanations based on the surrogate model. And finally Perturbation-based methods which perturb the input graph by removing or adding nodes, edges, or features, and observe the changes in the output to identify the influential components.

The most mainstream technique, GNNExplainer Ying et al. () achieves explanation by removing redundant edges from an input graph instance, maximising the mutual information between the distribution of subgraphs and the GNN’s prediction. It is able to provide an explanation both in terms of a subgraph of the input instance to explain, and a feature mask indicating the subset of input node features which is most responsible for the GNN’s prediction.

For model-based techniques, few methods come to mind Saha et al. (); Azzolin et al. (); Vu & Thai (); Wang et al. (); Xuanyuan et al. (); Yuan et al. (); Zhang et al. (). The most mainstream method seems to be XGNN Yuan et al. (). The authors of XGNN investigate the possible input characteristics used by a GNN for graph classification. But they formulate the problem as a reinforcement learning problem and generate graph patterns iteratively. Such an iterative approach is often intractable for large graphs. Moreover, it does not allow for both node classification and graph classification explanations, nor does it allow for an investigation of the learning process through the different layers of the GNN.

## 4 Datasets

All three datasets have the same setup: given a set of graphs $\left\{\mathcal{G}^{1},\mathcal{G}^{2},\ldots,\mathcal{G}^{N}\right\}$, predict the corresponding binary labels$\left\{y^{1},y^{2},\ldots,y^{N}\right\}$.

##### The Grid-House dataset

inspired by (Agarwal et al., ) is designed to evaluate the compositionality of Graph Neural Network (GNN). It features two concepts: a 3x3 grid and a house-shaped graph made of five nodes. The dataset consists of Barabási-Albert (BA) graphs (Barabási, ) with a normal distribution of the number of nodes. The negative class includes a BA graph connected to either a grid or a house, while the positive class contains a BA graph connected to both a grid and a house (see fig. 2). In order to ensure that the average number of nodes is the same between classes, the number of nodes is a uniformly distributed between
6 and 21 for the grid graphs, between 7 and 22 for the house graphs, and between 1 and 16 when both are present.
During generation, we ensure no test set leakage by removing isomorphisms.
On 2,000 graphs, we perform an 80/20 train/test split.

For accurate classification, models need to identify and combine simple patterns.
Recognizing isolated patterns or single node features is not sufficient.
The dataset helps investigate how GNNs combine multiple concepts and addresses the “laziness” phenomenon, where networks learn patterns characterising only one class and predict the other by default (Longa et al., ).

The dataset has been structured such that an optimal, linearly separable solution requires the combination of local properties, such as eigenvector centrality and betweenness centrality, or the identification of global structural motifs, like counting the number of squares (i.e., four-node cycles).
A random Barabási-Albert graph can’t contain any four-node cycles, while a grid subgraph will consistently exhibit four such cycles. A house subgraph contains exactly one four-nodes cycle and one three-nodes cycle. Therefore, a graph that contains both a grid and a house will have a total of five four-node cycles. The presence of a three-node cycle could help the diagnostic of one type of graph in the negative class but is not necessary nor sufficient for solving the classification problem. On the contrary, counting the number of four-node cycle is necessary and sufficient. Thus, distinguishing between the classes does not really necessitates leveraging centrality-based measures but only recognizing the presence of a specific number of four-node cycles, enabling the model to effectively differentiate between the positive and negative classes. Thus the interesting results of fig. 3.

Figure: Figure 2: Examples of the grid-house dataset. There are graphs with only a house and only a 3x3 grid, these are in one class. Graphs with both a house and a grid are in the other class.
Refer to caption: x2.png

##### ClinTox Molecular

contains molecular graphs representing compounds with binary labels indicating whether they are toxic or non-toxic. The dataset consists of 1,491 drug compounds with known chemical structures. Each molecule is represented as a graph where nodes correspond to atoms and edges to bonds, with node features representing atom types and edge features representing bond types. The task is to predict toxicity.

##### fMRI FC connectomes

consists of two parts. The *Autism Brain Imaging Data Exchange I* dataset contains 528 ASD patients and 571 typically developed (TD) individuals, the *REST-meta-MDD* dataset contains 848 MDD patients and 794 healthy controls. For both, the task is to classify these.
We use the datasets with functional connectivity (FC) graphs, as prepared by  Zheng et al. ().
We perform a 95/5% train-test random split.

In our paradigm, we hope probing functional connectivity matrices (FC) matrices (Farahani et al., ) of neurological disorders (ND) could help explore the link between stuctural properties of the brain’s functional connectivity and neurological disorders such as Autism Spectrum Disorders (ASD) and Major Depressive Disorders (MDD).
Similarly that probed graph properties in toxic molecules align with chemical knowledge.

## 5 Methodology

For each of the three datasets, we use a similar network architecture consisting of a number of GNN (GCN, GIN, or GAT) layers, followed by a pooling operation (mean- (Kipf & Welling, ), sum-(Xu et al., ), or max-pooling (Hamilton et al., )), and then a number of dense layers.
We optimize the hyperparameters to obtain good models for the binary classification task.

For the Grid-House dataset the hyperparameter information can be found in table 4. On this dataset we also compared different regularisation methods.
The explicit $L_{2}$ regularisation encourages the network to keep the weights small, and we expect that this would make the embeddings less sensitive to fluctuations in the input data and smoother. the latter would make them more linearly separable for our probing methods.
Dropout randomly disables a fraction of the neurons during each training iteration which forces the network to learn redundant representations, as any neuron could be dropped out.
These redundant representations might make it more difficult to linearly separate the graph embeddings. We ran each model 20 times and took the one with the best accuracy.

For the ClinTox Molecular dataset, we ranged the number of layers from 4 to 6 and hidden dimensions from 64 to 256. The final model architectures were selected based on optimal performance on the ClinTox dataset.

For fMRI FC connectomes the hyperparamter search space is described in table 15

### 5.1 Probing strategy

Probing is performed on the train and test sets, where train features $\{f_{\text{train}}^{(i)}\}$ and graph properties $\{z_{\text{train}}^{(i)}\}$ are paired for each graph (equally for the test set). Let’s define at least one example for the GCN model. Let $\mathcal{G}^{i}=(A^{i},X^{i})$ denote the $i$-th graph, where $A^{i}$ is the adjacency matrix and $X^{i}$ is the node feature matrix as previously defined. The GCN layers iteratively update the node features $H^{(l)}$ through graph convolutions defined previously as $H^{(l+1)}=\sigma(\hat{A}H^{(l)}W^{(l)})$, where $\hat{A}$ is the normalized adjacency matrix, $W^{(l)}$ are the trainable weights, and $\sigma$ is a non-linear activation function (ReLU). The node embeddings $H^{(l)}$ at each layer $l$ capture both local and global structural information by aggregating features from neighboring nodes. The final node embeddings $H^{(4)}$ are pooled using global max pooling to generate a graph-level embedding $H_{\text{global}}$, which is passed through three fully connected layers to produce the final prediction $\hat{y}$. We deinfe these post pooling operations as $H_{\text{global}}^{(5)}=\sigma(W_{1}H_{\text{global}}),H_{\text{global}}^{(6)}
=\sigma(W_{2}z_{1}),\hat{y}=\text{Softmax}(W_{3}z_{2})$. For probing purposes, we use $H^{(l)}$ at different layers to evaluate node-level properties, while $H_{\text{global}},H_{\text{global}}^{(5)},H_{\text{global}}^{(6)}$, and $\hat{y}$ are used to assess graph-level properties.

We aggregate node embeddings across all graphs to train a single probing classifier for each graph property. For each property, we construct a feature matrix by combining embeddings across all graphs, layer per layer. The classifier $g$ is then trained on this aggregated dataset to predict graph properties $z_{k}^{(i)}$, where $i$ denotes the $i$-th graph and $k$ represents the $k$-th graph property, as defined in table 3. This approach assumes that the relationships between node or graph embeddings and properties are consistent across graphs.

Probing pre-pooling layers to predict global graph properties presents challenges due to the varying numbers of nodes across graphs and the individual states for each node. To handle this, one approach would involve concatenating and flattening the embeddings into a matrix with dimensions $(\text{number of nodes},\text{number of features})$, padding with zeros if a graph has fewer nodes than the maximum in the dataset. However, flattening introduces issues because nodes do not have a canonical ordering; instead, they follow an arbitrary order based on their appearance in the dataset. This inconsistency can undermine permutation invariance, especially since a simple linear classifier applied to the flattened embeddings is not inherently permutation invariant.

To address this, we first sort the embeddings in descending order based on their norms before concatenating, which introduces permutation invariance. Sorting in this way ensures that any padding zeros align at the end of the sequence, enabling learnable representations for graphs with varying node counts. While sorting for permutation invariance is not widely discussed in the literature, it provides a practical solution by using the embeddings’ properties to enforce consistent ordering across graphs.

## 6 Results

### 6.1 Grid-House dataset

The models performed as anticipated thanks to their high expressiveness and the linearly separable nature of the classification problem as we can see in table 5.
The probing results on the Grid-House dataset demonstrate that the number of squares consistently yields the highest $R^{2}$ scores across all models in the global graph embeddings (after pooling aggregation has been applied). This aligns with our initial hypothesis.

In general, higher-layer embeddings filter out many other graph properties as they are less relevant for making the classification problem linearly separable. The GNN’s final layer focuses on the number of squares, effectively partitioning the graphs into two classes: those with $\#\text{squares}<5$ (indicating either the grid or house alone) and those with $\#\text{squares}=5$ (indicating the presence of both substructures). This reduction in feature space through the layers aligns with the model’s goal of optimizing the decision boundary for binary classification, where the number of squares becomes a clear and dominant factor for separability.
Further confirming expectations, density and average path length are also prominent
as the presence of both a house and a grid does slightly increase the average density and path length of graphs.
These findings confirm the correspondence between graph embeddings clustering and property hyperplane separation as shown in fig. 5.

Figure: Figure 3: T-SNE visualisation across different layers of our GIN architecture aligned with the probing $R^{2}$ scores plots (Grid House)
Refer to caption: x3.png

We further observe that, for both the GCN and GIN models, the application of $L_{2}$ regularization yields the expected behavior. The last layer of the GCN in fig. 7 shows a stronger dominance of the number of squares feature when $L_{2}$ regularization is applied compared to when it is not.
Similarly, in the GIN, both number of triangles and density become less detectable relative to the number of squares, by the probing classifier in the final layers under $L_{2}$ regularization, consistent with the anticipated effects on the feature representation.

We observe results consistent with our expectations for the models with dropout. The key property is less dominant, and multiple properties are represented in the final layers. Notably, in fig. 8 the last layer, the separability gap between the #square and the other properties is reduced, indicating a more distributed representation of features when dropout is applied.

When different architectures are compared, the results also align with what we expect from the expressivity of models. For GCN (control), the square detection is strong ($R\textsuperscript{2}=0.77$) in early layers, performance drops slightly in deeper layers, suggesting that GCN captures structural properties early without further refinement. There is less of a presence of #triangle in the control GCN than in the regularised one. The GIN (control) also consistently performs the best on squares ($R\textsuperscript{2}=0.93$) and shows the strong presence of the #triangle before filtering it out in the last layer.
The GIN in general is sharper in the aggregation of global graph properties has it shows results only for the three properties of interest (#square, #triangle, density) before filtering them out in the last layer.
It highlights that GIN excels at global feature detection and effectively isolates and leverages the most relevant structural property for the task, making it sharp in its ability to simplify complex graph data into essential information for decision-making. In other terms, its reliance on minimal yet critical features reflects its capacity for highly targeted feature extraction.

The GAT model stands out by capturing not only squares ($R^{2}=0.88$), but also performing well on other properties like triangles, cliques, and density. GAT assigns a weight to each neighbouring node based on a learned function of the node features, aggregating the neighbours’ information in a weighted manner. This feature-dependent mechanism introduces flexibility but also makes GAT’s performance contingent on the quality and richness of the node features. It seems that GAT’s broader capability compared to the GCN comes at the cost of focus, as GAT tends to incorporate multiple features, which may dilute its ability to pinpoint the most crucial property (in this case, the number of squares) for the classification task. This over-reliance on feature aggregation can lead to inefficiencies when simpler, more targeted properties suffice, as seen with GIN. These results also make sense with regard to the best score obtained by the GIN architecture as seen in table 5.

### 6.2 ClinTox Molecular

As expected GIN outperformed the other models. Based on scores and good properties (inferential mechanism with better expressivity) we focus on the GIN results for the results on other datasets. Detailed results can be found in the appendix, table 8.
When looking into the linear probing performance in table 1, we find that the highest scores are consistently yielded by the average degree, the spectral radius, the algebraic connectivity and the density, in that order.

**Table 1: Linear Probing $R^{2}$ Performance Across GIN Layers for Selected Graph Properties (ClinTox Dataset). Best Scores in Bold; Non-convergence indicated by —**
| GIN Layer | Avg. degree | Spectral radius | Alg. co. | Density | Avg. btw. cent. | Graph energy |
| --- | --- | --- | --- | --- | --- | --- |
| x_global | 0.81 | 0.74 | 0.67 | 0.58 | 0.48 | 0.44 |
| x6 (MLP) | 0.80 | 0.74 | 0.66 | 0.58 | 0.42 | 0.44 |
| x7 (MLP) | 0.75 | 0.71 | 0.56 | 0.50 | 0.47 | 0.46 |
| x8 (MLP) | — | 0.07 | 0.02 | 0.00 | 0.06 | 0.05 |

The average degree of atoms in a molecule provides a straightforward interpretation, as atoms with higher valencies are generally less stable and less biologically compatible. For instance, hydrogen with a valency of 1 and oxygen with a valency of 2 are more compatible with carbon-based molecules, whereas sulfur, with a valency of 6, is less favorable for biological systems  (Komarnisky et al., ). Therefore, the average degree serves as a useful indicator of molecular toxicity. Additionally, the spectral radius, often associated with molecular stability and reactivity, is another valuable graph property. Molecules with a lower spectral radius tend to be more stable, while those with a higher spectral radius may exhibit localized electron densities, increasing their reactivity. Using this property to predict molecular toxicity is a logical approach.
To the best of our knowledge, there is no comprehensive analysis exploring the role of spectral radius in the emergence of molecular toxicity, highlighting an opportunity for future research.

### 6.3 fMRI FC connectomes

In the detailed results, in the appendix table 14, GIN outperforms the other architectures in both parts of the dataset (reproducing the observation by Zheng et al. ().
Again, the strength of GIN lies in its injective aggregation mechanism.
The probing results on the ASD dataset reveal that the number of triangles consistently achieves high $R^{2}$ scores across all models, with particularly strong performance in GIN models. This property is followed by the spectral radius and the density.

As further detailed in the appendix tables 16, 18, 19 and 17 the number of edges is particularly well encoded in the representation of the GAT. This is a consequence of its reliance on feature-dependent weights and structure-free normalisation, which limit its ability to capture specific structural properties that do not directly depend on edges.
The GCN results are broadly comparable to those of the GIN, though they tend to be less precise and selective.

For the MDD results we also focus the GIN model. Detailed results are in the appendix tables 20, 23, 24 and 22.
The probing results MDD reveal that the number of triangles still consistently achieves high $R^{2}$ scores across all models while being less of a distinctive feature than in ASD.
This time, the spectral radius is dominated
by the density of the graph.
In general, the embeddings from the 7th layer of our GIN architecture exhibit higher $R^{2}$ scores for relevant graph properties, suggesting improved separability in the embedding space for MDD classification compared to ASD. This indicates that the learned representations at this depth capture more discriminative structural features, facilitating more effective class separation between MDD and healthy controls.

## 7 Discussion

### 7.1 Expectations

For Grid-House
we hypothesized that the GNN would benefit from leveraging both the local clustering coefficient and eigenvector centrality as node-level features. The first one would help characterize a house, the second a grid. However, neither feature alone is sufficient to render the problem linearly separable. We therefore expect either a combination of features or a single global property (e.g. number of squares) to be leveraged. If the tensor embeddings produced by the GNN can be used to predict these properties, this would indicate that the GNN is utilising them in solving the classification task.

For the ClinTox Molecular dataset, based on the literature Kengkanna & Ohue (); Chen et al. (); Jiang et al. () some few properties have been found to be link with toxicity such as the node degree (i.e. the valency), subgraph patterns (functional groups, chemical fragments), and the overall graph connectivity.

Based on existing literature on functional connectivity (FC) network properties in ASD and MDD, we hypothesized that specific properties will be critical in classifying brain networks for the fMRI FC connectomes dataset. For ASD, we expect betweenness centrality to play a significant role at the node level, reflecting local overconnectivity. At the graph level, we anticipate that clustering coefficient, characteristic path length, and small-worldness will be essential in capturing the local and global network disruptions seen in ASD, particularly the imbalance between local overconnectivity and long-range underconnectivity. For MDD, we hypothesise that increased clustering coefficients, modularity, number of triangles and number of squares will be key features for classification, as they could indicate of heightened local interconnectedness and disrupted global integration.

### 7.2 Findings

We first demonstrate the feasibility of our probing method through the Grid-House dataset. This acts as a proof of concept on probing classifiers plumbing the representations learned by the GNN. We made sure to choose a classification task which requires learning global structural properties of the graph, such as motifs like squares and triangles, and a long range dependencies between those. In line with our expectations, the results show significant dominance of number of squares in the post pooling layers of every model, while still highlighting the superiority of the GIN model in leveraging superior representations when trained to classify graphs with complex motifs like Grid and House. This shows the expressivity of the models and their ability to reduce the complexity of a graph related problem to known graph properties, making the problem linearly separable in the space of their embeddings. These results are consistently higher than those obtained from probing the models with randomized labeling, highlighting the relevance of the initial findings.

Using the ClinTox Molecular dataset to assess molecular toxicity, we explored how key graph properties, such as the average degree and spectral radius, are utilized by our GIN architecture. The average degree, closely linked to atomic valency, reflects a molecule’s potential for interactions. The spectral radius offers a complementary hypothesis, suggesting that the overall structural stability of a molecule, independent of specific atomic features, may also be a key factor in toxicity prediction.
These results suggest that despite the limitations of our approach, it still holds potential for assisting research in complex systems fields, such as neuroscience or social sciences, where emergent phenomena play a crucial role in understanding system dynamics.

Given the previous positive results, we explored a real life applications with the fMRI FC connectomes dataset.
Here, the results provide new insights that extended beyond our initial hypothesis. While we expected betweenness centrality, clustering coefficient, and characteristic path length and Small-worldness to be the most relevant for distinguishing ASD from healthy individuals, the prominence of the number of triangles highlighted the importance of local structural motifs.
This makes sense in the context of functional connectivity, where local overconnectivity in specific brain regions, such as sensory and association cortices, has been observed in individuals with ASD. The strong role of triangle motifs may reflect the tight, redundant local connections that characterize these regions, supporting the hypothesis that local overconnectivity is a key factor in ASD. The spectral radius and density and graph energy being particularly significant is also logical, as these properties are closely related to the overall connectivity strength and the compactness of connections within subnetworks.
The presence of SW is still quite significant in the post-pooling layers which is interesting. In ASD, where global integration is often reduced and local connectivity heightened, these metrics may provide an important reflection of the imbalance between short-range and long-range communication pathways in the brain.

The outcomes for both the ASD and the MDD datasets showed promising results that should be
discussed more deeply with neuroscientists. The results mainly suggest the importance of graph substructures
or spectral and small-world properties over more basic graph properties to explain how graph neural networks
predict these neurological disorders in the FC matrices of patients. The presence of algebraic connectivity
in the last two layers of the GIN MDD probing setting while being completely absent from the GIN ASD
probing one makes it a property of choice to examine for understanding better how this kind of neurological
disorders affect the brain’s connectivity. This is a result which, to the best of my knowledge, has never been
investigated before.

### 7.3 Comparison between datasets

Comparing the ClinTox and the fMRI datasets an interesting observation emerges: basic graph properties (such as the number of nodes, of edges or the average path length) are almost omnipresent in the early layers of the GIN trained on the ClinTox dataset. However, their presence is less pronounced in the GIN trained on the ASD or MDD datasets. This difference offers a clue in distinguishing the complexity of brain-related neurological disorders from the complexity of chemical qualities such as toxicity. This suggests that the emergent properties of the brain may not be as easily tied to simple, differentiable structural features as those seen in molecular systems.

As a confirmation, the types of global graph properties present in the post pooling layers of the GIN-clintox model are of less high level of abstraction than the ones in GIN-MDD or GIN-ASD.
The presence of the average degree, the spectral radius, the algebraic connectivity and the density as accurate explanations for the prediction of toxicity in molecules. The presence of the spectral radius in the last layer of the GIN makes it an even more interesting property to study for toxicity. On the other hand, the presence of motifs should be more investigated in the ASD and MDD datasets with eventually more complex motifs being probed (hexagons constituted of neighbored triangles, house, grid, etc).

### 7.4 Future work

Our methodology has several limitations. While we addressed dataset issues such as leakage and isomorphic graphs, a key challenge remains the lack of guarantees that GNNs find globally optimal solutions, despite their theoretical capacity as universal function approximators. This is particularly evident in fMRI data, where multiple layers of complexity—from MRI limitations and BOLD signal characteristics to Pearson correlation for functional connectivity—introduce noise and inaccuracies. Investigating additional graph properties like girth or complex motifs could be beneficial. Preliminary work on alternative architectures (e.g., GATv2, GraphSAGE, ChebNet, Set2Set, HO-Conv, DiffPool) has begun but is not yet complete.

## 8 Conclusion

We demonstrate the relevance of our model-agnostic explainability method for graph neural networks which probe for GNN graph theoretic representations on the Grid-House dataset. We anticipate any lazy learning bias. We manifest both the expressivity of different GNN architectures and their ability to solve a graph classification problem through optimal feature extraction. They render it linearly separable in the space of their embeddings through the computation of the number of squares in the graph.

That experiment prompted us to investigate both the Clintox Molecular dataset and fMRI FC connectomes dataset and anchored the possibility of formulating hypotheses on the emergent dependence of complex systems qualities to basic and more higher level structural properties.
This kind of higher order properties are very interesting in order to understand how the macroscopic behavior of complex systems emerges from the intricate interplay of their microscopic components.
For example how diseases spread in social networks, how neurons interact in the brain, or how information propagates through the Internet. There is a manifest emergence of molecular qualities like toxicity with regard to their structural properties like node degree (atom valency) and spectral radius (the molecule’s stability). But the complexity of complex systems like the brain makes blurrier the possibility of understanding what affects what as, for example, one could argue that behavioural therapies might influence the brain connectivity as the brain connectivity might influence behavioural qualities of one patient. Echoing this egg-chicken conundrum, the investigation of motifs and spectral properties’ role in neurological disorders like ASD and MDD could allow for promising avenues.

#### Author Contributions

Pelletreau-Duris Tom: Conceptualization (lead); methodology (lead); software (lead); formal analysis (lead); writing – original draft (lead); writing – review and editing (equal).
Ruud van Bakel: writing – review and editing (equal); Supervision
Michael Cochez:Writing – review and editing (equal); Supervision; Funding acquisition

#### Acknowledgments

RvB and MC are in part funded by the Graph-Massivizer project, the Horizon Europe research and innovation program of the European Union grant management number 101093202 [https://graph-massivizer.eu/](https://graph-massivizer.eu/)

## Appendix A R 2 Score

A good $R\textsuperscript{2}$ score gives a sense of how well the features at each layer can be separated linearly to predict the target labels. The second reason is that a more complex probe “bears the risk that the classifier infers features that are not actually used by the network” Hupkes et al. (). Of course, other non linear probes have been explored in the literature Belinkov (). If a few studies observed better performance with more complex probes, the logic remained the same $\operatorname{Perf}\left(g,f_{1},\mathcal{D}_{O},mathcal{D}_{P}\right)>
\operatorname{Perf}\left(g,f_{2},\mathcal{D}_{O},\mathcal{D}_{P}\right)$, of two representations $f_{1}(x)$ and $f_{2}(x)$, holds across different probes $g$. The important criteria is to compare the results obtained by the same measurement system. In general, if we can predict one property on one embedding for a given classification problem, then it means this properly is useful for the problem resolution.

From an information-theoretic perspective, training the probing classifier $g$ can be viewed as estimating the mutual information between the learned representations $f_{l}(x)$ and the property $z$. This mutual information is denoted as $\mathrm{I}(\mathbf{z};\mathbf{h})$, where $\mathbf{z}$ refers to the property and $\mathbf{h}$ represents the intermediate representations Belinkov ().

## Appendix B Local and global graph properties

**Table 2: Local Network Properties with definition and computational criteria**
| Local | Property | Visual Pattern & Definition | Computational Criteria |
| --- | --- | --- | --- |
| Degree | How many links a node has which is the simplest form of centrality | Count edges per node |  |
| Local clustering Coefficient | Are the neighbours of a node also connected together ? | Count triangles of neighbours / total possible triangles of neighbours |  |
| Betweenness Centrality | How much of a bridge between clusters is a node. Removing that node would break many shortest paths. Importance in information flow | Number of shortest paths through node |  |
| Closeness Centrality | Being in the middle of the network, the barycenter of the graph. | The average length of the geodesic distances to all the other nodes (inverse sum of shortest paths) |  |
| Eigenvector Centrality | Being connected to well connected nodes without necessarily having a large number of neighbours itself; influence based on connections | Recursive definition based on neighbours |  |
|  | PageRank | Nodes with important connections; web-inspired importance | Similar to Eigenvector but with random walk and teleportation |

**Table 3: Global Network Properties with definition and computational criteria**
| Global | Property | Visual Pattern & Definition | Computational Criteria |
| --- | --- | --- | --- |
| Number of Nodes | Graph size; total nodes in the network | Count vertices |  |
| Number of Edges | Graph density; total connections in the network | Count connections |  |
| Density | Overall graph connectivity; how densely connected | Ratio of actual to possible edges |  |
| Average Path Length | On average, how close are nodes to each other? Typical distance between node pairs | Average number of steps along the shortest paths for all possible pairs of nodes |  |
| Diameter | Graph span; longest of all shortest paths | Maximum shortest path |  |
| Radius | Graph core; minimum distance from central to farthest node | Minimum eccentricity |  |
| Transitivity | Triangle density; probability of connected node triplets | Ratio of triangles to triads |  |
| Assortativity | Node degree correlations; tendency of similar nodes to connect | Pearson correlation of degrees |  |
| Number of Cliques | Dense subgraphs; count of maximal fully connected subgraphs | Number of maximal complete subgraphs |  |
| Number of Triangles | Local density; fully connected 3-node subgraphs | Count 3-node cliques |  |
| Number of Squares | 4-node patterns; cycles in the graph | Count 4-node cycles |  |
| Largest Component Size | Main connected structure; size of biggest connected part | Largest set of connected nodes |  |
| Average Degree | Overall connectivity; average connections per node | Mean of all node degrees |  |
| Spectral Radius | Dominant graph structure; overall connectivity measure | Largest eigenvalue of adjacency matrix |  |
| Algebraic Connectivity | Graph cohesion; measure of how well-connected the graph is | Second smallest eigenvalue of Laplacian |  |
| Graph Energy | The eigenvalues capture deviations from regularity in the network. Complete graphs or highly connected networks tend to have higher energies due to the larger magnitude of their eigenvalues. In social networks, biology, and communication networks, graph energy can help assess robustness, synchronizability. | Sum of absolute Laplacian eigenvalues |  |
| Small World Coefficient | Balance of clustering and paths; small-world characteristics | Comparison to random graph |  |
|  | Small World Index | Refined small-world measure; comparison to random and lattice graphs | Comparison to random and lattice graphs |
|  | Betweenness Centralization | Central node dominance; degree of central bridging node | Variation in betweenness centrality across nodes |
|  | PageRank Centralization | Influence concentration; degree of dominant influential nodes | Variation in PageRank values across nodes |

### B.1 Artificial dataset

#### B.1.1 Grid House figures

.

Figure: Figure 4: Comparison of different centrality measures for the first graph in our Grid House dataset: (a) betweenness centrality, (b) eigenvector (PageRank) centrality, and (c) local clustering coefficients.
Refer to caption: x4.png

#### B.1.2 Grid House models

.

**Table 4: Range of Hyper-parameters and Final Specification for the Grid-House Dataset**
| Hyper-parameter | Range Examined | Final Specification |
| --- | --- | --- |
| Graph Encoder |  |  |
| #GNN Layers | {$[2,3,4,5]$} | 4 (GCN), 2 (GIN), 3 (GAT) |
| #MLP Layers | {$[2,3,4]$} | 3 (GCN), 2 (GIN), 2 (GAT) |
| Hidden Dimensions | {$[10,15,30,45,60,64,128,256]$} | 60 (GCN), 30 (GIN), 128 (GAT) |
| Attention Heads (GAT) | {$[4,8,16]$} | 8 heads, 32 dimensions per head |
| Learning Rate | {$[1\mathrm{e}-2,1\mathrm{e}-3,1\mathrm{e}-4]$} | $1\mathrm{e}-3$ |
| Batch Size | {$[32,64,128,256]$} | 64 |
| Weight Decay (when added) | {$[1\mathrm{e}-4,1\mathrm{e}-2]$} | $1\mathrm{e}-4$ (GCN), $1\mathrm{e}-2$ (GIN) |
| Batch Normalization | {$with,without$} | $without$ |
| Dropout (when added) | {[0.15, 0.5]} | 0.2 |
| Pooling Method | {$mean,sum,max$} | $max$ (GCN), $mean$ (GIN), $max$ (GAT) |

**Table 5: Performance of Different Models with Regularization on the Artificial Dataset (80%-20% Random Split). The highest performance is highlighted with boldface. All performances are reported under their best settings and rounded to 2 decimal places.**
| Method | Test Accuracy |
| --- | --- |
| GCN (control) | $0.90$ |
| GCN ($L_{2}$) | $0.97$ |
| GCN (dropout) | $0.93$ |
| GIN (control) | $\mathbf{1.00}$ |
| GIN ($L_{2}$) | $0.99$ |
| GIN (dropout) | $1.00$ |
| GAT | $0.97$ |

As expected the RGCN outperform the GCN on this node classification task.

#### B.1.3 Grid House RESULTS

#### B.1.4 Graph properties probing results

.

**Table 6: Linear Probing $R^{2}$ Performance Across models for Selected Graph Properties (GridHouse Dataset). Best Scores in Bold; Non-convergence indicated by —**
| Model | #nodes | #edges | density | avg path len | #cliques | #triangles | #squares | #Largest Component |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| GCN (control) |  |  |  |  |  |  |  |  |
| x_global | 0.36 | — | 0.66 | 0.33 | 0.02 | 0.31 | 0.77 | 0.36 |
| x5 | 0.33 | 0.22 | 0.64 | 0.29 | 0.27 | 0.39 | 0.77 | 0.33 |
| x6 | 0.19 | 0.08 | 0.56 | — | 0.07 | 0.06 | 0.74 | 0.19 |
| x7 | — | — | 0.45 | 0.13 | — | 0.03 | 0.72 | — |
| GCN ($L_{2}$) |  |  |  |  |  |  |  |  |
| x_global | 0.36 | 0.09 | 0.67 | 0.35 | 0.20 | 0.68 | 0.86 | 0.36 |
| x5 | 0.31 | 0.32 | 0.66 | 0.32 | 0.32 | 0.80 | 0.86 | 0.31 |
| x6 | 0.04 | — | 0.41 | 0.15 | 0.03 | 0.23 | 0.83 | 0.04 |
| x7 | — | — | 0.29 | 0.27 | — | 0.09 | 0.81 | — |
| GCN (dropout) |  |  |  |  |  |  |  |  |
| x_global | 0.21 | 0.07 | 0.67 | 0.33 | 0.07 | 0.63 | 0.72 | 0.22 |
| x5 | — | — | 0.59 | 0.26 | — | 0.66 | 0.74 | — |
| x6 | — | — | 0.42 | 0.21 | — | 0.49 | 0.65 | — |
| x7 | — | — | 0.35 | 0.10 | — | 0.26 | 0.51 | — |
| GIN (control) |  |  |  |  |  |  |  |  |
| x_global | 0.12 | 0.07 | 0.50 | 0.32 | 0.07 | 0.22 | 0.87 | 0.12 |
| x5 | — | — | 0.72 | 0.30 | — | 0.89 | 0.93 | — |
| x6 | — | — | — | 0.02 | — | 0.11 | 0.88 | — |
| GIN ($L_{2}$) |  |  |  |  |  |  |  |  |
| x_global | — | — | 0.49 | 0.30 | — | 0.18 | 0.85 | — |
| x5 | — | — | 0.51 | 0.15 | — | 0.52 | 0.89 | — |
| x6 | — | — | 0.40 | 0.12 | — | 0.10 | 0.80 | — |
| GIN (dropout) |  |  |  |  |  |  |  |  |
| x_global | — | — | 0.53 | 0.36 | — | 0.25 | 0.87 | — |
| x5 | — | — | 0.71 | 0.33 | — | 0.85 | 0.93 | — |
| x6 | — | — | — | 0.21 | — | 0.34 | 0.91 | — |
| GAT |  |  |  |  |  |  |  |  |
| x_global | 0.54 | 0.59 | — | 0.49 | 0.61 | 0.89 | 0.87 | 0.54 |
| x5 | — | — | 0.33 | 0.27 | — | 0.17 | 0.64 | — |
| x6 | — | — | 0.25 | 0.17 | — | 0.17 | 0.63 | — |

#### B.1.5 Graph properties probing plots

.

GCN

Figure: Figure 5: T-SNE visualisation across different layers of our GCN architecture aligned with the probing $R^{2}$ scores plots (Grid House)
Refer to caption: x7.png

Figure: Figure 6: Plot of the GCN (control) $R^{2}$ results across different layers probing for graph properties with post pooling layers only (Grid House)
Refer to caption: x8.png

Figure: Figure 7: Plot of the GCN ($L_{2}$) $R^{2}$ results across different layers probing for graph properties with post pooling layers only(Grid House)
Refer to caption: x9.png

Figure: Figure 8: Plot of the GCN (dropout) $R^{2}$ results across different layers probing for graph properties with post pooling layers only(Grid House)
Refer to caption: x10.png

GIN

Figure: Figure 9: Plot of the GIN (control) $R^{2}$ results across different layers probing for graph properties with post pooling layers only (Grid House)
Refer to caption: x11.png

Figure: Figure 10: Plot of the GIN ($L_{2}$) $R^{2}$ results across different layers probing for graph properties with post pooling layers only (Grid House)
Refer to caption: x12.png

Figure: Figure 11: Plot of the GIN (dropout) $R^{2}$ results across different layers probing for graph properties with post pooling layers only (Grid House)
Refer to caption: x13.png

GAT

Figure: Figure 12: Plot of the GAT $R^{2}$ results across different layers probing for graph properties with post pooling layers only (Grid House)
Refer to caption: x14.png

#### B.1.6 Grid House Node properties probing results

.

Using the probing method developed in the next section, we were not fully able to confirm our initial hypothesis.

**Table 7: Linear Probing $R^{2}$ Performance Across models for Selected Node Properties (GridHouse Dataset). Best Scores in Bold; Non-convergence indicated by —**
| GCN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.50 | 0.22 | 0.25 | 0.19 | 0.06 | 0.56 |
| x2 (GCN) | 0.54 | 0.32 | 0.28 | 0.24 | 0.09 | 0.57 |
| x3 (GCN) | 0.54 | 0.35 | 0.29 | 0.25 | 0.11 | 0.57 |
| x4 (GCN) | 0.55 | 0.37 | 0.28 | 0.30 | 0.17 | 0.57 |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | 0.55 | 0.18 | 0.24 | 0.22 | 0.05 | 0.56 |
| x2 (GIN) | 0.52 | 0.34 | 0.27 | 0.25 | 0.07 | 0.54 |
| GAT Layer |  |  |  |  |  |  |
| Layer 0 | 0.55 | 0.07 | 0.05 | 0.32 | 0.28 | 0.17 |
| Layer 1 | 0.52 | 0.48 | 0.08 | 0.31 | 0.30 | 0.14 |
| Layer 2 | 0.47 | 0.55 | — | 0.29 | 0.29 | — |
| Layer 3 | 0.41 | — | 0.14 | 0.19 | 0.26 | — |
| Layer 4 | 0.35 | 0.50 | 0.12 | 0.21 | 0.23 | — |

In these pre-pooling layers, we first observe the predominance of page rank and node degree in the early layers and in all the layers of the GCN and the GIN (which has only two of them). When considering the last layers of the GAT (unfortunately we should have have similar architecture with the GIN in order to fully test our hypothesis) it seems that closeness, node degree and clustering coefficient are the most significant. This aligns with our framing of the graph classification task, which is largely driven by the detection of squares and the fact that pre-pooling layers leading to this property detection should affect mostly these three properties. But this does not align with the use of node properties in a graph in order to do graph classification. This still makes a lot of sense. In general, contrary to the graph probing, and to the exception of the node degree, we see that there is not a single property clearly dominating others but that we go towards a combination of different properties just before the graph pooling method. We would have expect the GIN architecture to show similar results with four layers (as we already see an important increase with regard to the closeness between the first and second layer).

### B.2 Clintox dataset

.

#### B.2.1 model

.

**Table 8: Performance of Different Models on ClinTox with a 80%-20% Random Split. The highest performance is highlighted with boldface. All the performance of methods are reported under their best settings.**
| Method | ClinTox |
| --- | --- |
| GCN | $0.91$ |
| GAT | $0.92$ |
| GIN | $0.93$ |

#### B.2.2 Results

#### B.2.3 Graphs properties probing results

.

**Table 9: Linear Probing $R^{2}$ Performance across the GIN layers for basic graph properties (ClinTox dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GIN Layer | # Nodes | # Edges | Density | Avg. Path Length | Diameter | Radius |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GIN) | 1.00 | 1.00 | 0.66 | 0.76 | 0.55 | 0.60 |
| x2 (GIN) | 1.00 | 1.00 | 0.57 | 0.95 | 0.88** | 0.84 |
| x3 (GIN) | 1.00 | 1.00 | 0.62 | 0.97 | 0.93 | 0.89 |
| x4 (GIN) | 0.99 | 0.99 | 0.37 | 0.91 | 0.82 | 0.82 |
| x5 (GIN) | 0.99 | 0.99 | 0.29 | 0.90 | 0.82 | 0.82 |
| x_global | 0.41 | 0.44 | 0.58 | 0.20 | 0.20 | 0.20 |
| x6 (MLP) | 0.40 | 0.44 | 0.58 | 0.19 | 0.19 | 0.19 |
| x7 (MLP) | 0.42 | 0.46 | 0.50 | 0.27 | 0.23 | 0.25 |
| x8 (MLP) | 0.04 | 0.05 | 0.00 | 0.04 | 0.05 | 0.03 |

**Table 10: Linear Probing $R^{2}$ Performance across the GIN layers for clustering and centrality measures (ClinTox dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GIN Layer | Clustering coef. | Transitivity | Assortativity | Avg. clustering | Avg. btw. cent. | PageRank cent. |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GIN) | — | — | 0.32 | — | — | 0.18 |
| x2 (GIN) | — | — | 0.21 | — | — | — |
| x3 (GIN) | — | — | — | — | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — | — |
| x_global | — | — | 0.25 | — | 0.48 | 0.40 |
| x6 (MLP) | — | — | 0.27 | — | 0.42 | 0.39 |
| x7 (MLP) | — | — | — | — | 0.47 | — |
| x8 (MLP) | — | — | — | — | 0.06 | — |

**Table 11: Linear Probing $R^{2}$ Performance across the GIN layers for graph substructures (ClinTox dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GIN Layer | # Cliques | # Triangles | # Squares | Largest comp. size | Avg. degree | Graph energy |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GIN) | 0.99 | — | 0.00 | 0.99 | 0.53 | 1.00 |
| x2 (GIN) | 1.00 | — | 0.00 | 0.99 | 0.46 | 1.00 |
| x3 (GIN) | 1.00 | — | 0.00 | 0.99 | 0.53 | 1.00 |
| x4 (GIN) | 0.99 | — | 0.00 | 0.99 | 0.20 | 0.99 |
| x5 (GIN) | 0.99 | — | 0.00 | 0.99 | — | 0.99 |
| x_global | 0.43 | — | 0.00 | 0.40 | 0.81 | 0.44 |
| x6 (MLP) | 0.43 | — | 0.00 | 0.40 | 0.80 | 0.44 |
| x7 (MLP) | 0.46 | — | 0.00 | 0.42 | 0.75 | 0.46 |
| x8 (MLP) | 0.04 | — | 0.00 | 0.04 | — | 0.05 |

**Table 12: Linear Probing $R^{2}$ Performance across the GIN layers for spectral and small-world properties (ClinTox dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GIN Layer | Spectral rad. | Algebraic co. | Small world coef. | Small world idx | Avg. btw. cent. |
| --- | --- | --- | --- | --- | --- |
| x1 (GIN) | 0.70 | 0.78 | — | — | — |
| x2 (GIN) | 0.66 | 0.80 | — | — | — |
| x3 (GIN) | 0.61 | 0.80 | — | — | — |
| x4 (GIN) | 0.16 | 0.78 | — | — | — |
| x5 (GIN) | — | 0.69 | — | — | — |
| x_global | 0.74 | 0.67 | — | — | 0.48 |
| x6 (MLP) | 0.74 | 0.66 | — | — | 0.42 |
| x7 (MLP) | 0.71 | 0.56 | — | — | 0.47 |
| x8 (MLP) | 0.07 | 0.02 | — | — | 0.06 |

#### B.2.4 Plots

.

Figure: Figure 13: Plot of the GIN $R^{2}$ results across different layers probing for graph properties. ClinTox dataset (the negative $R^{2}$ values have been reduced to -0.05).
Refer to caption: x15.png

#### B.2.5 Node properties probing results

.

**Table 13: Linear Probing $R^{2}$ Performance across the GIN layers for various node properties (ClinTox dataset). Best Scores in Bold; Non-convergence indicated by —**
| GIN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| --- | --- | --- | --- | --- | --- | --- |
| x0 (GIN) | 0.99 | 0.06 | 0.57 | 0.30 | — | 0.16 |
| x1 (GIN) | 0.85 | 0.12 | 0.51 | 0.31 | 0.00 | 0.20 |
| x2 (GIN) | 0.89 | 0.11 | 0.59 | 0.29 | — | 0.26 |
| x3 (GIN) | 0.86 | 0.07 | 0.51 | 0.28 | — | 0.17 |
| x4 (GIN) | 0.85 | 0.09 | 0.49 | 0.32 | — | 0.14 |

Here again, the very strong presence of the node degree makes a lot of sense when we know this property prepares the aggregation of global properties in the post pooling layers. The interesting thing is the non negligible presence of the betweenness centrality in all the layers which suggests that the betweenness centrality of atoms is important in the aggregation of global molecule properties that help predict the toxicity of a molecule. This property is more than the closeness or the clustering coefficient. The irreplaceable nature of some atoms in the molecular graph, which is literally the meaning of having a high betweenness centrality, is an important feature which makes these atoms targets to be part of higher order molecular schemes and patterns.

Figure: Figure 14: Plot of the GIN $R^{2}$ results across different layers probing for node properties. ClinTox dataset (the negative $R^{2}$ values have been reduced to -0.05). (full results)
Refer to caption: x16.png

### B.3 fMRI datasets

#### B.3.1 models

.

**Table 14: Performance of Different Models on REST-meta-MDD and ABIDE with a 95%-5% Random Split. The highest performance is highlighted with boldface. All the performance of methods are reported under their best settings and round to the second decimal.**
| Method | ABIDE | REST-meta-MDD |
| --- | --- | --- |
| GCN | $0.56$ | $0.61$ |
| GIN | $0.69$ | $0.69$ |
| GAT | $0.62$ | $0.67$ |

**Table 15: Range of Hyper-parameters and Final Specification for FC datasets**
| Hyper-parameter | Range Examined | Final Specification |
| --- | --- | --- |
| Graph Encoder |  |  |
| #GNN Layers | {$[4,5,6]$} | 5 |
| #GIN Layers | {$[4,5,6]$} | 5 |
| #GAT Layers | {$[4,5,6]$} | 5 |
| #MLP Layers (for all models) | {$[2,3,4]$} | 2 |
| #GCN Hidden Dimensions | {$[64,128,256]$} | 128 |
| #GIN Hidden Dimensions | {$[64,128,256]$} | 128 |
| #GAT Hidden Dimensions | {$[64,128,256]$} | 128 |
| #GCN aggregation method | {$[mean,sum,max(pooling)]$} | max pooling |
| #GIN aggregation method | {$[mean,sum,max(pooling)]$} | mean pooling |
| #GAT aggregation method | {$[mean,sum,max(pooling)]$} | max pooling |
| GCN Learning Rate | {$[1\mathrm{e}-2,1\mathrm{e}-3,5\mathrm{e}-4,1\mathrm{e}-4]$} | $5\mathrm{e}-4$ |
| GIN Learning Rate | {$[1\mathrm{e}-2,1\mathrm{e}-3,5\mathrm{e}-4,1\mathrm{e}-4]$} | $5\mathrm{e}-4$ |
| GAT Learning Rate | {$[1\mathrm{e}-2,1\mathrm{e}-3,5\mathrm{e}-4,1\mathrm{e}-4]$} | $1\mathrm{e}-2$ |
| Batch Size (all models) | {$[32,64,128]$} | 32 |
| Weight Decay (alll models) | {$[1\mathrm{e}-3,1\mathrm{e}-4]$} | $1\mathrm{e}-4$ |
| batch normalisation | {$with,without$} | $without$ |
| dropout | {$with,without$} | $without$ |

#### B.3.2 RESULTS ABIDE (ASD) dataset

.

**Table 16: Linear Probing $R^{2}$ Performance across GNN layers for basic graph properties (ASD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | # Nodes | # Edges | Density | Avg. Path Length | Diameter | Radius |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | — | 0.90 | — | 0.21 | 0.13 | 0.07 |
| x2 (GCN) | — | 0.77 | — | 0.22 | 0.24 | — |
| x3 (GCN) | — | 0.62 | — | — | 0.31 | — |
| x4 (GCN) | — | 0.38 | — | — | 0.14 | — |
| x5 (GCN) | — | 0.02 | — | — | 0.09 | — |
| x_global | — | 0.58 | 0.56 | 0.48 | 0.36 | 0.37 |
| x6 (MLP) | — | 0.52 | 0.50 | 0.45 | 0.39 | 0.41 |
| x7 (MLP) | — | — | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | — | 0.94 | — | 0.41 | 0.47 | 0.45 |
| x2 (GIN) | — | 0.55 | — | 0.38 | 0.28 | 0.23 |
| x3 (GIN) | — | 0.25 | — | 0.25 | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | 0.18 | — | — | — | — |
| x_global | — | 0.56 | 0.58 | 0.11 | 0.07 | 0.00 |
| x6 (MLP) | — | 0.58 | 0.66 | 0.14 | 0.10 | 0.09 |
| x7 (MLP) | — | 0.36 | 0.37 | 0.09 | 0.11 | — |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | — | 0.93 | — | — | 0.16 | 0.04 |
| x2 (GAT) | — | 0.89 | — | 0.16 | 0.34 | 0.29 |
| x3 (GAT) | — | 0.84 | — | 0.30 | 0.39 | 0.31 |
| x4 (GAT) | — | 0.78 | — | 0.27 | 0.48 | 0.08 |
| x5 (GAT) | — | 0.67 | — | 0.52 | 0.44 | — |
| x_global | — | 0.74 | 0.70 | 0.60 | 0.29 | 0.40 |
| x6 (GAT) | — | 0.82 | 0.81 | 0.56 | 0.46 | 0.48 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 17: Linear probing performance (R^2 score) across GCN layers for clustering and centrality measures (ASD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | Clustering coe. | Transitivity | Assortativity | Avg. clustering | Avg. btw. cent. | PageRank cent. |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | — | — | — | — | — | — |
| x2 (GCN) | — | — | — | — | — | — |
| x3 (GCN) | — | — | — | — | — | — |
| x4 (GCN) | — | — | — | — | — | — |
| x5 (GCN) | — | — | — | — | — | — |
| x_global | 0.48 | 0.52 | 0.05 | 0.48 | 0.45 | 0.14 |
| x6 (MLP) | 0.33 | 0.30 | — | 0.33 | 0.41 | 0.06 |
| x7 (MLP) | — | — | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | — | — | — | — | — | — |
| x2 (GIN) | — | — | — | — | — | — |
| x3 (GIN) | — | — | — | — | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — | — |
| x_global | 0.19 | 0.04 | — | 0.19 | 0.12 | — |
| x6 (MLP) | 0.23 | 0.08 | — | 0.23 | — | — |
| x7 (MLP) | 0.04 | — | — | 0.09 | 0.11 | — |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | — | — | — | — | — | — |
| x2 (GAT) | — | — | — | — | — | — |
| x3 (GAT) | — | 0.02 | — | — | — | 0.02 |
| x4 (GAT) | — | — | — | — | — | — |
| x5 (GAT) | — | — | — | — | — | — |
| x_global | 0.44 | 0.08 | — | 0.41 | — | — |
| x6 (GAT) | 0.53 | 0.49 | 0.01 | 0.53 | — | 0.08 |
| x7 (GAT) | — | — | 0.00 | — | 0.00 | — |

**Table 18: Linear probing performance (R^2 score) across GCN layers for graph substructures (ASD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | # Cliques | # Triangles | # Squares | Largest comp. size | Avg. degree | Graph energy |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.51 | 0.88 | 0.54 | — | 0.85 | 0.90 |
| x2 (GCN) | 0.27 | 0.81 | 0.58 | — | 0.77 | 0.77 |
| x3 (GCN) | 0.06 | 0.73 | 0.40 | — | 0.64 | 0.62 |
| x4 (GCN) | — | 0.64 | 0.11 | — | 0.30 | 0.39 |
| x5 (GCN) | — | 0.61 | — | — | — | 0.04 |
| x_global | 0.46 | 0.42 | 0.61 | 0.19 | 0.57 | 0.58 |
| x6 (MLP) | 0.42 | 0.34 | 0.35 | 0.31 | 0.51 | 0.52 |
| x7 (MLP) | — | 0.00 | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | 0.58 | 0.95 | 0.69 | — | 0.94 | 0.95 |
| x2 (GIN) | — | 0.91 | 0.12 | — | 0.64 | 0.56 |
| x3 (GIN) | — | 0.74 | — | — | 0.22 | 0.25 |
| x4 (GIN) | — | 0.54 | — | — | — | — |
| x5 (GIN) | — | 0.75 | — | — | 0.23 | 0.17 |
| x_global | — | 0.86 | 0.14 | — | 0.57 | 0.56 |
| x6 (MLP) | — | 0.86 | 0.12 | — | 0.54 | 0.60 |
| x7 (MLP) | — | 0.59 | 0.00 | — | 0.37 | 0.36 |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | 0.58 | 0.86 | 0.66 | — | 0.93 | 0.93 |
| x2 (GAT) | 0.56 | 0.82 | 0.69 | — | 0.87 | 0.89 |
| x3 (GAT) | 0.54 | 0.80 | 0.70 | — | 0.80 | 0.84 |
| x4 (GAT) | 0.50 | 0.75 | 0.74 | — | 0.75 | 0.78 |
| x5 (GAT) | 0.24 | 0.72 | 0.59 | — | 0.71 | 0.67 |
| x_global | 0.32 | 0.56 | 0.40 | — | 0.73 | 0.74 |
| x6 (GAT) | 0.51 | 0.76 | 0.52 | 0.20 | 0.81 | 0.82 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 19: Linear probing performance (R^2 score) across GCN layers for spectral and small-world properties (ASD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | Spectral rad. | Algebraic co. | Small world coe. | Small world idx | Avg. btw. cent. |
| --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.72 | — | — | — | — |
| x2 (GCN) | 0.74 | — | — | — | — |
| x3 (GCN) | 0.56 | — | — | — | — |
| x4 (GCN) | 0.36 | — | — | — | — |
| x5 (GCN) | — | — | — | — | — |
| x_global | 0.46 | 0.43 | — | 0.48 | 0.45 |
| x6 (MLP) | 0.38 | 0.41 | — | 0.39 | 0.41 |
| x7 (MLP) | 0.00 | — | — | — | — |
| GIN Layer |  |  |  |  |  |
| x1 (GIN) | 0.88 | — | — | — | — |
| x2 (GIN) | 0.43 | — | — | — | — |
| x3 (GIN) | 0.25 | — | — | — | — |
| x4 (GIN) | — | — | — | — | — |
| x5 (GIN) | 0.24 | — | — | — | — |
| x_global | 0.76 | — | — | 0.40 | 0.12 |
| x6 (MLP) | 0.74 | — | — | 0.41 | — |
| x7 (MLP) | 0.18 | — | — | 0.23 | 0.11 |
| x8 (MLP) | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |
| x (GAT) | 0.79 | — | — | — | — |
| x2 (GAT) | 0.77 | — | — | — | — |
| x3 (GAT) | 0.02 | — | 0.02 | — | — |
| x4 (GAT) | 0.64 | — | — | — | — |
| x5 (GAT) | 0.49 | — | 0.09 | — | — |
| x_global | 0.58 | 0.20 | — | 0.38 | 0.56 |
| x6 (GAT) | 0.74 | 0.56 | 0.16 | 0.62 | 0.54 |
| x7 (GAT) | — | — | — | — | 0.00 |

#### B.3.3 RESULTS REST-meta-MDD (MDD) dataset

.

**Table 20: Linear probing performance (R^2 score) across GNN layers for basic graph properties (MDD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | # Nodes | # Edges | Density | Avg. Path Length | Diameter | Radius |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | — | 0.90 | — | — | — | — |
| x2 (GCN) | — | 0.85 | — | — | — | — |
| x3 (GCN) | — | 0.71 | — | — | — | — |
| x4 (GCN) | — | 0.64 | — | — | — | — |
| x5 (GCN) | — | 0.03 | — | — | — | — |
| x_global | 0.63 | 0.76 | 0.70 | 0.47 | 0.32 | 0.29 |
| x6 (MLP) | 0.60 | 0.67 | 0.60 | 0.33 | 0.23 | 0.18 |
| x7 (MLP) | — | — | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | — | 0.85 | — | 0.50 | — | — |
| x2 (GIN) | — | 0.67 | — | — | — | — |
| x3 (GIN) | — | — | — | — | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — | — |
| x_global | — | 0.55 | 0.89 | — | — | — |
| x6 (MLP) | — | 0.55 | 0.60 | — | — | — |
| x7 (MLP) | — | 0.74 | 0.77 | — | — | — |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | — | 0.94 | — | — | — | 0.04 |
| x2 (GAT) | — | 0.91 | — | — | — | — |
| x3 (GAT) | — | 0.86 | — | — | — | — |
| x4 (GAT) | — | 0.84 | — | — | — | — |
| x5 (GAT) | — | 0.73 | — | 0.20 | 0.16 | — |
| x_global | 0.52 | 0.80 | 0.74 | 0.29 | — | — |
| x6 (GAT) | 0.62 | 0.76 | 0.69 | 0.43 | 0.18 | 0.26 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 21: Linear probing performance (R^2 score) across GCN layers for clustering and centrality measures (MDD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | Clustering coe. | Transitivity | Assortativity | Avg. clustering | Avg. btw. cent. | PageRank cent. |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | — | — | — | — | — | — |
| x2 (GCN) | — | — | — | — | — | — |
| x3 (GCN) | — | — | — | — | — | — |
| x4 (GCN) | — | — | — | — | — | — |
| x5 (GCN) | — | — | — | — | — | — |
| x_global | 0.42 | 0.34 | — | 0.42 | 0.33 | — |
| x6 (MLP) | 0.35 | 0.33 | — | 0.35 | 0.41 | 0.11 |
| x7 (MLP) | — | — | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | — | — | — | — | — | — |
| x2 (GIN) | — | — | — | — | — | — |
| x3 (GIN) | — | — | — | — | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — | — |
| x_global | — | — | — | — | — | — |
| x6 (MLP) | 0.22 | — | — | 0.22 | — | — |
| x7 (MLP) | 0.43 | 0.33 | — | 0.43 | — | — |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | — | — | — | — | — | — |
| x2 (GAT) | — | — | — | — | — | — |
| x3 (GAT) | — | — | — | — | — | 0.02 |
| x4 (GAT) | — | — | — | — | — | — |
| x5 (GAT) | — | — | — | — | — | — |
| x_global | 0.45 | 0.59 | — | 0.45 | — | 0.24 |
| x6 (GAT) | 0.53 | 0.44 | — | 0.53 | — | 0.16 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 22: Linear probing performance (R^2 score) across GCN layers for clustering and centrality measures (MDD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | Clustering coe. | Transitivity | Assortativity | Avg. clustering | Avg. btw. cent. | PageRank cent. |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | — | — | — | — | — | — |
| x2 (GCN) | — | — | — | — | — | — |
| x3 (GCN) | — | — | — | — | — | — |
| x4 (GCN) | — | — | — | — | — | — |
| x5 (GCN) | — | — | — | — | — | — |
| x_global | 0.42 | 0.34 | — | 0.42 | 0.33 | — |
| x6 (MLP) | 0.35 | 0.33 | — | 0.35 | 0.41 | 0.11 |
| x7 (MLP) | — | — | — | — | — | — |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | — | — | — | — | — | — |
| x2 (GIN) | — | — | — | — | — | — |
| x3 (GIN) | — | — | — | — | — | — |
| x4 (GIN) | — | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — | — |
| x_global | — | — | — | — | — | — |
| x6 (MLP) | 0.22 | — | — | 0.22 | — | — |
| x7 (MLP) | 0.43 | 0.33 | — | 0.43 | — | — |
| x8 (MLP) | — | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | — | — | — | — | — | — |
| x2 (GAT) | — | — | — | — | — | — |
| x3 (GAT) | — | — | — | — | — | 0.02 |
| x4 (GAT) | — | — | — | — | — | — |
| x5 (GAT) | — | — | — | — | — | — |
| x_global | 0.45 | 0.59 | — | 0.45 | — | 0.24 |
| x6 (GAT) | 0.53 | 0.44 | — | 0.53 | — | 0.16 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 23: Linear probing performance (R^2 score) across GNN layers for graph substructures (MDD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | # Cliques | # Triangles | # Squares | Largest comp. size | Avg. degree | Graph energy |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.52 | 0.77 | 0.57 | — | 0.88 | 0.90 |
| x2 (GCN) | 0.58 | 0.84 | 0.69 | — | 0.83 | 0.85 |
| x3 (GCN) | 0.26 | 0.80 | 0.55 | — | 0.72 | 0.72 |
| x4 (GCN) | 0.04 | 0.79 | 0.51 | — | 0.52 | 0.64 |
| x5 (GCN) | — | 0.52 | — | — | 0.01 | 0.04 |
| x_global | 0.54 | 0.76 | 0.50 | 0.62 | 0.73 | 0.76 |
| x6 (MLP) | 0.55 | 0.66 | 0.44 | 0.62 | 0.63 | 0.67 |
| x7 (MLP) | 0.06 | 0.10 | — | — | 0.08 | 0.09 |
| GIN Layer |  |  |  |  |  |  |
| x1 (GIN) | 0.09 | 0.98 | 0.58 | — | 0.86 | 0.85 |
| x2 (GIN) | — | 0.97 | 0.45 | — | 0.48 | 0.67 |
| x3 (GIN) | — | 0.87 | — | — | 0.05 | — |
| x4 (GIN) | — | 0.65 | — | — | — | — |
| x5 (GIN) | — | 0.22 | — | — | — | — |
| x_global | — | 0.91 | — | — | 0.70 | 0.58 |
| x6 (MLP) | — | 0.85 | — | — | 0.67 | 0.54 |
| x7 (MLP) | 0.02 | 0.88 | 0.51 | — | 0.75 | 0.74 |
| x8 (MLP) | 0.02 | — | — | — | — | — |
| GAT Layer |  |  |  |  |  |  |
| x (GAT) | 0.67 | 0.82 | 0.70 | 0.07 | 0.93 | 0.94 |
| x2 (GAT) | 0.59 | 0.83 | 0.81 | — | 0.89 | 0.91 |
| x3 (GAT) | 0.56 | 0.82 | 0.83 | — | 0.82 | 0.86 |
| x4 (GAT) | 0.51 | 0.87 | 0.79 | — | 0.82 | 0.84 |
| x5 (GAT) | 0.45 | 0.83 | 0.26 | — | 0.67 | 0.74 |
| x_global | 0.56 | 0.79 | 0.68 | 0.56 | 0.78 | 0.80 |
| x6 (GAT) | 0.53 | 0.76 | 0.55 | 0.65 | 0.73 | 0.76 |
| x7 (GAT) | — | — | — | — | — | — |

**Table 24: Linear probing performance (R^2 score) across GNN layers for spectral and small-world properties (MDD dataset). Best Scores in Bold; Non-convergence indicated by —(full)**
| GCN Layer | Spectral rad. | Algebraic co. | Small world coe. | Small world idx | Avg. btw. cent. |
| --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.52 | — | — | — | — |
| x2 (GCN) | 0.60 | — | 0.20 | — | — |
| x3 (GCN) | 0.53 | — | — | — | — |
| x4 (GCN) | 0.47 | — | — | — | — |
| x5 (GCN) | 0.07 | — | — | — | — |
| x_global | 0.60 | 0.63 | 0.28 | 0.31 | 0.33 |
| x6 (MLP) | 0.51 | 0.58 | 0.23 | 0.41 | 0.16 |
| x7 (MLP) | 0.04 | 0.04 | 0.00 | 0.00 | — |
| GIN Layer |  |  |  |  |  |
| x1 (GIN) | 0.64 | — | — | — | — |
| x2 (GIN) | 0.52 | — | — | — | — |
| x3 (GIN) | 0.64 | — | — | — | — |
| x4 (GIN) | — | — | — | — | — |
| x5 (GIN) | — | — | — | — | — |
| x_global | 0.75 | 0.32 | 0.44 | 0.39 | — |
| x6 (MLP) | 0.73 | 0.20 | 0.43 | 0.40 | — |
| x7 (MLP) | 0.70 | 0.60 | 0.30 | 0.36 | — |
| x8 (MLP) | — | 0.03 | 0.01 | 0.00 | 0.01 |
| GAT Layer |  |  |  |  |  |
| x (GAT) | 0.70 | 0.02 | 0.00 | — | — |
| x2 (GAT) | 0.66 | — | 0.23 | — | — |
| x3 (GAT) | 0.68 | — | 0.26 | — | — |
| x4 (GAT) | 0.73 | — | 0.30 | — | — |
| x5 (GAT) | 0.66 | — | 0.12 | — | — |
| x_global | 0.68 | 0.63 | 0.18 | 0.52 | 0.04 |
| x6 (GAT) | 0.63 | 0.59 | 0.21 | 0.50 | 0.25 |
| x7 (GAT) | — | — | 0.00 | — | — |

#### B.3.4 plots

.

ASD .

Figure: Figure 15: Plot of the GIN $R^{2}$ results across post pooling layers probing for graph properties ($R^{2}$<0.1 have been hidden). (ABIDE dataset)
Refer to caption: x17.png

Figure: Figure 16: Plot of the GCN $R^{2}$ results across different layers probing for graph properties (ASD)
Refer to caption: x18.png

Figure: Figure 17: Plot of the GIN $R^{2}$ results across different layers probing for graph properties (ASD)
Refer to caption: x19.png

Figure: Figure 18: Plot of the GAT $R^{2}$ results across different layers probing for graph properties (ASD)
Refer to caption: x20.png

MDD .

Figure: Figure 19: Plot of the GIN $R^{2}$ results across different layers probing for graph properties ($R^{2}$<0.1 have been hidden). (REST-meta-MDD dataset).
Refer to caption: x21.png

Figure: Figure 20: Plot of the GCN $R^{2}$ results across different layers probing for graph properties (MDD)
Refer to caption: x22.png

Figure: Figure 21: Plot of the GIN $R^{2}$ results across different layers probing for graph properties (MDD)
Refer to caption: x23.png

Figure: Figure 22: Plot of the GAT $R^{2}$ results across different layers probing for graph properties (MDD)
Refer to caption: x24.png

#### B.3.5 Results ASD and MDD node properties

#### B.3.6 Node properties probing results

.

**Table 25: Linear probing performance ( $R^{2}$ score on the test set) across models for various node properties (ASD dataset). Best Scores in Bold; Non-convergence indicated by —**
| GCN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| --- | --- | --- | --- | --- | --- | --- |
| x1 (GCN) | 0.83 | 0.26 | — | 0.37 | 0.12 | — |
| x2 (GCN) | 0.73 | 0.29 | 0.02 | 0.37 | 0.16 | 0.43 |
| x3 (GCN) | 0.61 | 0.23 | 0.02 | 0.35 | 0.17 | 0.40 |
| x4 (GCN) | 0.53 | 0.19 | 0.03 | 0.31 | 0.17 | — |
| out (GCN) | 0.53 | 0.20 | — | 0.27 | 0.16 | — |
| GAT Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| x1 (GAT) | 0.55 | 0.07 | 0.05 | 0.32 | 0.28 | 0.17 |
| x2 (GAT) | 0.52 | 0.48 | 0.08 | 0.31 | 0.30 | 0.14 |
| x3 (GAT) | 0.47 | 0.55 | — | 0.29 | 0.29 | — |
| x4 (GAT) | 0.41 | — | 0.14 | 0.19 | 0.26 | — |
| out (GAT) | 0.35 | 0.50 | 0.12 | 0.21 | 0.23 | — |
| GIN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| x1 (GIN) | 0.90 | 0.38 | 0.05 | 0.42 | 0.14 | 0.57 |
| x2 (GIN) | 0.89 | 0.24 | 0.12 | 0.40 | 0.16 | 0.59 |
| x3 (GIN) | 0.80 | 0.35 | 0.12 | 0.38 | 0.13 | 0.51 |
| x4 (GIN) | 0.82 | 0.42 | 0.17 | 0.36 | 0.11 | 0.70 |
| out (GIN) | 0.83 | — | 0.13 | 0.30 | 0.13 | 0.70 |

For ASD results, the strong presence of Page Rank is interesting. Regardless of this, without surprise it’s the degree that is consistently the highest node property as it prepare for global properties to aggregate.

**Table 26: Linear probing performance ( $R^{2}$ score on the test set) across models for various node properties (MDD dataset). Best Scores in Bold; Non-convergence indicated by —**
| GCN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| --- | --- | --- | --- | --- | --- | --- |
| Layer 0 | 0.83 | 0.30 | 0.05 | 0.38 | 0.16 | 0.40 |
| Layer 1 | 0.74 | 0.26 | 0.04 | 0.38 | 0.25 | — |
| Layer 2 | 0.69 | 0.31 | 0.03 | 0.41 | 0.23 | — |
| Layer 3 | 0.61 | 0.32 | 0.04 | 0.37 | 0.22 | — |
| Layer 4 | 0.61 | 0.33 | — | 0.37 | 0.19 | — |
| GAT Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| Layer 0 | 0.54 | 0.34 | — | 0.33 | 0.34 | 0.00 |
| Layer 1 | 0.55 | 0.60 | — | — | — | — |
| Layer 2 | 0.48 | 0.40 | — | 0.33 | 0.30 | 0.15 |
| Layer 3 | 0.43 | 0.65 | — | 0.29 | 0.28 | — |
| Layer 4 | 0.39 | — | — | 0.23 | 0.27 | — |
| GIN Layer | degree | closeness | betweenness | eigenvector | clustering | pagerank |
| Layer 0 | 0.92 | 0.54 | 0.09 | 0.40 | 0.23 | 0.58 |
| Layer 1 | 0.82 | 0.53 | 0.06 | 0.29 | 0.16 | 0.45 |
| Layer 2 | 0.83 | 0.43 | 0.16 | 0.34 | 0.18 | 0.60 |
| Layer 3 | 0.73 | 0.37 | 0.13 | 0.34 | 0.16 | 0.47 |
| Layer 4 | 0.86 | 0.24 | 0.20 | 0.26 | 0.11 | 0.47 |

The MDD dataset shows similar results which are surely explained by the same arguments.

## Appendix C Brain Imaging and GNNs

Our brain is a network, more precisely a complex network of functionally interconnected regions specialised in specific cognitive tasks, sharing information with each other. In the last three decades, the field of biological neuroscience and computational cognitive neuroscience have provided and incredible amount of knowledge on the role, function and biological structure of such regions of interests, aiming at better understanding both the biological organisation of the brain (which we can refer to as the ’hardware implementation’), the representation embedded in this hardware and the computational strategy employed to treat this kind of representation Marr (). In other terms, we got better at understanding how each region independently organises itself and processes and forms information (cite Connecting network science and information theory). The main problem for modern computational neuroscience consists of understanding the brain’s plasticity (how regions change over time), the inter-individual differences (how regions specialise differently between people) and how the brain integrates the information (how regions communicate with regard to each other).

For example we understand very well more basic brain structures like the cerebellum due to its high inter-individual similarity but we have a lot more difficulties modelling the prefrontal cortex which is so different from an individual to the other Kanai & Rees (); Gu & Kanai (); Mills et al. (). In other terms, we do understand well the brain operating in segregation but not so much in integration Aine (). Functional segregation refers to the distinct specialisation of anatomical brain regions and functional integration refers to the possible temporal dependencies between the activity of anatomically separated regions of the brain.

Because the representation of a system composed by agents and interactions among them by a complex network is an effective way to extract information on the nature and topology of such interactions, it makes a lot of sense to study the integration of the brain network through its temporal dependencies. Understanding the mathematical properties of such a network with regard to some functional state of the brain network therefore helps understanding how the integration system of the brain and its architecture are linked to ways of processing information. Using Marr’s paradigm to reformulate : understanding the functional communicative structure of the brain network helps understanding its algorithmic footprint. In terms of information theory, we could say that it helps understanding the relationship between topology and dynamics.

One way of accessing the brain activity is to use fMRI imaging. With fMRI measurements at ultra-high-field (3 Tesla, 7 Tesla or even 11 Tesla), hydrogen nuclei present in water and fat molecules align with the scanner’s powerful magnetic field. When radio waves briefly disturb this alignment, the nuclei return to their initial alignment with the magnetic field, this is known as the resonance and causes local changes in the magnetic field. These changes are detected by receiver coils. The collected data from these interactions enable the precise determination of the 3D locations of these events, in the so-called voxels, which can then be visualised. This process underlies the BOLD (Blood Oxygen Level Dependent) response, which is crucial for functional Magnetic Resonance Imaging (fMRI) as it reflects changes in blood flow and oxygenation associated with neuronal activity. We use the magnetic response of blood flow as a proxy for brain activity.

Then, relying on fMRI, we have several ways to study the functional connectivity of the brain. Functional connectivity is defined as the temporal dependence of neuronal activity patterns of anatomically separated brain regions Aertsen et al. (); Friston et al. () and studies have shown that we could study functional connectivity between brain regions as the level of coactivation of functional MRI time-series Lowe et al. (; ). As a result, conceptualising the brain as an integrative network of functionally interacting brain regions offers a powerful framework for understanding large-scale neuronal communication. It provides a method for investigating how functional connectivity and information integration relates to human behaviour and how this organisation may be altered in neurodegenerative diseases Bullmore & Sporns (); Greicius et al. ().

To understand how a specific brain region interacts with others, researchers most often analyse its resting-state activity and use simple pearson correlation of time-series data of a region with the time-series data of all other brain regions, they create a functional connectivity map (fcMap), which visually represents the strength of these connections Biswal et al. (); Cordes et al. (). This is basically a matrix with value and we can understand it as a non relational data structure, in other terms, a graph.

More and more work in cognitive neurosciences explore the link between graph theory and connectomes (functional connectivity matrices) Farahani et al. (). By representing brain regions as nodes and their connections as edges, graph theory provides a powerful framework for analysing the structural and functional organisation of the brain. Notably, studies have begun to explore the link between structural properties of brain connectivity, as captured by connectomes, and the manifestation of neurological disorders such as Autism Spectrum Disorders (ASD) and Major Depressive Disorders (MDD). ASD, characterised by impairments in social communication and repetitive behaviours. MDD is characterised as a mood disorder marked by persistent sadness and loss of interest.

These findings highlight the potential of connectome analysis to elucidate the neurological underpinnings of NDs and pave the way for the development of novel diagnostic and therapeutic strategies. Studying the link between the brain’s Functional connectivity signature and behavioural quality of patients through probing learned embeddings of neural networks trained on classification tasks could thus be a promising avenue to help disentangle the gap between its segregational characteristics and the emergence Johnson (); Eccles (); Wang et al. (); Carroll & Parola () of higher level behavioural quality.

However, if NDs result in alterations in brain functional and structural connections, as well as local and global connections Seeley et al. (); Wang et al. (); Pasquini et al. (); Stam et al. (), traditional deep learning models such as CNN and LSTM are difficult to fit to the connectivity of the brain Zhang et al. (). These long range dependencies, though, are well captured by the relational models defined previously in this thesis : Graph Neural Networks.

Definition : Psychiatric diagnosis can be regarded as a graph classification task. Given an input graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with node feature matrix $X$, GNNs employ the message-passing paradigm to propagate and aggregate the representations of information along edges to generate a node representation $h_{v}$ for each node $v\in\mathcal{V}$ and then explore the modelled human brains using graph methods to extract abnormal brain networks, subnetworks, and local connections Palop et al. (); Thomas et al. ().

Similarly to Zheng et al. (), a GNN can be formally defined through an aggregation function A and a combine function C such that $h_{v}^{(k)}$ is the node embedding of node $v$ at the $k$-th layer and $\mathcal{N}(v)$ is the set of neighbour nodes of $v$:

$$ $\begin{gathered}a_{v}^{(k)}=\mathrm{A}^{(k)}\left(\left\{h_{u}^{(k-1)}:u\in \mathcal{N}(v)\right\}\right)\\ h_{v}^{(k)}=\mathrm{C}^{(k)}\left(h_{v}^{(k-1)},a_{v}^{(k)}\right)\end{gathered}$ $$

In the context of connectomes, many studies have focused on the relationship between general intellectual ability and small-world characteristics in intrinsic functional networks for describing individual differences in general intelligence van den Heuvel & Hulshoff Pol (); van den Heuvel et al. (); Langer et al. (); Hilger et al. (). Better intellectual performance was associated with shorter characteristic path length, the nodal centrality of hub regions in the salience network, as well as the efficiency of functional integration between the frontal and parietal areas Jung & Haier () In general, when connections between specialised brain regions are disrupted, even within localised areas, the result is often functional impairment. This impairment is linked to atypical integration of activity across distributed brain networks Ffytche & Catani (); Catani et al. (). Characterising this impairment through the use of GNN could be one application of our probing pipeline. So far, GNNs have achieved promising diagnostic accuracy on autism spectrum disorder (ASD) Rakhimberdina et al. (), schizophrenia Rakhimberdina & Murata (), bipolar disorder (BD) Yang et al. () and MDD Zheng et al. (). We’ll focus on ASD and MDD. But here as in other graph related fields, research has highlighted we were lacking Interpretability Zheng et al. ().

For ASD, The contribution of rs-fMRI studies based on graph theory for autism exploration is important Redcay et al. (); Rudie et al. (); Di Martino et al. (); Keown et al. (); Kazeminejad & Sotero (). Studies have found increased short-range connections in ASD, particularly within sensory and association cortices. This local overconnectivity may contribute to the sensory sensitivities and restricted interests often seen in ASD. Conversely, long-range connections between distant brain regions tend to be reduced in ASD. This underconnectivity affects integration of information across brain networks. Based on this literature Farahani et al. () we know that the modularity, clustering coefficient, and local efficiency are relatively reduced in ASD (i.e., inefficiency of information transmission in a particular module) while global communication efficiency is increased (shorter average path lengths). As another example, Redcay et al. () observed an increase in betweenness centrality and local connections by analysing the prefrontal brain areas in adolescents with ASD.

In the node property level, we would expect betweenness centrality to be one of the major properties linked with ASD. In the graph level level, we would thus expect the clustering coefficient, the degree to which connected nodes in the brain network are clustered together indicating increased local processing and functional segregation and over-connectivity in local brain regions. We would also expect the characteristic path length to be disrupted, the average shortest path length between all pairs of nodes in the network, suggesting differences in global information transfer efficiency. And small-worldness (SW) which quantifies the balance between local clustering and global integration. Atypical SW in ASD may reflect disrupted optimal network organisation imbalance between local and global processing. We would expect these properties to be critical in our GNNs embeddings trained on classification tasks.

For patients with MDD, several studies have reported topological changes in human brain connectome, including a loss of the small-world network Ye et al. (); Achard & Bullmore ()] and a significant reorganisation of the community structure Zhang et al. (); Leistedt et al. (); Lord et al. (). In general, MDD patients exhibit increased global and local clustering coefficients, indicating a higher degree of local interconnectedness and efficiency in information processing. Moreover, increased modularity in MDD patients indicated that there were relatively less inter-modular edges and more intra-modular edges, which may also be associated with the disruptions in emotion regulation by decreasing communications between the Default Mode Network (DMN) and the Cognitive Control Network (CCN) Ye et al. (). We would thus expect that classifying FC matrices with regard to MDD should use more clustering coefficient, clusterization properties and modularity measures than random (like the presence of motifs like the number small clusters, squares or triangles).

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. For the first sentence in the abstract — 'Graph learning models achieve state of the art performance on many tasks, but so at ever larger model sizes. Accordingly, the complexity of their representations increase.' — I found it hard to follow. As I understand it, this could be rephrased as: 'Graph learning models achieve state-of-the-art performance on many tasks, but this often requires increasingly larger model sizes, which in turn leads to more complex internal representations.'

2. For the images in Fig. 1, it would be better to improve the drawing styles and fonts. For example: the large black right bucket is better to be thiner and move down a little to vertically center alight with the GNN model and Dense layers.
Could the authors elaborate further on how their approach differs from other model-agnostic techniques, and clarify any specific novelty and contributions to the interpretability of intermediate layers?
soundness: 2
strengths: 1. Novelty: Proposes a new model-agnostic explainability approach using probing techniques for GNNs, which is less explored than instance-level explanations.

2. Comprehensive Dataset Testing: Validates the method across diverse datasets, showcasing its versatility.

3. Insightful Comparisons: Provides detailed comparative analysis of different GNN architectures (e.g., GCN, GAT, GIN), noting their strengths in capturing specific graph properties.
summary: This paper presents a model-agnostic approach to interpreting Graph Neural Network (GNN) embeddings by probing for encoded graph properties across various architectures and datasets. Using diagnostic classifiers, the study examines whether these embeddings represent structural properties such as clustering, density, and centrality. The method is tested on synthetic (Grid-House), molecular (ClinTox), and brain connectivity (fMRI) datasets, assessing each GNN’s capability to linearly separate embeddings based on these key properties.
weaknesses: 1. The padding and sorting steps for node embeddings, while intended to address non-uniform node counts, might introduce ordering biases and disrupt permutation invariance.

2. Some limitations in dataset representativity for real-world scenarios, as synthetic datasets like Grid-House may oversimplify the types of structures a GNN may encounter.

3. Fails to fully address the explainability of intermediate GNN layers, focusing more on final representations; probing deeper layers would provide more comprehensive insights into feature evolution.

4. This is my first time reviewing the ICLR paper, but I believe there is a large space to improve this manuscript, especially the presentation of the paper, which is the basic requirement for a good work.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: 1. The reviewer suggests that the findings section should encompass both the experimental settings and the derived conclusions. For instance, while the paper mentions an investigation into the effects of regularization techniques on learning, it fails to present the subsequent conclusions. A more comprehensive analysis, including both the setup and outcomes, would be beneficial.

2. Considering the expressiveness of GNNs, it is anticipated that stronger GNNs would more effectively capture graph properties. However, this may not hold true in practice due to generalization effects. It is recommended that the authors explore GNNs equivalent to 3-WL and consider examining graph transformers to bolster the paper's contributions.

3. Drawing parallels from computer vision research, where the supervision signal plays a pivotal role, it would be insightful to examine if self-supervised models in GNNs also capture structural information more effectively than supervised ones (like the results in CV).

4. While the extensive experiments contribute significantly to the field of Explainable AI (XAI), understanding the reasons behind varying performance levels is crucial. Questions such as why certain layers outperform others, why some substructures are more accurate, whether these findings correlate directly with the labels, and how the type of supervision influences the results, are all critical for a deeper comprehension of the study.
soundness: 3
strengths: 1. The investigation tackles a problem that is both intriguing and significant. To the best of the reviewer's knowledge, this is the inaugural paper to address it.

2. The proposed pipeline is streamlined yet yields effective results.

3. The experimental outcomes provide valuable insights that will benefit future research endeavors.
summary: The paper investigates the explainability of GNNs by proposing a model-agnostic pipeline that uses diagnostic classifiers to probe the embeddings of GNNs for known graph properties. It addresses the challenge of interpreting GNNs due to their complex, non-Euclidean structure and presents findings on how different GNN architectures and regularization techniques affect the representation of graph properties. The study demonstrates that GNN embeddings can indeed reflect graph-theoretic and domain-specific properties, offering insights into the models' internal states and potentially enhancing trust in their predictions. The approach is validated through applications to toxicity and fMRI datasets, confirming the alignment of probed properties with domain knowledge and uncovering previously unexplored structural properties.
weaknesses: This paper has several weaknesses. Addressing Weakness 2 & 3 & 4 will make a compelling case for 'accept.' Specifically:

1. The section designated for findings focuses more on articulating the experimental setting rather than the findings themselves. A clearer exposition of the results will fortify this section. (Question 1)

2. The use of GNNs, which are based on message passing and 1-WL, suggests that the inclusion of stronger GNNs could substantiate the experimental framework. (Question 2)

3. The authors have missed discussing the supervision signal, an important element that merits attention. (Question 3)

4. There is a lack of comprehensive explanations for the experimental results. Elaborating on these would enhance the paper's contribution to the field. (Question 4)

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: How does this work differ from the study "Probing Graph Representations" (Akhondzadeh et al., 2023)? In particular, I would like to focus on the limitations of the previous approach and the new insights provided by this current study.
soundness: 2
strengths: Linear probing is a straightforward yet powerful method, demonstrating strong results across the tested datasets. Its simplicity makes it easy to implement, while still providing valuable insights into the graph properties encoded by the GNNs.

Layer-wise analysis holds significant potential for identifying patterns in GNNs, much like the successes achieved in CNNs. This approach allows for a deeper understanding of how GNNs progressively capture and abstract graph properties across different layers.
summary: The paper explores the types of topological information encoded in node embeddings when solving downstream graph classification tasks. To reveal this information, the authors apply linear probing classifiers to predict various graph properties. The inputs are node embeddings computed at different depths across all graphs in the training data. Some of these properties can be inferred after the global pooling stage and are closely tied to the downstream task.
weaknesses: "We addressed this by sorting the embeddings in descending order and padding with zeros at the end". However, this solution for handling node features before the global pooling stage seems somewhat arbitrary. An immediate alternative approach could be to consider different pooling strategies and even their concatenation. 

Furthermore, Probing Graph Representations (Akhondzadeh et al., 2023) conducted a similar probing study, arriving at nearly identical conclusions. As a result, the amount of novel information introduced by this paper is limited. The authors also fail to compare their probing performance against that of Akhondzadeh et al., 2023, which would have been valuable.

Several key details are missing, such as the performance of the models on the main task and the correlation between task accuracy and probing accuracy. Although the referenced tables are provided in an extensive appendix, the main paper should be more self-contained to allow readers to verify the claims and key points easily.

Additionally, the text in Figures 2 and 3 is difficult to read, which hinders clarity.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. What is the mathematical definition of $R^2$?
2. How are the different local properties discussed in the experiments?
soundness: 2
strengths: 1. This paper discusses explainability from the model-agnostic aspect and explores the effectiveness of local graph properties.
2. This paper provides a comprehensive set of experiments on the local and global properties.
summary: This paper claims that existing explainability methods are mostly instance-based and not well-suited for GNNs. It presents a model-agnostic explainability pipeline using diagnostic classifiers to better understand and trust GNN representations across various architectures and datasets. The findings are interesting and provide unique insights in the XAI for GNNs field. However, the paper has a confusing organization.
weaknesses: 1. This paper is not well-organized. For example, in the introduction, the findings contain some figures and tables, but they can't be found in the main body. Since they are highly relevant to the findings, it would be better to present the results in the main body. The $R^2$ score is not well explained in the main body. In the appendix, there are some descriptions about the $R^2$ score; however, the mathematical formulation is not clear.
 2. Some key points are not clear about the experiments. For example, in line 133, it is confusing when discussing 'defining specific properties.' How to explore different local properties is not clear.
3. The figures are not well presented. For example, in Figure 3, the figures in the first row are hard to understand due to unclear captions. Additionally, it is not clear what $x$, $x_2$, and $x_3$ represent.

Meta-review: This paper proposes a model-agnostic pipeline for interpreting GNN embeddings by employing diagnostic classifiers to probe graph properties. The experiments demonstrate that the probed properties align well with domain knowledge.


However,all the reviewers tend to reject this paper and have pointed out some important weaknesses. First, the novelty of this paper is limited, as similar studies have already been conducted,  which arrive at comparable conclusions. Second, important details are missing. For example, some key points about the experiments are unclear, and there is a lack of comprehensive analysis of the results.


Therefore, I recommend the next version of this paper to incorporate more experiments and insights.

Decision: Reject