Preserving Motif-Relevant Structure in Graph Neural Networks: An Entropy-Aware Study of Aggregation and Depth

1. Introduction

Graphs are powerful tools for modeling real-world systems, especially in biology. Driven by growing interest in network biology, extensive cellular networks – such as transcriptional, metabolic, and protein-protein interaction maps – have already been characterized for numerous model organisms, most notably Escherichia coli and Saccharomyces cerevisiae [53]. In Protein-Protein-Interaction (PPI) networks nodes are proteins and edges correspond to physical bindings between them, forming what is known as the “interactome” [15,19,31,33,50,62]. To reconstruct these interactomes, researchers primarily combine literature-curated data, computational predictions, and systematic experimental techniques like yeast two-hybrid (Y2H) assays and affinity purification mass spectrometry (AP/MS) [41,48,62]. Unlike regulatory networks, the edges in PPIs are undirected, representing mutual, non-causal interactions between protein pairs [42,53]. The local structure of a network around a node can reveal important biological information. These localized wiring patterns are typically quantified by counting how frequently a node participates in so-called graphlets – small induced subgraphs [56]. A related concept, motifs, was introduced in 2002 as small subgraphs that occur more frequently than expected in randomized counterparts [15,31,38,45,47,57]. In contrast to graphlet-based approaches that rely on raw frequency vectors, motif analyses build on statistical significance profiles. Each subgraph is assigned a significance score based on its overrepresentation relative to a null model [37,38,45]. Functional modules such as protein complexes and signaling pathways are encoded as subgraphs within the interactome. These modules mediate enzymatic reactions, signal transduction and cellular decision-making [22,35,37]. However, identifying or counting subgraphs within large networks is computationally challenging: subgraph isomorphism is NP-indeterminate, and determining exact frequency becomes intractable as network size increases [3,20,43].

Motif discovery algorithms can be broadly divided into two main categories: subgraph census (exhaustive search) and subgraph sampling (approximate counting) [45,57]. The subgraph census approach aims to determine the exact frequency of all subgraphs of a certain size k that occur within the input network G. Due to the combinatorial nature of subgraph enumeration and the computational cost of subgraph isomorphism checking, this method becomes infeasible for large graphs. Popular tools implementing exhaustive search include: MAVisto, NeMoFinder and Kavosh [45,57]. Sampling-based algorithms estimate motif frequencies by drawing a large number of random subgraph samples of size k from the input network. Sampling methods are often used when motif detection is applied exploratorily, or when full enumeration is computationally infeasible. Notable implementations include: MFinder [4] and FANMOD [54]. Analytic motif counting methods aim to estimate motif frequencies mathematically. Two major techniques are matrix-based and decomposition-based counting. Matrix-based methods leverage combinatorial properties to derive linear relationships between subgraphs of different sizes [23,45]. The result is a system of linear equations that can be solved using matrix operations. Tools like ORCA [23] and JESSE [36] use such systems to compute graphlet counts efficiently without direct enumeration of all instances. Decompostion-based methods break down each motif into structural primitives, such as shared neighbors, wedge patterns, or overlapping triangles. By computing the frequency of these smaller building blocks, the algorithm infers the total count of larger subgraphs. However, for large biological networks, even state-of-the-art tools can require hours or even days to compute motif frequencies [23,45].

While graph neural networks (GNNs) are powerful tools for general node classification, their theoretical and empirical limitations in capturing local topological structures – such as subgraph isomorphism and motif counting – remain a critical challenge [13,59] Literature suggests that standard message-passing architectures are inherently limited by their equivalence to the Weisfeiler-Lehman (WL) graph isomorphism test, preventing them from reliably distinguishing specific graph patterns or counting induced subgraphs of size three or more [13]. This raises a fundamental question: Can GNNs generate distinct embeddings for nodes with different roles in a subgraph to accurately classify their membership, even if exact counting remains theoretically unreachable? To investigate this, we cast motif and graphlet detection as a multi-label node classification problem using GraphSAGE [21]. Our approach acknowledges the structural limitations of GNNs but seeks to determine if they can serve as reliable “hints” for the presence or absence of graphlets. By using predicted probabilities to guide a rule-based subgraph detector, we bridge the gap between representation learning and deterministic mining, evaluating whether GNNs can provide a useful complement to traditional search. Crucially, we move beyond the coarse homophily-heterophily dichotomy, analyzing the interplay between feature informativeness – ranging from constant scalars to spectral moments – and topological properties.

This work provides a comprehensive empirical investigation into motif-oriented graph representation learning through the following key contributions:

• Framing Motif Detection as a Multilabel Task: We utilize isomorphism-dependent downstream tasks to probe the expressiveness limits of message-passing architectures, testing the ability of GNNs to distinguish node roles within larger structures.
• Systematic Architectural Grid Search: Through a large-scale grid search of 396 configurations, we isolate the impacts of network depth, aggregation operators, and normalization. We demonstrate that performance degradation in deep models is driven by a complex interplay between feature informativeness and topological convergence, rather than oversmoothing alone.
• Fine-Grained Topological Analysis: We analyze how embedding separability and node homophily evolve across neighborhoods, offering a nuanced explanation for performance fluctuations that traditional “oversmoothing” narratives often overlook.
• Transductive vs. Inductive Generalization: By evaluating across both settings, we demonstrate that topology-aware feature engineering (such as spectral moments) enables GNNs to learn effectively even from minimal or dummy feature representations.
• Probability-Guided Hybrid Detectors: Recognizing that GNNs cannot perform exact subgraph matching, we introduce a post-hoc, probability-guided search. This approach leverages GNN outputs to prioritize candidates for rule-based detectors, investigating whether this hybrid pipeline can learn to locate structural patterns efficiently. While our analysis reveals that high label homophily in protein-protein interaction networks limits the reduction of computational complexity, the resulting precision highlights the potential of GNNs as structural indicators for motif mining.

