Predicting drug combination side effects based on a metapath-based heterogeneous graph neural network

Abstract

In recent years, combined drug screening has played a very important role in modern drug discovery. Generally, synergistic drug combinations are crucial in treatment for many diseases. However, the toxic side effects of drug combinations are probably increased with the increase of drugs numbers, so the accurate prediction of toxic side effects of drug combinations is equally important. In this paper, we built a Metapath-based Aggregated Embedding Model on Single Drug–Side Effect Heterogeneous Information Network (MAEM-SSHIN), which extracts feature from a heterogeneous information network of single drug side effects, and a Graph Convolutional Network on Combinatorial drugs and Side effect Heterogeneous Information Network (GCN-CSHIN), which transforms the complex task of predicting multiple side effects between drug pairs into the more manageable prediction of relationships between combinatorial drugs and individual side effects. MAEM-SSHIN and GCN-CSHIN provided a united novel framework for predicting potential side effects in combinatorial drug therapies. This integration enhances prediction accuracy, efficiency, and scalability. Our experimental results demonstrate that this combined framework outperforms existing methodologies in predicting side effects, and marks a significant advancement in pharmaceutical research.

Introduction

The use of combinatorial drug therapy or polypharmacy has become a primary clinical approach in modern disease treatment, with most hospital visits resulting in prescriptions for more than two types of drugs [1, 2]. The absorption and metabolism of different drugs within the body can lead to various negative interactions, such as enhancing or diminishing clinical efficacy and increasing the risk of drug interactions, and potential adverse drug reactions [3,4,5,6]. Therefore, predicting the side effects of multiple drugs and adequately preparing for them is critically important.

However, as the number of drugs increases, the potential combinations sometimes reach an explosive level, making precise prediction challenging, akin to finding a needle in a haystack. Mathematical models can effectively address this complexity by systematically analyzing and prioritizing treatment options, thereby aiding clinicians in making effective decisions and devising suitable treatment strategies [7, 8].

To address these challenges, pioneering computational methods have been developed to predict drug side effects. For instance, Atias et al. [9] first proposed a computational approach focusing on individual drugs. Building on this foundation, the field has seen advancements with models like the Decagon [10], which integrates multimodal networks including protein and drug-protein interactions to predict multiple drug side effects using a graph convolutional network and a tensor decomposition decoder. This model significantly outperforms traditional methods by 69% in accuracy. Subsequent enhancements by Xu et al. [11] with the Tripartite Information Propagation model, and Carletti et al. [12] through the incorporation of a graph attention network, have further refined the accuracy of these predictions by emphasizing different weights in side effect relationships within heterogeneous graphs. Liu et al. [13] developed SC-DDIS, introducing additional drug features into the model to optimize predictions further. Extended applications like those by Brandon et al. [14] utilize Decagon data to create drug-protein knowledge graphs employing learning frameworks such as DistMult and KBlrn for improved prediction outcomes. New methodologies like the TriVec [15] and SimVec [16] models respectively, focus on enhancing node initialization in knowledge graphs, illustrating ongoing innovation in the field. Most recently, Jinwoo et al. [17] presented a unified embedding prediction model based on knowledge graphs, setting a new direction for future drug research by improving the accuracy of multi-drug side effect predictions.

In response to these advancements, this study proposes two distinct heterogeneous graph neural network models tailored for single and combination drugs within heterogeneous information networks. These models are designed to enhance the learning of feature vector representations of drugs and side effect nodes, optimizing the prediction of combinatorial side effects. An in-depth evaluation of these models demonstrates their predictive effectiveness, underscoring the potential of heterogeneous graph embedding in drug research.

Heterogeneous graph embedding, essential for projecting nodes into a low-dimensional vector space, has evolved significantly with two main categories of models. The first category, including shallow embedding algorithms like those based on the skip-gram model, captures semantic similarities through proximity in low-dimensional space, assuming nodes with similar neighbors should be closely positioned [18,19,20,21,22,23,24, 57]. The second category encompasses Graph Neural Networks (GNNs), which employ complex encoders to model network structure and node features more deeply [25,26,27,28,29,30,31,32,33 ]. This category includes significant innovations such as the application of GNNs to heterogeneous graphs by models, each contributing uniquely to the field by handling different types of data and network dynamics. For our Single Drug–Side Effect Heterogeneous Information Network, integrating metapaths based on proteins and side effects, as well as their features, is crucial to better represent drug node features and enhance prediction accuracy.

In this paper, we introduce two innovative models: the Metapath-based Aggregated Embedding Model for Single Drug-Side Effect Heterogeneous Information Networks (MAEM-SSHIN) and the Graph Convolutional Network for Combinatorial Drugs and Side Effect Heterogeneous Information Networks (GCN-CSHIN). MAEM-SSHIN effectively extracts features from a heterogeneous information network, focusing on the side effects associated with individual drugs, while GCN-CSHIN simplifies the complex task of predicting multiple side effects between drug pairs by forecasting interactions between combinatorial drugs and individual side effects. The integration of these models into a unified framework significantly enhances the accuracy, efficiency, and scalability of predicting potential side effects in combinatorial drug therapies. Our experimental results demonstrate that this combined approach substantially outperforms existing methodologies in side effect prediction, representing a major advancement in pharmaceutical research. This achievement highlights the transformative potential of advanced computational models in enhancing drug safety and efficacy.

