EpitopeGNN: A Graph Neural Network for Influenza A Virus Hemagglutinin Subtype Classification Based on 3D Structure

1. Introduction

Influenza A virus causes seasonal epidemics and pandemics, posing a serious threat to public health. Hemagglutinin (HA) plays a key role in its identification—a protein responsible for attachment to the host cell and defining the virus subtype (H1–H18). Accurate subtype determination is essential for epidemiological surveillance, vaccine strain selection, and pandemic risk assessment.

Existing classification methods (ClassyFlu [1], INFINITy [2]) are based on comparing nucleotide or amino acid sequences using hidden Markov models (HMMs) or alignment-free machine learning approaches [3,4]. However, they ignore the three-dimensional structure of the protein, which can be critical for functional differences even when sequence homology is high.

In recent years, graph neural networks (GNNs) have demonstrated high efficiency in structural bioinformatics tasks, enabling the use of residue interaction networks (RINs) [5].

We propose EpitopeGNN – a model that predicts the HA subtype from its RIN graph using informative node features and an architecture based on an attention mechanism (Graph Attention Network, GAT). This work is the first to address influenza subtype classification using three-dimensional structures and GNNs, allowing the incorporation of topological and physicochemical features of protein spatial folding.

While existing sequence-based methods, including subtype classifiers (ClassyFlu, INFINITy) and even epitope predictors (e.g., CLBTope [6]]), rely solely on primary structure, they cannot capture the three-dimensional conformational features that determine antigenicity. The ability to classify subtypes from 3D structure, as proposed here, provides a complementary level of information that is essential for understanding immune recognition and for rational vaccine design.

2. Materials and Methods

2.1. Dataset

The list of hemagglutinin PDB IDs was obtained by searching the RCSB PDB portal [7] using the keyword “hemagglutinin”. After filtering out structures with fewer than 50 residues, the final dataset consisted of 249 structures for which subtype annotations were present in the .cif.gz files. The subtype annotations extracted from the PDB files were cross-referenced with data from the GISAID database [8] for quality control.

The distribution of subtypes is shown in Table 1.

2.2. Construction of the RIN Graph

Let a set of residues R={r₁, r₂, …, r_n} be defined for each hemagglutinin structure, where n is the number of amino acid residues satisfying the selection criteria (standard amino acids, presence of a C_α atom). Each residue r_i is associated with a point ${\mathbf{c}}_i\in {\mathbb{R}}^3$—the coordinates of its C_α- atom, extracted from the .cif.gz file.

Graph Definition

An undirected graph $G=(V,E)$ is constructed, where:

$V=\{\mathrm{1,2},\dots ,n\}$ – is the set of vertices corresponding to residues;

$E\subseteq V\times V$ – is the set of edges defined by a threshold condition:

$$\left(i,j\right)\in E\hspace{0.25em} ⟺\hspace{0.25em} i\ne j \mathrm{and} \parallel {\mathbf{c}}_i-{\mathbf{c}}_j{\parallel }_2<\delta ,$$

where $\delta =8$ Å – is the chosen spatial proximity threshold. The choice of $\delta =8$ Å is based on literature data [9]: this distance is sufficient to cover most non-covalent interactions (hydrogen bonds, van der Waals contacts, ionic interactions) while excluding overly distant contacts that do not affect the local structure.

Adjacency matrix

The graph $G$ is represented by a symmetric adjacency matrix $\mathbf{A}\in \{\mathrm{0,1}{\}}^{n\times n}$, where:

$$A_{ij}=A_{ji}=\{\begin{array}{cc}1,& \mathrm{if} (i,j)\in E,\\ 0,& \mathrm{otherwise}.\end{array}$$

Representation for a graph neural network

For subsequent use in architectures such as GCN, GAT, or GraphSAGE, the graph is stored as a tuple $\left(\mathbf{X}\right|\mathbf{E})$, where $\mathbf{E}$ is an edge list (index scheme) that enables directed traversal during message aggregation. In an undirected graph, each edge $\left(i\right|j)$ is added twice: $\left(i\right|j) and \left(j\right|i)$.

2.3. Node Features

Each vertex $i\in V$is associated with a feature vector ${\mathbf{x}}_i\in {\mathbb{R}}^d$, $d=26$, which combines:

• One hot encoding of the amino acid type (20 dimensions). Each of the 20 canonical amino acids is represented by a binary vector, where a 1 is placed in the position corresponding to that residue.
• Physicochemical properties (4 dimensions): charge, hydrophobicity (according to the Kyte–Doolittle scale [10]), volume, polarity (according to the Grantham scale [11]).

These values are tabulated constants for each amino acid:

• Charge – the electric charge of the side chain at physiological pH (≈7.4). For example, Lys (+1), Arg (+1), Asp (−1), Glu (−1), others 0.
• Hydrophobicity – the index according to the Kyte-Doolittle scale, reflecting the hydrophobicity/hydrophilicity of the side chain. Values range from −4.5 (Arg) to +4.5 (Ile).
• Volume – the van der Waals volume of the side chain (in ų). For example, Gly ~48 ų, Trp ~163 ų.
• Polarity – the polarity measure according to Grantham, characterizing the ability of the residue to participate in polar interactions
• Topological characteristics of the graph (2 dimensions):
• Vertex degree $\mathrm{d}\mathrm{e}\mathrm{g}\left(i\right)={\sum }_j{A}_{ij}$–the number of neighbors in the graph.

• Betweenness centrality BC(i) – a measure of how often a node lies on shortest paths between other nodes. It was computed approximately for large graphs using Brandes’ algorithm [12].

The feature matrix of the entire graph $\mathbf{X}\in {\mathbb{R}}^{n\times d}$ is formed as

$$\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_n{]}^{\mathrm{\top }}.$$

2.4. EpitopeGNN Model Architecture

For classifying hemagglutinin subtypes based on RIN graphs, a graph neural network based on the attention mechanism – Graph Attention Network (GAT) [13] – was developed. The model consists of three consecutive GAT layers, followed by a global mean pooling operation and a fully connected layer with a softmax activation function for multi-class classification.

Input data. Each graph $G=(V,E)$ with n nodes is fed into the model as a set of node features $\mathbf{X}\in {\mathbb{R}}^{n\times d}$, where d = 26 is the dimensionality of the feature description (see Section 2.3), and a list of edges E that defines the connections between nodes.

GAT layer. Each GAT layer transforms node features using a self-attention mechanism. For node i at layer l, the updated representation ${\mathbf{h}}_i^{\left(l\right|1)}$is computed as follows:

$${\mathbf{h}}_i^{\left(l\right|1)}=\sigma \left(\sum _{j\in \mathcal{N}\left(i\right)}{\alpha }_{ij}^{\left(l\right)}{\mathbf{W}}^{\left(l\right)}{\mathbf{h}}_j^{\left(l\right)}\right),$$

where:

• ${\mathbf{h}}_i^{\left(l\right)}\in {\mathbb{R}}^F$ – feature of node i at the input of layer l (for the first layer, ${\mathbf{h}}_i^{\left(0\right)}={\mathbf{x}}_i$);
• $\mathcal{N}\left(i\right)$ – set of neighbors of node i (including itself if self-attention is used);
• ${\mathbf{W}}^{\left(l\right)}\in {\mathbb{R}}^{F^{\mathrm{\text{'}}}\times F}$ – trainable linear transformation matrix applied to all nodes;
• ${\alpha }_{ij}^{\left(l\right)}$ – attention coefficient reflecting the importance of node j’s features for updating the representation of node $i$;
• $\sigma$ – nonlinear activation function (ELU is used in this work).

The attention coefficients are computed using an additive attention mechanism:

$${\alpha }_{ij}^{\left(l\right)}={\displaystyle \frac{\exp \left(\mathrm{LeakyReLU}\left({\mathbf{a}}^{\left(l\right)\mathrm{\top }}\left[{\mathbf{W}}^{\left(l\right)}{\mathbf{h}}_i^{\left(l\right)}\parallel {\mathbf{W}}^{\left(l\right)}{\mathbf{h}}_j^{\left(l\right)}\right]\right)\right)}{\sum _{k\in \mathcal{N}\left(i\right)}\exp \left(\mathrm{LeakyReLU}\left({\mathbf{a}}^{\left(l\right)\mathrm{\top }}\left[{\mathbf{W}}^{\left(l\right)}{\mathbf{h}}_i^{\left(l\right)}\parallel {\mathbf{W}}^{\left(l\right)}{\mathbf{h}}_k^{\left(l\right)}\right]\right)\right)}}.$$

Here:

• $\parallel$ denotes vector concatenation;
• ${\mathbf{a}}^{\left(l\right)}\in {\mathbb{R}}^{2F^{\mathrm{\text{'}}}}$ – trainable attention parameter vector;
• LeakyReLU is used with a negative slope $\alpha =0.2$.

To stabilize training and enrich representations, each GAT layer employs multi-head attention. In this model, 4 independent attention heads are used, and their outputs are concatenated:

$${\mathbf{h}}_i^{\left(l\right|1)}={\parallel }_{k=1}^4\hspace{0.25em} \sigma \left(\sum _{j\in \mathcal{N}\left(i\right)}{\alpha }_{ij}^{\left(l\right|k)}{\mathbf{W}}^{\left(l\right|k)}{\mathbf{h}}_j^{\left(l\right)}\right).$$

Thus, the dimensionality of the hidden representation after each layer is $4\times F^{\mathrm{\text{'}}}$. In our model, $F^{\mathrm{\text{'}}}=32$, is chosen for all layers, so after concatenation the dimensionality becomes 128. After three GAT layers, we obtain node representations ${\mathbf{h}}_i^{\left(3\right)}\in {\mathbb{R}}^{128}$.

Global mean pooling. To obtain a vector representation (embedding) of the whole graph, a global mean pooling over nodes is applied:

$${\mathbf{h}}_G={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\mathbf{h}}_i^{\left(3\right)}.$$

The resulting vector ${\mathbf{h}}_G\in {\mathbb{R}}^{128}$ is a compact description of the hemagglutinin structure.

Classification head. The graph embedding ${\mathbf{h}}_G$ is fed into a fully connected layer with a softmax activation function, which transforms it into a probability distribution over $C$ classes (subtypes):

$$\mathbf{p}=\mathrm{softmax}\left({\mathbf{W}}_{\mathrm{cls}}{\mathbf{h}}_G\right|{\mathbf{b}}_{\mathrm{cls}}),$$

where ${\mathbf{W}}_{\mathrm{cls}}\in {\mathbb{R}}^{C\times 128}$ and ${\mathbf{b}}_{\mathrm{cls}}\in {\mathbb{R}}^C$ are trainable parameters. The predicted subtype corresponds to the class with the highest probability.

2.5. Training

The dataset was randomly split into training (70%) and test (30%) sets while preserving class proportions (stratification). Training was performed using the cross-entropy loss function:

$$\mathcal{L}=-{\displaystyle \frac{1}{N}}\sum _{i=1}^N\mathrm{l}\mathrm{o}\mathrm{g}\left({\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(z_{i,y_i}\right)}{\sum _{c=1}^C\mathrm{e}\mathrm{x}\mathrm{p}(z_{i,c})}}\right)$$

The optimizer was Adam (learning rate = 0.001), weight decay = 5e-4, batch size = 4, and the number of epochs was 100. To assess stability, five-fold cross-validation was performed.

3. Results

3.1. Classification Performance

On the test set, the model achieved 100% accuracy for all subtypes. Five-fold cross-validation yielded an average accuracy of 97.6% (σ=2.3%). Table 2 presents the classification report for the test set.

3.2. Embedding Visualization

To visualize the data structure, 128-dimensional embeddings extracted from the penultimate layer of the EpitopeGNN model were projected onto a plane using the t-SNE algorithm (Figure 1). The projection includes all 249 structures belonging to eight hemagglutinin subtypes: H1N1 (84), H2N2 (11), H3N2 (80), H5N1 (48), H7N3 (4), H7N7 (4), H7N9 (15), and H9N2 (3).

Visualization. In Figure 1, seven compact clusters are clearly visible. Five of them correspond to the abundant subtypes H1N1, H2N2, H3N2, H5N1, and H7N9. The clusters of the rare subtypes H7N3 and H7N7 (4 structures each) almost completely overlap with the H7N9 cluster on the plane, which is not unexpected since all these subtypes share H7 HA molecules (distances between them in the t-SNE space are only 0.7–1.5). Despite the small number of samples (3 structures), the H9N2 subtype appears as a separate cluster located near H2N2 (distance between centroids 2.21), which is consistent with the phylogenetic relationship of these subtypes.

Quantitative analysis. To objectively assess cluster separability, centroid coordinates in the t-SNE space (Table 3) and pairwise Euclidean distances between them (Table 4) were calculated. In addition, the average silhouette coefficient was computed for each subtype, characterizing cluster compactness and separation from neighboring clusters (Table 5). The overall silhouette score for all 249 structures was 0.495, indicating good separation quality (values above 0.5 are considered excellent). It should be noted that for subtypes with a small number of samples (H7N3, H7N7, H9N2), silhouette estimates may be less stable due to high sensitivity to outliers.

High silhouette values for H1N1, H2N2, H3N2, H5N1, and especially for H9N2 (0.96) indicate exceptional compactness and isolation of these clusters in the original 128-dimensional space. The negative silhouette for H7N3 (–0.025) is explained by the small number of points (4) and their proximity to the larger H7N9 cluster (distance between centroids 0.74). The visual merging of the H7N3, H7N7, and H7N9 clusters in the t-SNE projection does not contradict the classification results: in the high-dimensional space, these subtypes remain distinguishable (100% classification accuracy), but the two-dimensional projection, with its inevitable loss of information, cannot capture the subtle differences between such close clusters. Thus, the visualization results are fully consistent with the performance metrics and confirm the effectiveness of the proposed approach.

3.3. Hierarchical Clustering

To further verify the structural similarity between subtypes, we performed Ward’s hierarchical clustering on the same 128-dimensional embeddings. The resulting dendrogram (Figure 2) shows a clear separation into seven main branches, which completely matches the visual pattern of the t-SNE plot and confirms the classification results.

The branches corresponding to the H1N1, H2N2, H3N2, H5N1, and H7N9 subtypes separate at high hierarchy levels, indicating their substantial separation in the embedding space. The rare subtypes H7N3, H7N7, and H9N2 also form distinct branches, albeit at lower heights, which is consistent with their small numbers and proximity to larger clusters (e.g., H7N3 and H7N7 are nested within the H7N9 branch but retain internal clustering). Importantly, no branch overlaps with another, further confirming the model’s ability to capture subtle structural differences between hemagglutinin subtypes.

3.4. Comparison of Structural and Genetic Distances

To assess the biological relevance of the obtained structural embeddings, we quantitatively compared the distance matrices derived from them with classical phylogenetic distances calculated from amino acid sequences. For all 249 structures in the dataset, amino acid sequences were extracted from the corresponding PDB files. Multiple sequence alignment was performed using MUSCLE [14], after which a genetic distance matrix was calculated based on the p-distance (proportion of mismatched positions). The structural distance matrix was obtained by computing pairwise Euclidean distances between the 128-dimensional embeddings extracted from the penultimate layer of the EpitopeGNN model.

Mantel correlation. To evaluate the similarity between the two distance matrices, the Mantel test was applied with 999 permutations. Pearson’s correlation coefficient between the matrices was r = 0.484 with a high level of statistical significance (p = 0.001). This value indicates a moderate but significant positive correlation. The result demonstrates that the structural embeddings, trained without any evolutionary information (based solely on three-dimensional structure), nonetheless capture a phylogenetic signal. The squared correlation coefficient (r²) indicates the proportion of variance in one variable that is explained by a linear relationship with the other variable:

$$r^2=(0.4842{)}^2\approx 0.234 (\mathrm{e.g.}~23.4).$$

Thus, the proportion of variance not explained by a linear relationship is $1-r^2=1-0.234=0.766$, or ~76.6%, round² ≈ 0.23) suggests that the structural embeddings reflect not only evolutionary differences encoded by the sequences but also additional aspects of spatial organization. This implies that such embeddings may carry information important for the analysis of functional and antigenic properties, which are directly linked to the three-dimensional protein structure. Further studies are needed to test this hypothesis. In more detail, the unexplained variance is not “empty” but rather comprises several meaningful factors:

Nonlinear relationships. The relationship between genetic and structural distances may not be strictly linear. For example, small changes in sequence can lead to disproportionately large structural rearrangements (or vice versa). If nonlinear measures (e.g., rank-based distances) were used, the proportion of explained variation might change.

Noise and measurement errors. Genetic distances were calculated from sequences, which may contain uncertainties (especially in alignment). Structural embeddings, in turn, were obtained from a GNN model trained on a limited dataset and may include some noise or artifacts.

Additional information. The three-dimensional structure of a protein carries information that is not completely determined by the amino acid sequence. In particular, it has been shown [15] that when sequence identity falls below 25–30%, reliably predicting structural similarity from sequence alone is impossible. This fundamentally justifies the use of structural information as an additional source. For instance, conformational variability within a single subtype (different H1N1 strains may have somewhat different structures, which is reflected in the embeddings) represents important structural signal that is not reducible to ordinary amino acid sequence. Furthermore, side-chain packing, surface charge distribution, etc., may also influence the structural embeddings but are poorly captured by simple p-distance. This part is often referred to as “additional information” and is potentially important for protein function.

Thus, given that the GNN model was trained exclusively on three-dimensional structures and physicochemical properties (without directly using sequences), it is reasonable to assume that a significant portion of the unexplained variation reflects structural features that are not reducible to sequence.

3.5. Distance Matrix Heatmaps

For a more detailed analysis, Figure 3 shows heatmaps of the genetic and structural distance matrices. For clarity, the rows and columns of both matrices are ordered according to the hierarchical clustering of the structural embeddings. In these plots, each cell represents the distance between two hemagglutinin structures:

• The horizontal and vertical axes correspond to all 249 structures (in the order obtained from hierarchical clustering).
• The color of each cell represents the distance magnitude:

○

Dark blue → structures are very close (small distance).

○

Green/yellow → intermediate distance.

Red/orange → large distance (strongly different).

Because the rows and columns have been reordered (clustered), similar structures appear near each other. Along the diagonal (from the top left to the bottom right), structures are compared with themselves – the distance is always 0, so the diagonal is dark blue. The blocks along the diagonal represent groups of structures that are close to each other. To the left of each heatmap, a color annotation of subtypes is added. This subtype annotation is a vertical bar to the left of the heatmap, where each structure is represented by a vertical strip whose color encodes its subtype. Continuous regions of the same color correspond to compact groups of structures belonging to the same subtype, as formed by hierarchical clustering. The structural matrix shows clearer boundaries between some clusters (e.g., H1N1 and H3N2), which is consistent with the high classification accuracy.

3.6. Outlier Analysis

Ordering the data by structural similarity revealed potential structural anomalies. Five structures have a mean Euclidean distance to all other objects exceeding the threshold of two standard deviations (Figure 3, Table 6). These structures can be considered potential structural outliers. Interestingly, four of them belong to subtypes H1N1 (6Z2T, 6GOL) and H2N2 (8TP4, 2WRF). Since the classification accuracy for these subtypes is 100%, these structures are not misclassified; however, their spatial conformation deviates significantly from typical representatives of their subtype. This may be due to features of the specific strain (e.g., mutations in epitopes), structure quality, or crystallization conditions. Structure 5DUT (H5N1) also appears in the outlier list, making it an interesting target for further structural and functional analysis.

Thus, the analysis confirms that the proposed embeddings are not merely an abstract feature space for classification, but a compact and informative representation reflecting both evolutionary relationships and unique structural features of hemagglutinin.

4. Discussion

The achieved 100% classification accuracy on the test set and 97.6% with five-fold cross-validation convincingly demonstrate that the three-dimensional structure of hemagglutinin (HA) contains all the necessary information for discriminating influenza A virus subtypes. In contrast to existing methods (ClassyFlu [1], INFINITy [2]) that work exclusively with DNA or amino acid sequences, the proposed EpitopeGNN approach is the first to use a graph representation of the protein spatial fold and graph neural network mechanisms. This allows it to account not only for the local sequence but also for the global topology of residue interactions and their physicochemical properties, which is critical for understanding antigenic differences.

Quantitative analysis of embedding clustering confirmed high separation quality: the overall silhouette coefficient was 0.495, and for most subtypes (H1N1, H2N2, H3N2, H5N1, H9N2) the average silhouette scores exceeded 0.46, reaching 0.96 for H9N2. The proximity of the centroids of H7N3, H7N7 and H7N9 (distances 0.7–1.5) and the partial overlap of these clusters in t-SNE are explained by the small number of samples (4 and 15) and inevitable distortions upon dimensionality reduction. At the same time, the high classification accuracy and the absence of errors on the test set indicate that in the original 128-dimensional embedding space these subtypes remain completely separable. Furthermore, a significant correlation was quantitatively demonstrated between structural embeddings and phylogenetic distances, indirectly confirming the presence of an evolutionary signal in the spatial protein fold.

Comparison with existing services highlights the fundamental novelty of this work: none of them use three-dimensional structures. Thus, EpitopeGNN does not replace but complements traditional genetic classifiers by providing a fundamentally different level of information – structural. This opens new possibilities for monitoring virus evolution, enabling tracking not only of genetic changes but also of rearrangements in protein spatial folding that may affect antigenic properties and pandemic potential.

Comparison with sequence-based methods. Traditional approaches for influenza subtype classification, such as ClassyFlu [1], INFINITy [2], operate exclusively on nucleotide or amino acid sequences. Although they are fast and widely applicable, they cannot account for the spatial folding of the protein, which is directly linked to antigenic properties. A recent method, CLBTope [6], also uses only sequence information to predict B-cell epitopes, but its accuracy for hemagglutinin is naturally limited because conformational epitopes are defined by spatial proximity rather than linear sequence. Our structural embedding-based approach overcomes this limitation by explicitly modeling residue interactions and surface accessibility, enabling not only accurate subtype classification but also the identification of potential antigenic determinants. This structural perspective is critical for tracking antigenic drift and for guiding vaccine strain selection.

Limitations and prospects

The model requires an experimentally determined or reliably predicted three-dimensional structure, which is not always available for new strains. However, the rapid development of protein structure prediction methods such as AlphaFold will soon remove this limitation.

Rare subtypes (H7N3, H7N7, H9N2) are represented by a small number of samples, which may affect clustering stability. Enlarging the dataset with new experimental structures or high-quality predicted models will improve classification reliability for these subtypes.

Future work will include enriching the node feature representation with evolutionary information (PSSM), as well as exploring the use of the obtained embeddings for finding structural analogs and predicting other functional properties of HA, such as receptor specificity (human α2,6 vs. avian α2,3).

5. Conclusions

1. An original graph neural network (EpitopeGNN) was developed for classifying hemagglutinin subtypes of influenza A virus, utilizing the three-dimensional protein structure and physicochemical features of amino acid residues.

2. On a test set of 249 structures (8 subtypes), the model achieved 100% accuracy; five-fold cross-validation confirmed the stability of the result (97.6 ± 2.3%).

3. The proposed method complements traditional sequence-based classifiers and opens up opportunities for structural monitoring of virus evolution, searching for structural analogs, and predicting functional properties of hemagglutinin (e.g., receptor specificity).

4. Further development of the model involves incorporating evolutionary information (PSSM), expanding the dataset with predicted structures (AlphaFold, ESMFold, etc.), and applying the embeddings to related tasks.

5. Moreover, the structural embeddings generated by EpitopeGNN may serve as a foundation for predicting functional properties such as receptor specificity and epitope immunogenicity, tasks that are challenging for sequence-only methods like CLBTope.