Motif mining typically involves an additional significance testing step by comparing raw counts against randomized versions. In this work, we bypass explicit significance testing, investigating whether GNNs can directly learn representations to guide subgraph search. The remainder of this paper is organized as follows. Section 2 establishes the theoretical framework, detailing graphlet degree distributions, message-passing mechanics, and node-level feature representations, followed by our multi-label classification and homophily formalisms. Section 3 outlines the experimental setup, describing the dataset properties, label binarization scheme, and full hyperparameter search space. Section 4 presents the empirical findings across both transductive and inductive settings, examining the interplay between feature quality, aggregation operators, network depth, and oversmoothing, alongside the performance of our post-hoc probability-guided motif search. Section 5 discusses the broader implications of embedding separability, topological limitations, and computational trade-offs and concludes the study with directions for future research.

2. Theoretical Framework and Modeling Approach

2.1. The Rudiments of Graph Theory

Motif finding inherently requires detecting whether target structures exist within a larger network via subgraph isomorphism. In this section we first establish the foundational graph constructs that govern these matching tasks. A graph $G=(V,E)$ models entities as vertices $V\left(G\right)$ and their pairwise relationships as edges $E\left(G\right)\subseteq \left({\displaystyle \genfrac{}{}{0pt}{}{V}{2}}\right)$, where two vertices $u,v\in V$ are adjacent ($u\sim v$) if $\{u,v\}\in E$ [2]. At the local level, local message passing in GNNs operates over the open neighborhood $N\left(v\right)=\{u\in V\mid \{u,v\}\in E\}$, whose size determines the vertex degree $\mathrm{deg}\left(v\right)=\left|N\right(v\left)\right|$ [2]. To enable matrix-based GNN operations, G is compactly represented via its adjacency matrix $A\in {\mathbb{R}}^{\left|V\right|\times \left|V\right|}$ [2,9], where:

$$A(i,j)=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}{v}_i\sim v_j\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

Building upon this matrix representation, motif discovery fundamentally relies on solving the NP-complete subgraph isomorphism problem [7,45] to identify whether a target motif H appears within G. Formally, $H=\left(V\right(H),E(H\left)\right)$ is a subgraph of $G=\left(V\right(G),E(G\left)\right)$ if $V\left(H\right)\subseteq V\left(G\right)$ and $E\left(H\right)\subseteq E\left(G\right)$. In the literature, motif mining is defined in subtly different ways depending on whether the target patterns are treated as general (non-induced) or induced subgraphs. Formally, given a subset of vertices $W\subseteq V\left(G\right)$, the induced subgraph $G\left[W\right]$ is defined as:

$$E\left(G\left[W\right]\right)=\left\{\{u,v\}\in E\left(G\right)\mid u\in W,v\in W\right\}.$$

To differentiate fine-grained structural positions within these induced motifs, vertex roles are classified via graph symmetries. An automorphism of a graph $G=(V,E)$ is a self-isomorphism, defined as a bijection $\sigma :V\to V$ that preserves edge connectivity:

$$\{u,v\}\in E\iff \left\{\sigma \right(u),\sigma (v\left)\right\}\in E.$$

The set of all such bijections forms the automorphism group $\mathrm{Aut}\left(G\right)$. This group induces an equivalence relation over V, partitioning the vertex set into equivalence classes referred to as orbits [2,45].

In the remainder of this work, we use the terms vertex (plural: vertices) and node interchangeably.

2.2. Graphlet Degree Distribution

Formally, a k-graphlet denotes a small, connected, and non-isomorphic induced subgraph g comprising a subset of k interconnected nodes in G [34,45]. To account for structural symmetries within these subgraphs, automorphism orbits partition the nodes of a graphlet into equivalence classes, grouping together topologically equivalent positions that can be interchanged without modifying the overall graphlet structure [2,45,56]. As the subgraph size k grows, the combinatorial diversity of distinct undirected graphlets alongside their respective automorphism orbits grows exponentially. For instance, while a size of $k=3$ yields merely 2 unique graphlets and 3 orbits, increasing the order to $k=5$ expands this space to 21 graphlets with 58 distinct orbits, eventually scaling to more than two million orbits when considering graphlets of order $k=9$ [23,43].

The graphlet degree vector (GDV) of a node encodes how many times a node touches each graphlet at a particular automorphism orbit [23,43,45,56]. In this context, the node degree can be interpreted as a special case – namely, orbit o(0) [45]. For any given node, the total count for a specific graphlet is the sum of the orbit counts associated with that graphlet [43,56]. These vectors can be aggregated across the graph to form the graphlet degree distribution matrix ${\mathrm{GDD}}_G$, where entry $(i,j)$ indicates the number of nodes that appear i times in orbit j [45,62]. The computational complexity of counting all graphlets of size k is $\mathcal{O}\left(n{d}^{k-1}\right)$, where d is the maximum degree of the graph G. Due to the combinatorial explosion, the computation for graphlets is usually restricted to $k\le 5$ nodes. [11,34,60]. Figure 1 shows graphlets of size up to five nodes.

2.3. Deep Learning on Graphs and Information Diffusion Mechanism

Graph Neural Networks (GNNs) leverage both the graph structure and node features ${\mathbf{X}}_v$ to learn a representation vector for individual nodes ${\mathbf{h}}_v$ or entire graphs ${\mathbf{h}}_{\mathcal{G}}$. Modern GNNs commonly follow a neighborhood aggregation scheme. After k iterations, the node embedding ${\mathbf{h}}_v^{\left(k\right)}$ captures the structure of its k-hop neighborhood [18,20,28,58].

Formally, the general update rule at the k-th layer of a GNN is given by:

$$h_u^{(k+1)}={\mathrm{UPDATE}}^{\left(k\right)}\left(h_u^{\left(k\right)},{\mathrm{AGGREGATE}}^{\left(k\right)}\left(\{h_v^{\left(k\right)}\mid v\in \mathcal{N}\left(u\right)\}\right)\right)$$

or equivalently:

$$h_u^{(k+1)}={\mathrm{UPDATE}}^{\left(k\right)}\left(h_u^{\left(k\right)},m_{\mathcal{N}\left(u\right)}^{\left(k\right)}\right)$$

