Predicting Biomedical Interactions with Higher-Order Graph Convolutional Networks

Biomedical interaction networks have incredible potential to be useful in the prediction of biologically meaningful interactions, identification of network biomarkers of disease, and the discovery of putative drug targets. Recently, graph neural networks have been proposed to effectively learn representations for biomedical entities and achieved state-of-the-art results in biomedical interaction prediction. These methods only consider information from immediate neighbors but cannot learn a general mixing of features from neighbors at various distances. In this paper, we present a higher-order graph convolutional network (HOGCN) to aggregate information from the higher-order neighborhood for biomedical interaction prediction. Specifically, HOGCN collects feature representations of neighbors at various distances and learns their linear mixing to obtain informative representations of biomedical entities. Experiments on four interaction networks, including protein-protein, drug-drug, drug-target, and gene-disease interactions, show that HOGCN achieves more accurate and calibrated predictions. HOGCN performs well on noisy, sparse interaction networks when feature representations of neighbors at various distances are considered. Moreover, a set of novel interaction predictions are validated by literature-based case studies.


INTRODUCTION
A Biological system is a complex network of various molecular entities such as genes, proteins, and other biological molecules linked together by the interactions between these entities. The complex interplay between various molecular entities can be represented as interaction networks with molecular entities as nodes and their interactions as edges. Such a representation of a biological system as a network provides a conceptual and intuitive framework to investigate and understand direct or indirect interactions between different molecular entities in a biological system. Study of such networks lead to systemlevel understanding of biology [1] and discovery of novel interactions including protein-protein interactions (PPIs) [2], drug-drug interactions (DDIs) [3], drug-target interactions (DTIs) [4] and gene-disease associations (GDIs) [5].
Recently, the generalization of deep learning to the network-structured data [6] has shown great promise across various domains such as social networks [7], recommendation systems [8], chemistry [9], citation networks [10]. These approaches are under the umbrella of graph convolutional networks (GCNs). GCNs repeatedly aggregate feature representations of immediate neighbors to learn the informative representation of the nodes for link prediction. Although GCN based methods show great success in biomedical interaction prediction [3], [11], the issue with such methods is that they only consider information from immediate neighbors. SkipGNN [12] leverages skip graph to aggregate feature representations from direct and • K. KC second-order neighbors and demonstrated improvements over GCN methods in biomedical interaction prediction. However, SkipGNN cannot be applied to aggregate information from higher-order neighbors and thus fail to capture information that resides farther away from a particular interaction [13].
To address the challenge, we propose an end-to-end deep graph representation learning framework named higher-order graph convolutional networks (HOGCN) for predicting interactions between pairs of biomedical entities. HOGCN learns a representation for every biomedical entity using an interaction network structure G and/or features X. In particular, we define a higher-order graph convolution (HOGC) layer where the feature representations from korder neighbors are considered to obtain the representation of biomedical entities. The layer can thus learn to mix feature representations of neighbors at various distances for interaction prediction. Furthermore, we define a bilinear decoder to reconstruct the edges in the input interaction network G by relying on feature representations produced by HOGC layers. The encoder-decoder approach makes HOGCN an end-to-end trainable model for interaction prediction.
We compare HOGCN's performance with that of stateof-the-art network similarity-based methods [14], network embedding methods [15], [16], and graph convolutionbased methods [10], [12], [17] for biomedical link prediction. We experiment with various interaction datasets and show that our method makes accurate and calibrated predictions. HOGCN outperforms alternative methods based on network embedding by up to 30%. Furthermore, HOGCN outperforms graph convolution-based methods by up to 6%, alluding to the benefits of aggregating information from higher-order neighbors.
We perform a case study on the DDI network and arXiv:2010.08516v1 [cs.
LG] 16 Oct 2020 observe that aggregating information from higher-order neighborhood allows HOGCN to learn meaningful representation for drugs. Moreover, literature-based case studies illustrate that the novel predictions are supported by evidence, suggesting that HOGCN can identify interactions that are highly likely to be a true positive. In summary, our study demonstrates the ability of HOGCN to identify potential interactions between biomedical entities and opens up the opportunities to use the biological and physicochemical properties of biomedical entities for a follow-up analysis of these interactions.

