SGCL-DPI: structure-guided curriculum learning for drug-protein interaction prediction

Dataset preprocessing

This study involved multiple preprocessing steps to prepare a high-quality dataset suitable for binary classification of drug–protein interactions. These steps included filtering affinity data based on IC50 thresholds, labeling interaction pairs, and splitting the dataset for model training and evaluation. Additional preprocessing was carried out to construct a hard evaluation set based on strict compound and protein separation criteria, following benchmark protocols.

Binding database (BindingDB)

Given its extensive coverage, BindingDB allows researchers to explore a diverse set of interactions across multiple therapeutic areas. The dataset includes both approved drugs and experimental compounds from various chemical classes, facilitating the study of diverse protein-drug interactions and the identification of potential therapeutic candidates (Liu et al., 2025). Moreover, BindingDB’s open-access nature promotes global collaboration and reproducibility, while its high-quality, carefully curated data ensures reliability for computational modeling, virtual screening, and structure-activity relationship analysis. The original BindingDB dataset is available at: https://www.bindingdb.org/rwd/bind/index.jsp.

The BindingDB dataset, as utilized in the study, contains affinity data for 2,286,319 drug-protein pairs, encompassing 8,536 proteins and 989,383 drugs. Gao et al. (2018) refined this dataset by selecting data with IC50 values and converting these into binary labels: interactions (IC50 < 100 nM) were labeled as 1, and no interactions (IC50 > 10,000 nM) were labeled as 0 (Gao et al., 2018; Wu et al., 2022). This resulted in a binary classification dataset comprising 33,777 positive samples and 27,493 negative samples.

The data preparation process also involved splitting the dataset into training and evaluation sets. A holdout setting was employed due to the complexity of the training process and the substantial volume of training data. We used 80% of the dataset for training the models to uncover patterns and relationships between drugs and proteins, while the remaining 20% was reserved for evaluating model performance and generalization capabilities.

Benchmark dataset and hard split construction

In the present work, a partitioning protocol presented by Watanabe, Ohnuki & Sakakibara (2021) was employed to evaluate model performance under challenging conditions. In their study, the authors downloaded protein–compound, protein–protein, and compound–compound interactions from the STITCH (Szklarczyk et al., 2016) and STRING (Szklarczyk et al., 2025) databases, compiling 22,881 protein–compound, 175,452 protein–protein, and 69,231 compound–compound interactions. From this comprehensive collection, they constructed three distinct cross-validation schemes: the baseline dataset, the unseen compound-test dataset, and the hard dataset.

The hard dataset, used here for rigorous evaluation, was constructed by partitioning the data so that neither proteins nor compounds appearing in any test fold were included in the corresponding training folds. This demanding split requires accurate prediction of interactions for completely unseen proteins and compounds, closely mimicking real-world scenarios where both entities are novel.

It is important to note that the aforementioned splitting methodology was directly adopted from Watanabe, Ohnuki & Sakakibara (2021). This ensures that evaluations are performed on a benchmark setting validated in previous studies. In particular, the hard dataset serves as a stringent test of generalization capability and provides a robust basis for comparing these results with other state-of-the-art methods.

Overview of SGCL-DPI architecture

SGCL-DPI is presented (see Fig. 1) as an innovative single-modality deep learning framework for predicting drug–protein interactions, distinguished by a unified architecture that integrates both drug and protein information. Although it incorporates both GNN and RF components, it produces a single, coherent prediction from a single model. Unlike ensemble methods that aggregate outputs from multiple models at inference time, SGCL-DPI uses the RF solely during training as a teacher guiding the learning process via curriculum strategies and indirectly shaping the final prediction. This auxiliary role for the RF is central to the curriculum learning strategy, where its guidance is progressively diminished as training advances. The method first encodes drugs and proteins into feature vectors and constructs similarity graphs to encapsulate their structural relationships. These graphs are subsequently processed by a graph neural network that enriches the embeddings, which are then fused in a final prediction network. This fusion integrates the distilled knowledge from the RF, ensuring that the overall prediction pipeline remains fundamentally structure-guided and curriculum-trained rather than dependent on an ensemble of disparate predictors.

Feature generation and molecular representation

Drug features were generated using a two-component approach: extended-connectivity fingerprints (ECFPs) and chemical descriptors (Soufan et al., 2018). ECFPs were computed using RDKit’s Morgan fingerprint algorithm with radius 2 and 1,024 bits. Chemical descriptors included 15 physicochemical properties: molecular weight, LogP, hydrogen bond donors/acceptors, topological polar surface area, ring count, fraction of SP3 carbons, heavy atom count, rotatable bonds, and partial charge statistics. These descriptors were normalized and combined with ECFPs to create comprehensive molecular representations.

Protein features were generated using a novel combination of sequence-based descriptors. The approach included amino acid composition (AAC), physicochemical properties, and a hashed one-hot encoding scheme. AAC features captured the normalized frequency of each standard amino acid. Physicochemical properties incorporated hydrophobicity, molecular weight, and isoelectric point metrics for each residue, scaled using MinMaxScaler. A fixed-length representation was achieved through a hashing trick applied to one-hot encoded sequences, mapping to a 500-dimensional space while preserving sequence information.

Similarity network construction and integration

Drug-drug (DDS) and protein-protein (PPS) similarity networks were constructed using a weighted combination of feature-based similarities. For drugs, the similarity scores combined Tanimoto similarity of ECFPs (weight = 0.5) with normalized Euclidean distance of chemical descriptors (weight = 0.5). In particular, for two drugs i and j, we compute: (1) Fingerprint similarity as the Tanimoto coefficient between their 1,024-bit ECFPs, and (2) descriptor similarity by taking the Euclidean distance between their descriptor vectors and converting it to a similarity score via a normalized inverse-distance function $si{m}_{desc}=1/\left(1+dist\right)$. We then take a weighted average (we use equal weights 0.5 each by default) of these two measures to obtain a unified similarity score. Using this score, each drug is connected to its k nearest neighbors (highest similarity) in the dataset. The edge weight in the graph is the computed similarity value, reflecting the strength of the relationship. Protein similarities were computed analogously using their sequence-derived features. To manage computational complexity while preserving network information, similarities were computed using a k-nearest neighbor approach (k = 5), with edges established only for similarity scores exceeding 0.5. The analysis generated 14,971,621 drug-drug similarity links and 309,057 protein-protein similarity links.

