We constructed the PPI network of Escherichia coli based on the STRING database version 11.5 (Szklarczyk et al., 2020). The interactions with a score less than 0.7 were ignored. The ultimate network included 59, 496 interactions involving 4, 020 proteins.

We collected drug–target interactions from previous literature reports or the DrugBank database (Wishart et al., 2017). Since we used the data from the in vitro antimicrobial test, proteins from bacteria were considered and proteins of Homo sapiens were ignored in this study.

In this study, we focused on pairwise antibiotic combinations by recent experimental data of the Escherichia coli strain MG1655 (Chandrasekaran et al., 2016). Interactions were quantified based on the α-score, and the following three types were used: synergy (α-score ≤ −0.25), additive (−0.25< α-score < 1), and antagonism (α-score ≥ 1) (Cokol et al., 2011). In this study, we only considered antibiotics with known targets (protein or RNA). In total, 91 pairwise antibiotic combinations involving 14 antibiotics were retained.

Targets of the aforementioned antibiotics were fed into a network propagation model (Vanunu et al., 2010) to simulate the propagation of drug-action effects on the PPI network: F t + 1 = β A ′ F t + ( 1 − β ) F 0 , (1) where β is a parameter ( 0 ≤ β ≤ 1 ). Vanunu et al., (2010) confirmed that the algorithm is not sensitive to the choice of β as long as it is above 0.5, so we set β to 0.7. A ' = D − 1 / 2 A D − 1 / 2 in which D is a diagonal matrix, and the values of the diagonal element (d _ii ) are equal to the degree of the vertexes (k _i ), and A is an adjacency matrix. F ₀ denotes a matrix, in which “1” indicates that the drug can bind to this target, and if the drug can bind multiple targets (n), the values were set as 1/n. During each iteration, nodes can not only receive the information from their neighbors (first term) but can also retain their initial information (second term), resulting in a DAPM. Next, let us show that formula (1) converges. The general term formula of formula (1) is F t = ( β A ′ ) t F 0 + ( 1 − β ) ∑ i = 0 t − 1 ( β A ′ ) i F 0 (2) Since 0 < β < 1 and the eigenvalues of A ' are in the range of −1 to 1 (according to the Perron–Frobenius theorem), lim t → ∞ ( β A ′ ) t = 0   and   lim t → ∞ ∑ i = 0 t − 1 ( β A ′ ) i = ( I − β A ′ ) − 1 (3)

Hence, lim t → ∞ F t = ( I − β A ′ ) − 1 (4)

Subsequently, the Jaccard index (Eq. 5) and network proximity model (Eq. 6) (Menche et al., 2015; Cheng et al., 2019) were used to quantify the relationships of each DAPM: J A B = | S A ∩ S B | | S A ∪ S B | (5) where S _A and S _B are the nodes of drug A and drug B in their DAPMs, respectively. S A B ≡ < d A B > − < d A A > + < d B B > 2 (6) where < d A A > and < d B B > are the mean shortest distances between each pair of nodes in the DAPM. < d A B > is the mean shortest distance between each pair of nodes between the DAPM of drug A and the DAPM of drug B: < d A B > = 1 | | A | | + | | B | | ∑ y ∈ A m i n x ∈ B d ( x , y ) (7) where A and B are the DAPMs of drug A and drug B, respectively. d ( x , y ) is the shortest distance between node x and node y. In the next section, we demonstrated how to build the affinity matrix and graph regularization model based on network proximity.