The AGGREGATE function collects information from the set of neighboring nodes $\mathcal{N}\left(u\right)$, resulting in a combined message $m_{\mathcal{N}\left(u\right)}^{\left(k\right)}$ (Figure 2). Then, this message is passed to the UPDATE function, which combines it with the current embedding of node u to produce the new representation [18,20,28,58,59].

The initial embeddings at $k=0$ are set to the input features for each node:

$$h_u^{\left(0\right)}=x_u\forall u\in V$$

After running K iterations of the message passing process, the final node embedding is obtained as:

$$z_u=h_u^{\left(K\right)}\forall u\in V$$

The neighborhood aggregation operation is fundamentally a set function. We are given a set of neighbor embeddings $\{h_v\mid \forall v\in \mathcal{N}\left(u\right)\}$ and must map this set to a single vector $m_{\mathcal{N}\left(u\right)}$. A GNN must distinguish different multisets of node features – otherwise, it will map structurally different neighborhoods to the same embedding. The more discriminative the aggregation function is with respect to different multisets, the higher the representational power of the GNN [58].

Figure 3 illustrates this issue: although the graph structures around node v differ, simple aggregation schemes like mean or max assign the same embeddings to v and $v^{\prime }$. These functions fail to capture multiplicities in the multiset of neighbor features – they reduce the multiset either to its average or to its maximum element. This leads to an inability to distinguish graphs that differ only in neighborhood cardinality or composition.

2.4. Node-Level Features

Instead of computing embeddings for each node individually, GraphSAGE learns a function that can generate node embeddings by aggregating information from local neighborhoods [14]. This function can generalize to new nodes by leveraging node features. The method can also be applied to graphs without explicit node features, relying on structural signals instead. Node attributes are typically represented as a feature matrix $X\in {\mathbb{R}}^{\left|V\right|\times m}$, where $\left|V\right|$ is the number of nodes and m is the dimensionality of the feature vectors. Basic statistics such as the degree of a node are frequently used as features in traditional node classification tasks. Further informative descriptions can be derived from the centrality measures, which quantify different aspects of a node’s importance within the graph structure [20,27].

Centrality measures estimate how important a node or edge is for maintaining connectivity or facilitating information flow [20,27]. The concept of centrality depends on the specific type of influence under investigation. Three widely used measures include:

• Degree centrality: describes the number of direct connections a node has.
• Betweenness centrality: describes how often a node lies on the shortest paths between other nodes.
• Closeness centrality: defined as the inverse of the average shortest path from a node to all others.

The eigenvector centrality is defined recursively as:

$$e_u={\displaystyle \frac{1}{\lambda }}\sum _{v\in V}A[u,v]·e_v\forall u\in V$$

where $\lambda$ is the leading eigenvalue of the adjacency matrix A, and $\mathbf{e}$ is the corresponding eigenvector representing node centralities [20,27].

Clustering captures the tendency of a node’s neighbors to also be connected to each other, forming triangles [3,8,20,52]. The local clustering coefficient $c_i$ for a node i with degree $k_i$ is defined as the ratio between the number $E_i$ of existing edges between the $k_i$ neighbors and the possible total number ${\displaystyle \frac{k_i(k_i-1)}{2}}$:

$$c_i={\displaystyle \frac{2·E_i}{k_i(k_i-1)}}$$

2.5. Cospectrality and Eigen Decomposition Pruning (EDP)

Two graphs are cospectral if they share the exact same adjacency eigenvalues. While all isomorphic graphs are inherently cospectral, the converse does not necessarily hold [2,55]. Spectral similarity is based on matrix similarity: two matrices $A_1$ and $A_2$ are similar if there exists an invertible matrix P such that:

$$A_2=P^{-1}{A}_{1}P.$$

If $A_1$ and $A_2$ are similar, then they share the same eigenvalues [2,9,16,55].

Verifying cospectrality does not require full eigenvalue decomposition. Instead, it is sufficient to demonstrate that the traces of powers of their adjacency matrices are equal across all walk lengths::

$$\mathrm{tr}\left(A_1^k\right)=\mathrm{tr}\left(A_2^k\right)\forall 1\le k\le n.$$

Here, $\mathrm{tr}\left(A_1^k\right)$ corresponds to the total number of closed walks of length k in the graph. Thus, $G_1$ and $G_2$ are cospectral if and only if for each $k\in \{1,2,\dots ,n\}$, the total number of closed walks of length k in $G_1$ equals that in $G_2$ [2].

Building on spectral properties, Moayed et al. [40] introduced Eigen Decomposition Pruning (EDP) to pinpoint candidate subgraphs in large networks that potentially match a given query motif. The method is based on the following key concepts:

• Neighborhood Pattern $p(u,k)$: is defined as the structural signature of the k-hop neighborhood around a node u.
• Neighborhood Subgraph ${\mathrm{NS}}_G^k\left(v\right)$: the induced subgraph containing all nodes and edges within a distance k from a node v.
• Local Laplacian Matrix (LLM):

  $$L_L\left({\mathrm{NS}}_G^k\left(v\right)\right)=D\left({\mathrm{NS}}_G^k\left(v\right)\right)-A\left({\mathrm{NS}}_G^k\left(v\right)\right),$$
  where D and A denote the degree and adjacency matrices of the induced subgraph, respectively.

Inspired by the EDP method, we approximate the Laplacian spectrum of each node’s 2-hop egonetwork using the ApproxSpectralMoment algorithm proposed by [16,26]. This method estimates the eigenvalue distribution by simulating 10.000 random walks of length 20 and computing the proportion of closed walks, which reflect the spectral properties of the transition matrix. For each egonetwork, we repeat the estimation 20 times and average the results. From the resulting spectra, we extract the first six spectral moments per node. Since the first moment is zero we use moments 2 - 6 as features. These moments are isomorphism-invariant, meaning they remain unchanged for structurally equivalent graphs [16]. We use the intersection of the set of spectral features of the query graph and subgraphs as a clue of similarity.

