Relation-aware graph structure embedding with co-contrastive learning for drug-drug interaction prediction

We construct a multi-relational DDI graph $\mathcal{G}=(\mathcal{D},\mathcal{E},\mathcal{R})$, where nodes represent drugs, and edges denote labeled interactions. Herein, $\mathcal{D}$ is the set of all drugs (including known and new drugs), and $|\mathcal{D}|=N$ represents the number of drugs. Let ${\bf X}=\left\{{\bf x}_{i}\right\}_{i=1}^{N}\in\mathbb{R}^{N\times d}$ be the initial feature matrix of drugs, derived from Section 2.3. $\mathcal{R}$ is the set of interaction event types, and $|\mathcal{R}|=R$ is the number of interaction type. $\mathcal{E}=\left\{\mathcal{E}_{r}\right\}_{r=1}^{R}$ is the set of known interactions, and $\mathcal{E}_{r}$ represents the set of interactions under interaction type $r$. Let $\mathcal{A}=\left\{{\bf A}_{r}\right\}_{r=1}^{R}\in\mathbb{R}^{N\times N\times R}$ be the multi-relational adjacency tensor, where ${\bf A}_{r}\in\mathbb{R}^{N\times N}$ represents the adjacency matrix under interaction relation $r$. Let $A_{r(i,j)}$ $(i,j=1,\dots,N)$ be the element of $\bf{A}_{r}$, $A_{r(i,j)}=1$ if $(i,j)\in\mathcal{E}_{r}$ and $A_{r(i,j)}=0$ if $(i,j)\notin\mathcal{E}_{r}$.

Subsequently, we use the R-GCN layer [34] to learn the RaGSEs of drugs from the multi-relational DDI graph. The forward propagation function is defined as follows:

where $\mathbf{h}_{i}\in\mathbb{R}^{1\times d^{\prime}}$ represents the RaGSEs of drug $i$, $\mathbf{x}_{i}\in\mathbb{R}^{1\times d}$ are the initial feature vector of drug $i$, $i\in\mathcal{D}$. Herein, $\mathbf{W}_{r}\in\mathbb{R}^{d\times d^{\prime}}$ represents the relation-specific weight matrix and the adoption of a set of $\mathbf{W}_{r}$ $(r=1,\dots,R)$ supports multiple edge types. On the contrary, $\mathbf{W}_{0}\in\mathbb{R}^{d\times d^{\prime}}$ is a single weight matrix regardless of relations. Here, $\mathcal{N}_{i}^{r}$ represents the set of drugs connected to drug $i$ under relation $r$. $\sigma(\cdot)$ is an element-wise activation function: $\operatorname{ReLU}(\cdot)=\max(0,\cdot)$. To normalize the incoming messages of each drug, we use $\left\{R_{i}\right\}_{i=1}^{N}$ to define a set of normalized constants. Here, $R_{i}$ equals the number of interaction types in which drug $i$ is involved. $\hat{A}_{r(i,j)}$ is the aggregation weight assigned to the neighboring drug $j$. Specifically, $\hat{A}_{r(i,j)}$ denotes the element in $\hat{\bf A}_{r}$. The calculation of $\hat{\bf A}_{r}$ is based on a classic graph-based normalization method [19].

where ${\bf D}_{r}$ is the degree matrix of ${\bf A}_{r}$ and a diagonal matrix.

Herein, $\mathbf{h}_{i}$ can represent relation-aware graph structure information of drug $i$ if drug $i$ is known. Nevertheless, if drug $i$ is new, $\mathbf{h}_{i}$ would be influenced only by drug $i$ itself. Consequently, directly using the output of Eq. (2) as the drug embeddings for DDI prediction could potentially result in severe over-fitting issues for Tasks 2 and 3.

