In this section, we would download known microbe-disease associations from public databases including HMDAD and Disbiome separately, among them, HMDAD (Ma et al., 2017) was compiled by Ma et al., in 2017, and after eliminating duplicate entries, we downloaded 450 distinct association pairs involving 39 diseases and 292 microbes. Besides, Disbiome (Janssens et al., 2018) was compiled by Janssens Y et al., in 2018, and after eliminating duplicate entries, we extracted 5,573 established associations between 240 diseases and 1,098 microbes. Detailed statistical information was shown in Table 1.

The incidence matrix E ∈ R n d × n m is a square matrix used to represent a bipartite graph where one set of vertices represents diseases ( n d ) and the other set represents microbes ( n m ). The matrix is structured such that the rows correspond to diseases and the columns correspond to microbes. Each entry e i j in the matrix E indicates the presence or absence of an interaction between disease d i and microbe m j . Specifically, if E i j = 1 , it means there’s a relationship between disease d i and microbe m j . If E i j = 0 , there’s no known relationship Equation 1 shows how it works: E i j = 1 , i f   d i   a s s o c i a t e s   w i t h   m j 0 , o t h e r w i s e (1)

The similarity D G I P d i , d j ∈ R n d × n d between diseases d i and d j , as measured by the Gaussian kernel, can be determined using Equation 2: D G I P = exp − γ d E d i − E d j 2 (2)

In the Gaussian kernel similarity, E d i − E d j refers to the Euclidean distance between two diseases. The parameter γ d as shown in Equation 3, is key in controlling how the distance affects the similarity measure: γ d = 1 / 1 n d ∑ i = 1 n d E d i 2 (3)

Absolutely, the Gaussian kernel similarity M G I P m i , m j ∈ R n m × n m can be applied to measure the similarity between microbes as well as shown in Equations 4, 5: M G I P = exp − γ m E m i − E m j 2 (4) γ m = 1 / 1 n m ∑ i = 1 n m E m i 2 (5)

Under the premise that diseases with similar characteristics are likely to interact with analogous genes (Wei and Liu, 2020; Xu and Li, 2006), we proceeded to calculate the functional similarity of diseases based on the functional associations among genes implicated in these diseases. The recently unveiled HumanNet v2.0 database serves as a potent tool for efficiently accessing gene interactions (Hwang et al., 2019; Long et al., 2021), with each interaction being accompanied by a log-likelihood score (LLS). This LLS quantifies the likelihood of a functional connection existing between genes. For a pair of diseases, denoted as d i and d j , we initially extracted their respective associated gene sets, denoted as G i = g i 1 , g i 2 , … , g i m and G j = g j 1 , g j 2 , … , g j m . Here, m represents the count of genes within set G i , and n represents the count of genes within set G j . We then defined the functional association between a single gene g and a gene set G = g 1 , g 2 , … , g k using Equation 6: P G g = max P S S g , g i (6) where g i ∈ G , and the functional similarity score between genes, represented by P S S , is defined as shown in Equation 7: P S S g , g i = L L S ′ g i , g j , i f   i ≠ j 1 , i f   i = j (7) where L L S ′ represents the normalized L L S of genes, which is defined as shown in Equation 8: L L S ′ g i , g j = L L S g i , g j − L L S min L L S max − L L S min (8) where L L S max and L L S min denote the maximum and minimum L L S in the HumanNet database, respectively. Ultimately, we articulated the disease functional similarity using Equation 10: D F U N d i , d j = ∑ g t ∈ G d i P G d i g t + ∑ g t ∈ G d j P G d j g t m + n (9)

In terms of microbe functional similarity, we employed the methodology advanced by Kamneva. (2017). To ascertain the functional similarity among microbes. To meticulously determine the functional similarity for any given pair of microbes, we initially sourced the protein-protein functional association network from the STRING v11 database (Szklarczyk et al., 2019). Utilizing the similarity scores derived therefrom, we constructed a microbe functional similarity matrix, M F U N , wherein each entry M F U N m i , m j denotes the degree of similarity between microbe m i and microbe m j .

By fusing the microbe-disease adjacency matrix with the disease Gaussian kernel similarity matrix and the microbe Gaussian kernel similarity matrix, as shown in Equation 10, we have crafted a heterogeneous network: B = D G I P E E T M G I P (10) where E T represents s transposition.

Bilinear Attention Networks (BAN), introduced by Kim in 2018, are composed of a central component known as the bilinear attention mechanism, which is designed to learn the distribution of attention by considering the bilinear interactions between input channels. This network employs two critical techniques to enhance feature interaction and manage complex data relationships: bilinear transformation and attention mechanisms. Bilinear transformation, which uses a weight matrix and an additive bias to process input features, is adept at revealing nuanced relationships within complex datasets, providing a robust framework for analyzing interactions. Its formula can be expressed as: x = a T H a + b (11)

In Equation 11, a represents the input vector to the BAN, H is a trainable weight matrix, b denotes the bias term, and x is the resulting output vector from the BAN. The forward propagation process of the BAN can be described as follows: y R e L U = R e L U H 1 x + b 1 (12) y = H 2 y R e L U + b 2 (13) y = H 2 R e L U H 1 x + b 1 + b 2 (14)

