In 2018, Rajput et al. introduced the aBiofilm (http://bioinfo.imtech.res.in/manojk/abiofilm/) dataset, which is of great significance for the development of the bioinformatics and graph neural network fields (Rajput et al., 2018). Over the last three decades, many anti-biofilm agents have been experimentally verified to disrupt biofilms. aBiofilm organizes these data, which contain a database, a predictor, and a data visualization module. The database contains biological, chemical, and structural details of 5,027 anti-biofilm agents (1720 different ones) reported from 1988 to 2017. After eliminating redundant associations among them, a total of 2,884 known interaction associations of 1720 drugs and 140 microbes were finally obtained.

MDAD (https://github.com/Sun-Yazhou/MDAD/) is also a valuable microbe-drug association dataset, which was proposed by Sun et al. based on a variety of drug-related databases as well as a large amount of literature (Sun et al., 2018). Specifically, MDAD contains 5,505 associations between 180 microbes and 1,388 drugs collected from 993 documentation. After filtering out redundant information, a total of 2,470 microbe-drug associations were obtained, involving 173 microbes and 1,373 drugs.

DrugVirus (https://drugvirus.info/) compiles interactions involving 118 virus-targeting drugs and 83 human viruses, encompassing SARS-CoV-2 (2019-nCoV) (Andersen et al., 2020). Building upon this foundation, Lond et al. systematically extracted and curated 57 drug-virus associations from pertinent drug databases and scholarly publications, which involved 76 unique drugs and 12 distinct viruses. Ultimately, they assembled a dataset comprising 175 drugs and 95 viruses, yielding a total of 933 documented drug-virus interaction records.

First, considering that the functions of drugs are determined by their microstructures, and drugs with similar structures have similar chemical properties. So, the SIMCOMP2 tool based on the maximum common substructure between drugs is used in this paper to calculate the drug structure similarity (Hattori et al., 2010). For two drugs d _i and d _j respectively, their structure-based similarity can be expressed as DSS(d _i , d _j ). After calculating all the similarities between all drug pairs, an N _d × N _d matrix D S S ∈ R N d × N d can be obtained to represent the chemical structure similarities between N _d different drugs.

Next, for any two given drugs or microbes, the Gaussian interaction profile kernel similarity between them is calculated herein by utilizing a Gaussian kernel function based on known microbe disease associations as shown in Eq. 2: D G S d i , d j = exp − γ d A i , : − A j , : 2 (2) where A (i, :) and A (j, :) denote the ith and jth rows of the adjacency matrix A, respectively, and γ _d denotes the drug-normalized kernel bandwidth, which can be calculated by Eq. 3. γ d = 1 1 N d ∑ i = 1 N d A i , : 2 (3)

Also, this paper measures microbe similarity in two ways. The first one is the functional similarity of microbe proposed by Kamneva (2017). This computational method is mainly based on the microbial gene family information kernel protein-protein interaction association network. The second similarity between microbes is the Gaussian interaction profile kernel similarity MGS. similar to the drug similarity based on the Gaussian interaction profile kernel, for any given microbe pair m _i and m _j , it is computed using the Gaussian kernel function based on the known microbe drug associations as shown in Eq. 4. M G S m i , m j = exp − γ m A : , i − A : , j 2 (4) where A (:, i) and A (:, j) denote the ith and jth columns of the adjacency matrix A, respectively, and γ _m denotes the microbe normalized kernel bandwidth that can be computed according to Eq. 5. γ m = 1 1 N m ∑ i = 1 N m A : , i 2 (5)

Considering that not all drugs have their structures retrieved from databases, it is not possible to obtain all chemical structure similarities between drugs lacking structural information and other drugs. Therefore, in this paper, a comprehensive similarity is constructed to estimate the similarity between drugs and microbes by integrating Gaussian interaction profile nuclear similarity, microbe functional similarity, and drug chemical structure similarity. Specifically, for any two given drugs d _i and d _j , the integrated similarity between them is calculated as shown in Eq. 6: D S d i , d j = 1 2 D S S d i , d j + D G S d i , d j if  D S S d i , d j ≠ 0 , D G S d i , d j  otherwise (6) In addition, for any given microbe pair m _i and m _j , the combined similarity between them is calculated as shown in Eq. 7: M S m i , m j = 1 2 M F S m i , m j + M G S m i , m j if  M F S m i , m j ≠ 0 , M G S m i , m j otherwise (7) Then, the heterogeneous network H ∈ R ( N d + N m ) × ( N d + N m ) , shown in Eq. 8, can be constructed by combining the above integrated microbe similarity network D S ∈ R N d × N d , the integrated disease similarity network M S ∈ R N m × N m and the known drug-microbe association network A ∈ R N d × N m . H = D S A A T M S (8) Next, the model uses above newly constructed heterogeneous network H as an input to the GNN-based encoder to learn the low dimensional embedding representations of the drugs and microbes.

OGNNMDA is a graph neural network that directly processes the graph as input, effectively utilizing both node information and structural characteristics. Graph neural networks have gained significant popularity in link prediction tasks (Zhang and Chen, 2018), showcasing their widespread adoption. By leveraging the adjacency matrix H obtained earlier, Eq. 9 defines the specific formulation of the GNN. h v l = γ h v l − 1 , □ u ∈ N v , ϕ h v l − 1 , h u l − 1 , H (9) Here, l ∈ 1 … L conv , h v l ∈ R 1 × k is the embedding feature of the layer l, N v denotes the set of neighboring nodes for the node v, L _conv corresponds to the number of layers in the GNN network and the number of message-passing rounds. The dimension of the node’s embedding feature is denoted by k. In this study, the final embedding dimension is set to match the embedding dimensions used across the GNN layers. H is the microbe-drug heterogeneous network graph defined in Eq. 8, which is processed for embedding and provides edge information for the GNN. The node representation h 0 ∈ R N d + N m × k is obtained by a two-layer MLP defined by Eq. 10 and 11. The trainable variables W f c 1 , W f c 2 ∈ R N d + N m × k and B f c 1 , B f c 2 ∈ R k are involved in this process. Additionally, H init ∈ R N d + N m × N d + N m represents the initial node representation, and σ denotes the ReLU activation function. h 0 = σ W f c 2 σ W f c 1 H init + B f c 1 + B f c 2 (10) H init = 0 A A T 0 (11) The function ϕ calculates the messages transmitted between nodes, where the edge attribute is directly used as the message. The symbol □ represents the message aggregation function, and in this paper, the mean method is employed (Huan et al., 2021). This means that messages received from multiple neighboring nodes are aggregated by taking their average, resulting in message characteristics used for updating node representations. Finally, γ represents the node representation update function, which implements the ordered message-passing mechanism discussed in this paper.

In the message-passing process of a single-level GNN, a node only exchanges messages with its immediate neighbors. This pattern of neighbor message transmission at different orders aligns with the structure of the node root tree in a multi-layer GNN (Liu et al., 2020). As illustrated in Figure 2, for a node v, N v l represents the neighbor information of node v at the lth layer, and the nesting relationship of its neighbor messages at each layer can be described using Eq. 12. N v 1 ⊆ N v 2 ⊆ ⋯ ⊆ N v L conv (12)

In single-layer message passing, direct-neighbor node messages and higher-order neighbor node messages are differentially encoded to ensure orderly message delivery. Specifically, the neuron rows are aligned with the node root tree at each layer, enabling the acquisition of node feature representations with consistent nesting relationships. To implement this alignment encoding method, the neurons can be ordered by linearly arranging the neurons of each layer and considering a segmentation point, denoted as s. The information of the neighbors of the current node v, at order one or higher, can be encoded as s v l (Song et al., 2023). The segmentation point s corresponds to the nested nature of node v, and its size relationship is determined by Eq. 13. s v 1 ≤ s v 2 ≤ ⋯ ≤ s v L conv (13)

Next, we describe the node feature update function γ, which is exemplified below for a specific node v. The function can be divided into three distinct steps.

1. Compute the aggregated message representation m l ∈ R N d + N m × k for layer l.

After the previous rounds of the ordered message passing process, the final node embedding representation h L conv ∈ R N d + N m × k is obtained. This representation can be considered as the concatenation of the final embedding features of the drugs, h d ∈ R N d × k , and the microbes, h m ∈ R N m × k . In this paper, the final embedding features h _d and h _m are obtained separately using the matrix splicing approach defined in Eq. 18. h d h m = h L conv (18)

To reconstruct the adjacency matrix A′ representing possible microbe-disease associations, the bilinear decoder is employed. It is a structural component employed for predicting the probability of potential edges or links based on node embedding vectors. These decoders commonly integrate the embedding vectors of node pairs within a graph to generate a score function that assesses the likelihood of a link between two nodes. The key characteristic of bilinear decoders lies in their utilization of bilinear transformations to capture the interaction effects among nodes. Specifically, for a drug node and microbe node pair (u, v) with their respective embedding vectors h _d (u) and h _m (v), a bilinear decoder might compute the score by Eq. 19. s c o r e h d u , h m v = h d u T W h m v (19) Where W is a learnable weight matrix. This score can be interpreted as the probability of link occurrence after a nonlinear activation function transformation, so that A′ can be obtained by the bilinear decoder as shown in Eq. 20. A ′ = σ h d W B h m T (20) In the above formula, where W B ∈ R k × k represents a trainable matrix and σ x = 1 / 1 + e − x is the sigmoid function. Overall, the complete computational flow of OGNNMDA can be seen in Algorithm 1.

Algorithm 1

OGNNMDA.

Require: Known associations matrix A ∈ R N d × N m , drug similarity matrix D S ∈ R N d × N d , microbe similarity matrix M S ∈ R N m × N m and α = 600 is the number of iterations for OGNNMDA

During the experiment, positive samples were the drug-microbe pairs with known associations, while negative samples were the drug-microbe pairs without known associations. These sets of positive and negative samples are denoted as Ω⁺ and Ω⁻, respectively, for ease of description. It is important to note that the number of pairs with known associations in both the aBiofilm dataset and the MDAD dataset is significantly smaller than the number of pairs without known associations. Therefore, when training OGNNMDA, the loss function incorporates a weighted cross-entropy loss, as defined in Eq. 21. L = − 1 N d × N m λ ∑ i , j ∈ Ω + log a i , j ′ + ∑ i , j ∈ Ω − log 1 − a i , j ′ (21) In the above formula, (i, j) represents a pair of the drug d _i and microbe m _j . λ is introduced as a balancing factor, calculated as the ratio of the number of samples in Ω⁻ to the number of samples in Ω⁺. This factor helps attenuate the impact of data imbalance and emphasizes the reinforcement of known correlation information.

In this paper, the Xavier initialization method (Duong et al., 2019) is employed to initialize the trainable parameter matrices in various components of the model. These include the 2-layer fully connected layer, the ordered message-passing graph neural network layer, the bilinear decoder, and others, denoted as W f c l , B f c l | W f c l ∈ R N d + N m × k , B f c l ∈ R k , 1 ≤ l ≤ K f c , W g l , B g l | W g l ∈ R 2 * k × c s , B g l ∈ R c s , 1 ≤ l ≤ K conv , and the bias matrix W B ∈ R k × k . Furthermore, the Adam optimizer (Wang et al., 2023) is utilized to minimize the loss function. Adam combines the benefits of momentum optimization and adaptive learning rate, enabling quick convergence and adaptation to different parameter learning rates during the training process. This optimization technique enhances the training effectiveness of the deep learning model.

To prevent overfitting, the paper introduces node dropout (Piotrowski et al., 2020) and regularized dropout (Berg et al., 2017) schemes in the graph convolution layer. Node dropout can be seen as training multiple models on various sub-nodes, and the combination of these sub-nodes is used to predict unknown microbe-drug pairs (Tan et al., 2020).

In this paper, all experimental evaluations are conducted within a five-fold cross-validation setup. To ensure statistical robustness, each method is executed ten independent times for every experiment, thereby enabling the calculation of the mean value for each performance metric across these repetitions. In detail, this involves dividing all known associations in the dataset equally into 5 parts, denoted as t e s t p = t p 1 , t p 2 , t p 3 , t p 4 , t p 5 . Additionally, a subset of the same size as the known associations is randomly selected from the unknown association set. This subset is divided equally into 5 parts, denoted as t e s t n = t n 1 , t n 2 , t n 3 , t n 4 , t n 5 .

During the i − th (1 ≤ i ≤ 5) cross-validation iteration, the training set is defined as t r a i n i = t e s t p − t p i , and the test set is defined as t e s t i = t p i ∪ t n i . The final test result of the 5-fold cross-validation experiment is calculated based on the combined test set, test = test _p ∪ test _n .

Based on the previous description of the model structure, OGNNMDA incorporates several hyperparameters, including the dimension size (k) of embedded features, the number of fully-connected layers (L _fc ), the number of ordered message-passing GNN layers (L _conv ), the initial learning rate (r) of Adam’s optimizer, the total training period (α), the node dropout metrics (β), and the regularized dropout parameter (γ).

To establish initial values for these parameters, we set L _fc = 2, r = 0.008, α = 600, β = 0.6, and γ = 0.4. Subsequently, we examine the effects of different values for parameters k and L _conv through experimental analysis.

To investigate the impact of different hyperparameter values on the model, this paper performed 5-fold cross-validation (5 cv) experiments on the aBiofilm and MDAD datasets. The results for the AUROC were plotted in Figure 3, showcasing the outcomes for various combinations of the parameters L _conv and k.

From Figures 3A, B, it is evident that the optimal combination of L _conv and k is L _conv = 12 and k = 512. Therefore, this parameter setting will be utilized for OGNNMDA in subsequent experiments.

In this study, we replicate the code and data based on publicly accessible resources of these six methodologies, with all competing methods’ parameter configurations set according to their optimal values as reported in their respective publications. The 6 methods we compared OGNNMDA with are HMDAKATZ (Zhu et al., 2019a), GCNNMDA (Long et al., 2020), GSAMDA (Tan et al., 2022), SCSMDA (Tian et al., 2023), LAGCN (Yu et al., 2021), and NTSHMDA (Luo and Long, 2018), which are widely used in linkage prediction problems across various bioinformatics domains. However, due to GSAMDA not having performed experiments on DrugVirus in their paper nor specifying the construction process for the microbe-disease associations and drug-disease associations used to derive disease-based microbial and drug-Hamming similarities, comparative evaluations on DrugVirus are limited to the remaining five competing approaches.

To train and evaluate these methods, a 5-fold cross-validation experimental framework was employed. Performance evaluation was based on metrics such as AUC, AUPR, accuracy, and F1 score, chosen for their effectiveness in assessing performance. The experimental results, including the performance metrics, are presented in Tables 2–4. Additionally, ROC curves (see Figure 4A, 5A, 6A) and PR curves (see Figure 4B, 5B, 6B) were plotted to facilitate comparison among the different methods on the respective datasets.

Based on the experimental results from Table 2, it is evident that OGNNMDA achieves the highest AUC values on the aBiofilm dataset, with an average AUC of 0.9693 ± 0.0008. This is 0.65% higher than the next highest AUC value of 0.9628 ± 0.0021 obtained by SCSMDA. OGNNMDA also outperforms other methods in terms of AUPR, Accuracy, and F1-Score, with values of 0.9690 ± 0.0009, 0.9141 ± 0.0031, and 0.9151 ± 0.0026, respectively.

Similarly, in Table 3, which presents the results on the MDAD dataset, OGNNMDA exhibits superior performance across all four evaluation metrics. The comparison between the two tables suggests that OGNNMDA performs better on the aBiofilm dataset compared to MDAD. This disparity can be attributed to the sparser nature of the data in MDAD, resulting in a smaller ratio of positive to negative samples and a more pronounced sample imbalance issue.

Finally, we examine the results from Table 4, which presents the performance of all methods on the DrugVirus dataset. OGNNMDA achieved the highest AUPR score with a mean value of 0.8633 ± 0.0078; however, SCS-MDA outperformed others in terms of the AUC (0.8810 ± 0.0053), Accuracy (0.8098 ± 0.0071), and F1-score (0.8201 ± 0.0038). Notably, OGNNMDA did not maintain its leading position on the DrugVirus dataset as it did on the aBiofilm and MDAD datasets. This relative underperformance may be attributed to the smaller scale of the DrugVirus dataset compared to aBiofilm and MDAD, potentially limiting OGNNMDA’s ability to effectively train its more complex weighting parameters for optimal prediction.

To evaluate the efficacy of the ordered message-passing mechanism, this section presents ablation experiments, the results of which are presented in Table 5. In this context, GNN refers to a simple graph neural network model utilizing a mean aggregator as an encoder, while OGNN represents an enhanced ordered message-passing graph neural network model based on GNN, specifically the model proposed in this paper, OGNNMDA. The evaluation entails 5-fold cross-validation experiments on the aBiofilm and MDAD datasets, with specific parameter settings described in previous sections.

Based on the data presented in Table 5, the underlying GNN encoder exhibits poor performance on both datasets, showing a significant gap in all metrics compared to the OGNNMDA model utilizing OGNN as the encoder. Therefore, it is reasonable to conclude that the ordered message-passing mechanism effectively enhances the embedding performance of GNN, leading to improved prediction results in microbe drug association prediction.

To validate the prediction performance of OGNNMDA, case study experiments were conducted using two popular drugs and two microbes as targets. First, OGNNMDA was trained on the complete aBiofilm dataset to obtain the predicted association information neighbor matrix. Then, the top 20 most relevant objects for each target microbe and drug were filtered out. Finally, the relevant published PubMed literature was searched to validate the predicted microbe-drug association pairs against existing references. The first drug selected for the case study was ciprofloxacin, a fluorinated quinolone antibiotic, which has been extensively studied and shown to be associated with a wide range of human microbiome (Yayehrad et al., 2022). For instance, Rehman et al. (2019) demonstrated the effectiveness of amphotericin-B and 5% ciprofloxacin in blocking the growth mechanisms of Pseudomonas aeruginosa and Candida albicans. Ciprofloxacin has also shown susceptibility against Staphylococcus aureus, Staphylococcus epidermidis, Mycobacterium subspecies, Escherichia coli, and Mycobacterium tuberculosis (Smirnova and Oktyabrsky, 2018). The second drug chosen for the case study is moxifloxacin, a fluoroquinolone antibiotic (Rodríguez-López et al., 2020), known to be associated with antibiotic-resistant bacteria (ARB) (Loyola-Rodriguez et al., 2018) and Listeria monocytogenes (Rodríguez-López et al., 2020). The specific experimental results for the two drugs are presented in Tables 6, 7, respectively. These tables provide supporting literature information for the top 20 predicted microbes associated with ciprofloxacin and moxifloxacin. Upon observing Tables 6, 7, it is evident that 20 and 17 out of the top 20 predicted microbes associated with ciprofloxacin and moxifloxacin, respectively, have been validated by the available literature.

Furthermore, the first microbe selected for the case study was Aggregate Actinobacteria Accompanying Bacteria, a Gram-negative bacterium belonging to the family Pasteuriaceae (Krueger and Brown, 2020). It is primarily found in the oral cavity and is associated with various oral diseases and systemic infections (Jensen et al., 2019). In terms of its impact on human health, aggregates of Actinobacillus companionis are commonly linked to periodontal diseases, particularly aggressive forms of periodontitis. This bacterium has the ability to invade and colonize periodontal tissues, leading to inflammation, destruction of the periodontal ligament, and bone loss. Consequently, it is often found at a higher rate in individuals with severe periodontal disease. Sol et al. demonstrated that sub-killer concentrations of LL-37, Cathelicidin, and Scrambled LL-37 inhibit the biofilm formation of Actinobacillus actinomycetemcomitans and act as conditioning agents and lectins, greatly enhancing clearance by neutrophils and macrophages (Sol et al., 2013). Basavaraju et al. found that AHL lactonase hydrolyzes the lactone ring in the high serine portion of AHL, without affecting the rest of the signaling molecular structure. This inhibitory effect of AHL lactonase on group sensing of actinomycete aggregates has been observed (Basavaraju et al., 2016). The second microbe chosen for the case study was Clostridium nucleatum, a bacterium known for causing opportunistic infections and recently associated with colorectal cancer (Brennan and Garrett, 2019). In this study, Tables 8, 9 present the top 20 predicted drugs that are most relevant to Aggregate Actinobacteria Accompanying Bacteria and Clostridium nucleatum, respectively. Based on the information in the tables, 17 out of the top 20 predicted drugs for Aggregate Actinobacteria Accompanying Bacteria and 18 out of the top 20 predicted drugs for Clostridium nucleatum have been validated in the existing literature. Therefore, it can be concluded that OGNNMDA achieves satisfactory predictive performance in both microbe and drug case studies.