RELATED WORKS
With the increasing coverage of the interactome, various network-based approaches have been proposed to exploit already available interactions to predict missing interactions [14], [18]- [20]. These methods can be broadly classified into (1) network similarity-based methods (2) network embedding methods (3) graph convolution-based methods. We next summarize these categories of methods for biomedical interaction prediction.
Given a network of known interactions, various similarity metrics are used to measure the similarity between the biomedical entities [18] with an assumption that higher similarity indicates interaction. Triadic closure principle (TCP) has been explored in biomedical interaction prediction with the hypothesis that biomedical entities with common interaction partners are likely to interact with each other [19]. TCP relies on a common neighbor algorithm to count the number of shared neighbors between the nodes and is quantified by A 2 where A is the adjacency matrix. Recently, L3 heuristic [14] shows the common neighbor hypothesis fails for most protein pairs in PPI prediction and proposes to consider nodes that are similar to the neighbors of the nodes and can be quantified by A 3 . This indicates that higher-order neighbors are important for interaction prediction.
Next, network embedding methods embed the existing networks to low-dimensional space that preserves the structural proximity such that the nodes in the original network can be represented as low-dimensional vectors. Deepwalk [15] is a popular approach that generates the truncated random walks in the network and defines a neighborhood for each node as a set of nodes within a window of size k in each random walk. Similarly, node2vec performs a biased random walk by balancing the breadth-first and depth-first search in the network. The random walks generated by these methods can be considered as a combination of nodes from various order of neighborhoods such as 1-hop to khop neighborhood. In other words, DeepWalk and node2vec learn the embeddings for the nodes in the network from the combination of A 1 , A 2 , A 3 , . . . A k where A i is the i th power of the adjacency matrix. These embeddings can then be fed into a classifier to predict the interaction probability. These methods are only limited to the structure of the biomedical networks and cannot incorporate additional information about the biomedical entities. Also, they cannot learn the feature difference between nodes at various distances.
Furthermore, graph convolution-based methods use a message-passing mechanism to receive and aggregate information from neighbors to generate representations for the nodes in the network. Graph convolutional networks (GCNs) [17] and variational graph convolutional autoencoder (VGAE) [10] aggregate feature representation from immediate neighbors to learn the representation of biomedical entities in an end-to-end manner using link prediction objective. These methods are only limited to the average pooling of the neighborhood features [13]. SkipGNN [12] therefore proposes to use skip similarity between the biomedical entities to aggregate information from secondorder neighbors. However, these methods cannot aggregate feature representations from higher-order neighbors and also cannot learn feature differences between neighbors at various distances.

PRELIMINARIES
A biomedical network is defined as G = (V, E, X) where V denotes the set of nodes representing biomedical entities (e.g. proteins, genes, drugs, diseases) and |V| denotes the number of nodes. E ⊆ (V × V) denotes a set of interactions between biomedical entities. X ∈ R |V|×F is the features of biomedical entities where F is the dimension of features.
Let A denote the adjacency matrix of G, where A ij indicates an edge between nodes v i and v j . We assume the case of binary adjacency matrix A ij ∈ {0, 1} n×n where A ij represents the existence of edge between the nodes v i and v j , indicating the presence of the experimental evidence for their interaction (i.e. A ij = 1) or the absence of the experimental evidence for their interaction (i.e. A ij = 0). Note that the same notation of adjacency matrix can be used to represent weighted graphs such that A ij = [0, 1]. Table 1 shows the notations and their definitions used in the paper.
Ground-truth interaction between nodes i and j p ij ∈ [0, 1] Probability of interaction between nodes i and j Z ∈ R |V|×d * Final node embeddings W