In Equation 12, H₁ denotes the weight matrix from the input layer to the hidden layer, b₁ represents the bias vector of the hidden layer, and x, defined in Equation 11, corresponds to the input vector. In Equation 13, H₂ and b₂ are the weight matrix from the hidden layer to the output layer and the output layer’s bias vector, respectively. By substituting Equation 12 into Equation 13, we derive the final output y and a streamlined forward propagation formula, Equation 14, which explicitly formalizes the computation process. The activation function used within the network is ReLU, as defined in formula 15, which introduces non-linearity to the model and helps in learning complex patterns. The feature vector that undergoes processing by this ReLU activation function is referred to as y R e L U . r e L U z = z , z > 0 0 , o t h e r w i s e (15)

By feeding the heterogeneous network B into the BAN, a low-dimensional matrix C = C d C m ∈ R n d + n m × l is produced, where indices C d and C m denote the disease nodes and microbial nodes, respectively. During the computation, the cross-entropy function is utilized for optimization purposes.

To effectively capture both the local and global topological intrinsic features of nodes, we have further implemented an enhanced version of Random Walk with Restart (RWR) on the D F U N . The RWR is defined as shown in Equation 16: r i l + 1 = φ X r i l + 1 − φ θ i (16)

In Equation 16, φ denotes the restart probability. X signifies the transition probability matrix, and θ i represents the initial probability vector for node i . The definition of the initial probability vector is as shown in Equation 17: θ i j = 1 , i f   i = j 0 , o t h e r w i s e (17)

Following the aforementioned RWR process, it becomes evident that we can derive a new matrix D R W R . Similarly, after applying the improved RWR on M F U N , we can obtain a new matrix M R W R .

Consequently, by amalgamating all the matrices E, D F U N , and D R W R , as shown in Equation 18: we can construct a new disease matrix D N : D N = E ; D F U N ; D R W R (18)

Similarly, by integrating E,   M F U N , and M R W R a new microbial matrix M N can be obtained as shown in Equation 19: M N = E T ; M F U N ; M R W R (19)

Then, use the above two matrices as inputs to the sparse autoencoder (SAE). SAE excel in feature extraction and dimensionality reduction, enabling them to distill crucial features from intricate microbial data and reduce its complexity, which is particularly valuable for managing high-dimensional biomedical data. Additionally, the incorporation of sparsity penalties within SAE helps to constrain the activation of neurons in the hidden layer, thereby enhancing the model’s feature extraction capabilities. This sparse representation not only boosts predictive accuracy but also contributes to the interpretability of the model. SAE consists of the following steps:

Encoding process: Input data x is converted into a hidden layer representation h through an encoder, and R e L U is used as a non-linear activation function after linear transformation. The specific formula is as shown in Equation 20: h = R e L U W e n c o d e r x + b e n c o d e r (20)

Among them, W e n c o d e r is the weight matrix of the encoder, and b e n c o d e r is the bias term.

Decoding process: The hidden layer representation h is reconstructed back to the original data x ′ through the decoder. The definition of x′ is shown in Equation 21. This process is also a linear transformation followed by a nonlinear activation function R e L U : x ′ = R e L U W d e c o d e r x + b d e c o d e r (21)

Among them, W d e c o d e r is the weight matrix of the decoder, and b d e c o d e r is the bias term.

Refactoring loss: Refactoring error is an indicator that measures the difference between the reconstructed data x ′ and the original data x , represented by the binary cross entropy (BCE) loss function. The specific form is shown in Equation 22: L r e c o n = B C E x , x ′ = − ∑ i x i ⁡ log x i ′ + 1 − x i log 1 − x i ′ (22)

Sparsity loss: To introduce sparsity, SAE adds a sparsity penalty term to the loss function, which is typically based on L1 regularization. The specific form is shown in Equation 23: L s p a r s i t y = λ ∑ j h j (23) where λ is the regularization coefficient and h j is the activation value of the hidden layer.

Total loss function: The total loss function, which is the sum of the reconstruction loss and the sparsity loss, serves as the objective function for optimization during the training process. It can be articulated as shown in Equation 24: L t o t a l = L r e c o n + L s p a r s i t y (24)

Consequently, by feeding the disease matrix D N and the microbe matrix M N into the SAE individually, we can derive matrices D S and M S , respectively.

Based on the processing results of BAN and SAE models, by integrating the disease matrix C d , D S , D R W R , D F U N and the adjacency matrix E, inspired by Xuan et al. (2020) (Long et al., 2021), we can construct a new disease feature matrix D as shown in Equation 25: D = C d ;   D S ;   D R W R ;   E ;   D F U N ;   E   (25)

Similarly, integrating the microbial matrix C m , M S , M R W R , M F U N and the adjacency matrix E, we can construct a new microbe feature matrix M as shown in Equation 26: M = C m ;   M S ;   E T ;   M R W R ;   E T ;   D F U N   (26)