The resulting similarity networks were integrated into the deep learning architecture through edge-weighted graph convolutional layers. Edge weights, derived from the computed similarity scores, were normalized using L1 normalization and constrained to [0, 1]. These weighted networks were processed through multiple GCN layers, each incorporating layer normalization and Gaussian error linear units (GELU) activation functions. To maintain numerical stability during message passing, safe normalization was implemented with a minimum norm threshold (ε = 1e−12).

This comprehensive feature generation and network construction approach enables effective capture of both local molecular properties and global similarity patterns, while the sparse graph representation achieved through similarity thresholding provides computational efficiency.

To test robustness, additional graphs were generated with (i) cosine similarity at $k=3$ and $k=7$, (ii) Euclidean similarity ( $k=5$), and (iii) Pearson-correlation similarity ( $k=5$). All other settings were kept identical. Refer to Table S1 for details.

Graph representations and data encoding

The proposed framework models drugs and proteins as nodes within similarity networks (Meng et al., 2024; Xu et al., 2024), capturing global relationships derived from their chemical and biological descriptors. In the drug–drug similarity graph, each drug is represented as a node. Formally, the graph is defined as

$$G_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}=\left(V_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}},E_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}},H_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}\right),$$where $V_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}$ is the set of drugs, and $H_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}\in {\mathbb{R}}^{\left|V_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}\right|\times d}$ is a node feature matrix that encodes chemical descriptors such as molecular weight, LogP, and other properties. The edges $E_{\mathrm{d}\mathrm{r}\mathrm{u}\mathrm{g}}$ connect drugs based on similarity metrics, with Tanimoto coefficients computed from molecular fingerprints serving as the basis for quantifying structural similarity.

Proteins are similarly modeled in a protein–protein similarity graph:

$$G_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}=\left(V_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}},E_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}},H_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}\right)$$

Here, $V_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}$ denotes the set of proteins, and $H_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}\in {\mathbb{R}}^{\left|V_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}\right|\times d^{\prime }}$ is the node feature matrix that contains descriptors derived from sequence or structural data. Edges $E_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}}$ are defined based on similarities computed, for example, from normalized Euclidean distances between physicochemical descriptors.

In both graphs, edge weights are assigned to reflect similarity. For drugs, these weights are derived from Tanimoto coefficients computed on molecular fingerprints. For proteins, a dual approach is employed: sequence alignment scores capture evolutionary and structural relationships, and normalized Euclidean distances of relevant physicochemical descriptors quantify differences in protein properties. This combination produces a robust and comprehensive measure of protein similarity. For numerical stability, the weights are normalized and clipped according to

$$w_{ij}^{\prime }=clip\left({\Vert w_{ij}\Vert }_1,ϵ,1-ϵ\right),$$with $\mathrm{ϵ}={10}^{-12}$.

Graph neural network architecture and message passing

A GNN is employed to learn enriched embeddings for drugs and proteins from the similarity graphs (Hao et al., 2025; Wu et al., 2020). The GNN processes the drug graph (i.e., DDS-GNN-see Fig. 1) and protein graph (i.e., PPS-GNN–see Fig. 1) separately (i.e., a separate GNN encoder for each similarity graph), producing a latent representation for each drug and each protein that accounts for its neighbors. Both graph encoders share a similar architecture: a multi-layer graph convolutional network with normalization and attention mechanisms to ensure stable and informative embeddings.

Each encoder is a stacked GCN with $L=3$ layers (graph convolution operations). For the drug encoder, the input feature size is 1,039 (fingerprint + descriptors); for the protein encoder, it is the length of the protein feature vector (i.e., 523 using a 500-dim hash + 23-dim descriptors). We first project the input features to a lower-dimensional hidden space using a linear layer of 128 dimensions followed by LayerNorm and a GELU nonlinear activation function. This initial projection ensures a manageable embedding size for graph convolution. At each GNN layer $l$, the feature vector of node $i$ is updated by aggregating information from its neighbors according to

$$h_i^{\left(l+1\right)}=\mathrm{\sigma }\left({\displaystyle \frac{1}{\sqrt{d_i{d}_j}}\sum _{j\in \mathcal{N}\left(i\right)}{\mathrm{\alpha }}_{ij}{W}^{\left(l\right)}{h}_j^{\left(l\right)}}\right).$$

In this expression, $h_i^{\left(l\right)}$ is the feature vector at layer $l$, $W^{\left(l\right)}$ is the learnable weight matrix, $d_i$ represents the node degree, and ${\mathrm{\alpha }}_{ij}$ are scaling factors computed from edge features. The non-linear activation $\mathrm{\sigma }$ (implemented as GELU) introduces necessary non-linearity, while normalization by $\sqrt{d_i{d}_j}$ prevents nodes with high connectivity from overwhelming the aggregation process.

Additionally, GCN layers refine these node embeddings. The GCN formulation, incorporating self-loops, is expressed as

$$H^{\left(l+1\right)}=\sigma \left({\tilde{D}}^{-1/2}\tilde{A}{\tilde{D}}^{-1/2}{H}^{\left(l\right)}{W}^{\left(l\right)}\right),$$where $\tilde{A}=A+I$ is the augmented adjacency matrix and $\tilde{D}$ is its corresponding degree matrix. This formulation effectively aggregates both the local neighborhood and the node’s own features, resulting in robust and discriminative embeddings for subsequent analysis.

Multi-modal integration and feature fusion

The framework integrates features from GNNs with complementary predictions from a RF model that uses traditional molecular descriptors, thereby leveraging both deep and classical machine learning insights.

In our approach, drugs and proteins are processed through separate GNN pipelines to extract nuanced representations from their respective similarity graphs. A GCN layer takes as input a node feature matrix $H$ and an adjacency matrix $A$ that captures the pairwise similarities between nodes. Conceptually, the code implements an operation analogous to