This spectral approach is further validated by Jin and Zafarani [25], who established formal equivalences between low-order spectral moments and fundamental node attributes (e.g., degree and clustering coefficient) as well as motif patterns such as edges, wedges, triangles, and cycles.

2.6. Multilabel Classification

Multilabel classification is a generalization of binary classification to multiple independent labels [44]. Each instance is associated with a subset of labels $S\subseteq \mathcal{L}$ so the dataset $\mathcal{D}$ becomes:

$$\mathcal{D}=\{(x_1,S_1),(x_2,S_2),\dots ,(x_n,S_n)\}$$

The binary relevance method trains $\left|L\right|$ binary classifiers $C_1,C_2,\dots ,C_{\left|L\right|}$, where each classifier is responsible for predicting the presence or absence of its corresponding label [44,63]. In the multilabel setting, the following F1 variants are commonly applied: micro, macro and instance-based averaging [30,44]. In practice, each label can be assigned a separate decision threshold. Here, we assume a fixed probability cutoff (e.g., 0.5) to make the approach more reproduceble.

2.7. Graph Homophily

Label homophily quantifies the similarity among neighboring nodes based on their label distributions. Graphs are typically categorized as either homophilic or heterophilic [18,39,63].

The node homophily ${\mathcal{H}}_{node}$ is defined as the average proportion of edge-label consistency of all the nodes [18,39,63]. Formally, we have:

$${\mathcal{H}}_{node}={\displaystyle \frac{1}{\left|\mathcal{V}\right|}}\sum _{v\in \mathcal{V}}{\displaystyle \frac{\left|\right\{u\in \mathcal{N}\left(v\right):y_v=y_u\left\}\right|}{\left|\mathcal{N}\right(v\left)\right|}},$$

(1)

where $\mathcal{V}$ is the node set and $\mathcal{N}\left(v\right)$ is the neighbor set of node v. $y_v$, $y_u$ denote the labels of nodes v and u, respectively. Gong et al. [18] extend this definition to the k-hop neighbors to measure high-order homophily:

$${\mathcal{H}}_{high-order}={\displaystyle \frac{\left|\right\{u\in {\mathcal{N}}^{\left(k\right)}\left(v\right):y_v=y_u\left\}\right|}{|{\mathcal{N}}^{\left(k\right)}\left(v\right)|}}$$

(2)

where ${\mathcal{N}}^{\left(k\right)}\left(v\right)=\{u:1\le \mathrm{ShortestPath}(u,v)\le k\}$ represents the k-hop neighbors of v. Both node and edge homophily range from [0,1], and high homophily indicates low heterophily, and vice versa [18].

Zhao et al. [63] argue the basic homophily definition can not be directly used in the multi-label graph datasets, as each node can have more than one label. Further, they define multi-label label homophily as follows:

Given a multi-label graph G, the label homophily h of G is defined as the average of the Jaccard similarity of the label set of all connected nodes in the graph:

$$h={\displaystyle \frac{1}{\left|\mathcal{E}\right|}}\sum _{(i,j)\in \mathcal{E}}{\displaystyle \frac{\left|l\right(i)\cap l(j\left)\right|}{\left|l\right(i)\cup l(j\left)\right|}}.$$

(3)

2.8. Oversmoothing in Graph Neural Networks

When striving for generalization, it is important to account for oversmoothing – a phenomenon where deep GNNs tend to produce almost constant representation, independent of the original node-feature distribution, which can ultimately lead to a loss of discriminative power [5,6,10,12,17,32,46]. Formally, let $H^{(\ell )}\in {\mathbb{R}}^{n\times d}$ be the matrix of node embeddings after ℓ message-passing layers in a GNN, where n is the number of nodes and d is the hidden dimension. Consider a similarity function $\pi :\mathcal{H}\to \mathbb{R}$, where $\mathcal{H}={\mathbb{R}}^{n\times d}$ is the space of all possible embedding matrices. We say that a deep GNN experiences oversmoothing if

$$\underset{\ell \to \infty }{lim}\pi \left(H^{(\ell )}\right)=c,$$

where c is a constant indicating a collapse of embeddings [6,39]. This creates an almost constant representation, independent of the original node-feature distribution [5,6,10,12,17,32,46].

In the literature, drop in test accuracy in deep GNNs is frequently blamed on oversmoothing, following the hypothesis that once embeddings collapse into a trivial space, nodes become hard to separate, leading to a decline in performance [6,10]. To pin this down, recent studies have worked on more precise definitions of oversmoothing using metrics like graph Dirichlet, cosine similarity, or mean-average distance (MAD) [6,46]. To determine how far information needs to propagate, Alon and Yahav [5] define a ’problem radius’ r, which is usually approximated by tuning the number of layers K. In practice, since most nodes in a connected component are within a few hops of each other, GNNs with just two or three layers are often sufficient to achieve competitive results [5,32].

While Chen et al. [12] attribute oversmoothing to the over-mixing of signals – where inter-class interactions introduce harmful noise and impair representation separability – the nature of our task introduces a unique constraint. Because our target labels are inherently tied to the network topology, a high local label density is practically unavoidable. Consequently, we hypothesize that this task requires only a short receptive field, as the underlying graphlet structures span a maximum neighborhood of two hops. Nonetheless, to examine how GNN performance scales with network depth, we evaluate architectures ranging from 2 to 12 layers – extending one layer beyond the diameter of the interactome. We systematically investigate how label homophily varies across different graphlets as network depth increases, and how these structural changes impact performance. In the literature, various architectural techniques – such as residual connections, normalization layers, or gating mechanisms – are commonly used to mitigate oversmoothing and preserve local distinctions [6,17,46]. Building on this, we specifically incorporate batch normalization and dropout into our models to stabilize training and counteract embedding collapse across these varying depths.

3. Experimental Setup

3.1. Dataset Description