Message Passing
Given a biomedical network G, message passing algorithms learn the representation of biomedical entities in the network by aggregating information from immediate neighbors [9]. Additional information about biomedical entities Given a biomedical interaction network G with initial features X for biomedical entities, the encoder mixes the feature representation of neighbors at various distances and learn final representation Z. The decoder takes the representation z i and z j of nodes v i and v j to learn the representation e ij for the edge (denoted by ?) and predict probability p ij of its existence.
can be used to initialize the feature matrix X. These algorithms involve the message passing step in which each biomedical entity sends its current representation to, and aggregates incoming messages from its immediate neighbors. Representation for each biomedical entity can be obtained after L steps of message passing and feature aggregation. However, such message passing operation is limited to average pooling of features from immediate neighbors and thus is unable to learn feature differences among neighbors at various distances [13]. Neighborhood nodes at various distances provide network structure information at different granularities [21]- [25]. Taking k-hop neighborhoods into consideration, we aim at aggregating information from various distances at every message passing step. Different powers of adjacency matrices such as A 1 , A 2 , A 3 , . . . , A k provide information about the network structure at different scales. Higherorder message passing operations can therefore learn to mix their representations using various powers of the adjacency matrix at each message passing step.

Graph Convolutional Networks (GCNs)
Graph convolutional networks (GCNs) are the generalization of convolution operation from regular grids such as images or texts to graph structured data [6], [26]. The key idea of GCNs is to learn the function to generate the node's representation by repeatedly aggregating information from immediate neighbors. The graph convolutional layer is defined as: where H (l−1) and H (l) are the input and output activations, W (l) is a trainable weight matrix of the layer l, σ is the element-wise activation, andÂ is a symmetrically normalized adjacency matrix with self-connectionŝ A GCN model with L layers is then defined as: and H (L) can be used to predict the probability of interactions between biomedical entities.

HIGHER-ORDER GRAPH CONVOLUTION NET-WORK (HOGCN)
In this work, we develop a higher-order graph convolutional network (HOGCN) that takes an interaction network G as input and reconstruct the edges in the interaction network ( Fig. 1). HOGCN has two main components: • Encoder: a higher-order graph convolution encoder that operates on an interaction graph G and produces representations for biomedical entities by aggregating features from the neighborhood at various distances and • Decoder: a bilinear decoder that relies on these representations to reconstruct the interactions in G.

Higher-Order Graph Encoder
We first describe the higher-order graph encoder, that operates on an undirected interaction graph G = (V, E, X) and learns the representations for biomedical entities.
We develop an encoder with higher-order Graph Convolution (HOGC) layer to mix feature representations from neighbors at k-distances. Specifically, HOGC layer considers the neighborhood information at different granularities and is defined as: where P is a set of integer adjacency powers, A j denotes the adjacency matrix A multiplied j times, and denotes column-wise concatenation [13]. Graph convolutional network [17] only considers the 1 st power of adjacency matrix and can be exactly recovered by setting P = {1} in Equation (2). Similarly, SkipGNN [12] considers direct and skip similarity and can be recovered by setting P = {1, 2} in Equation (2). Fig. 2 shows a HOGC layer with P = {0, 1, . . . , k} where k is maximum order of neighborhood considered in

The feature representation H (l) is a linear combination of the neighbors
j represents feature representation of neighbors at j distances for layer l.
each HOGC layer. If k = 0, HOGC layer only considers the features of the biomedical entities and can capture the feature similarity between various biological entities. This is equivalent to a fully connected network with features of biomedical entities as input. For the HOGC layer, A 0 is the identity matrix I |V| where |V| is the number of nodes in the network. This allows the HOGC layer to learn the transformation of node features separately and mix it with feature representations from neighbors.
The maximum order of neighborhood k and the number of trainable weight matrices |P |, one per each adjacency power, can vary across layers. However, we set the same k for neighborhood aggregation and the same dimension d for all the weight matrices across all layers.
Neighborhood features from different adjacency powers j ∈ {0, 1, . . . , k} at layer (l − 1) are column-wise concatenated to obtain feature representation H (l−1) . As shown in (·) at layer l can learn the arbitrary linear combination of the concatenated features to obtain H (l) . Specifically, the layer can assign different coefficients to different columns in the concatenated matrix. For instance, the layer can assign positive coefficients to the columns produced by certain power of A and assign negative coefficients to other powers. This allows the model to learn feature differences among neighbors at various distances. We apply L HOGC layers to learn the latent representation Z ∈ R |V|×d * for biomedical entities in the network, where d * = d×|P | and d is the dimension of node's representation for each adjacency power.