Given three drugs (drug A, drug B, and drug C), if drug A–drug B is a synergistic antibiotic combination and drug A and drug C are pharmacologically similar, then drug C–drug B will likely be a synergistic antibiotic combination. Therefore, we can define a loss function: E ( Y ) = 1 2 ( ∑ i , j = 1 l + u W i j | | Y i d i − Y j d j | | 2 ) + μ ∑ i = 1 l | | Y i − Y ^ i | | 2 (8) where Y and Y ^ are the entire and the known drug combination matrix, respectively. μ ( μ > 0 ) is the regularization parameter. d _i is the degree of node i. The key to this model rests on the construction of the affinity matrix W, which is calculated by Eq. 9. W i j = { 1 ,   − 1 < S i j < 0 ; 0 ,   o t h e r w i s e (9) where S _ij is the network proximity (Eq. 6) between drug i and drug j. The classifying model is as follows. Y ∗ = a r g m i n   E ( Y ) (10)

We then take the derivation of E ( Y ) with respect to Y to solve the optimization problem. ∂ E ( Y ) ∂ Y | Y = Y ∗ = Y ∗ − Q Y ∗ + μ ( Y ∗ − Y ^ ) = 0 (11) where Q = D − 1 / 2 W D − 1 / 2 . The detailed derivation of Eq. 11 can be found in the supporting information, and the analytical solution of Eq. 11 is Y ∗ = δ ( I − γ Q ) − 1 Y ^ (12) where I is the identity matrix, γ = 1 1 + μ , and δ = μ 1 + μ .

The performance of the graph regularization model was estimated using the precision (Eq. 13), recall (Eq. 14), accuracy (Eq. 15), and F1 (Eq. 16), and these evaluation metrics can be defended as follows: p r e c i s i o n = T P T P + F P (13) r e c a l l = T P T P + F N (14) a c c u r a c y = T P + T N T P + F P + T N + F N (15) F 1 = 2 × p r e c i s i o n × r e c a l l p r e c i s i o n + r e c a l l (16) where TP, FP, FN, and TN are true positive, false positive, false negative, and true negative, respectively.

We used previously reported antibiotic combinations involving 14 antibiotics (Chandrasekaran et al., 2016) listed in Table 1. These antibiotics range over various mechanisms of action, including protein biosynthesis, DNA and RNA replication, folate metabolism, and cell wall biosynthesis. Since we concentrated on the subtle differences among synergy, additive, and antagonism, all 91 pairwise combinations fall into three categories, according to the α-score (Supplementary Table S1). Targets of these antibiotics were collected from previous literature studies (Pongs et al., 1973; Shen and Pernet, 1985; Buck and Cooperman, 1990; Pan et al., 1996; Onodera et al., 2002; Aracena et al., 2014; Kocaoglu and Carlson, 2015; Wekselman et al., 2017; Lin et al., 2018; Salehi et al., 2020; Wróbel et al., 2020). Because some antibiotics are RNA-targeted small molecules, ribosomal proteins that affect antibiotic binding are considered targets of antibiotics. For example, 30S ribosomal proteins S7 (rpsG) and S14 (rpsN) were shown to be the most important for tetracycline binding (Buck and Cooperman, 1990). Mutations of 50S ribosomal proteins L22 (rplV) and L4 (rplD) will lead to macrolide (erythromycin, etc.) resistance (Wekselman et al., 2017).

Network analysis showed that the shortest distance between targets of antibiotic combinations ranged from 0 to 5 (Supplementary Figure S1). Most antibiotic combinations (92.3%) did not share the same targets. Approximately thirty percent of antibiotic combinations were adjacent, and almost half of synergistic antibiotic combinations (57.1%) were included (Supplementary Figure S1). However, a considerable portion of antagonistic or additive antibiotic combinations have adjacent targets, but they are not synergistic (Supplementary Figure S1). Therefore, mere knowledge of the network structure may not be sufficient to explain the intricate interactions among antibiotic combinations and their targets. To investigate the network-based relationship between antibiotic combinations and their targets, we introduced network propagation (Vanunu et al., 2010) to predict the effect of antibiotics and antibiotic combinations on the PPI network.

Network propagation has been used to quantify the influence of mutations in colorectal tumorigenesis (Shin et al., 2017). When a mutation arises in a node, perturbation spreads out along the protein–protein interaction (PPI) network and eventually forms a mutation-propagating module. Similar to mutation, if a drug acts on a node, neighboring nodes are also affected (Figure 1B) (Saraswathi et al., 2009). Predictably, the impact is greatest in its neighbors, whereas nodes far away from targets receive attenuated influences. Therefore, we can generate a subnetwork with drug targets as hubs, and the nodes ( F i * ≥ 0.0065 ) will be incorporated into the subnetwork.

Based on the network propagation method (Eq. 1), these antibiotics and antibiotic combinations were mapped to the PPI network to investigate the potential relationships of these subnetworks (Figure 1). On average, DAPMs include approximately 13 nodes, although almost all drugs only have 1 to 2 targets. Because of the high threshold, each DAPM consisted almost exclusively of nearest neighbors. Interestingly, we found that there are areas of overlap for some antibiotic combinations and that antibiotic combinations are associated with the synergy effect (Figure 1C). Hence, we inferred that the structure of DAPKs can be used to quantify interactions between drug pairs and thereby predict synergetic antibiotic combinations.

Network proximity was used to investigate FDA-approved drug combinations (Cheng et al., 2019) and herb combinations in traditional Chinese medicine (Wang Y. et al., 2021; Zhang et al., 2021). Compared with random herd pairs, herd pairs in traditional Chinese medicine formulas tend to have smaller network proximity (Wang Y. et al., 2021). To probe whether it could also be used to distinguish synergy, additive, and antagonism, we used the Jaccard index (Eq. 5) and network proximity (Eq. 6) to quantify DAMP–DAMP interactions. We found that all possible antibiotic combinations can be divided into three topologically distinct categories: a) overlap: two DAMPs overlap but do not equate (Figure 2A); b) separation: two DAMPs are topologically separated (Figure 2B); and c) identical: two DAMPs are completely consistent (Figure 2C).