In this study, we used the human proteome-scale PPI network published by Thul et al. [50], accessible via the Human Cell Map. The network was modeled as an undirected graph, and for all subsequent analyses, we extracted its greatest connected component (GCC) containing 7185 nodes and 59968 edges. For training and validation in the inductive setting, we employed the BioSNAP human PPI dataset [1], which contains 24 biologically meaningful networks.

3.2. Label Distribution

ORCA provides a 73-component graphlet GDV for each node in the graph. To establish our ground truth, we convert the GDV into a count vector and binarize it to indicate simply whether a node participates in a given graphlet or not. To quantify class prevalence, we define the normalized label frequency by dividing the count of active nodes by the total number of nodes $\left|V\right|$ for each of the 29 labels. Figure 4 illustrates the normalized label frequency for graphlets 1 through 29 (orbit 0 and graphlet 0, corresponding to a simple edge, are omitted). The normalized label frequency exhibits pronounced class imbalance, an inherent consequence of explicitly searching for induced subgraphs.

3.3. Model Architecture and Hyperparameter Space

Each GNN layer in this framework can be defined as:

$$h_v^{(k+1)}=\mathrm{AGG}\left(\left\{\mathrm{ACT}\left(\mathrm{DROPOUT}\left(\mathrm{BN}\left(W^{\left(k\right)}{h}_u^{\left(k\right)}+b^{\left(k\right)}\right)\right)\right):u\in \mathcal{N}\left(v\right)\right\}\right)$$

where $h_v^{\left(k\right)}$ is the k-th layer embedding of node v, $W^{\left(k\right)}$, $b^{\left(k\right)}$ are trainable weights, and $\mathcal{N}\left(v\right)$ is the local neighborhood of v. We stack multiple layers and define the full model configuration as the Cartesian product of the following intra-layer design options:

• Aggregation function: $\{\mathrm{mean},\max ,\mathrm{sum}\}$
• Activation function: $\left\{\mathrm{ReLU}\right\}$
• Batch normalization: $\{\mathrm{True},\mathrm{False}\}$
• Dropout rate: $\{0.0,0.1,0.2\}$
• Hidden dimension: $\{128,256\}$

We evaluate architectures ranging from 2 to 12 layers. This results in a total of $11\times 3\times 1\times 2\times 3\times 2=396$ unique configurations, each of which was trained for 1000 epochs per feature set. To ensure robustness and account for initialization variance, every configuration was evaluated using five different random seeds. All models were trained using the Adam optimizer with a fixed learning rate of 0.001, complemented by early stopping and weight decay. We used GraphSAGE as the base model for its inductive capabilities [14,21]. The simulation environment was implemented in Python, using widely adopted graph learning libraries such as PyTorch Geometric and NetworkX, together with standard packages like NumPy, Matplotlib, and Pandas. For the multilabel classification task, we used the BCEWithLogitsLoss function, which combines a Sigmoid layer and the binary cross-entropy loss in a single class. This enables the model to make independent binary predictions for each graphlet membership per node.

4. Results

4.1. Results Transductive Setting

4.1.1. Grid Search and Effect of Varying Feature Quality

To analyze the mechanics of model depth and feature representations, we conducted a full grid search across 396 architectural configurations over 5 random seeds. Figure 6 reports the mean Macro-F1 performance across all feature sets, where error bars denote the standard deviation across configurations sharing the same design choices. As illustrated in Figure 5, the decline in performance with increasing layer depth is strongly dependent on the choice of the aggregation operator and the utilized feature set. Specifically, our results demonstrate that the max, mean, and sum strategies react differently to increased depth, suggesting that the observed performance degradation is not merely a consequence of the number of layers (i.e., ’excessive smoothing’), but rather a result of the specific mechanism of information aggregation and its interaction with the input features.

In alignment with Figure 5, performance fluctuates heavily with network depth. This variance is reflected in the large standard deviations observed in Figure 6, particularly for deeper architectures where representation collapse occurs. As the network depth increases from 2 to 12 message-passing layers, models consistently suffer from degradation due to overfitting and embedding collapse. However, the severity of this performance drop varies significantly depending on the choice of aggregation operator and feature informativeness:

• Sum Aggregation and Expressiveness Trade-offs: While sum aggregation enhances expressiveness for uninformative or dummy features (e.g., constant scalar) by preserving multiset cardinalities, it becomes inherently unstable with increasing depth. Because the receptive field grows exponentially with node degree, summing unnormalized neighbor features leads to an uncontrolled accumulation of ones, causing training loss to surge drastically and remain stuck at high levels.
• Hidden Dimensions and Feature Similarity: Model performance remains virtually on par across different hidden dimensions (128 vs. 256).
• Dropout Impact: By randomly zeroing out features during message passing, dropout disturbs consistent neighborhood aggregation. Under sum aggregation, this introduces violent amplitude fluctuations, whereas under max aggregation, it severely exacerbates information loss.
• Batch Normalization: The impact of Batch Normalization fluctuates strongly based on the aggregation function. While it acts as a crucial stabilizer that prevents gradient explosion under sum aggregation and mitigates over-smoothing, its interaction with max and mean aggregations shows varying dynamics across hyperparameter choices (see Appendix B for details).

Figure 6. Mean macro F1 across different feature sets. Each subpanel reports the results of a full-grid search over architectural design dimensions: aggregation function, dropout rate, batch normalization, and hidden dimension size. Error bars represent standard deviations across all configurations sharing the same design choice. (a): constant scalar feature set, (b): spectral moments, (c): centralities and clustering coefficient.

Figure 6. Mean macro F1 across different feature sets. Each subpanel reports the results of a full-grid search over architectural design dimensions: aggregation function, dropout rate, batch normalization, and hidden dimension size. Error bars represent standard deviations across all configurations sharing the same design choice. (a): constant scalar feature set, (b): spectral moments, (c): centralities and clustering coefficient.