Interaction Decoder
We introduced the encoder based on HOGC layers that learns feature representation Z for biomedical entities by mixing neighborhood information at multiple distances. Next, we discuss the decoder that reconstructs the interactions in G based on the representation Z.
We adopt a bilinear layer to fuse the representation of biomedical entities v i and v j and learn the edge representation e ij . More precisely, we define a simple bilinear layer that takes the representation z i and z j as input: where W b ∈ R d×d * ×d * represents the learnable fusion matrix, e ij is the representation of edge e ij between nodes v i and v j , b denotes the bias of the bilinear layer. ELU is nonlinearity.
The edge representation e ij is then fed into 2-layered fully connected neural network to predict probability p ij for edge e ij : where FC 1 (e ij ) = W 1 ·e ij +b 1 denotes fully connected layer with weight W 1 and bias b 1 , p ij represents the probability that biomedical entities v i and v j interact. So far, we have discussed the encoder and decoder of our proposed approach. Next, we describe the training procedure of our proposed HOGCN model. In particular, we explain how to optimize the trainable neural network weights in an end-to-end manner.

Training HOGCN
During HOGCN training, we employ binary cross entropy loss to optimize the model parameters and encourage the model to assign higher probability to observed interactions (v i , v j ) than to randomly selected non-interactions. p ij is the predicted interaction probability between v i and v j and A ij denotes the ground-truth interaction label between these nodes. The final loss function considering all interactions is We follow an end-to-end approach to jointly optimize over all trainable parameters and backpropagate the gradients through encoder and decoder of HOGCN.

Algorithm
HOGCN leverages biomedical network structure A along with additional information about biomedical entities as the initial feature representation X. In this paper, we initialize the initial features X to be one-hot encoding i.e. I |V| .
The feature matrix X can be initialized with properties of biomedical entities or pre-trained embeddings from other network-based approaches such as DeepWalk, node2vec. Given an adjacency matrix A and the initial node representations X, higher-order neighborhood indicated by the higher power of the adjacency matrix is iteratively computed that makes the model more efficient. By adopting right-to-left multiplication, for instance, we can calculatê j learned for the neighborhood at j distances are concatenated to obtain the representation H (l) as shown in Fig. 2 (Line 11 in Algorithm 1). After passing through Algorithm 1 Training of HOGCN for biomedical interaction prediction 1: Inputs:Â, X, k 2: H (0) = X 3: for t = 1 to T do 4: Sample mini-batch of training edges and their interaction labels 5: for l = 1 to L do 6: B := H (l−1)

7:
for j = 1 to k do p ij := Interaction Decoder(Z) 15: Compute loss in (6) 16: Update model parameters via gradient descent 17: end for L HOGC layers, we obtain the final representation Z for biomedical entities. With the final representations Z and the mini-batch of training edges, we retrieve the embeddings for the nodes in training edges and feed them into the interaction decoder to compute their interaction probabilities.
The parameters of HOGCN are optimized with a binary cross-entropy loss (Equation (6)) in an end-to-end manner. Given two biomedical entities v i and v j , the trained model can predict the probability of their interactions.

EXPERIMENTAL DESIGN
We view the problem of biomedical interaction prediction as solving a link prediction task on an interaction network. We consider various interaction datasets and compare our proposed method with the state-of-the-art methods.

Datasets
We conduct interaction prediction experiments on four publicly-available biomedical network datasets: • BioSNAP-DTI [27]: DTI network contains 15,139 drugtarget interactions between 5,018 drugs and 2,325 proteins. • BioSNAP-DDI [27]: DDI network contains 48,514 drugdrug interactions between 1,514 drugs extracted from drug labels and biomedical literature. • HuRI-PPI [2]: HI-III human PPI network contains 5,604 proteins and 23,322 interactions generated by multple orthogonal high-throughput yeast two-hybrid screens. • DisGeNET-GDI [28]: GDI network consists of 81,746 interactions between 9,413 genes and 10,370 diseases curated from GWAS studies, animal models and scientific literature. Table 2 provides summary of datasets used in our experiments. We provide the number of interactions used for training, validation, and testing for each interaction datasets. Also, the table includes the average number of interactions for each biomedical entity which can be computed as 2|E| |V| .