Methods

Notations and definitions

Heterogeneous Information Graph: Let G = (V, E,$\mathcal{A}$ , $\mathcal{R}$) be a graph, where V is the set of nodes of the graph G, E is the set of edges of the graph G, is the set of node types, each representing a distinct category of nodes based on their features or roles within the graph, and $\mathcal{R}$ is the set of edge types, each signifying a distinct kind of relationship or interaction between the nodes. If |$\mathcal{A}$| +| $\mathcal{R}$ |> 2, then the graph G is called a heterogeneous graph, which can also be referred to as a Heterogeneous Information Network (HIN).

Metapath: Given a heterogeneous graph G = (V, E, , $\mathcal{R}$), a metapath $P$ is a sequence composed of node types and edge types, represented as $A_{1} \to ^{{R_{1} }} A_{2} \to ^{{R_{2} }} \cdots \to ^{{R_{l} }} A_{l + 1}$, which can be abbreviated as ${A}_{1}{A}_{2}\cdots {A}_{l+1}$. Here, ${A}_{1},{A}_{2},\dots ,{A}_{l+1}\in \mathcal{A},{R}_{1},{R}_{2},\dots ,{R}_{l}\in \mathcal{R}$, describing a composite relationship $R={R}_{1}\circ {R}_{2}\circ \cdots \circ {R}_{l}$, where ∘ is a composite relation operator.

Metapath Instance: Given a metapath $P={A}_{1}{A}_{2}\cdots {A}_{l+1}$ in a heterogeneous graph $G=\left(V,E,\mathcal{A},\mathcal{R}\right)$, an instance of metapath P refers to a sequence of nodes in graph G that follows the pattern of metapath P, denoted as ${v}_{{1}_{i}}{v}_{{2}_{j}}\dots {v}_{{\left(l+1\right)}_{k}}$. Here, ${v}_{{1}_{i}},{v}_{{2}_{j}},\dots ,{v}_{{\left(l+1\right)}_{k}}\in V$, with ${v}_{{1}_{i}}$ being of type ${A}_{1}$, ${v}_{{2}_{j}}$ of type ${A}_{2}$, and ${v}_{{\left(l+1\right)}_{i}}$ of type ${A}_{l+1}$.

Metapath-based Neighbors: Let $G=\left(V,E,\mathcal{A},\mathcal{R}\right)$ be a heterogeneous graph and $P$ be a metapath. Define the set of neighbors ${\mathcal{N}}_{\mathcalligra{v}}^{{\mathcal{P}}}$ for a node $v$ based on metapath $P$ as:

$${\mathcal{N}}_{\mathcalligra{v}}^{\mathcal{P}}=\{u\in V:\text{there exists a path from }v\text{ to }u\text{ following instances of }P\}$$

Furthermore, if $P$ is symmetric, then include $v$ in ${\mathcal{N}}_{\mathcalligra{v}}^{{\mathcal{P}}}$, i.e.,$v\in {\mathcal{N}}_{\mathcalligra{v}}^{{\mathcal{P}}}$.

Metapath-based Graph: Given a heterogeneous graph $G=\left(V,E,\mathcal{A},\mathcal{R}\right)$ and a metapath $P$.

The metapath-based graph ${G}^{P}$ is defined as ${G}^{P}=\left(V,{E}^{P}\right)$, where $V$ is the set of all nodes in the original graph $G$, ${E}^{P}$ is the set of edges in the graph ${G}^{P}$, specifically defined as:

$${E}^{P}=\{\left(u,w\right)\in V\times V:u,w\in {\mathcal{N}}_{\mathcal{v}}^{\mathcal{P}}\text{ for some} \, v\in V\text{ and there exists a path following }P\text{ between }u\text{ and }w\}$$

If $P$ is symmetric, then ${G}^{P}$ is defined as a homogeneous graph, implying that all edges $\left(u,w\right)$ conform uniformly to the relations specified by $P$.

Heterogeneous Graph Embedding: Let G = (V, E, A, $\mathcal{R}$) be a heterogeneous graph. For each node type ${A}_{i}\in \mathcal{A}$, let ${V}_{{A}_{i}}\subseteq V$ denote the set of nodes of type ${A}_{i}$, and let ${X}_{{A}_{i}}\in {\mathbb{R}}^{\left|{V}_{{A}_{i}}\right|\times {d}_{{A}_{i}}}$ be the feature matrix associated with nodes of type ${A}_{i}$, where ${d}_{{A}_{i}}$ represents the dimensionality of the feature space for nodes of type ${A}_{i}$. Here, ${\mathbb{R}}$ denotes the real numbers.

Heterogeneous graph embedding aims to learn a $d$-dimensional vector representation for each node $v \in V$, denoted as $v\in {R}^{d}$. This representation should satisfy the condition that $d$ is significantly smaller than $\left|V\right|$, the total number of nodes in the graph.

Heterogeneous information network construction

In this study, we initially constructed six network structures using data from the bioinformatics database at Stanford University BioSNAP Datasets [34]. The sources and methods of construction for these networks are detailed as follows:

Combinatorial drugs and side effect network