To analyze the underlying mechanics of oversmoothing, we examine how feature distributions evolve across network layers. Formally, a GNN reaches a state of complete representation collapse ($\mu \left(X^{\left(l\right)}\right)=0$) when its covariance matrix contracts to constant entries, $K_{\alpha \beta }^{\left(l\right)}=k^{\left(l\right)}$ [17]. Under this condition, nodes receive identical pre-activations, ultimately collapsing into uniform embeddings ($x_{\alpha }^{\left(l\right)}=x_{\beta }^{\left(l\right)}$). As demonstrated in prior work, this uniform state acts as an invariant fixed point of the network dynamics [17]. Driven by strong layer-wise contraction properties, this fixed point behaves as a powerful attractor: even when starting from distinct initial features ($d(x_{\alpha }^{\left(0\right)},x_{\beta }^{\left(0\right)})\ne 0$), node representations rapidly converge toward zero feature distance ($\mu \left(X^{\left(l\right)}\right)\to 0$) as network depth increases [17].

To empirically evaluate this behavior, we recorded the final node embeddings across all evaluated network depths ($2\to 12$) and computed their pairwise cosine similarities (Figure 7). While several configurations exhibit a steep increase in similarity toward higher depths – indicating severe embedding collapse – our findings show that this effect is highly sensitive to model architecture. In particular, Batch Normalization acts as a vital stabilizer across all aggregation schemes: by re-scaling and re-centering features at each step, it continuously reintroduces variance into the feature space, thereby suppressing oversmoothing and maintaining layer-wise feature separability.

Building upon this theoretical framework, our grid search reveals an intriguing interplay between feature informativeness and aggregation methods:

• Constant Scalar: homogeneous initializations, such as a constant value of 1 assigned to every node, places the network almost directly onto the fixed point, maximally accelerating the collapse of representations.
• Spectral Moments & Centralities: Spectral moments derived from 2-hop ego-networks exhibit highly overlapping neighborhoods due to the interactome’s high clustering coefficient.
• Centralities and Clustering Coefficient: centrality measures paired with the clustering coefficient provide node-specific attributes, offering the highest degree of distinguishability and performing best.
• Aggregation Function:max-pooling proves to be the weakest and most unstable aggregation function, consistently characterized by aggressive information loss and a rapid collapse of representations. sum aggregation acts as a strong contraction driving the network toward the attractive fixed point, though it can theoretically prevent collapse when paired with proper weight variance tuning (${\sigma }_w^2$) [17]. In contrast, mean aggregation excels when paired with informative input features, acting as a stabilizing damper through constant degree normalization.

Finally, introducing dropout modulates this behavior significantly depending on the chosen aggregation scheme. Dropout acts as a regularizer but, in deep GNNs, leads to stochastic feature dropping that accelerates embedding collapse:

• Under sum aggregation, which scales neighbor messages directly by the node degree, dropout breaks the consistent accumulation of signals, resulting in violent amplitude fluctuations and a rapid drift into the attractive fixed point.
• max aggregation, inherently preserves only the dominant neighborhood signal while discarding all alternative features. Consequently, random feature dropping severely compounds this information bottleneck, depriving the network of sufficient variance and causing rapid representation degeneration.
• Conversely, mean aggregation remains robust against dropout disturbances. By dividing by the node degree, it compensates for missing signals and prevents the abrupt breakdown of the information flow, as further evidenced in Figure 7.

For a detailed breakdown of results across varying dropout rates and hidden dimensions, we refer to Appendix B.

Figure 8 reports the graphlet precisions averaged over 5 train/val/test splits, highlighting clear performance variations across both feature sets and aggregation mechanisms. The anticipated failure cases are prominently reflected in the constant scalar setting, where cycles (e.g., $G_2,G_5,G_{15}$) and clique-graphs ($G_8,G_{29}$) exhibit vanishing precision. Theoretically, these non-isomorphic topologies cannot be distinguished through local feature distributions alone (Figure 3), as they differ strictly in neighborhood cardinality and composition.

Under uninformative inputs, less expressive, non-injective aggregators (e.g., mean or max pooling) map these distinct local neighborhoods onto identical vectors, inducing severe embedding collapse. In contrast, sum aggregation preserves multiset injectivity by accumulating unscaled neighbor signals, thereby retaining neighborhood cardinality and enabling the model to disambiguate subtle structural differences. Beyond these expressiveness limitations, a comparison with the label distribution (Figure 4) reveals that class imbalance further amplifies underperformance, as graphlets with sharply declining precision values are heavily underrepresented in the interactome.

4.1.2. Effect of Label Homophily

A complementary explanation for the precision fluctuations in Figure 8 is oversmoothing, caused by information overmixing across layers. To analyze this, we extend our evaluation to higher-order homophily based on Gong et al. [18], adapting their formulation to our multi-label setting by computing the metric independently for each label across k-hop neighborhoods ($k\in \{1,\dots ,11\}$). As shown in Figure 9, this reveals distinct topological patterns: graphlets like $G_2$, $G_5$, and $G_{15}$ show low local homophily ($k=1$) that decays rapidly over higher hops. In these regions, nodes aggregate signals from neighbors with different labels, blurring feature boundaries. For graphlet classification, this is worsened by non-induced neighbor connections: nodes within an induced graphlet often connect to non-induced variants of the same motif, pulling in signal from “forbidden” structures. Combined with their low representation in the dataset (Figure 4), this lack of neighborhood consistency and structural noise prevents the model from separating classes effectively, causing sharp precision drops for minority graphlets in deeper networks.

4.2. Results Inductive Setting

4.3. Grid Search and Effect of Varying Feature Quality

To evaluate the model’s ability to generalize to unseen graph structures, we adapted the experimental design for the inductive setting by constraining the network depth to 2 – 4 layers. Table 1 shows the condensed configuration space for the inductive setting. The dataset comprises a heterogeneous collection of BioSNAP and cell-compartment graphs, varying significantly in size and density (see Table A1). To ensure robust generalization, we stratified graphs into three size categories (small, moderate, large) and performed a random 17/8/8 train/val/test split.

Figure 10 presents the hyperparameter landscape across the configuration space. Compared to the transductive setting, the performance profiles appear markedly more stable, as the constraints on network depth effectively prevent the inclusion of severely oversmoothed configurations.