The dot product of two vectors serves as an effective mechanism for modeling interactions, highlighting the shared aspects of these interactions while diminishing the distinct information they might carry. Consequently, for any given disease d i and microbe m j , we can determine their predicted association scores by computing the inner product of their feature representations, as shown in Equation 27: R i j = S i g m o i d D d i · M m j T (27)

In this part, we began by conducting a sensitivity analysis of crucial parameters to improve the model’s effectiveness. Next, we chose six state-of-the-art techniques to benchmark against BANSMDA. Additionally, to confirm the model’s reliability, we selected two exemplary microbes and diseases for evaluation.

Considering the actual situation of the model, we identified and analyzed four parameters that have a significant impact on the final prediction results. These include the L2 regularization parameter λ ,which is named l 1 ,in the BAN model, φ in the RWR of formula (14), the learning rate l 2 and the average activation ρ in the SEA model.

In this section, our objective is to determine the optimal settings while maintaining the separation of the training and testing datasets. Specifically, The range of values for l 1 is set to 0.0001 , 0.0005 , 0.001 , 0.005 , 0.01 , 0.05 , 0.1 . The range of values for φ is 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 . The range of values for l 2 is 0.0001 , 0.0005 , 0.001 , 0.005 , 0.01 , 0.05 , 0.1 . The range of values for ρ is {0.001,0.005,0.01,0.05,0.1}.Subsequently, we employed a 5 - fold cross - validation (CV) method to evaluate the area under the receiver operating characteristic curve (AUC) and the area under the precision - recall curve for the parameter configurations. In the parameter validation experiment, we first set one parameter, e.g., l ₂_lamda, to a fixed value. Then, in each epoch, we changed the value of one of the other parameters. After that, we collected all AUC and AUPR values under such circumstances. Finally, we calculated the average of the obtained AUC and AUPR values respectively, and used them as the final results for that particular parameter setting. In the 5-fold CV experiment, we first randomly assigned 80% of the dataset, including both identified and unidentified associations, to the training set, with the remaining 20% reserved as the independent test set. We then divided the training set into five equal-sized subsets to perform 5-fold cross-validation. Using the HMDAD dataset, we independently conducted five cross-validation runs. After the cross-validation was completed, the model’s performance was evaluated on different subsets of the training set. Finally, we used the pre-defined independent test set to assess the model’s final performance. As shown in Figure 2, the model achieves the best performance when the parameter value are configured as follows: l 1 = 0.01 , φ = 0.4 , l 2 = 0.005 , ρ = 0.005 .

To rigorously assess the predictive performance of the BANSMDA model, we conducted case study validation using two prevalent diseases (Asthma and Colorectal carcinoma) and two clinically significant microbial genera (Escherichia and Bacteroides).

Asthma is a common chronic inflammatory disease originating in the lower airways, characterized by persistent airway inflammation (Polyxeni et al., 2021). Clinical manifestations include recurrent episodes of wheezing, coughing, chest tightness, and shortness of breath, typically exacerbated during nocturnal or early morning periods (James, 2015). A growing body of evidence from multiple references underscores a significant correlation between the pathogenesis of asthma and specific microbiota, such as Helicobacter pylori (Zhi et al., 2021), Proteobacteria (Kian, 2017), and Bacteroidetes (Chie, 2023). Based on the predictive scores, microbes associated with asthma were ranked in descending order according to their respective scores. As illustrated in Table 3, among the top 20 predicted microbes associated with Asthma, 19 have been confirmed by existing research indexed in PubMed.

Colorectal carcinoma ranks as the third most common cancer globally (Chie, 2023; Inés et al., 2017). The gut microbiota is intricately involved in its development, with ecological imbalances capable of inducing colorectal carcinoma through chronic inflammatory pathways. Key bacterial taxa implicated in this multifaceted process include Clostridium (Hui et al., 2023), Bacteroides (Yasutoshi et al., 2024), and Enterobacteriaceae (Rashmi et al., 2016). As illustrated in Table 4, all of the top 20 predicted microbes associated with Colorectal carcinoma have been confirmed by existing studies in PubMed.

Escherichia is a bacterium that embodies a dual identity, capable of functioning as both a symbiotic microbe and a pathogenic agent within the host’s body (Olivier et al., 2010). Recent research has demonstrated that specific strains of Escherichia are capable of causing a range of intestinal infections, including diarrhea and enteritis (James, 2005). Moreover, Escherichia can extend its pathogenicity beyond the gut to cause extraintestinal infections through mechanisms like fecal contamination or hematogenous dissemination (Kevin et al., 2019). As illustrated in Table 5, among the top 20 predicted diseases associated with Escherichia, 19 have been confirmed by existing research indexed in PubMed.

Bacteroidetes are significant clinical pathogens that, when they breach the intestinal barrier, can induce severe pathology. This includes bacteremia and the formation of abscesses in various parts of the body (Hannah, 2007). As illustrated in Table 6, among the top 20 predicted diseases associated with Bacteroidetes, 19 have been confirmed by existing research indexed in PubMed.

In summary, these case studies provide additional evidence for the ability of the BANSMDA model to predict potential associations between microbes and diseases.