Baselines
We compare our proposed model with the following network-based baselines for interaction prediction: • network similarity-based methods -L3 [14] counts the number of paths with length-3 normalized by the degree for all the node pairs.
• Network embedding methods -DeepWalk [15] performs truncated random walk exploring the network neighborhood of nodes and applies skip-gram model to learn the d-dimensional embedding for each node in the network. Node features are concatenated to form edge representation and train a logistic regression classifier. -node2vec [16] extends DeepWalk by running biased random walk based on breadth/depth-first search to capture both local and global network structure.
• Graph convolution-based methods -VGAE [10] uses graph convolutional encoder with two GCN layers to learn representation for each node in the network and adopts inner product decoder to reconstruct adjacency matrix. -GCN [17] uses normalized adjacency matrix to learn node representations. The representation for nodes are concatenated to form feature representation for the edges and fully connected layer use these edge representation to reconstruct edges, similar to HOGCN. Setting P = {1} in our proposed HOGCN is equivalent to GCN. -SkipGNN [12] learns the node embeddings by combining direct and skip similarity between nodes. Setting P = {1, 2} in our proposed HOGCN is equivalent to SkipGNN.

Experimental setup
We split the interaction dataset into training, validation, and testing interactions in a ratio of 7:1:2 as shown in Table 2.
Since the available interactions are positive samples, the negative samples are generated by randomly sampling from the complement set of positive examples. Five independent runs of the experiments with different random splits of the dataset are conducted to report the prediction performance. We use (1) area under the precision-recall curve (AUPRC) and (2) area under the receiver operating characteristics (AUROC) as the evaluation metrics. With these evaluation metrics, we expect positive interactions to have higher interaction probability compared to negative interactions. So, the higher value of AUPRC and AUROC indicates better performance.
We implement HOGCN using PyTorch [29] and perform all experiments on a single NVIDIA GeForce RTX 2080Ti GPU. We construct a 2-layered HOGC network with k = 3 for each layer. At each HOGC layer, the node mixes the feature representations from neighbors at distances P = {0, 1, 2 and 3}. The dimension of all weight matrices in HOGC layers is set to d = 32. All the weight matrices are initialized using Xavier initialization [30]. We train our model using mini-batch gradient descent with Adam optimizer [31] for a maximum of 50 epochs, with a fixed learning rate of 5 × 10 −4 . We set the mini-batch size to 256 and the dropout probability [32] to 0.1 for all layers. Early stopping is adopted to stop training if validation performance does not improve for 10 epochs. The dimension of the edge feature e ij from the bilinear layer is 64 followed by linear layers to project the edge features to edge probabilities. For baseline methods, we follow the same experimental settings discussed in [12].

RESULTS
In this section, we investigate the performance and flexibility of HOGCN on interaction prediction using four different datasets. We further explore the robustness of HOGCN to sparse networks. Finally, we demonstrate the ability of HOGCN to make novel predictions with literature-based case studies.

Biomedical interaction prediction
We compare HOGCN against various baselines on biomedical interaction prediction tasks using four different types of interaction datasets including protein-protein interactions (PPIs), drug-target interactions (DTIs), drug-drug interactions (DDIs) and gene-disease associations (GDIs). We randomly mask 20% of interactions from the network as a test set and 10% as a validation set. We train all models with 70% of interactions and evaluate their performances on test sets. The best set of hyperparameters is selected based on their performances on the validation dataset. Finally, the experiment is repeated for five independent random splits of the interaction dataset and the results with ± one standard deviation are reported in Table 3. All of our models used for the reported results are of same capacity (i.e. P = {0, 1, 2, 3} and d = 32). Table 3 shows that HOGCN achieves huge improvement over network embedding methods such as DeepWalk and node2vec across all datasets. Specifically, HOGCN outperforms Deepwalk on AUPRC by 24.44% in DTI, 28.51% in DDI, 30.07% in PPI, and 13.79% in GDI. Although node2vec achieves better performance compared to DeepWalk by adopting a biased random walk, HOGCN still outperforms node2vec by a significant margin. DeepWalk and node2vec consider different orders of neighborhood defined by the window size and learns similar representations for the nodes in that window. In contrast, HOGCN learns feature differences between neighbors at various distances to obtain feature representation for the node and thus achieves superior performance. The improved performance suggests that feature differences between different order neighbors provide important information for interaction prediction.
A network similarity-based method, L3 [14] outperforms network embedding methods across four datasets but is limited to a single aspect of network similarity i.e. the number of paths of length 3 connecting two nodes. So, L3 cannot be applied when other similarities between nodes such as similarity in features and common neighbors at various distances need to be considered. HOGCN overcomes these limitations and outperforms L3 across all interaction datasets with huge gain. In particular, HOGCN gains 3.5% AUPRC and 7.09% AUROC on PPI over L3 [14], which recently outperformed 20 network science methods in the PPI prediction problem.
Graph convolution-based methods such as GCN and VGAE achieves significant improvements over network embedding approaches but achieves comparable performance with L3. SkipGNN shows improvement over all other methods by incorporating skip similarity to aggregate information from second-order neighbors. Moreover, HOGCN with k = 3 achieves an improvement over all graph convolution-   Table 3 demonstrate that our approach with higherorder neighborhood mixing outperforms the state-of-the-art methods on real interaction datasets.