Based on the assumption that similar drugs may interact with the same drugs [25, 5], we propose a strategy that propagates the RaGSEs of known drugs to new drugs with similar features. This strategy can effectively overcome the over-fitting issues that arise in Tasks 2 and 3. To facilitate the propagation of embeddings, we construct a multi-attribute DDS graph $\mathcal{G}_{s}=(\mathcal{D},\mathcal{E},\mathcal{S})$, where nodes represent drugs, and edges denote Jaccard similarities between drugs under various attributes. Herein, we still employ these three biological attributes (i.e., chemical substructure, enzyme, and target). In the DDS graph, $\mathcal{D}$ is the set of nodes, with ${\bf H}=\left\{{\bf h}_{i}\right\}_{i=1}^{N}$ as the feature matrix of nodes, which is the output of Eq. (2). $\mathcal{E}$ represents the set of edges. $\mathcal{S}$ is the set of similarity types, and $|\mathcal{S}|=3$ is the number of similarity types. Let $\mathcal{A}=\left\{{\bf A}_{s},{\bf A}_{e},{\bf A}_{t}\right\}\in\mathbb{R}^{N\times N\times 3}$ be the adjacency tensor, where ${\bf A}_{s}$, ${\bf A}_{e}$, and ${\bf A}_{t}$ represent the substructure-based, enzyme-based, and target-based adjacency matrix, respectively.

To normalize the incoming messages, these three adjacency matrices also need to be normalized. The normalization method is referred to as Eq. (3). Let $\hat{\bf A}_{s}$, $\hat{\bf A}_{e}$, and $\hat{\bf A}_{t}$ be the normalized adjacency matrices. We use these normalized adjacency matrices to propagate the learned RaGSEs to the strongly associated drugs within the DDS graph. The propagation procedure can be expressed as follows:

where $\hat{\bf A}^{n}_{s}\bf{H}$, $\hat{\bf A}^{n}_{e}\bf{H}$, and $\hat{\bf A}^{n}_{t}\bf{H}$ represent graph convolutions on the DDS graph. $n$ is a hyperparameter, and $\hat{\bf A}_{s}^{n}$ is $\hat{\bf A}_{s}$ to the power of $n$. The same definition is extended to $\hat{\bf A}_{e}^{n}$ and $\hat{\bf A}_{t}^{n}$. This strategy allows RaGSECo to capture higher-order neighbor information flexibly. $\mathbf{W}_{s}$, $\mathbf{W}_{e}$, $\mathbf{W}_{t}\in\mathbb{R}^{d^{\prime}\times d^{\prime}}$ are the attribute-specific trainable weight matrices. $\mathbf{H}_{s}$, $\mathbf{H}_{e}$, and $\mathbf{H}_{t}$ are the propagation results along with different pathways. $\mathbf{H}_{embed}=\left\{{\bf h}_{embed,i}\right\}_{i=1}^{N}\in\mathbb{R}^{N\times d^{\prime}}$ is the sum of the three propagation results and represents the final embeddings of drugs. Particularly, ${\bf h}_{embed,i}$ can represent effective RaGSEs of drug $i$, whether drug $i$ is known or new.

To increase data diversity, facilitate gradient propagation, and generate discriminative DP representations, we employ multiple drug features containing RaGSEs, initial features, and Simplified-Molecular-Input-Line-Entry-System (SMILES) string.

SMILES string is a line notation that uses a predefined set of rules to describe the structure of compounds, and each drug has its own SMILES string [4]. The characters of a SMILES string represent chemical atoms or chemical bonds, and the lengths of the SMILES strings are unconstrained [31]. In this work, we convert each SMILES string to a $p\times q$ dimensional feature matrix for simplicity, where $p$ denotes the number of characters and $q$ is the unified length of the SMILES string [16]. Accordingly, the columns of the SMILES-based feature matrix are one-hot vectors. Since the initial lengths of the SMILES strings are flexible, we will cut off the extra part if the actual length of the SMILES string is greater than $q$, and we will pad zeros if the actual length is less than $q$. The $p$ and $q$ are set as $64$ and $100$ according to Huang et al. [16].