$$H^{\prime }=ClampValue\left(SafeNormalize\left(GCNConv\left(H,A\right)\right)\right),$$where $SafeNormalize\left(x\right)={\displaystyle \frac{x}{max\left(||x|{|}_2,ϵ\right)}}$ ensures numerical stability and $ClampValue\left(x\right)=min\left(max\left(x,-1\right),1\right)$ bounds the feature values. By normalizing activations, we prevent the embeddings from growing unbounded or becoming too small to be informative.

To obtain fixed-size graph representations, node features are passed through a multi-layer perceptron:

$$MLP\left(H\right)=W_2\cdot GELU\left(LayerNorm\left(W_{1}H\right)\right),$$after which attention weights are computed via

$$\mathrm{\alpha }=Softmax\left(Mask\left(MLP\left(H\right)\right)+\mathrm{ϵ}\right).$$The final graph representation $g=\alpha H$ is computed as a weighted sum of the node features $H$, where the weights are given by the attention vector $\alpha$. Therefore, nodes with higher attention weights contribute more to the overall graph representation.

To integrate the GNN-derived representation with the RF model’s output, two complementary mechanisms are used. First, direct feature fusion projects both $g_{\mathrm{G}\mathrm{N}\mathrm{N}}$ and the RF predictions $p_{\mathrm{R}\mathrm{F}}$ into a shared latent space via linear transformations and nonlinear activations as follows:

$$g_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}=LayerNorm\left(GELU\left(W_g{g}_{\mathrm{G}\mathrm{N}\mathrm{N}}\right)\right)$$

$$p_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}=LayerNorm\left(GELU\left(W_p{p}_{\mathrm{R}\mathrm{F}}\right)\right).$$

The resulting embeddings are then concatenated or combined through a learnable fusion function, creating a unified representation that benefits from both the GNN’s structural insights and the RF model’s robust, descriptor-based predictions:

$$f_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}=Fusion\left(\left[g_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}},p_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}\right]\right),$$which encapsulates both the detailed local structure from the GNN and the robust global patterns from the RF model. Additionally, RF outputs further guide training via a knowledge distillation loss, aligning the GNN’s predictions with those of the RF model.

Curriculum learning strategy

The model is trained using a two-phase curriculum learning (Hacohen & Weinshall, 2019; Soviany et al., 2022) strategy designed to gradually shift emphasis from global patterns captured by RF predictions to the detailed local structural features extracted by the GNN. In the initial phase, known as RF-guided learning, the loss function is configured to heavily weight the RF predictions typically assigning a weight of $\alpha =0.8$ to this component while the structural component receives a lower weight, such as $\beta =0.2$. This approach allows the model to quickly learn reliable global patterns before focusing on finer details.

As training proceeds, the weight distribution is smoothly adjusted using an exponential schedule:

$${\alpha }_t={\alpha }_0\cdot \exp \left(-{\displaystyle \frac{t}{\tau }}\right),{\beta }_t={\beta }_0\cdot \left(1-\exp \left(-{\displaystyle \frac{t}{\tau }}\right)\right),$$where $\tau$ is a hyperparameter that determines the rate of transition. This gradual re-weighting minimizes abrupt changes in the learning dynamics, ensuring a smooth transition from RF-guided to structure-guided learning.

Multi-component loss function

The training objective is defined by a composite loss function that simultaneously addresses interaction prediction, RF-guided knowledge distillation, and structural consistency. The primary loss is a weighted binary cross-entropy loss formulated as:

$${\mathcal{L}}_{BCE}=-\sum _i{w}_i\left[y_i\log \left(\widehat{y_i}\right)+\left(1-y_i\right)\log \left(1-\widehat{y_i}\right)\right],$$where the adaptive weights are computed by

$$w_i=1+\exp \left(-\mathrm{\gamma }\left|\widehat{y_i}-y_i\right|\right)\cdot c_i,$$Here, $c_i$ quantifies the confidence of the RF predictions, and $\gamma$ is a scaling hyperparameter that accentuates the focus on samples with higher prediction error.

To ensure the model’s predictions align with those from the RF model, a knowledge distillation loss is introduced:

$${\mathcal{L}}_{\mathrm{K}\mathrm{D}}={\lambda }_{\mathrm{m}\mathrm{e}\mathrm{d}}\mathrm{M}\mathrm{S}\mathrm{E}\left(\hat{y},p_{\mathrm{R}\mathrm{F}}\right)\cdot I\{p_{\mathrm{R}\mathrm{F}}>0.4\}+{\lambda }_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}\mathrm{M}\mathrm{S}\mathrm{E}\left(\hat{y},p_{\mathrm{R}\mathrm{F}}\right)\cdot I\{p_{\mathrm{R}\mathrm{F}}>0.8\},$$with the dynamic weights defined as

$${\lambda }_{\mathrm{m}\mathrm{e}\mathrm{d}}=\exp \left(-\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}{\mathrm{n}}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}}\right),{\lambda }_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=\exp \left(-\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}{\mathrm{l}}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}}\right).$$

This dual-threshold strategy ensures that high-confidence RF predictions have a stronger influence on the training process while still incorporating moderately confident signals. Additionally, a structural consistency loss ${\mathcal{L}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}}$ is incorporated to maintain the fidelity of the learned representations with respect to established biochemical and biophysical similarity metrics. This term is calibrated based on empirical similarity distributions.

The overall loss function is then given by

$${\mathcal{L}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}={\lambda }_1{\mathcal{L}}_{\mathrm{B}\mathrm{C}\mathrm{E}}+{\lambda }_2{\mathcal{L}}_{\mathrm{K}\mathrm{D}}+{\lambda }_3{\mathcal{L}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}},$$