Exploration of HOGCN's drug representations
Next, we evaluate if HOGCN learns meaningful representation when feature representations of higher-order neighbors are aggregated. To this aim, we train GCN, SkipGNN, and HOGCN models on the DDI network to obtain the drug representations Z. The learned drug representations are mapped to 2D space using t-SNE [33] and visualize them in Fig. 3.
Drugbank [34] provides information about drugs and their categories based on different characteristics such as involved metabolic enzymes, class of drugs, side effects of drugs, and the like. For this experiment, we collect drug categories from Drugbank and limit the selection of drug categories such that the training set doesn't contain any interactions between the drugs in the same category. The selected drug categories are ACE Inhibitors and Diuretics (DBCAT002175), Penicillins (DBCAT000702), and Antineoplastic Agents (DBCAT000592) with 10, 24, and 16 drugs respectively. Although these drugs don't have direct interactions in the training set, we assume that these drugs share neighborhoods at various distances and can be explored accordingly with HOGCN. Fig. 3 shows the clustering structure in drugs' representations as neighborhood information at multiple distances are considered. Examining the figure, we observe that drugs in the same category are embedded close to each other in the 2D space when the model aggregates information from farther neighbors. For example, 24 drugs in the Penicillins (DBCAT000702) category (marked with blue triangles in Fig. 3) are scattered in the representation space learned by GCN that only considers feature aggregation from immediate neighbors (Fig. 3a). Note that these drugs don't have any direct interaction between themselves in the training set. Since GCN-based models can only average the representation from immediate neighbors, these drugs are mapped relatively farther to each other and closer to other interacting drugs. SkipGNN considers skip similarity to aggregate features from second-order neighbors and show relatively compact clusters compared to GCNs (Fig. 3b). On the other hand, HOGCN considers the higher-order neighborhood and learns similar representations for drugs that belong to the same category demonstrated by compact clustering structure in Fig. 3c even though no information about categorical similarity is provided to the model. This analysis demonstrates that HOGCN learns meaningful representation for drugs by aggregating feature representations from the neighborhood at various distances.
Next, we test if the clustering pattern in Fig. 3 holds across many drug categories. With this aim, we consider all drug categories in DrugBank and compute the average Euclidean distance between each drug's representation and representations of other drugs within the same drug category. We then perform 2-sample Kolmogorov-Smirnov tests and found that HOGCN learns significantly more similar representations of drugs than expected by chance (p-value = 4.93e − 106), GCNs (p-value = 5.05e − 56) and SkipGNN (1.29e − 12). Thus, this analysis indicates that HOGCN learns meaningful representations for drugs by aggregating neighborhood information at various distances.