While the qualitative preferences – such as the inherent instability of max aggregation and the sensitivity to feature informativeness – persist, the inductive performance highlights two key observations:

• Optimal Depth: all three feature sets maintain their strong preference for shallow, 2-layer architectures, confirming that depth optimization is an intrinsic property of the feature set rather than a setting-specific artifact.
• Aggregation Dynamics: While sum and mean aggregations perform nearly on par when utilizing informative Spectral Moments and Centralities, the sum operator consistently outperforms mean for Constant Scalar inputs.

Ultimately, the inductive results corroborate the findings from our transductive analysis: by limiting message-passing to a moderate depth, the GNN avoids the performance degradation induced by the attractive fixed point of vanishing feature distance.

The label-wise precision plot in Figure 11 reveals a marked reduction in the performance drops previously observed in the transductive setting. Benefiting from a substantially larger and more diverse pool of training data points, the GNN achieves smoother and more consistent precision profiles across label classes. Nevertheless, the Constant Scalar features continue to exhibit characteristic performance drops across specific failure cases, heavily dependent on the chosen aggregation method. A potential explanation for the observed weaknesses of Spectral Moments comes from the local nature of their computation. Because Spectral Moments are derived from bounded subgraphs (e.g., 2-hop ego-networks), adjacent nodes share a substantial portion of their local neighborhood. Consequently, neighboring nodes receive highly correlated initial feature representations. This spatial overlap renders the representations particularly vulnerable to accelerated oversmoothing during message passing.

4.3.1. Effect of Label Homophily

Table 2 outlines the structural characteristics and Macro-F1 performances across the test graphs, highlighting three core insights into inductive generalization:

• Stable Generalization under High Homophilie: The BioSNAP graphs exhibit exceptionally high label homophily ($\approx 0.96--0.98$) alongside very low Jensen-Shannon divergences (JSD). Under such extreme homophily, message passing acts primarily as a smoothing operator over aligned labels, mitigating misclassifications from oversmoothing. The minimal Jensen–Shannon divergence confirms that local neighborhood features remain topologically consistent This structural consistency explains why Macro-F1 scores remain constantly high ($\approx 0.97$).
• Impact of Graph Scale and Distribution: While larger BioSNAP networks display near-perfect homophily and homogeneous label distributions, smaller or more structurally diverse graphs – such as Golgi ($\left|V\right|=192$) and Vesicles ($\left|V\right|=585$) – exhibit a noticeable drop in homophily and JSD values. This structural divergence correlates directly with lower Macro-F1 scores, particularly when paired with uninformative Constant Scalar features.
• Robustness of Informative Features: Even in structurally challenging environments, such as the Golgi compartment with a homophily of $0.81$, Spectral Moments and raw node features maintain remarkable stability. They significantly outperform the Constant Scalar baseline, underscoring that informative feature-engineering effectively mitigates performance losses induced by low homophily.

4.4. Graphlet-Matching Guided by Predicted Probabilities

To complement the GNN predictions, we implement a post-hoc motif detection strategy that explicitly counts motif instances. While classification models indicate whether a node participates in a specific motif, they do not quantify the exact frequency or local context of its occurrence. To address this limitation, we developed custom detectors for each of the 29 graphlets. Obtaining a reliable ground-truth count requires aggregating node-specific information; to this end, we leverage ORCA, which provides the frequency of a node v appearing in a given orbit $o\left(x\right)$. Through combinatorial transformations, these orbit counts are converted into graph-wide totals: one can either sum the asymmetric orbit counts across all nodes or utilize symmetric orbit counts adjusted by the number of unique symmetric orbits within the respective graphlet. We restrict our evaluation to $k\in \{3,4\}$ graphlets, as detectors repeatedly exceeded time limits for $k=5$.

Detection is simulated for each feature set using its best-performing configuration $\{\mathrm{Layers},\mathrm{Aggregation}.,$$\mathrm{Batch}\mathrm{Norm}.,\mathrm{Dropout},\mathrm{Hidden}\mathrm{Dim}.\}$ based on Macro-F1 score:

• Constant Scalar: $\{2,\mathrm{sum},\mathrm{False},0.0,128\}$
• Spectral Moments: $\{2,\mathrm{sum},\mathrm{False},0.1,256\}$
• Centralities and Clustering Coefficient: $\{4,\mathrm{mean},\mathrm{True},0.1,256\}$

For each graphlet, target nodes exceeding a classification probability threshold of $0.5$ are assigned an anchor role (orbit) based on their local neighborhood. We then extract their 1-hop neighborhood and filter out candidate nodes that violate structural constraints by containing “deal-breaker” edges – connections explicitly prohibited within the target graphlet’s induced subgraph definition (see Figure A4). For a detailed overview of the absolute errors across feature sets and test graphs, see Tables A3–A5 in the Appendix.

Table 3 presents the ground-truth counts computed via ORCA alongside the absolute deviation from our detected counts. Consistent with Figure 8, the absolute differences indicate that graphlets with lower precision were omitted more often; however, these deviations are several orders of magnitude below the true counts.

To limit computational complexity, our workflow replaces exhaustive search with a targeted orbit-matching strategy, identifying positive nodes only when positioned on detector-suited orbits. Integrating an analytical counting algorithm, a GNN, and heuristic detectors, we search for graphlets across 21 test graphs. Missed graphlet instances primarily stem from misclassified orbits that serve as anchor nodes for the detector (see Figure A4).

Although we initially intended to use GNN-predicted probabilities to bypass exhaustive mining, high local homophily in PPIs causes probability-guided search to closely mimic exhaustive search in practice. Nevertheless, because our detector selectively targets nodes lacking prohibited “deal-breaker” edges, subtle runtime gains may still exist. A formal time-complexity comparison remains beyond the scope of this work. Furthermore, formal statistical significance testing to isolate true motifs is not viable under this setup, as the remaining predictions preserve inherent uncertainty. These directions are left for future research.

5. Conclusion