For overlap (Figures 2A,D), these antibiotic combinations are probably synergetic (87.5%, p − v a l u e = 0.118 , permutation test). From the perspective of network pharmacology, if DAMPs of two drugs overlap, it indicates that the two drugs are pharmacologically similar (Cheng et al., 2019). For example, chloramphenicol and erythromycin both target the 50S ribosome, and their binding sites are the peptidyl transferase center (PTC) and the nascent peptide exit tunnel (NPET) on the 50S subunit, respectively (Lin et al., 2018). They can inhibit protein synthesis in a synergistic way (Figure 3B) (Chandrasekaran et al., 2016). As shown in Figure 3A, DAMPs of chloramphenicol and erythromycin have common nodes. Hence, the network proximity of the two DAMPs is negative, S T E T − C H L = − 0.97 . To verify this idea, we performed virtual screening for nodes in the DAMP of trimethoprim (a dihydrofolate reductase inhibitor). Eventually, we identified a dihydropteroate synthase inhibitor—sulfamethoxazole. The DAMPs between sulfamethoxazole and trimethoprim overlap ( S T R I − S U L = − 0.12 , Supplementary Figure S2A). Previous studies have shown that a combination of trimethoprim and sulfamethoxazole not only interferes with folic acid synthesis synergistically (Yeh et al., 2006) but also reduces the risk of bacterial resistance (Pappas et al., 2009). In summary, synergistic drug combinations tend to act on the same biological pathways.

For separation (Figures 2B,E), these antibiotic combinations were almost not synergetic (90.1%, p − v a l u e < 10 − 4 , permutation test, see more from SI). In other words, the two drugs are pharmacologically distinct in this case. For example, nalidixic acid (an inhibitor of DNA gyrase) and chloramphenicol (an inhibitor of protein synthesis) take effect in different biological processes, so their DAMPs are topologically separated ( S N A L − C H L = 0.92 , Supplementary Figure S2B), and nalidixic acid and chloramphenicol do not show the synergy effect. Levofloxacin not only inhibits DNA gyrase but also inhibits DNA topoisomerase (Table 1). Hence, the DAMPs of levofloxacin and nalidixic acid overlap, resulting in the synergy effect. In addition, DAMPs of some synergistic drug combinations are topologically separated. This may result from the following reasons: a) experimental data itself: the correlation coefficient of the α-score between two replicates is only 0.81, which leads to a random error; b) some drugs have unknown targets: recent evidence suggests that gentamicin has a second binding site around H69 of the 23S rRNA of the 50S ribosome (Serio et al., 2018). This may be the reason for the synergy between gentamicin and tetracycline.

For identical (Figures 2C,F), these antibiotic combinations showed a definite additive effect (100%). For example, clarithromycin and erythromycin not only act on the same targets (Table 1) but also have similar chemical structures (98.1%, Tanimoto similarity; more details can be found in SI). Hence, we consider the two drugs to be pharmacologically identical, which leads to an additive effect.

To demonstrate the usefulness of the PPI network, an ablation test was performed where the PPI network was randomized. Supplementary Figure S3 shows that the randomized PPI network produces worse results, so an accurate PPI network is crucial for our model.

Graph regularization is a useful model to predict different relationships between various types of biological entities (Luo et al., 2018; Ding et al., 2020). Through the aforementioned analysis, we found that if two drugs are pharmacologically similar, then the drug pair is probably a synergistic antibiotic combination (Figure 2D). Therefore, we can define a loss function of Y (Eq. 8). However, if two drugs are pharmacologically identical (S _AB = −1), then the drug pair shows an additive effect (Figure 2C). Therefore, we set W _ij of these drug pairs to 0 (Eq. 9). Next, we used the aforementioned 14 antibiotics (Table 1) for the training set to predict interactions with the following three antibiotics (Table 2) by Eq. 12. The entire predicted scores are listed in Table 3. A larger predicted score of drug pairs suggests that they would probably be the synergistic antibiotic combinations. In Supplementary Tables S2–S6, we confirmed that the algorithm is not sensitive to the choice of γ , so it was simply fixed at 0.7. In Supplementary Figure S4, we demonstrated that the impact of the threshold changes on the performance of our method. When the threshold increases from 0.1 to 0.5, the precision increases, and the recall and accuracy decrease. When the threshold is larger than 0.2, the F1 decreases. Therefore, we set the threshold to be 0.2 in our model. Eight potential synergistic antibiotic combinations were found: TET-ROX, ROX-CLA, OXA-PNG, CEF-PNG, ROX-ERY, ROX-CHL, PNG-TET, and PNG-TRI. In the experiments conducted by Mason et al. (Mason et al., 2017), TET-ROX, ROX-CLA, OXA-PNG, CEF-PNG, and PNG-TET were identified as synergistic antibiotic combinations, and ROX-ERY, ROX-CHL, and PNG-TRI were additive. However, as alluded to above, there are random errors in experimental measurements, which might have some impact on the classification of antibiotic combinations. As expected, we found that ROX-ERY and ROX-CHL were identified as synergistic antibiotic combinations in the experiments by Yilancioglu (Yilancioglu, 2019). This means that our model has good stability for the prediction of synergistic antibiotic combinations.

Previously, there have been studies to predict synergistic antibiotic combinations through computational methods. In this section, we compared the performance of our model (Eqs 13–16) with other methods, such as CosynE (Mason et al., 2017) and INDIGO (Chandrasekaran et al., 2016) on the benchmark dataset. As shown in Table 4, our model achieved better performance in terms of various metrics.