This composite loss function is key to balancing the multiple objectives and ultimately achieving robust performance in drug–protein interaction prediction. ${\mathcal{L}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ is optimized using the AdamW optimizer (Adam with weight decay) for robust convergence. Throughout training, the contribution of each loss component is monitored to ensure that none of them dominates in an unhealthy way (for instance, if ${\mathcal{L}}_{\mathrm{K}\mathrm{D}}$ remains large when $\mathrm{\alpha }$ is supposed to be low, that indicates the model is still trying to chase the RF–we avoided this by tuning the decay schedule).

Theoretical rationale for the descriptor-to-structure curriculum

Drug–protein interaction (DPI) prediction involves two information sources for each compound–target pair $x=\left(x_g,x_s\right)$:

• Global descriptors $x_g$: fixed-length physicochemical fingerprints that encode coarse, low-noise signals (e.g., logP, topological indices).

• Structural graphs $x_s$: high-dimensional, task-specific atom–bond graphs that capture fine-grained three-dimensional complementarity but are harder to model.

Let

• $T\left(x_g\right)$ be a teacher hypothesis learned by a Random Forest (RF) on $x_g$.

• $S\left(x_s;\theta \right)$ be the student GNN we ultimately care about.

• $y\in \left\{0,1\right\}$ be the ground-truth interaction label.

At training epoch $t$ $\left(0\le t<T_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$ we minimise the following curriculum objective:

$$\mathcal{L}\left(t\right)=\left(1-\alpha \left(t\right)\right)\cdot BCE\left(S\left(x_s;\theta \right),y\right)+\alpha \left(t\right)\cdot KD\left(S\left(x_s;\theta \right),T\left(x_g\right)\right)+{\lambda }_{struct}\cdot {\mathcal{L}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}}$$where

• $BCE\left(S\left(x_s;\theta \right),y\right)$ is the label-driven term

• $KD\left(S\left(x_s;\theta \right),T\left(x_g\right)\right)$ is the teacher-driven term

• $0\le \alpha \left(t\right)\le 1$ is a monotonically decreasing curriculum weight with

  $\alpha \left(0\right)=1$ and $alpha\left(T_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=0$.

• KD is any differentiable divergence (e.g., KL, soft-MSE) that pushes the student logits toward the teacher’s.

No specific numerical schedule is assumed and any smooth decay that satisfies the boundary conditions can implement the curriculum.

This formulation is equivalent to the earlier formulation ${\mathcal{L}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}={\lambda }_1{\mathcal{L}}_{\mathrm{B}\mathrm{C}\mathrm{E}}+{\lambda }_2{\mathcal{L}}_{\mathrm{K}\mathrm{D}}+{\lambda }_3{\mathcal{L}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}}$, with the identification ${\lambda }_1=1-\alpha \left(t\right)$, ${\lambda }_2=\alpha \left(t\right)$, ${\lambda }_3={\lambda }_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}}$. Thus λ₁ and λ₂ are not tuned independently; they are deterministically governed by the curriculum schedule α(t).

Early epochs are dominated by $T\left(x_g\right)$, whose estimation variance is low because RF ensembles average over descriptors with limited dimension. This regularises the GNN, preventing it from over-fitting the sparse, high-variance structural space. As $\mathrm{\alpha }\left(\mathrm{t}\right)\downarrow$, the student gradually reduces bias by exploiting structural nuances inaccessible to the teacher. Formally, under the classic bias–variance trade-off, a convex combination of low-variance and low-bias estimators yields an expected risk no worse than either extreme (Hacohen & Weinshall, 2019).

Curriculum learning suggests that presenting samples in order of increasing task difficulty accelerates convergence and improves generalization (Bengio et al., 2009). In DPI, mapping ${\mathrm{x}}_{\mathrm{g}}\to \mathrm{y}$ is an easier sub-task (lower VC-dimension) than mapping $\left(x_g,x_s\right)\to y$. The schedule $\alpha \left(t\right)\downarrow$ therefore adheres to the self-paced-learning principle that easier hypotheses guide the learner toward regions of parameter space associated with flatter minima.

Model optimization and training protocol

Robust optimization techniques are critical to ensuring the stable convergence of deep graph networks. In this work, several strategies are employed to maintain training stability and efficient convergence.

A cosine annealing schedule with warm restarts is utilized to modulate the learning rate over the course of training. Specifically, the learning rate at time $t$ is defined as

$$\mathrm{l}{\mathrm{r}}_t=\mathrm{l}{\mathrm{r}}_{min}+0.5\left(\mathrm{l}{\mathrm{r}}_{max}-\mathrm{l}{\mathrm{r}}_{min}\right)\left(1+\cos \left({\displaystyle \frac{t_{\mathrm{c}\mathrm{u}\mathrm{r}}}{T_i}\mathrm{\pi }}\right)\right),$$ $\mathrm{l}{\mathrm{r}}_{min}$ and $\mathrm{l}{\mathrm{r}}_{max}$ denote the lower and upper bounds of the learning rate, $t_{\mathrm{c}\mathrm{u}\mathrm{r}}$ is the number of epochs since the last restart, and $T_i$ is the period of the current cycle. This schedule facilitates periodic exploration of the parameter space, which is beneficial for escaping local minima and promoting a thorough search of the weight landscape. The periodic warm restarts allow the learning rate to reset to a higher value, thereby reintroducing diversity into the optimization process.

To further stabilize the training process, a gradient clipping strategy is applied. This strategy constrains the norm of the gradients during backpropagation to prevent gradient explosion, which is particularly crucial in deep architectures such as those involving graph convolutions.

Formally, gradients are clipped such that $\left|\mathrm{\nabla }\mathrm{\theta }\right|\le \mathrm{c}\mathrm{l}\mathrm{i}\mathrm{p}\mathrm{\_}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}$, with $\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{p}\mathrm{\_}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}=1.0$. This constraint ensures that the parameter updates remain within a reasonable range, preserving both the stability and the efficiency of the optimization process.

In addition to these measures, a batch size of 512 is chosen to balance computational efficiency with the statistical reliability of gradient estimates. Consistent application of layer normalization and safe normalization techniques further ensures that feature distributions remain stable across mini-batches. These normalization techniques help mitigate internal covariate shifts, thereby improving the overall training dynamics.

This methodological framework presents a sophisticated multi-modal approach to drug–protein interaction prediction that goes beyond a simple ensemble of RF and GNN predictions. Rather than merely combining the outputs of two independent models, RF predictions serve as a critical guiding signal throughout the training process. Initially, the RF model—known for its robustness and interpretability based on traditional molecular descriptors—provides reliable global predictions that help the network quickly capture broad interaction patterns. This RF guidance is then gradually phased out in favor of the more detailed local structural features learned by the GNN, achieved via a carefully designed curriculum learning strategy.