Robustness to network sparsity
We next explore the robustness of network-based interaction prediction models to network sparsity. To this aim, we evaluate the performance with respect to the percentage of training edges varying from 10% to 70%. We make predictions on the rest of the interactions. We further use 10% of test edges for validation to select the best hyperparameter settings. For a fair comparison, we compare with graph convolution-based methods that aggregate information from direct and/or second-order neighbors.   4 shows the robustness of HOGCN to network sparsity. HOGCN achieves strong performance in all tasks with different network sparsity. The performance of HOGCN steadily improves with the increase in training edges. The mixing of features from a higher-order neighborhood in HOGCN and SkipGNN shows improvement over GCN and VGAE that only consider direct neighbors. Since HOGCN can learn the linear combination of features from a 3-hop neighborhood for this experiment, it shows improvement over SkipGNN in almost all cases. This demonstrates that features from farther distances are informative for interaction prediction in sparse networks.

Calibrating model's prediction
All graph convolution-based model proposed for biomedical link prediction predicts the confidence estimate p ij for interaction between two biomedical entities v i and v j . We thus test if a predicted confidence p ij represents the likelihood of being true interaction. In other words, we expect the confidence estimate p ij to be calibrated, i.e. p ij represents true interaction probability [35]. For example, given 100 predictions, each with the confidence of 0.9, we expect that 90 interactions should be correctly classified as true interactions.
To evaluate the calibration performance, we use reliability diagrams [35] and Brier score [36]. In particular, the reliability diagram provides a visual representation of model calibration. These diagrams plot the expected fraction of positives as a function of predicted probability [35]. A model with perfectly calibrated predictions is represented by a diagonal in Fig. 5. In addition to reliability diagrams, it is more convenient to have a scalar summary statistics of calibration. Brier score [36] is a proper scoring rule for measuring the accuracy of predicted probabilities. Lower Brier score indicates better calibration of a set of predictions. It is computed as the mean squared error of a predicted probability p ij and the ground-truth interaction label A ij . Mathematically, the Brier score can be computed as: where |E | denotes the number of test edges. Fig. 5 shows the calibration plots for GCN, SkipGNN and HOGCN (k = 3). For DTI dataset, SkipGNN show better calibration compared to GCN and HOGCN (Fig. 5a), indicating that second-order neighborhood information is appropriate and aggregating features from farther away makes model overconfident. For other datasets, GCNs are relatively overconfident for all predicted confidence. For example, approximately 20% − 30% of interactions are true positives among the interactions with high predicted confidence 0.8 in PPI (Fig. 5c) and GDI dataset (Fig. 5d). In contrast, HOGCN achieves a lower Brier score in comparison to the GCN and SkipGNN across DDI, PPI, and GDI datasets, alluding to the benefits of aggregating higher-order neighborhood features for calibrated prediction. This analysis demonstrates that HOGCN with higher-order neighborhood mixing makes accurate and calibrated predictions for biomedical interaction.

Impact of higher-order neighborhood mixing
In Section 6.3, we contrast HOGCN's performance with that of alternative graph convolution-based methods in varying fraction of edges. In this experiment, we aim to observe the performance of HOGCN when the order k is increased to allow the model to aggregate neighborhood information from farther away. We follow a similar setup as discussed in 6.3.  show that HOGCN's performances are not sensitive to the hyperparameter settings of k for all datasets since for settings k = {3, 4, 5}, we achieve comparable performances across the datasets. This analysis indicates that the 3-hop neighborhood provides sufficient information for interaction prediction across all datasets and the performance remains stable even with a large value for k.