Let $Q$ be the number of the known DDIs. Considering computational burdens, we adopt a batch-wise scheme to train the RaGSECo. Let $K$ be the number of DPs in a batch. Given DP $k$ $(k=1,\dots\,K)$ that involves two drugs $i$ and $j$ $(i,j=1,\dots\,N;i\neq j)$. Let ${\bf{S}}_{i}$, ${\bf{S}}_{j}\in\mathbb{R}^{p\times q}$ be the SMILES-based feature matrices of drugs $i$ and $j$. ${\bf{h}}_{embed,i}$ and ${\bf{h}}_{embed,j}$ are the RaGSEs of two drugs. ${\bf x}_{i}$ and ${\bf x}_{j}$ are initial feature vectors. Subsequently, we employ two feedforward neural networks (FNNs) to process the RaGSEs and initial features of DPs, respectively. In addition, a CNN encodes SMILES strings of DPs. The employed CNN is a multi-layer 1-D CNN followed by a global max pooling layer, and with reference to [16]. Accordingly, three features of DP $k$ are encoded as follows:

where the symbols $\Theta_{\mathrm{FNN1}}$, $\Theta_{\mathrm{FNN2}}$, and $\Theta_{\mathrm{CNN}}$ represent the trainable weights involved in FNN1, FNN2, and CNN, respectively.

The initial feature vectors of drugs are denoted by Jaccard similarity scores between drugs. Thus, ${\bf p}_{k}^{initi}\in\mathbb{R}^{1\times d^{FNN}}$ can be considered as the DP feature representations under the similarity view. Meanwhile, the drug RaGSEs primarily concentrate on interaction information between drugs. As such, ${\bf p}_{k}^{embed}\in\mathbb{R}^{1\times d^{FNN}}$ can be perceived as the DP representations under the interaction view.

The interaction view mainly focuses on known interaction relationships between nodes, while the similarity view can infer potential therapeutic effects, adverse reactions, or drug interactions of new drugs by discovering underlying drug correlations. Therefore, these two views are complementary and mutually reinforcing. Herein, we define positive and negative samples for contrastive learning. In this paper, given the feature representations of a DP under the interaction view, we can simply define its feature representations under the similarity view as the positive sample. Nevertheless, our analysis of the employed datasets suggests that there may be underlying correlations between DPs. Consequently, we propose a novel positive selection strategy.

This paper utilizes two real-world, multi-relational DDI datasets: Dataset 1 and Dataset 2. In Dataset 1, the known DDIs are categorized into 65 types of DDI events, while the number of interaction event types in which each drug is involved ranges from 1 to 20. The median of the event type counts is 10, significantly lower than 65. Similar findings are also observed in Dataset 2. These findings suggest that a particular drug is likely to associate with specific interaction event types strongly. Thus, the interaction event types a drug participates in can be considered significant characteristics of the drug. We refer to these significant characteristics as interaction characteristics of drugs. Let ${{\bf t}_{i}}\in{\mathbb{R}}^{1\times R}$ represent the interaction characteristics of drug $i$, where $R$ is the number of interaction event types. The vector ${\bf t}_{i}$ is binary, and its nonzero elements represent the interaction event types the drug $i$ participates in. Given a DP $k$ that involves two drugs $i$ and $j$, the interaction characteristics of DP $k$ are defined as ${\bf c}_{k}={\bf t}_{i}+{\bf t}_{j}$. Based on this, we propose a new positive selection strategy: if two DPs exhibit similar interaction characteristics, they are designated as positive samples. One advantage of such a strategy is that the selected positive samples may reflect the underlying interaction tendency of the target DP. For DPs $k$ and $q$, we calculate the cosine similarity of their interaction characteristics as follows:

Given the threshold values $t_{pos}$ and $t_{neg}$, where $t_{pos}>t_{neg}$. DP $q$ is deemed a positive sample of DP $k$ if ${S}_{k,q}\geq t_{pos}$, and conversely, is considered a negative sample if ${S}_{k,q}\leq t_{neg}$. Let $\mathcal{P}$ represent the set of positive pairs of one batch, and $\mathcal{N}$ be the set of negative pairs. We then model a self-supervised learning task using the standard binary cross-entropy loss. With this, we present the co-contrastive learning loss function as follows:

where the contrastive discriminator $\Psi(\cdot,\cdot)$ is instantiated as $\sigma({\bf p}_{k}^{embed}{\bf W}{{\bf p}_{q}^{initi}}^{T})$. ${\bf W}$ is a learnable weight matrix. The activation function $\sigma$ is the sigmoid function, which produces a score representing the probability of being a positive sample. The co-contrastive learning mechanism maximizes the mutual information of positive pairs under the interaction and similarity views. To stabilize the co-contrastive learning process and facilitate two views supervising each other, we present another co-contrastive learning loss function:

In conclusion, to fully utilize all significant information and learn more discriminative DP representations, we define the final DP representations ${\bf p}_{k}$ as follows:

where the symbol $||$ denotes the concatenation.

The hyperparameter $n$ represents the power of adjacency matrices in the DDS graph and determines the distance of RaGSE propagation. It plays a crucial role in RaGSE propagation. To understand the influence of $n$ on the prediction performance, we conduct experiments on three tasks of Dataset 1 and evaluate six metric scores across different $n$ values. Referring to Fig. 4 (a), in Task 1, the metric scores of RaGSECo do not change drastically as $n$ increases, and RaGSECo can consistently achieve good prediction performances. On the other hand, in Tasks 2 and 3, we observed that the prediction performance of RaGSECo is sensitive to the value of $n$. Specifically, when $n$ is set to 0, the prediction performance of our RaGSECo is at its lowest. The prediction performance significantly improves as the value of $n$ increases. This is because the test DDIs in Tasks 2 and 3 involve new drugs that do not possess relation-aware interaction information when $n$ is 0. As a result, the model tends to be over-fitting. As $n$ increases, new drugs can aggregate effective relation-aware interaction information from similar known drugs, which significantly improves the generalization ability of RaGSECo. Based on the experimental results, we determined that setting $n$ as 0 in Task 1 and 3 in Tasks 2 and 3 leads to desirable prediction performances.

The output dimension $d^{FNN}$ of FNN1 and FNN2 is an important hyperparameter for our RaGSECo model. Increasing $d^{FNN}$ can enhance the generalization ability of RaGSECo to some extent. In order to investigate the impact of $d^{FNN}$ on the prediction performance of RaGSECo, we conducted experiments on three tasks of Dataset 1 and examined the changes in metric scores as $d^{FNN}$ varies. From Fig. 5, it can be observed that the prediction performances of RaGSECo gradually improve as $d^{FNN}$ increase among all tasks. This indicates the robustness of our model. To strike a balance between accuracy and efficiency, we set $d^{FNN}$ to 1000 for Task 1 and 1500 for Tasks 2 and 3.

The threshold values $t_{pos}$ and $t_{neg}$ are important parameters for selecting positive and negative samples of DPs. In our experiments, we set $t_{pos}$ to the range of $[0.8,0.85,\dots,1]$ and $t_{neg}$ to the range of $[0,0.1,\dots,0.4]$. Fig. 6 presents the comparisons of three representative metric scores (ACC, AURR, and F1) for different values of $t_{pos}$ and $t_{neg}$. It can be observed that the results are not significantly affected by variations in $t_{pos}$ and $t_{neg}$. This indicates the stability of our co-contrastive learning strategy. Herein, the feature representations from the two views of each DP serve as positive samples for each other when $t_{pos}$ equals 1. As observed, the model produces superior prediction results when $t_{pos}$ is set to 0.95, as opposed to when it is set to 1. This observation indicates the effectiveness of our positive sample selection strategy. Finally, we set $t_{pos}=0.95$ and $t_{neg}=0.1$ for subsequent experiments.