By integrating RF predictions through both direct feature fusion and knowledge distillation losses, the approach not only leverages the complementary strengths of classical machine learning and deep learning but also orchestrates a seamless transition from coarse global insights to fine-grained local understanding. This dynamic integration is more than an ensemble strategy; it underpins the training process itself, ensuring that the model converges to a robust representation that is both interpretable and highly predictive.

Performance evaluation metrics

To rigorously evaluate the model’s performance, several standard classification metrics were employed, namely accuracy, precision, recall, F1-score, area under the receiver operating characteristic curve (AUC-ROC), and area under the precision-recall curve (AUC-PR). These metrics provide a comprehensive assessment of the model’s ability to correctly predict drug–protein interactions while accounting for the trade-offs between precision and recall, especially under class imbalance.

Accuracy represents the proportion of correct predictions (both positive and negative) among all predictions made:

$$\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}.}$$

Precision is defined as the proportion of correctly predicted positive interactions among all interactions predicted as positive:

$$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}={\displaystyle \frac{TP}{TP+FP}.}$$

Recall quantifies the proportion of true positive interactions that were correctly identified:

$$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}={\displaystyle \frac{TP}{TP+FN}.}$$

The F1-score, which balances precision and recall, is computed as the harmonic mean of the two:

$$\mathrm{F}1\text{-}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=2\times {\displaystyle \frac{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\times \mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}.}$$

To further characterize overall discriminative performance, the AUC-ROC metric was employed. The AUC-ROC measures the model’s ability to distinguish between positive and negative classes across varying thresholds. Additionally, the AUC-PR metric was computed, as it is more sensitive to performance on imbalanced datasets. The AUC-PR evaluates the trade-off between precision and recall across different threshold settings and provides a more informative summary in settings where true interactions are relatively rare.

Here, TP (true positives) refers to correctly predicted positive drug–protein interactions, FP (false positives) to incorrectly predicted positives, and FN (false negatives) to positive interactions missed by the model. Together, these metrics enable a robust and nuanced evaluation of model effectiveness, balancing both the ability to recover true interactions and to minimize false discoveries.

Model selection rationale

The selection of a graph-based deep learning framework was motivated by the inherent relational structure of drug–protein interaction data, which naturally lends itself to graph representation. By modeling drugs and proteins as nodes in similarity networks, the architecture can effectively capture both molecular structure and relational context. GNNs were therefore selected to enable localized message passing and representation learning across these networks.

The RF classifier was incorporated as a teacher model during training to guide the learning process through knowledge distillation. This choice was based on its strong baseline performance using traditional molecular descriptors, as well as its interpretability and robustness. The integration of RF guidance allowed the model to benefit from handcrafted feature knowledge without requiring a full multi-modal ensemble.

The decision to maintain a single-modality architecture was made, with the goal of evaluating the effectiveness of a focused, structure-guided model without the added complexity of multi-modal integration. This design simplifies implementation and enhances interpretability while still benefiting from molecular feature knowledge through the distillation process. By embedding guidance from the RF teacher during training, the model is able to leverage informative patterns without requiring direct access to raw molecular descriptors at inference time.

Interpretability analysis with integrated gradients

To examine which molecular substructures and protein regions drive SGCL-DPI predictions, we applied integrated gradients (IG), a gradient-based attribution method. IG was computed on the attention-weighted graph embeddings, highlighting the relative contribution of each atom and residue to the final output. The visualization pipeline converts SMILES strings into molecular diagrams and color-codes atomic contributions (red for strong positive attribution, blue for negative). This allows direct inspection of how the model propagates information across the drug and protein graphs, providing insight into its decision-making process. Details are provided in Article S1.

Evaluation methods

To ensure a thorough and transparent assessment, the following evaluation strategies were applied in this study: (1) Baseline comparison using a RF model trained on traditional molecular descriptors, (2) generalization testing through the “hard split” protocol introduced by Watanabe, Ohnuki & Sakakibara (2021), (3) comparative analysis against single- and multi-modality baselines reported in prior literature, (4) ablation studies to examine the contribution of individual loss components and architectural elements, and (5) curriculum learning sensitivity analysis to explore the effect of varying the knowledge distillation weight (Alpha).

The performance and generalization capability of the proposed SGCL-DPI framework were evaluated through a comprehensive set of experimental strategies. As an initial baseline, a RF model trained on traditional molecular descriptors was employed. This model, built using established cheminformatics features, served as a classical benchmark and was assessed across 10 random splits of the BindingDB dataset to provide robust statistical averages for comparison.

To test model robustness under realistic and challenging conditions, the evaluation incorporated the “hard split” protocol introduced by Watanabe, Ohnuki & Sakakibara (2021). In this setting, no drugs or proteins appearing in the training folds are allowed in the test folds, thereby simulating scenarios involving entirely novel molecular entities. This approach is widely recognized for its rigor and relevance to drug discovery tasks.

Comparative analysis was conducted using baseline models reported in the literature, particularly those from Watanabe, Ohnuki & Sakakibara (2021). These included single-modality models utilizing only molecular features or only network-based embeddings, as well as a more complex integrated model that combines both data types. This comparison established a contextual benchmark for assessing the SGCL-DPI model’s standing among contemporary deep learning approaches.

To isolate the contribution of each architectural and training component, several ablation studies were performed. These experiments examined model variants trained using only individual loss components including weighted binary cross-entropy (BCE), knowledge distillation (KD) from the RF teacher, or structure consistency loss, as well as a version that excluded the graph encoder. The resulting performance differentials provided insight into the importance of each component.

The effect of curriculum learning parameters was also investigated by varying the weighting factor (Alpha) associated with the KD loss during training. Different Alpha values were tested to determine the optimal balance between guidance from the RF teacher and autonomous learning within the graph-based architecture.

All model variants were evaluated on held-out validation sets using standard classification metrics, including accuracy, precision, recall, F1-score, AUC-ROC, and AUC-PR. These metrics offered a comprehensive view of each model’s predictive behavior across different evaluation scenarios.