The TWOSIDES database [35] revealed 4,649,441 combinatorial drugs and side effect relationships, involving 63,473 drug pairs, 645 drugs, and 1,317 side effects. Here, the term 'drug' includes both combinatorial and single drugs unless specified otherwise. The frequency distribution of combinatorial drugs with side effects is shown in Fig. 1, with 964 side effects occurring in at least 500 drug combinations, forming our final dataset for experimentation.

Fig. 1

Distribution of side effects across drug combinations and single drugs. A Distribution of side effects across different drug ranges. B Distribution of side effects across different combinatorial drug ranges. C Distribution of drugs across different ranges of side effects. D Distribution of combinatorial drugs across different ranges of side effects

Drug-side effect network

Individual drug side effect information was derived from the SIDER [36] and OFFSIDES [35] databases. SIDER data, includes 286,399 drug-side effect relationships involving 1556 drugs and 5868 side effects. OFFSIDES data, collected from adverse event reporting systems, consists of 487,530 drug-side effect pairs involving 1332 drugs and 10,097 side effects. Combining these databases, we compiled 174,977 drug-side effect relationships involving 639 drugs and 10,184 side effects. The frequency distribution of single drugs with side effects generally ranges from 100 to 300, up to over 1000 (Fig. 1). In contrast, single drug side effects are more specific, with 9271 side effects observed in a maximum of only 50 drugs.

These four charts illustrate the distribution of drugs and their combinations across different ranges of side effect counts, providing an in-depth analysis of the relationship between side effect frequency and the number of drug combinations.

The first chart (A) shows the distribution of side effects across different drug count ranges. It is evident that the number of side effects is concentrated in a small range of drugs, particularly in the [1, 9] range, where the number of side effects is significantly higher than in other ranges. This suggests that a small number of drug combinations are responsible for generating many side effects.

The second chart (B) displays the distribution of side effects across different drug combination count ranges. Similar to the first chart, the number of side effects is largely concentrated in the lower drug combination ranges (such as [0, 100] and [100,500]), indicating that these combinations are more likely to produce a higher number of side effects.

The third chart (C) focuses on the distribution of drugs within different ranges of side effect counts. The data shows that most drugs have a small range of side effects (e.g., [100,200] and [200,300]), while drugs with a high number of side effects are relatively rare.

The final chart (D) illustrates the distribution of drug combinations across different ranges of side effects. The results indicate that the majority of drug combinations generate fewer side effects (e.g., [0, 50] and [50, 100]), while the number of drug combinations sharply decreases as the number of side effects increases.

Drug-drug similarity network

Utilizing chemical structural information, we obtained SMILES notation for each drug from DrugBank and transformed it into 2048-dimensional binary vectors using RDKit [37, 38]. Let A and B be vectors in ${Z}_{2}^{2048}$, representing the corresponding binary vectors for drugs a and b. Here, ${Z}_{2}$ refers to the set of binary values $\left\{0, 1\right\}$, where each vector component is either 0 or 1. The similarity of their chemical structures is calculated using the following formula (1). Subsequently, a threshold $c$ is selected. If the chemical structure similarity ${\text{SimChem}}$ is greater than or equal to $c$, an edge is established between the nodes of drug $a$ and drug $b$; otherwise, no edge is formed.

$${\text{SimChem}}\text{=}\frac{2{A}^{T}B}{\left|A\right|+\left|B\right|}$$

(1)

Building on the combinatorial drugs and side effect data from the previous database, we analyzed the similarity distributions among 639 drugs, as summarized in Table 1, finding that most similarities were below 0.3. To more effectively leverage chemical structural similarities in our network construction, we established a similarity threshold of 0.2. Drugs exhibiting a similarity above this threshold were linked, resulting in the formation of 52,557 edges in the network.

Table 1 Statistics of chemical structure similarity of 582 drugs

Drug-protein network

Relationships between drugs and proteins were obtained from the STITCH database [39], focusing solely on experimentally verified interactions between small molecular compounds and target proteins. This dataset includes 8,083,600 drug-protein interactions, involving 8934 proteins and 519,022 compounds. In the resulting network, each protein and each compound are represented as nodes. Edges between nodes represent the verified interactions from the STITCH database.

Protein–protein interaction network

We utilized the human PPI network compiled by Menche et al. [40] and Chatr-Aryamontri et al. [41], combined with additional PPI data from STRING [42] and Rolland et al. [43]. This network includes physically interacting proteins as recorded in human experiments, unweighted and undirected, comprising 19,085 proteins and 719,402 interactions.

Side effect-protein network

While no database directly records side effect-protein interactions, indirect relationships can be established through associated diseases. We extracted side effect-protein relationships from disease-protein data sourced from the DisGeNET database [44], which includes 21,357 relationships between 519 diseases and 7294 proteins. In this network, nodes represent both proteins and side effects. Edges are established based on a documented association between a protein and a disease that manifests a given side effect.

In this study, we focused on predicting side effects related to combinatorial drug therapies. We selected 964 side effects that appeared in at least 500 drug combinations as our experimental dataset. Detailed dataset and network information is in Table 2.

Table 2 Experimental data for predicting combination drug side effects

Utilizing the above network information, we integrated and constructed two heterogeneous information networks as follows:

The first network is the Single Drug–Side Effect Heterogeneous Information Network (SSHIN), depicted in Fig. 2. This network is composed of five sub-networks: the drug-drug chemical structure similarity network, the drug-side effect network, the drug-protein network, the side effect-protein network, and the protein–protein interaction network. It is formally represented as the heterogeneous graph ${G}_{SSHIN}=\left({V}_{1}, {E}_{1}, \left\{D, P, S\right\}, \left\{D-D, D-P, P-P, D-S, S-P\right\}\right)$. Here, D (Drug) denotes the set of drug nodes, P (Protein) denotes the set of protein nodes, and S (Side effect) denotes the set of side effect nodes. To analyze and interpret this complex network, we utilized a metapath-based graph neural network model, which was specifically designed to learn feature vector representations of drugs, thereby laying the groundwork for subsequent predictive tasks.

Fig. 2

Single drug–side effect heterogeneouss information network

The second network we constructed is the Combinatorial Drugs and Side Effect Heterogeneous Information Network (CSHIN), as illustrated in Fig. 3. This network establishes a bipartite graph structure between combinations of drugs and their corresponding side effects. An edge in this network indicates that a specific combination of drugs is associated with a particular side effect. However, it is important to note that the edges between drug combinations and side effects in Fig. 3 do not imply that individual drugs necessarily cause these side effects (i.e., this result is not derived directly from the information in Fig. 2); rather, it is the combination of different individual drugs that may lead to these side effects (independent of Fig. 2's information). It is formally represented as the heterogeneous graph ${G}_{CSHIN} = \left({V}_{2}, {E}_{2}, \left\{C, S\right\}, \left\{CD-S\right\}\right)$, where C (Combinatorial drugs) represents nodes of drug combinations, and S (Side effect) represents nodes of side effects. For analysis and prediction within this network, we utilized advanced graph convolutional neural network models specifically designed for edge prediction tasks, effectively identifying potential adverse effects arising from drug combinations.

Fig. 3

Combinatorial drugs and side effect heterogeneous information network

Computational models

Inspired by the metapath-based models [33, 45], we propose the Metapath-based Aggregated Embedding Model on Single Drug–Side Effect Heterogeneous Information Network (MAEM-SSHIN), to learn embeddings of drug and side effect nodes in SSHIN.

In this paper, within the Combinatorial Drugs and Side Effect Heterogeneous Information Network, we transform the relationships between combinatorial drug pairs and side effects into an edge prediction problem. To address this, we employ a two-layer graph convolutional neural network model, which is specifically designed for predicting new connections. This model utilizes feature vectors of drug nodes and side effect nodes, which are learned from the SSHIN (Single Drug–Side Effect Heterogeneous Information Network). We denote this model as GCN-CSHIN (Graph Convolutional Network Model on Combinatorial Drugs and Side Effect Heterogeneous Information Network), emphasizing its tailored application for intricate network dynamics and interactions.

MAEM-SSHIN model

This section outlines the Metapath-based Aggregated Embedding Model in the Single Drug–Side Effect Heterogeneous Information Network (MAEM-SSHIN). The model includes the following components: metapath selection, node feature transformation, intra-metapath information aggregation, metapath inter-information aggregation, and node feature vector embedding learning.

Metapath selection

In the construction of the single drug and side effect network, our objective is to elucidate the relationships among drugs in terms of chemical structure, target protein, and side effects, as well as the connections between side effects, proteins, and individual drugs. To achieve this, we begin by setting the maximum length of metapaths at 3, allowing for comprehensive exploration within the network. This constraint includes two types of nodes: firstly, the first-order neighboring nodes of drug nodes (or side effect nodes); secondly, the nodes at the second and third orders of adjacency from drug nodes (or side effect nodes), essentially limiting the selection to metapaths that start and end at the same type of node. This strategic selection results in 11 types of metapaths for drug nodes $\left(DD,DP,DS,DPD,DSD,DDPD,DDSD, DPDD,DSDD,DPPD,DPSD\right)$, and 8 types for side effect nodes $\left(SD,SP,SDS,SPS,SDDS,SPPS, SPDS,SDPS\right)$, effectively mapping the intricate interconnections within the network. Conversely, longer metapaths (e.g., length of 4) tend to introduce excessive noise, potentially diminishing the accuracy of predictions. Our approach to restrict the path length aims to balance computational complexity and the need for predictive accuracy, ensuring the model remains efficient without being overly complex.

Node feature transformation

In heterogeneous graphs, which contain nodes with diverse features, different node types may have feature vectors—numerical representations capturing the attributes or characteristics of nodes—of varying dimensions, meaning the number of elements in each vector can differ. Even when these vectors share the same dimensionality, they may still originate from distinct feature spaces, which are mathematical spaces defined by the possible values and relationships of the features. To effectively process entities like drugs, target proteins, and side effects within a unified analytical framework, it is essential to map the initial feature vectors of various node types into a common latent feature space. This mapping is accomplished through type-specific linear transformations, ensuring that disparate data types are harmonized for subsequent analysis and modeling [32].

In our model, the initial feature vectors of all nodes are encoded as one-hot vectors based on their type, and subsequently projected into the same feature space using a node type-specific transformation matrix ${W}_{l}$. The transformation is represented by the equation:

$$h^\prime_{v} = W_{l} \cdot x_{v} ,v \in V$$

(2)

where $h^{\prime}_{v}$ and ${x}_{v}$ are the post-transformation and original feature vectors of node $v$, respectively, with $l\in \left\{D,P,S\right\}$ denoting the node's type.

Intra-metapath information aggregation

This section is dedicated to elucidating the process of information aggregation. It uses metapaths originating from drug nodes as a primary example to demonstrate how this process unfolds.

Given a metapath related to a drug node ${P}_{m}$, where $m$ belongs to the set $\{\text{1,2},\dots ,11\}$, and considering ${p}_{m}\left(v,u\right),u\in {N}_{{P}_{m}}^{v}$ as the metapath instance for node $v$, the feature information of node $v$ is transmitted along the metapath instance ${p}_{m}\left(v,u\right)$ as follows:

$${h}_{{p}_{m}\left(v,u\right)}={f}_{\uptheta }\left({h}^{\prime}_{t},\right), \forall t\in {V}_{{p}_{m}}\left(v,u\right)$$

(3)

In this equation, ${h}_{{p}_{m}}\left(v,u\right)\in {\mathbb{R}}^{{d}^{\prime}}$ is a vector of dimension ${d}^{\prime}$, where ${V}_{{p}_{m}}\left(v,u\right)$ signifies all nodes present in the metapath instance ${p}_{m}\left(v,u\right)$, and ${f}_{\uptheta }\left(\cdot \right)$ represents an information transfer function connected with the parameter $\theta$. The metapath instances are encoded into feature vectors, following which a graph attention mechanism is utilized to assign weights to diverse metapath instances for the target protein node $v$. The information from all metapath instances is aggregated through a weighted summation approach.

$$\begin{array}{*{20}c} {e_{{p_{m} \left( {v,u} \right)}} = {\text{LeakyReLU}}\left( {a_{{P_{m} }}^{T} \cdot \left[ {h^{\prime}_{v} \parallel h_{{p_{m} \left( {v,u} \right)}} } \right]} \right)} \\ {\alpha_{{p_{m} \left( {v,u} \right)}} = \frac{{\exp \left( {e_{{p_{m} \left( {v,u} \right)}} } \right)}}{{\mathop \sum \nolimits_{{s \in {\mathcal{N}}_{v}^{{P_{m} }} }} \exp \left( {e_{{p_{m} \left( {v,s} \right)}} } \right)}}} \\ {h_{v}^{{P_{m} }} = {\upsigma }\left( {\mathop \sum \limits_{{u \in {\mathcal{N}}_{v}^{{P_{m} }} }} \alpha_{{p_{m} \left( {v,u} \right)}} \cdot h_{{p_{m} \left( {v,u} \right)}} } \right)} \\ \end{array}$$

(4)

Here, ${\varvec{a}}_{{Pm}} \in \mathbb{R}^{{2d^{\prime}}}$ represents the attention parameter vector for metapath ${P}_{m}$, with '$\parallel$' denoting the concatenation of two vectors. The term ${e}_{{p}_{m}\left(v,s\right)}$ signifies the importance of the metapath instance ${p}_{m}\left(v,u\right)$ for node $v$. This importance is then normalized using the softmax function to obtain the weight $\alpha_{{p_{m} \left( {v,u} \right)}}$ for the metapath instance ${p}_{m}\left(v,u\right)$. Consequently, all metapath instances are aggregated to derive the feature vector ${h}_{v}^{{P}_{m}}$ for node $v$ on metapath ${P}_{m}$. The symbol $\sigma \left(\cdot \right)$ denotes the activation function. To enhance the model's expressive power and stabilize the learning process, we employ a multi-head attention strategy, executing $K$ independent attention mechanisms and subsequently concatenating their outputs as the final result.

$${h}_{v}^{{P}_{m}}={\parallel }_{k=1}^{K}\sigma \left(\sum_{u\in {\mathcal{N}}_{v}^{{P}_{m}}}{\left[{\alpha }_{{p}_{m}\left(v,u\right)}\right]}_{k}\cdot {h}_{{p}_{m}\left(v,u\right)} \right)$$

(5)

Here, ${\left[{\alpha }_{{p}_{m}\left(v,u\right)}\right]}_{k}$ denotes the weight of the metapath instance ${p}_{m}\left(v,u\right)$ under the ${k}^{th}$ attention mechanism.

The methodology for metapaths starting from side effects is analogous and will not be elaborated on further.

Metapath inter-information aggregation

After aggregating node and edge data within each metapath, we further aggregate the features across different metapaths to achieve a comprehensive feature expression for all nodes across all metapaths. Here is an example for drug nodes, which are representations of pharmaceutical compounds in the network, each characterized by their vector feature expressions. While this example focuses on drug nodes, the same aggregation process is applied to side effect nodes, which also involve assembling vector feature expressions to capture the comprehensive impact of potential adverse effects associated with the drugs.

For a node $v$ in a graph $G$, after the steps previously described, we obtain feature vector expressions for a drug node $v$ across metapath ${P}_{m}$, represented as ${h}_{v}^{{P}_{m}}$, where $m \in \left\{1, 2, \cdots , M\right\}$.

Next, we assign a different weight to each metapath to aggregate all metapath information comprehensively. Specifically, we calculate the average of the feature vectors for drug node $v$, represented in ${V}_{D}$, across each metapath ${P}_{m}$, where $m \in \left\{1, 2, \cdots , M\right\}$. This average reflects the consolidated information from all instances of node $v$ on metapath ${P}_{m}$, capturing the integrated effects of various relational paths in the network.

$${\varvec{s}}_{{Pm}} = \frac{1}{{\left| {V_{D} } \right|}}\sum\limits_{{v \in V_{D} }} {\tan \;h} \left( {M_{D} \cdot h_{v}^{{Pm}} + b_{D} } \right)$$

(6)

where, $M_{D} \in {\mathbb{R}}^{{d_{D} \times d^{\prime}}}$ and $b_{D} \in {\mathbb{R}}^{{d_{D} }}$ are the parameters to be learned.

Then, an attention mechanism is employed to aggregate the node vectors of node v across all metapaths.

$$\begin{gathered} e_{{Pm}} = {\varvec{q}}_{D}^{T} \cdot {\varvec{s}}_{{Pm}} \hfill \\ \beta _{{Pm}} = \frac{{\exp \left( {e_{{Pm}} } \right)}}{{\sum\nolimits_{{i = 1}}^{M} ~ \exp \left( {e_{{Pi}} } \right)}} \hfill \\ h_{v} = \sum\limits_{{i = 1}}^{M} {\beta _{{Pm}} } \cdot h_{v}^{{Pm}} \hfill \\ \end{gathered}$$

(7)

where, ${\text{q}}_{\text{D}} \in {\mathbb{R}}^{{\text{d}}_{\text{D}}}$ is the attention parameter for drug nodes to be learned, and ${h}_{v}$ is the final aggregated embedding vector for node $v$.

Node feature vector embedding learning

To predict drug-drug and drug-side effect linkages, we employ the Binary Cross Entropy loss (BCE) as the loss function to assess the discrepancy between the predicted values and the actual outcomes. The formula is:

$$\mathcal{L}=-\frac{1}{B}{\sum }_{b=1}^{B}\left[{y}_{b}\text{log}\left({y}^{\prime}_{b}\right)+\left(1-{y}_{b}\right)\text{log}\left(1-{y}^{\prime}_{b}\right)\right]$$

(8)

Here, $B$ represents the number of training data per batch, ${y}_{b}$ denotes the actual value of the ${b}^{th}$ edge prediction sample, where the presence of an edge is indicated by 1, and absence by 0. ${y}^{\prime}_{b}=F\left({h}_{u},{h}_{v}\right)$ signifies the predicted probability of an edge existing between node pair $\left(u,v\right)$ in the model, where $F\left(\cdot \right)$ is a function calculating the similarity between two feature vectors. Ultimately, when the loss function falls below a set threshold, the resulting ${h}_{v}$ is the final embedding vector for node $v$.

With this, the complete process of the MAEM-SSHIN model is thoroughly described. The specific flowchart is presented in Fig. 4.

Fig. 4

Flowchart of MAEM-SSHIN

GCN-CSHIN model

The goal of GCN-CSHIN model is to predict whether an edge exists between combinatorial drugs and side effects in a heterogeneous graph. An edge signifies the emergence of a particular side effect due to the combinatorial drugs. We employ a Graph Convolutional Network (GCN) model to aggregate structural information of combinatorial drugs and side effect nodes for this edge prediction task.

Leveraging the MAEM-SSHIN model, we obtain feature vector representations for each drug node. By combining the feature vectors of two drug nodes, we derive the initial feature ${f}_{0}^{c}={h}_{{v}_{i}}+{h}_{{v}_{j}}$ for the combinatorial drugs node in the heterogeneous graph ${G}_{CSHIN}$, where $c\in {V}_{C}$ and ${v}_{i},{v}_{j}\in {V}_{D}$. This additive approach is chosen to mitigate the influence of the drug combination's sequence on prediction outcomes. The initial feature of side effect nodes is as obtained from the MAEM-SSHIN model, i.e., ${f}_{0}^{s}={h}_{{v}_{s}}$, where ${v}_{s}\in {V}_{S}$.

In the GCN model, the node information propagation formula for each layer (illustrated here for combinatorial drugs nodes) is:

$${f}_{k+1}^{c}=\phi \left({\sum }_{u\in N\left(c\right)}{\gamma }_{cu}{W}_{k}{f}_{k}^{u}+{f}_{k}^{c}\right)$$

(9)

Here, ${f}_{k}^{c}$ and ${f}_{k+1}^{c}$ denote the feature representations of node $c$ at the ${k}^{th}$ and $(k+{1)}^{th}$ layer, respectively, $N\left(c\right)$ is the set of neighboring nodes of $c$, ${\gamma }_{cu}=1/\left|N\left(c\right)\right|$, ${W}_{k}$ is the parameter matrix, and $\phi \left(\cdot \right)$ is the activation function. We use a two-layer GCN model to ultimately derive the feature representations ${f}_{c}$ and ${f}_{s}$ for combinatorial drugs and side effect nodes in the aggregated CSHIN graph, where $c\in {V}_{C} and {v}_{s}\in {V}_{S}$. The node pair combines as $f={f}_{c}{||f}_{s}$ and inputs into a Multilayer Perceptron (MLP) model to yield the prediction result:

$$y=\sigma \left({W}_{f}f\right)$$

(10)

Here, $W_{f}$ represents the trainable parameters, and $\sigma \left( \cdot \right)$ is the activation function. The Binary Cross-Entropy loss is used to evaluate the difference between predicted and actual results:

$$\mathcal{L}=-\frac{1}{B}{\sum }_{b=1}^{B}\left[{y}_{b}\text{log}\left({y}_{b}{\prime}\right)+\left(1-{y}_{b}\right)\text{log}\left(1-{y}_{b}{\prime}\right)\right]$$

(11)

In this formula, $B$ is the number of training data per batch, ${y}_{b}=1$ indicates the presence of an edge between the ${b}^{th}$ pair of combinatorial drugs and side effect nodes, and ${y}_{b}{\prime}$ denotes the model's predicted probability of an edge existing between them. For the flowchart of the GCN-CSHIN process, refer to Fig. 5.

Fig. 5

Flowchart of GCN-CSHIN

Results

Evaluation metrics

To assess the performance of our model, we utilized three commonly employed metrics in the field of multi-drug side effect prediction:

• AUROC (Area under the curve ROC) This is the area under the Receiver Operating Characteristic (ROC) curve. The ROC curve's vertical axis represents the True Positive Rate (TPR), defined as $TPR=Recall =\frac{TP}{TP+FN}$, where the horizontal axis is the False Positive Rate (FPR), defined as $FPR=\frac{FP}{FP+TN}$.

• AUPRC (area under the curve PRC) This is the area under the Precision-Recall Curve (PRC). The PRC curve's vertical axis is Precision, defined as $Precision=\frac{TP}{TP+TN}$, and the horizontal axis is Recall, defined as $Recall =\frac{TP}{TP+FN}$.

Parameter settings

Parameter settings were determined using a grid search approach for evaluating candidate values, while ensuring the consistency of the training, testing, and validation nodes. The key settings for the MAEM-SSHIN model include employing a dual attention mechanism and a two-layer neural network architecture, with each layer having a dimension of 128. Negative sample selection is based on random walk sampling, and the loss function is set to binary cross-entropy. Similarly, the GCN-CSHIN model utilizes a two-layer graph neural network, where each layer is also 128-dimensional.

We conduct a fivefold cross-validation, ensuring that the training, validation, and test nodes are consistent across both MAEM-SSHIN and GCN-CSHIN. During each fold of cross-validation, we keep the node division fixed to ensure consistency of the train, validation, and test nodes. For each fold, 80% of the positive and negative samples are randomly selected as the training set, and the remaining samples are used as the test set. All models are trained and evaluated on the same partitioned datasets. The learning rate is a dynamic learning rate, where the initial rate is 1e-3 and decreases by 0.05 every 100 steps. The batch sizes for training and testing differ, with the training batch size being 1024 and the test batch size being 4096. The smaller batch size during training allows for more frequent updates and better memory management, while the larger batch size during testing enables faster and more efficient evaluation of the model's performance.

Case study

To validate whether our model is capable of predicting side effects for drug combinations that currently lack recorded side effect information, we designed a case study specifically to test the model's predictive capabilities.

We conducted experiments using a dataset containing 915,583 edges. Given the large size of the dataset, we adopted a batch processing approach, predicting 4098 positive edges per batch. Through multiple experiments, the model generated over 200 top-10 prediction results. Among all the predicted results, the side effect of "anaemia" was predicted the most frequently, appearing in 25 combinations. We then conducted an in-depth study on all the drug combinations related to anaemia (Table 3).

Table 3 The predicted drug combinations associated with anaemia side effect and their corresponding literature validation

During the testing phase, all drug combinations were treated as having no direct side effects, allowing the model to predict potential side effects from scratch. The 25 combinations related to "anaemia" were the result of predictions made on the entire test set, not specifically pre-selected combinations without direct side effects.

Although the side effects of these 25 drug combinations are not fully documented in the literature, upon validation, we found that 9 of these combinations have indeed been reported to be associated with anaemia. For the remaining combinations, literature has documented individual drugs that may cause similar side effects. For example, the study [46] reported that some patients treated with Taxol for recurrent anaplastic astrocytomas experienced various toxicities, including anaemia. Similarly, the study by reference [47] explored haemolytic anaemia induced by Taxol. These literature validations support the model's predictions and further demonstrate its potential in identifying previously unreported side effects of drug combinations.

Model comparison and analysis

In our model comparison, we selected commonly used multi-relational prediction models, the classical Decagon model, and its improved versions. The specific models are as follows:

• Concatenated Drug Features [10]: This approach constructs feature vectors for each drug based on PCA representations of drug-target interaction matrices and individual drug side effects. Combinatorial drugs features are obtained by concatenating corresponding single-drug feature vectors, followed by using a gradient boosting tree classifier to predict the side effects of combinatorial drugs.

• DeepWalk Neural Embeddings [48, 49]: This method learns d-dimensional feature representations of drug nodes through a biased random walk process. Combinatorial drug features are derived by concatenating individual drug vectors. A logistic regression classifier is applied to predict the occurrence of specific side effects for each drug combination.

• DEDICOM Tensor Decomposition [50]: A tensor decomposition method suitable for sparse data. Given a drug-drug feature vector matrix ${X}_{i}$, it is decomposed as: ${X}_{r}= A{U}_{r}T{U}_{r}{A}^{T}$. The association prediction between drugs i and j for r is ${a}_{i}{U}_{r}T{U}_{r}{a}_{j}$.

• RESCAL Tensor Decomposition [51]: A tensor decomposition method that considers multi-relational structures. Given a drug-drug feature vector matrix ${X}_{i}$, based on side effect r, it is decomposed as: ${X}_{r}= A{T}_{r}{A}^{T}$, r $\in${1, 2, …, 963}. The association prediction between drugs i and j under side effect r is: ${a}_{i}T{a}_{j}$.

• Decagon [10]: Constructs a multimodal network of protein-drug multi-drug side effects, where in the drug-drug interaction network, different edge types are established based on different side effects. It develops a combination model of a new graph convolutional neural network and a tensor decomposition decoder for multi-relational link prediction in multimodal networks.

• TIP [11]: A Tripartite Information Propagation (TIP) model using the multimodal network of protein-drug multi-drug side effects from the Decagon method. It builds a message passing neural network (MPNN) framework through GCN (P-P graph embedding model), Decagon (graph-to-graph information propagation model), and R-GCN (D-D graph embedding model), propagating pharmacological information from P-P to D-D graphs via P-D graphs, obtaining drug feature representations, and using a neural network as the decoder.

• DistMult [14]: Constructs a drug-protein knowledge graph, using the classic knowledge graph embedding method—DistMult—and an end-to-end knowledge graph representation learning framework—KBlrn—to predict the side effects of combinatorial drugs.

• SC-DDITS [13]: A novel heterogeneous signature network model that combines more drug features, dividing multi-modal relationships into positive or negative categories, and builds a graph neural network model based on extended structural balance theory in social networks to learn drug feature representations.

In our study, we ensured that the data partitioning process was consistent across all phases. The test data were isolated at the beginning and never used during training or tuning. This partitioning was strictly followed from preprocessing to model evaluation, ensuring that no test data were inadvertently exposed during model development.

The test dataset remained unseen throughout the framework and was solely used for the final evaluation to assess the model's generalization performance on unseen data. This ensures that the results reported are unbiased and adhere to best practices in machine learning.

We utilized fivefold cross-validation to ensure robust evaluation across all models. Each fold involved using one subset of the data as the test set, with the remaining used for training, ensuring comprehensive testing across all data points. We also conducted ablation studies on the GCN-CSHIN model, replacing the original GCN + MLP with a simple MLP or other graph neural network models (such as GraphSAGE) + MLP. The results showed that neither using MLP directly nor replacing GCN with other graph models performed as well as the current approach. These ablation study results further emphasize the critical role of incorporating graph structure information in enhancing overall model performance. Table 4 shows our model outperforms other models on three metrics: AUROC and AUPRC. Unlike Decagon and other models, ours integrates drug chemical structure and side effect data, simplifying multi-relational prediction within CSHIN. This allows extending to higher-order drug combination side effect prediction and incorporating more drug feature information in learning feature vectors via the MAEM-SSHIN model (Table 5).

Table 4 Comparison of AUROC and AUPRC between our proposed model and the classical models

Table 5 Comparison of AUROC and AUPRC under different SSHIN networks

To better illustrate the model's performance and the importance of drug features, we modified the drug-side effect heterogeneous feature network, omitting some relational features. We denote this adapted network as MAEM-SSHIN_ddpps, which includes drug-drug chemical structure, drug-protein, drug-side effect, protein–protein interaction, and side effect-protein information. Similarly, we define MAEM-SSHIN_dds(only containing drug-drug chemical structure and drug-side effect information), MAEM-SSHIN_dpp(only containing drug-protein and protein–protein interaction information), MAEM-SSHIN_dpps(with drug-protein, drug-side effect, protein–protein interaction, and side effect-protein information), and MAEM-SSHIN_ddpp(with drug-drug chemical structure, drug-protein, drug-side effect, and protein–protein interaction information). MAEM-SSHIN_ddpps performed best, while MAEM-SSHIN_dds performed worst, highlighting the significance of drug-drug chemical structure and drug-side effect information in predicting combinatorial drug side effects.

Discussion and conclusions

In this study, we have innovatively proposed a predictive framework for combinatorial drugs side effects based on heterogeneous graph neural networks. The framework comprises two core components: a metapath-based heterogeneous graph embedding model (MAEM-SSHIN) for single drug side effects, and a graph convolutional network prediction model (GCN-CSHIN) for combinatorial drugs side effects. The integration of these models not only enriches methodology for predicting multiple drug side effects, but also opens new avenues for practical clinical applications. The MAEM-SSHIN model extends beyond traditional drug-protein heterogeneous network analyses by incorporating drug chemical structure and individual drug side effect data, resulting in a more comprehensive single drug side effect heterogeneous network. The GCN-CSHIN model presents a novel approach for combinatorial drugs side effect predictions, transforming the problem of predicting multiple side effect relations for drug pairs into a network relation prediction problem between combinatorial drugs and side effects, utilizing the potent learning capabilities of GCN for more accurate and efficient predictions.

Our experimental results demonstrate that the combined MAEM-SSHIN and GCN-CSHIN models outperform existing models in predicting combinatorial drugs side effects. This validates the effectiveness of the proposed methods and indicates potential application of the framework in addressing real-world issues. Future research directions include multidimensional data integration, model interpretability, drug repositioning, personalized medicine, and design of multi-drug combination therapies. Through these studies, we aim to provide a solid scientific foundation for safety assessment and efficacy prediction of drug combinations, ultimately benefiting a wide range of patients.