In this study, we investigated the capabilities and limitations of GNNs in solving subgraph-related tasks. We demonstrated that while GNNs exhibit strong performance in common tasks such as node classification or link prediction, their expressiveness falls short when it comes to more complex structural reasoning tasks, such as subgraph isomorphism or graphlet counting. Ultimately, learning such induced substructures is governed by a multifaceted interplay rather than a single bottleneck. By leveraging foundational metrics like homophily and initial embedding separability, we map out these performance variations through the following key takeaways.

• Oversmoothing and Feature Separability. Oversmoothing remains a heavily debated topic in GNN literature; however, low accuracy in deep GNNs cannot be attributed to representation collapse alone. Instead, initial node embedding separability plays a decisive role alongside training pathologies such as vanishing gradients or overfitting [6]. Remedies like skip connections, normalization, or gating are standard architectural tools designed to preserve local distinctions [6,17,46]. In our experiments, oversmoothing emerges as an intricate interplay between feature informativeness, network depth, and the configuration space. Notably, higher initial feature similarity drastically accelerates the onset of embedding collapse, shifting optimal performance toward shallower architectures.
  Interplay of Depth, Topology, and Homophily We evaluate a non-trivial task where feature distributions across classes partially or completely overlap, placing structural topology at the center of the learning objective. In a short-range task characterized by high label homophily, we target induced structures within local neighborhoods, yielding correlated label distributions. By extending maximum network depth beyond the network diameter, we explicitly analyze performance decay under deep aggregation. Due to the small-radius nature of graphlet detection, shallow networks (2-3 layers) consistently achieve peak performance. However, classifying strictly induced structures introduces unique topological constraints: mutual exclusivity does not protect the model from information overmixing as depth increases. While classes with high local homophily build robust representations, low-homophily graphlets suffer from severe overmixing across inter-class and non-induced edges. As a result, the model aggregates signals from “forbidden” neighbor contexts, causing rapid oversmoothing and degraded precision for minority classes in deeper layers.
  Impact of Feature Quality and Network Models. Utilizing domain-agnostic or dummy feature sets (e.g., a constant scalar) can enhance model generalizability, but it tightly couples performance to graph topology. In our inductive evaluations, this topological dependency became particularly evident: smaller, more densely clustered networks with higher Jensen-Shannon divergence consistently yielded significantly weaker results than medium and large graphs. Remarkably, in these medium-to-large networks, the constant scalar feature performed almost on par with more expressive feature sets, as neighborhood summation effectively allows the model to recover structural proxies like node degree. However, relying purely on constant features limits inductive transferability and leaves scalability questionable, as unnormalized summation over expanding receptive fields introduces size and magnitude effects in larger graphs. Furthermore, synthetic network models – such as Erdős-Rényi or Watts-Strogatz graphs – possess fundamentally different degree distributions, clustering behaviors, and graphlet counts compared to biological networks [3,8]. Consequently, topology-dependent features may struggle to generalize across distinct graph families.
  Aggregation and Regularization Mechanics. We propose a full-batch approach operating without neighborhood sampling that maintains strong discriminative capacity even under dummy feature regimes. While GraphSAGE’s sampling strategy could potentially mitigate size-related effects and enhance generalizability in an inductive setting, feature informativeness introduces a non-negligible time-performance trade-off. A crucial interplay emerges between network topology, feature expressiveness, and the choice of aggregation function (mean, max, and sum). Because our models operate on very small feature sets, traditional regularization techniques exhibit counterproductive dynamics. Specifically, applying dropout on low-dimensional input vectors disturbs neighborhood aggregation consistency and severely accelerates oversmoothing. Conversely, Batch Normalization acts as a vital stabilizer: by mitigating feature alignment and managing amplitude explosions under sum aggregation, it counteracts the initial high feature similarity and preserves embedding separability across layers.
  Substructure Localization and Post-Hoc Search. While GNNs excel at similarity-based tasks like node classification, they struggle with explicit substructure localization and counting. Existing methods typically reduce subgraph matching to binary node classification, labeling query nodes as +1 and remaining nodes as -1, which localizes substructures but fails to resolve exact node-to-node correspondences [29]. We address this limitation using post-hoc, probability-guided detectors. A key advantage of our proposed pipeline is its dynamic design. PPIs are inherently dynamic, continually evolving – either through newly discovered interactions or by merging sources (e.g., HPRD, BioGRID, STRING) to complete the interactome. At the same time, handling PPI data requires caution: techniques like yeast two-hybrid screening carry high false-positive rates. To ensure reliable downstream hypothesis testing, researchers must often filter out these spurious edges, further altering the graph topology [48]. To handle these topological changes, we selected GraphSAGE, an inductive GNN framework. Dynamic node additions from database integration or removals from quality filtering are easily accommodated, as the model generates new predictions on the updated graph without requiring retraining. To avoid an exhaustive search across the entire network, predictions can be restricted to a local ego-network around the affected nodes. However, high label homophily in real-world interactomes causes probability-guided searches to closely mimic exhaustive searches in practice, offering limited relief at high label densities. Although our hybrid detector does not yield a general-purpose complexity breakthrough for arbitrary graph types, it successfully transforms brute-force motif mining into a targeted refinement process with high precision in inductive settings. Future work should evaluate synthetic benchmarks with controlled label densities to formally delineate the time-complexity boundary between probability-guided and exhaustive search strategies.
  Expressivity Trade-offs and Open Challenges. Importantly, we showed that current GNN architectures still rely heavily on heuristic supervision or ground-truth annotations for tasks like subgraph matching. No automatic matching strategy has yet emerged that can reliably recover isomorphic structures in a general setting. Our inclusion of spectral, centrality-based, and constant-scalar features provided a thorough feature-based benchmarking, underscoring the trade-offs between model depth, expressivity, and task-specific performance. Our experiments revealed that even with minimal or noisy input features, GNNs are capable of learning useful representations. However, challenges such as over-smoothing, limited receptive field adaptation, and the inability to produce fine-grained node-to-node mappings persist.

In conclusion, while GNNs remain a powerful tool in graph representation learning, their application to motif-centric tasks requires further architectural and theoretical advances.