Performance and scalability

Graph construction represents the primary computational bottleneck in SGCL-DPI, as it requires exhaustive pairwise similarity calculations across all drug–drug and protein–protein combinations. To mitigate this cost, precomputed similarities are cached and stored using efficient data structures; RDKit handles fingerprint generation, while NumPy is used for descriptor distance computations. Parallel execution on an NVIDIA A100 GPU reduces this preprocessing to a one-time operation, after which persisted graphs enable rapid training. On a dataset with tens of thousands of interactions, SGCL-DPI completes over 100 epochs using a cosine learning-rate schedule in approximately 16 min per run on the same GPU. Evaluation on a standard laptop (16 GB RAM, quad-core Intel Core i7) confirms that, with precomputed graphs, both training (with a reduced number of epochs) and inference are feasible without specialized hardware. The PyTorch/PyG implementation supports fully vectorized GNN operations and dynamic computation of the multi-component loss, maintaining scalability and efficiency in typical research settings. In a production setting, graph construction can be incremental, only the similarity rows and columns for newly added compounds or targets need to be computed, so the preprocessing cost is spread over time and typically falls to a few seconds per update.

Random forest baseline performance: revealing the potential of traditional molecular descriptors

As shown in Table 1, the RF classifier using traditional molecular descriptors achieved exceptionally high performance on the BindingDB benchmark. In particular, the RF model attained an average AUC of 99.18% and AUPR of 99.14% across 10 random splits, along with a Precision of ~96.8% and Recall of ~95.7%, yielding an F1-score of ~96.2%. This outperforms several deep learning models reported in the literature on the same dataset. For example, GraphDTA (which employs graph neural networks on molecular graphs and CNNs on protein sequences) achieved an AUC of about 93.6% and AUPR of 93.4%. Similarly, the TransformerCPI model (a transformer-based sequence-only approach) reached roughly 95.7% AUC and 95.8% AUPR. Even BridgeDPI reported an AUC/AUPR of 97.5%/97.3%, which, while high, is still lower than the RF’s performance in our experiments.

Notably, these results were achieved without using any graph neural networks, protein sequence models, or other complex multi-modal frameworks. The RF model relied solely on well-established molecular descriptors (traditional chemical features), yet it matched or exceeded state-of-the-art performance. This highlights the surprising strength of carefully engineered classical descriptors in this context. They appear to capture the key determinants of binding interactions effectively enough that a simple ensemble classifier can leverage them to rival deep learning models. In other words, for the BindingDB dataset, the added complexity of deep neural networks (GNNs, transformers, etc.) did not translate into better performance than the descriptor-based approach, underlining how informative the descriptors are for this task.

Comparison with Random Forest baseline–hard dataset

To better assess robustness under more realistic and demanding conditions, we evaluated SGCL-DPI and baseline methods on a harder test set. On the hard test set (see Materials and Methods), the proposed graph-based model achieves an accuracy of ~76.3%, slightly higher than the RF baseline’s ~75.3% (Table 2). More importantly, our model exhibits a much better balance between precision and recall. It attains recall = 73.00% (±6.68), meaning it correctly recovers a large fraction of true drug–protein interactions. This recall is over 2.2× higher than that of the RF baseline (32.64% ± 9.53), which struggled to identify positive interactions. The RF model’s precision is indeed higher (82.52% vs. 62.24%), indicating it makes fewer false-positive predictions, but this comes at the cost of missing most true interactions. The proposed model’s more moderate Precision (62.24% ± 3.25) coupled with its strong recall yields a substantially better F1-score (67.06% ±3.50) compared to the RF’s F1 (46.20% ± 10.35). In other words, our method achieves a more favorable precision–recall trade-off, improving F1 by ~21 points. This suggests that the graph-based approach is more effective at finding true interactions without being overly conservative. We note that the AUC-ROC values of the two models are practically equivalent (≈82.3% vs. 82.1%), indicating that overall ranking ability (in terms of true-positive rate vs. false-positive rate across thresholds) is comparable. However, as Watanabe, Ohnuki & Sakakibara (2021) emphasize, the AUC-PR is a more informative metric on class-imbalanced problems like this (where true interactions are relatively scarce) (Watanabe, Ohnuki & Sakakibara, 2021). By this measure, our model’s performance (AUC-PR ≈71.7%) is on par with the RF (72.6%), suggesting that both models achieve a similar area under the precision–recall curve despite their different operating points. The key distinction is that our model operates at a higher-recall point on that curve, capturing far more positives. This balanced performance is critical in practical settings where missing potential interactions (false negatives) is as problematic as having some false positives.

Critically, we acknowledge that the improved recall of the proposed model comes with somewhat lower precision than the RF. This indicates a higher false-positive rate, which is a trade-off to gain more true hits. Depending on the application, one might adjust the decision threshold to tune this balance; however, the substantially higher F1 of our method demonstrates that overall it achieves a better compromise between precision and recall than the RF baseline. A threshold sweep (see Fig. S1) shows that while both models reach a similar peak F1 of ≈ 0.70 at a narrow threshold (~0.25), the Random Forest baseline’s F1 drops sharply thereafter, whereas SGCL-DPI maintains F1 ≥ 0.50 across a broad threshold range (0.25–0.70), demonstrating that the balanced precision–recall performance is not tied to a single threshold. In summary, compared to a traditional machine learning baseline using the same input features, the proposed graph-based model delivers more robust and comprehensive predictions, identifying many more true interactions while maintaining reasonable precision.

Comparative analysis with baseline models and Watanabe, Ohnuki & Sakakibara’s (2021) approaches

The work of Watanabe, Ohnuki & Sakakibara (2021) was selected as a benchmark due to its rigorously defined hard dataset and comprehensive evaluation of deep learning models. While our RF-based approach outperformed several recent GNN methods on BindingDB, Watanabe’s framework presented a more challenging and widely recognized generalization test. Although reproducing their experiments was not feasible due to difficulties running their published code, their dataset and methodology remain influential and were instrumental in guiding the comparative evaluation in this study.