Investigation of novel predictions
Next, we perform the literature-based validation of novel predictions. Our goal is to evaluate the quality of HOGCN's novel predictions compared to that of GCN and SkipGNN and show that HOGCN predicts novel interactions with higher confidence. We consider GDIs and DDIs for this evaluation.
We first evaluate the potential of the HOGCN to make novel GDI predictions. We collect 1,134,942 GDIs and their scores from DisGeNET [28]. The score corresponds to the number and types of sources and the number of publications supporting the associations. With the score threshold of 0.18, we obtain 17,893 new GDIs that are not in the training set. We make predictions on these 17,893 GDIs with GCN, SkipGNN, and HOGCN (k = 3). Out of 17,893 GDIs, HOGCN predicts a higher probability than GCN for 17,356 (96.99%) GDIs and than SkipGNN for 11,418 (63.8%) GDIs. Table 4 shows the top 5 GDIs with a significant increase in interaction probabilities when higher-order neighborhood mixing is considered. We also provide the number of evidence from DisGenNet [28] to support these predictions. Improvement in predicted probabilities by HOGCN models shows that aggregating feature representations from higherorder neighbors make HOGCN more confident about the potential interactions as discussed in Section 6.4. We select two predicted GDIs with a large number of supporting evidence and investigate the reason for the improvement in predicted confidence with HOGCN. Specifically, we choose gene-disease pairs (a) ABO and Pancreatic carcinoma (26 pieces of evidence) and (b) GPC3 and Hepatoblastoma ((17 pieces of evidence). To explain the prediction, the subnetworks containing all shortest paths between these pairs are selected. In particular, there are 49 shortest paths of length 3 between ABO and Pancreatic carcinoma including 6 diseases and 15 genes (Fig. 7a). Similarly, there are 20 shortest paths of length 3 between GPC3 and Hepatoblastoma including 6 diseases and 9 genes (Fig. 7b). Since these nodes are 3-hop away from each other and GCNs can only consider immediate neighbors, GCNs assign low confidence to these interactions.
Examining the subnetwork in Fig. 7a, we found that most of the diseases are related to a cancerous tumor in the pancreas and the prostate. Furthermore, pancreatic carcinoma is associated with other diseases such as Pancreatic neoplasm, malignant neoplasm of pancreas, and malignant neoplasm of prostate [28]. Since ABO is linked with diseases that are related to pancreatic carcinoma and other genes are related to these diseases as well, HOGCN captures such association (Fig 7a) even though they are farther away in the network. Similarly, HOGCN predicts association for GPC3 and Hepatoblastoma (Fig. 7b).
Next, we perform a similar case study for DDIs and evaluate the predictions against DrugBank [34]. For this experiment, we make a prediction for every drug pair with GCN, SkipGNN and HOGCN and exclude the interactions that are already in the training set. Table 5 shows the top 5 interactions with an increase in interaction probabilities when highers orders of the neighborhood are considered. As discussed in Section 6.4, HOGCN makes predictions with higher confidence compared to GCN and SkipGNN for the interactions that are likely to be a true positive.
Moreover, we validate the false positive DDI predictions of GCNs and investigate the subnetwork for these drugs in DDI networks to reason the predictions. Table 6 shows the top 5 interactions with a significant decrease in predicted confidence compared to GCN-based models. Since these DDIs are false positives [34], GCN-based models make overconfident predictions for such DDIs. In contrast, HOGCN   [41] significantly reduces the predicted confidence for these DDIs to be true positive, indicating that the higher-order neighborhood allows HOGCN to identify false positive predictions. In particular, HOGCN can identity false positive DDI between Belimumab and Estazolam even though they are 3-hop away from each other. We select subnetwork involving the drugs to investigate the reason for such predictions. Fig. 8 shows the subnetwork with all shortest paths between the drugs in Table 6. Examining the figure, we observe that the drugs in these false positive DDIs have common immediate neighbors for all cases. GCN makes wrong predictions for these DDIs with high confidence. However, SkipGNN becomes less confident about the interaction being true positive by considering the skip similarity. HOGCN further reduces the predicted confidence for Tranylcypromine and Melphalan to 0.065, indicating that there is no association between these drugs. These case studies show that HOGCN with higher-order neighborhood mixing not only provide information for the identification of novel interactions but also help HOGCN to reduce false positive predictions.

CONCLUSION
We present a novel deep graph convolutional network (HOGCN) for biomedical interaction prediction. Our proposed model adopts a higher-order graph convolutional layer to learn to mix the feature representation of neighbors at various scales. Experimental results on four interaction datasets demonstrate the superior and robust performance of the proposed model. Furthermore, we show that HOGCN makes accurate and calibrated predictions by considering higher-order neighborhood information.
There are several directions for future study. Our approach only considers the known interactions to flag potential interactions. There are other sources of biomedical information such as various physicochemical and biological properties of biomedical entities that can provide additional information about the interaction and we plan to investigate the integration of such features into the model. As HOGCN aggregates the neighborhood information at various distances and can flag novel interactions, it would be interesting to provide interpretable explanations for the predictions in the form of a small subgraph of the input interaction network G that are most influential for the predictions [42].