To assess the competitiveness of our RaGSECo, we compared it with several state-of-the-art DDI prediction methods, namely MCFF-MTDDI, MDF-SA-DDI, RANEDDI, DDIMDL, DeepDDI, and DNN. Table 2 presents the metric scores achieved by these methods across the three tasks of Dataset 1. As observed, in most cases, the comparison results demonstrate that our RaGSECo outperforms the competitors in terms of performance metrics in three tasks. In Task 1, although RANEDDI also considers relation-aware graph structure information to learn drug embedding and achieve outstanding prediction performance, our RaGSECo achieves better results than RANEDDI. Specifically, the improvements of RaGSECo over RANEDDI are $2.33\%$, $1.81\%$, $3.74\%$, $2.83\%$, and $3.33\%$ in terms of ACC, AUPR, Precision, Recall, and F1, respectively. The reasons causing this phenomenon are manifold. On the one hand, RANEDDI ignores multiple drug-related attributes, such as targets, enzymes, and chemical substructures. This makes RANEDDI lack the ability to capture relationships beyond interactions between drugs. On the other hand, simply concatenating embeddings of two drugs to obtain drug-pair features also limits the generalization ability of RANEDDI. MDF-SA-DDI considers multiple attributes representing drugs and employs multiple drug fusion methods to learn drug-pair features. Nevertheless, our RaGSECo still has $1.63\%$, $1.01\%$, $2.83$, and $1.72\%$ improvements over MDF-SA-DDI with respect to ACC, AUPR, Recall, and F1, respectively. The main reason is that MDF-SA-DDI does not take advantage of specific interaction information between drugs.

For further insight, we investigate the performances of our RaGSECo and four competitive baselines for each event. Fig. 10 displays the AUPR and AUC scores of the five prediction models for each event on Task 1. As observed, RaGSECo produces greater AUPR and AUC scores than other methods in most event types. In most cases, RaGSECo can achieve a satisfactory result. All unsatisfactory prediction results are observed in low-frequency event types, such as $\#39$, $\#52$, and $\#64$, with frequencies of only $98$, $24$, and $10$ samples, respectively. The limited availability of training samples for these low-frequency event types may contribute to relatively poorer performance.

In Tasks 2 and 3, we compare our RaGSECo with five competitive prediction methods, i.e., MCFF-MTDDI, MDF-SA-DDI, DDIMDL, DeepDDI, and DNN. Since RANEDDI only focuses on the embedding of known drugs, RANEDDI can not be applied to Task 2 and Task 3. In Task 2 and Task 3, the test DDIs include new drugs, and the lack of interaction information on the new drugs weakens the generalization ability of models. Therefore, the prediction accuracy of models performed on Task 2 and Task 3 is lower than that on Task 1. Nevertheless, RaGSECo outperforms other competitors in most cases. The reasons are three folds: 1) RaGSECo enables all drugs, including new drugs, to capture effective relation-aware interaction information. 2) RaGSECo inventively combines initial features, SMILES information, and drug RaGSEs, enhancing the information diversity. 3) RaGSECo captures the underlying correlation between DPs and enables two views to supervise each other by co-contrastive learning.

Table 3 presents the performance metrics of the proposed RaGSECo and other outstanding approaches on three tasks of Dataset 2. Dataset 2 exhibits greater diversity compared to Dataset 1. It has more drugs, more DDIs, and a wider range of DDI event types. As observed, our RaGSECo can achieve better prediction performances than other competitors in most cases. In Task 1, compared with MCFF-MTDDI, the proposed RaGSECo has $0.73\%$, $2.25\%$, $1.47\%$, and $2.11\%$ performance improvements in terms of ACC, Precision, Recall, and F1, respectively. As observed in Tasks 2 and 3, our RaGSECo still acquires the best results compared with five competitive prediction models. The experiment results demonstrate the effectiveness and robustness of RaGSECo.