The performance of the proposed graph-based model was compared against the benchmarks reported by Watanabe, Ohnuki & Sakakibara (2021) on the hard split derived from the STITCH database (see Fig. 2). Watanabe, Ohnuki & Sakakibara (2021) presented three variants in their work: an integrated model (combining molecular sequence/structure features with network context), a single-modality model based on molecular features alone, and a network-only model. Notably, the present work focuses on a single-modality approach, wherein the drug–protein interaction graph serves as the sole data source; however, it benefits indirectly from molecular features via a RF module integrated through curriculum learning. Unlike traditional ensemble methods, the RF component in this approach is not used for independent prediction; instead, its output is used for knowledge distillation and feature fusion during training, thereby informing the graph-based predictor without requiring the simultaneous optimization of heterogeneous data sources.

Performance against single-modality network baseline

The proposed graph-based model markedly outperforms Watanabe, Ohnuki & Sakakibara’s (2021) single-modality network approach in most metrics. In particular, it achieves a higher F1-score (67.06% vs. 60.1%; an improvement of about 11.6%) and AUPR (71.69% vs. 70.6%) than the purely network-driven model of Watanabe, Ohnuki & Sakakibara (2021) despite a slightly lower overall accuracy (76.26% vs. 78.4%) (see Table 2). This improvement indicates that our method’s incorporation of molecular information provides an advantage over an interactome-only strategy. By integrating knowledge from molecular descriptors, the proposed method detects more true drug–protein interactions (recall 73.00% vs. an implied lower recall for the network baseline) while maintaining a reasonable precision, thus yielding a better balance between precision and recall. In contrast, Watanabe’s network-only model, which relies solely on similarity networks of proteins and compounds, appears to miss many interactions (lower recall), resulting in a diminished F1-score. These gains indicate that the graph-based formulation, even when operating as a single-modality approach, is able to capture more nuanced patterns of interactions—most likely by exploiting local graph topology and higher-order connectivity—than a pure network embedding method.

Performance against single-modality molecular baseline

The SGCL-DPI framework demonstrated competitive performance when compared with the single-modality molecular baseline reported by Watanabe, Ohnuki & Sakakibara (2021) which leveraged protein sequence and compound structure features without incorporating interaction networks. On the STITCH hard split, SGCL-DPI achieved an F1-score of 67.1%, slightly surpassing the molecular baseline’s 66.2%, representing an improvement of approximately 1.4%. Although the molecular baseline reported higher AUROC (85.1% vs. 82.3%) and AUPR (77.0% vs. 71.7%), the higher F1-score of SGCL-DPI highlights a more favorable balance between precision and recall. This advantage is particularly noteworthy considering that SGCL-DPI operates without direct access to raw molecular sequences or structures, relying instead on similarity networks and knowledge distillation from a RF teacher model. These findings underscore the effectiveness of SGCL-DPI’s graph-based design and training strategy in capturing predictive interaction patterns with minimal feature complexity.

SGCL-DPI vs. integrated multi-modal models

Watanabe, Ohnuki & Sakakibara’s (2021) integrated model, which fuses molecular sequence/structure features with network context, achieved an AUROC of 88.2% (±3.5), an AUPR of 83.4% (±4.1), and an F1-score of 71.4% (±6.4). Although the integrated model outperforms the proposed approach in absolute terms, it should be noted that the present study deliberately adheres to a single-modality design. The goal is to demonstrate that even without the complexity of integrating multiple data modalities (i.e., protein sequences, chemical structure features, and extensive interactome data) a focused graph-based method, when augmented via curriculum learning with an RF-based molecular guidance, can yield competitive performance. In fact, the proposed model’s F1-score (67.06%) attains approximately 94% of the value achieved by the integrated model, and its AUROC is within 6 percentage points. These findings highlight that a streamlined, single-modality approach can capture a significant portion of the predictive signal, while preserving simplicity and reducing computational and architectural complexity.

While the obvious extension of integrating complementary data types (e.g., through fully ensemble methods) remains a promising avenue for future research, the current work demonstrates that a graph-based method using only the drug–protein interaction network enhanced indirectly with molecular feature guidance via a curriculum learning process achieves a robust balance between precision and recall. The strength of the method lies in its ability to substantially outperform conventional network-only methods and to deliver results competitive with more complex multi-modal systems, despite not utilizing a full ensemble strategy. This single-modality approach offers practical advantages including reduced model complexity, easier interpretability, and greater ease of deployment in resource-constrained settings.

Loss function ablation and model variants

To evaluate the contribution of each loss component and the graph encoder, we conducted ablation experiments on the hard split of the STITCH dataset. Table 3 summarizes key performance metrics (accuracy, precision, recall, F1-score, AUC-ROC, AUC-PR) for the full multi-component model (using weighted binary cross-entropy (BCE), knowledge distillation (KD) from the Random Forest, and structure consistency loss) compared to several ablated variants using only weighted BCE loss, only KD loss, only structure consistency loss, and a no-graph variant combining BCE + KD losses without the graph encoder (no GNN-based structure learning).

From these results, the full model clearly achieves the strongest overall performance, with the highest F1-score (67.06%) and AUC-ROC (82.33%). This indicates that all loss components together produce a synergistic effect, improving both precision and recall compared to any single-loss model. In contrast, using only one loss type led to significant drops in performance. For example, a model trained only with weighted BCE (i.e. using graph-based learning on labels alone, without KD or structure loss) reached an F1 of 57.82%, notably lower than the full model.

The KD-only model (trained solely by distilling the Random Forest’s predictions, without direct supervision on the interaction labels) achieved perfect recall (100%) at the expense of lower precision (33.22%). By predicting a large number of positives, it successfully captured nearly all true interactions, but at the cost of introducing many false positives. This behavior suggests that while knowledge distillation can effectively transfer general decision patterns, relying solely on the teacher’s outputs without ground-truth labels can cause the student model to overgeneralize. These findings highlight the importance of complementing distillation with task-specific supervision to achieve more balanced and reliable predictions. Notably, the KD-only model’s AUC-ROC (82.32%) and AUC-PR (73.37%) were very close to those of the full model, indicating that the teacher’s knowledge remains valuable, particularly for ranking predictions even when not directly trained on the interaction labels.

Including the structure consistency loss on its own (without BCE or KD) resulted in moderate predictive performance (Table 3). This variant achieved a relatively high accuracy (77.87% ± 1.84) and precision (76.62% ± 6.03), but much lower recall (48.73% ± 3.7), indicating a tendency to make confident but conservative predictions. The F1-score (59.39% ± 2.75) reflects this imbalance. While the model performed reasonably well in ranking (AUC-ROC: 81.93% ± 1.83, AUC-PR: 71.93% ± 3.86), its lack of direct supervision likely limited its ability to generalize across diverse interaction types. This suggests that while the structure loss offers useful regularization, it is insufficient on its own to drive robust interaction prediction. These findings reinforce the role of the structure loss as a complementary signal, best used in conjunction with a primary supervised objective such as BCE or knowledge distillation.

Crucially, removing the graph encoder degraded performance significantly. The No GNN variant (which uses both BCE and KD losses but no graph-based similarity learning) achieved an F1 of 43.74% and AUC-PR of 37.59%, considerably lower than the full model. Both precision and recall dropped relative to the full model, indicating that structural learning via the GNN is important for generalization. Without the graph, the model cannot leverage similarities between drugs or between proteins; effectively it must treat each entity in isolation (aside from what the RF teacher provides). This leads to a notable loss in recall on the hard split, because the model struggles with novel drug–protein combinations that were not seen during training. In contrast, the full graph-based model can propagate interaction signals across similar drugs and proteins, enabling it to catch interactions involving new compounds by analogy to known ones. For example, if a new drug has a close structural analog in the training set that interacts with a given protein, the GNN can transfer that relational signal through the drug similarity graph, something the no-graph model cannot do. The higher AUC-PR for the full model compared to the no-GNN model (71.69 vs. 37.59) reflects this advantage, especially in identifying the minority positive class instances. We also ablated the knowledge-distillation term while keeping the GNN intact (see row ‘BCE only’ in Table 3). Eliminating RF guidance lowered F1 on STITCH-hard from 67.06 ± 3.5% to 57.82 ± 5.8% and reduced Recall by 14%, confirming that a descriptor-first curriculum is critical for capturing difficult, cross-domain interactions.

Overall, these ablations demonstrate that each component of the proposed approach contributes meaningfully to performance. The weighted BCE loss serves as the core learning signal for interaction prediction, the KD loss transfers useful guidance from the RF teacher to enhance both precision and recall, and the structure consistency loss adds a regularization effect that supports more stable learning. Together, these components enable the model to strike a better balance between false positives and false negatives than any individual part alone. By combining teacher knowledge with graph-based structural learning, the full model achieves improved generalization and a more effective precision–recall trade-off on this challenging prediction task.

Importantly, the proposed model attains this performance with a single-modality architecture. It does not incorporate the chemical descriptors or protein features directly into a multi-branch neural network, nor does it ensemble the RF with the GNN at inference. Instead, molecular feature knowledge is infused through distillation during training, allowing the final model to remain focused on graph-based representations. This design means that at test time the model uses only the similarity graphs (and no external features or ensemble voting), yet it performs on par with or better than approaches that explicitly fuse multiple data types. The strong F1 and AUC-PR achieved without late-stage feature fusion highlight the efficiency of this strategy. We gain the informative signal of descriptors in an interpretable way (via the teacher’s guidance) while avoiding the complexity of full multi-modal integration.

Impact of curriculum learning parameters on model performance

We observed differential performance across three curriculum learning setups (Alpha = 20, 35, 50) on drug–protein interaction prediction F1-scores, highlighting the influence of curriculum parameter tuning (see Fig. 3). In this context, Alpha denotes the weighting factor applied to the knowledge distillation loss from the RF model during the initial training phase. Larger Alpha values place stronger emphasis on aligning the graph neural network’s predictions with those of the RF teacher, thus providing more global guidance at early stages of learning. These values were chosen based on initial experimentation, where smaller Alpha values failed to provide adequate RF supervision to guide early learning.

Among these, the intermediate setting (Alpha = 35) yielded the highest single-fold F1-score (70.73% on Fold 2) as well as the highest overall mean performance (~66.3%, vs. ~63–64% for the other setups). These results suggest that a moderately paced curriculum may be optimal for this task. However, the Alpha = 35 configuration also showed greater variability between folds (e.g., 64.45% on Fold 1 vs. 70.73% on Fold 2) relative to the other configurations (with ~6–7 percentage-point F1 swings across folds, compared to ~3–5 for the others), indicating a sensitivity to training data splits or initial conditions. In comparison, the more conservative curriculum (Alpha = 20) produced more consistent F1-scores across folds (65.26%, 62.94%, 61.76%) but did not achieve the same peak performance, while the more aggressive curriculum (Alpha = 50) had intermediate outcomes (60.64–66.46% across folds). Notably, each fold’s highest score was achieved with a different Alpha value (Alpha=20 in Fold 1, 35 in Fold 2, 50 in Fold 3), underscoring the non-monotonic influence of the curriculum parameter on performance. This sensitivity and the fact that only the appropriately tuned value (Alpha = 35 in this experiment) delivered a pronounced performance gain underscore the importance of careful hyperparameter optimization. Although these results are preliminary, they demonstrate that tuning curriculum learning parameters can yield meaningful gains in predictive performance. Given the observed sensitivity, a more exhaustive hyperparameter search could be pursued in the future to further refine these gains and validate the robustness of the approach.

Qualitative model interpretation

Atom-level attributions were next examined using IG (see Article S1). In a true-positive example, the attributions appeared sharply localized with high intensity over biologically meaningful functional groups, including a phosphate moiety, substituted aromatic rings, and polyhydroxy regions. This pattern suggests that the model bases correct predictions on chemically relevant substructures. In contrast, a false-positive case revealed diffuse, low-contrast attributions scattered across non-specific hydroxyl chains, lacking any dominant substructure to support the model’s high predicted probability which is indicative of over-generalization from widely occurring motifs. Collectively, these examples demonstrate how interpretability analyses can differentiate between well-grounded and unreliable predictions in SGCL-DPI.