HerbSyner_Finder: a network community-based model for identifying synergistic combinations from herbal medicines and complex systems

Herbal medicines have been used for thousands of years for disease treatment and are often formulated by combining several herbs^[1]. Compared with single-drug therapies for complex diseases, herbal medicines show more comprehensive therapeutic effects and lower drug resistance^[2−4]. However, the multiple ingredients in herbal medicines can yield hundreds or thousands of possible drug combinations, making it challenging to explore all synergistic interactions exhaustively. Consequently, the mechanisms underlying their synergistic effects remain unclear despite the evident therapeutic effects of herbal formulas. Therefore, it is crucial to understand the multicomponent synergy in herbal medicines to unlock their potential for new drug discovery^[5,6].

Recently, network-based approaches have shown promise in accelerating drug discovery, offering systematic tools to unravel complex biological interactions, such as neural network methods^[7−9] and network proximity. In particular, target-based network proximity can quantitatively assess functional closeness among components of complex biological systems, such as compounds, therapeutic targets, and disease modules. Based on network proximity, many advanced frameworks have been developed for drug discovery^[10], including the prediction of novel drug combinations^[11] and drug repurposing^[12,13]. Among the diversity of complex biological systems, the network-based strategies are particularly valuable for investigating the synergistic effects of herbal medicines, as their theraputic effects largely depend on the underlying interactions between ingredients and herbs^[14,15]. Gan et al. established a network medicine framework to map disease symptoms and herbal targets onto the human protein interactome^[16], which was later applied to repurpose a PD-L1 inhibitor for cancer treatment^[17], and to identify synergistic ingredients for Alzheimer's disease^[18,19].

Although network-distance-based models have been successfully utilized for drug combinations, it is worth noting that synergistic effects underlying herbal medicines may vary when the same TCM formula is practiced for different diseases. Thus, a deeper analysis also requires evaluating interactions across multiple frameworks. To address this issue, we propose a network community-based framework, named HerbSyner_Finder, to prioritize synergistic ingredients or herbs by generating a combinational interaction landscape to capture their combinational characteristics (Fig. 1).

Figure 1. 

Overview of our network models. (a) Systematic comparison and key target investigation of CVA-related herbal prescriptions. (b) Integrating network proximity and community analysis for the synergistic combination discovery. (c) Pharmacodynamic validation and mechanistic exploration of the synergistic ingredient combinations.

Cough-variant asthma is a type of bronchial asthma characterized by bronchial hyperresponsiveness and chronic airway inflammation^[20−23]. Inhaled corticosteroids (ICSs) are typically considered first-line treatment for CVA due to their bronchodilation and anti-inflammatory effects^[24]. However, prolonged and excessive use of ICSs can adversely affect cardiovascular, endocrine, metabolic, and skeletal systems^[25]. Herbal medicines, particularly traditional Chinese medicine (TCM), have demonstrated clear clinical efficacy in treating CVA^[26−30]. However, the complex, synergistic molecular mechanisms of action (MOAs) within these TCM prescriptions remain underexplored. Thus, we use CVA as a case study to systematically investigate synergistic MOAs across five standard herbal formulas. With the HerbSyner_Finder model, we successfully prioritized two synergistic pairs, namely kaempferol-quercetin and berberine-luteolin, from thousands of candidate ingredients. In summary, our HerbSyner_Finder model provides novel insights into combination strategies and ultimately facilitates the discovery of combinational drugs from herbal sources.

We conducted a PubMed search using the terms 'CVA' and 'Prescription' to identify prescriptions for CVA with demonstrated therapeutic effects. In total, five prescriptions were selected for further study: Suhuang Zhike Capsule (SHZKJN, SH), Huanglong Zhike Granule (HLZKKL, HL), Suzi Jiangqi Decoction (SZJQT, SZ), Shegan Mahuang Decoction (SGMHT, SG), and Jiawei Dingchuan Decoction (JWDCT, JW). Herb ingredients of herbal formulae were collected from seven databases, including ETCM (version: 2019)^[31], SymMap2.0 (version: 2022)^[32], TCM-ID (version: 2006)^[33], TCMSP2.3 (version: 2014)^[34], TCMID 2.0 (version: 2018)^[35], TCM-Mesh (version: 2017)^[36], and TM-MC (version: 2015)^[37]. The quality of this data was manually checked, and the datasets were harmonized to integrate them. Our main rules are as follows: (1) firstly, only ingredients with structure information were extracted; (2) different types of structure information, such as ingredient name, SMILES, PubChem ID, etc., were standardized with InChIKey IDs; (3) for those duplicated ingredients with the same InChIKey IDs from the same herb, only one InChIKey ID was kept. Oral bioavailability (OB) and drug-likeness (DL) are critical parameters for evaluating the potential pharmacological effects and drug properties of compounds. Only ingredients with OB ≥ 30%, and DL ≥ 0.18 were selected for further analysis after clustering. The targets of these ingredients were obtained from the STITCH database with a Score > 700^[38]. All targets were standardized to official gene symbols for the 'Homo sapiens' species using the AnnotationDbi package in R. Disease-related genes were retrieved from our previous study^[39], and the overlap among targets was calculated using the ggVennDiagram package. Differentially expressed genes (DEGs) were analyzed using the DESeq2 package based on RNA sequencing with |log₂FoldChange| > 1, and P-adj < 0.05. Subsequently, only genes with |log₂FoldChange| > 2, and P-adj < 0.01 were adopted as CVA-associated genes for further network proximity calculations.

To capture the general characteristics of CVA prescriptions, only the disease-related genes common to all five prescriptions were selected for further investigation of the mechanisms of CVA treatment. The protein-protein interaction (PPI) networks among these shared genes were retrieved from the STRING database^[40], with the species set to 'Homo sapiens'. The critical genes in the PPI network were identified using the maximum clique centrality (MCC) algorithm in the cytoHubba app for Cytoscape. MCC prioritizes hub genes based on network topology and node centrality analysis. To elucidate the main biological processes involved in these prescriptions, we performed Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, using the clusterProfiler package in R. Only pathways or terms with P-value < 0.05 were considered significant.

Construction of a combinational network for herbal medicines between herbs, ingredients, and disease, using network proximity

Network proximity methods were used to further explore the combinatorial relationships among prescriptions, calculating ingredient-ingredient, herb-herb, ingredient-disease, and herb-disease interactions.

Ingredient-ingredient synergistic score

For all potential ingredient pairs within a single formula, we used a closest-distance algorithm based on their target networks to calculate their combined distance. The Closest distance algorithm is defined as follows:

In function (1), A and B are sets of targets separately. d_(a,b) is the shortest path length between target a in network A, and target b in network B. $ \| A\| $ and $ \| B\| $ denote the number of nodes in networks A and B separately. For each node a in network A, its shortest path length to each node in herb B is calculated, and the minimum is kept. The same process is applied to network B. Finally, the minimum values of each node in both networks A and B are summed and averaged, indicating how closely these two networks interact. Here, A represents the targets of one ingredient, and B represents the targets of another ingredient. Here, smaller d_(X,Y) indicate strong interactions among PPI networks.

Herb-herb interaction score

The distance between herb-herb combinations was calculated from their ingredient-ingredient network, with the shortest distance used as the edge weight. The distance between two herbs is the average shortest distance of their center ingredients:

where, B is the subnetwork covering all the ingredients in one herb. d_(b,u) represents the shortest path of each pairwise ingredient within a herb. The central ingredient is the one with the minimum sum of distances to the other ingredients. The shortest algorithm for calculating the distance between two ingredients is defined as follows:

where, a' and b' are the targets from ingredients a and b, separately. The $ \| a\| $ and $ \| b\| $ represents the number of targets for ingredients a and b. Finally, the shortest (a, b) is the sum of all shortest path lengths between the two target sets, averaged to provide a measure of network proximity.

Ingredient-disease association score

To measure the association between ingredients and disease, we applied the Z-score method. This network module method evaluates the relationship between the drug and the disease. The shortest path length between drug X and disease Y is defined as:

where, d_{(x, y)} is the shortest path length between target node x in drug X and target node y in disease Y. The minimum values across all nodes in disease Y are summed and averaged, similar to the closest distance algorithm, but in one direction from Y to X.

Herb-disease association score

Similar to ingredient-disease associations, we combined all ingredient targets into a single herb target set. Denoting that C_hA = (C₁, C₂, ..., C_i) is a set ingredient in herb A. C_hB = (C₁, C₂, ..., C_i) is a set ingredient in herb B. T_C = (T₁, T₂, ..., T_i) is a set of targets in one compound. Then, the targets for one herb can be denoted as:

where, T_H is the union of targets from ingredients of this herb X. We extracted DEGs from RNA-seq analysis of CVA models as disease-related genes Y. The herb-disease proximity was also calculated.

Community analysis for detection of CVA-related ingredients, and herbs, using the Louvain algorithm

To further identify synergistic ingredient-ingredient and herb-herb associations with CVA disease, we employed a community analysis based on network distance to prioritize hub modules that cluster with the disease. The Louvain algorithm detects the community structure in a combination network by optimizing a modularity metric. Modularity is a value between −1 and 1 that measures the density of links within communities. We hypothesize that herbs and ingredients clustered into a single community with disease indicate functional association with disease. Additionally, the herb-herb, and ingredient-ingredient interactions in this key module may be synergistic.

For a given network division into communities, it provides a quantitative measure of the division's quality. The modularity Q of a partition of a network is defined as:

where, A is the adjacency matrix of the network, A_ij represents the edge weight between nodes i and j, k_i is the sum of the weights of the edges attached to node i, m is the sum of all of the edge weights in the network, c_i is the community of node i, δ is the Kronecker delta function, which is '1' if c_i = c_j, and '0' otherwise.

Initially, each node in the network is assigned to its own community. For each node i, consider each neighbor j and evaluate the gain in modularity if i was removed from its community and placed in the community of j. The change in modularity Q can be approximated as:

where, $ \sum \nolimits_{in}^{} $ is the sum weight of edges inside the community to which i is moved, k_i,in is the sum weight of edges from i to other nodes in the community to which i is moved, and $ \sum \nolimits_{tot}^{} $ is the sum weight of edges to nodes in the community to which i is moved.

The node i is then placed in the community for which this gain is maximum (if positive). Once the modularity optimization is complete, a new network is built. Now, each community from the previous step is represented as a single node, and the weights of the edges between the new nodes are given by the sum of the weights of the edges between nodes in the corresponding two communities. Steps are repeated until Q is non-positive, indicating that no further improvement in modularity is possible by this process. A modularity value greater than 0.3 generally indicates a strong community structure. This means the network has well-defined, tightly connected communities with relatively sparse connections between them. In this way, the number of intra-community links in the resulting communities is maximized, while the number of inter-community links is minimized relative to a random link distribution.

In this study, the community analysis is performed on an integrated network comprising ingredient-ingredient, disease-ingredient, disease-herb, and herb-herb pairs, with edge weights defined by network distance. To identify the optimal modularity (Q value), we systematically evaluated cutoff thresholds for network distances in ingredient-ingredient and herb-herb pairs, while retaining all disease-ingredient and disease-herb associations. Ingredients and herbs that co-clustered into a single module, with modularity scores Q > 0.3, were interpreted as potential synergistic components contributing to therapeutic effects.

Experimental validation of potential synergistic ingredients with RAW264.7 cells, and airway smooth muscle cells (ASMCs)

Reagents and chemicals

Berberine (Cat.# B414323, purity ≥ 99%), luteolin (Cat.# L107329, purity ≥ 99%), kaempferol (Cat.# K107144, purity ≥ 98%), quercetin (Cat.# Q111273, purity ≥ 98.5%), and N-Nonylvanylamide (Cat.# V107237, purity ≥ 97%) was purchased by Aladdin (Shanghai, China). Dexamethasone (Dex, Cat.# 200126) was purchased from Zhejiang Xianju Pharmaceutical Co., Ltd (Hangzhou, China). Lipopolysaccharide (LPS, Cat.# 297-473-0), 3-(4,5-Dimethyl-2-thiazolyl)-2,5-diphenyltetrazolium bromide (MTT, 98%, Cat.# M2128), Aluminum hydroxide (Al[OH]3, Cat.# 239186), Dimethyl sulfoxide (DMSO, Cat.# 472301), Ovalbumin (OVA, Cat.# A5503), and Griess reagent (Cat.# MAK367) were purchased from Sigma-Aldrich (St. Louis, MO, USA). Antibodies for NLRP3 and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were supplied by Proteintech (Chicago, IL, USA), or Abcam (Cambridge, MA, USA). Antibodies for p65, phosphorylated p65 (p-p65), IκBα, and phosphorylated IκBα (p-IκBα) were purchased from Cell Signaling Technology (Beverly, MA, USA).

Preparation and assay procedures

The standard solution, Griess reagent solutions A and B, and the procedures for cell revival, culture, passaging, and cryopreservation were performed as described in our previous work^[26].

Cell viability assay

RAW264.7 cells, which had undergone at least three passages and demonstrated robust growth, were seeded into 96-well plates at a density of 1 × 10⁴ cells per well. Upon reaching 60% confluence, the original culture medium was discarded. The cells were exposed to varying concentrations of berberine, luteolin, and their combination (0–80 μmol·L⁻¹), as well as kaempferol, quercetin, and their combination (0–64 μmol·L⁻¹) for 24 h, with five replicates per group. After treatment, 20 µL of MTT (5 mg·mL⁻¹) was added to each well, and the cells were incubated for 4 h. After incubation, the culture medium was removed, and 150 µL of DMSO was added to each well. The plate was then agitated for 5 min, and the absorbance was measured at 490 nm using a microplate reader.

ASMCs were seeded at a density of 5,000 cells/well in a 96-well culture plate and cultured in RPMI-1640 medium, supplemented with 12% fetal bovine serum. Then, they were treated with LPS, and with ingredient pairs for 24 h. Next, as described above, the detection was performed after adding MTT.

Nitric oxide (NO) release inhibition assay

To evaluate the inhibition of NO release, RAW264.7 cells, after at least three passages, were seeded in 96-well plates at a density of 1 × 10⁴ cells per well. Once the cells reached 60% confluence, the original medium was replaced, and the cells were stimulated with LPS at 1 μg·mL⁻¹. The experiment was divided into the following groups: (1) blank group, to which 200 μL of serum-free culture medium was added; (2) the model group received 100 μL of serum-free culture medium and 100 μL of LPS solution; (3) 100 μL of LPS solution plus 100 μL of drug solution, with final drug concentrations of 0–5 μmol·L⁻¹; (4) 100 μL of LPS solution plus 100 μL of Dex solution. An appropriate concentration of Dex is used as a positive control drug to compare pharmacological effects. Each group had five replicates. Following 24 h of incubation, 100 μL of the supernatant from each well was transferred to a new 96-well plate. Griess reagent solutions A and B were prepared and mixed in a light-protected environment. Then, 100 μL of the mixed Griess reagent was added to each well, and the plate was shaken for 5 min. Optical density (OD) was measured at 540 nm using a microplate reader. The inhibition rate of NO release was calculated using the following formula:

Synergistic score evaluation

The Synergy Finder tool in R was used to calculate synergy scores for each drug combination using models such as Zero Interaction Potency (ZIP), Loewe Additivity (Loewe), Highest Single Agent (HSA), and Bliss Independence (Bliss). These models compare the observed effects of drug combinations with the expected effects of the drugs acting independently. The combination pairs were classified as synergistic (Score > 10), additive (–10 < Score < 10), and antagonistic (Score < –10).

Real-time quantitative PCR (qPCR) analysis

RAW264.7 cells were plated in 6-well plates at a density of 2 × 10⁶ cells, pretreated with synergistic ingredients, and then stimulated with 1 μg·mL⁻¹ LPS for 4 h. Total RNA was isolated from cells using Trizol reagent (Takara, Otsu, Japan) according to the manufacturer's protocol. RNAs were reverse-transcribed into cDNAs for quantification using the previously described method^[41]. The GAPDH gene was used as an internal control.

Western blotting (WB) analysis

After RAW264.7 cells were prepared according to the above method, the total proteins were collected using radioimmunoprecipitation assay buffer containing phenylmethanesulfonylfluoride. The following methods are based on our group's previous articles^[42]. Western blot analysis used specific antibodies, including anti-IκBα, anti-p-IκBα, anti-p-p65, anti-p65, anti-NLRP3, and anti-GAPDH. Detection was performed by using a Tanon 5200 Multi chemiluminescent substrate system (Tanon, Shanghai, China).

Experimental validation of potential synergistic ingredients on the CVA animal model

Animals and experimental protocols

Male Wistar rats (220–260 g, 8–10 weeks old) were sourced from the Laboratory Animal Center of Nanjing Qinglongshan (Nanjing, Jiangsu, China). The animals were maintained in a controlled environment at 25 ± 1 °C, with 55%–60% relative humidity and a 12 h light/dark cycle. They had ad libitum access to a standard rodent diet, and water. The CVA rat model was established as described in our previous study^[26]. In brief, rats were sensitized on days 0 and 7 with a mixture of 1 mg OVA, and 100 mg AL(OH)₃. Starting on day 14, the rats were stimulated with a 1% OVA saline solution for 10 d. A total of six groups, each consisting of five rats (n = 5), were formed based on the drug administration protocol: control, model, luteolin (20 mg·kg⁻¹), berberine (40 mg·kg⁻¹), co-administration of luteolin (20 mg·kg⁻¹), berberine (40 mg·kg⁻¹), and Dex (0.5 mg·kg⁻¹). The optimal dosages and ratios of berberine and luteolin were determined through preliminary experiments. OVA aerosol and the respective drug treatments were administered as a single daily dose from day 14 to day 24. All experimental procedures involving animals were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals, and were approved by the Ethics Committee of China Pharmaceutical University (Ethical clearance number: 2025-01-021).

Hematoxylin-eosin (H&E) staining, Masson's trichrome staining, and periodic acid-Schiff (PAS) staining

The rat lung tissue was dissected and fixed in 4% paraformaldehyde. After 24 h, the tissue was subjected to graded dehydration, then cleared, embedded in paraffin, and sectioned at 5 μm. Finally, the slices were dried on slides and stained with H&E, Masson's trichrome, and PAS for routine histological examination. The detailed staining scoring criteria have been thoroughly described in our previous work^[43].

Explore the underlying synergistic mechanisms of berberine and luteolin for CVA treatment

The gene expression profiles associated with ingredients were obtained from HERB42, which manually derived the herb/ingredient-centric classification of these data by defining a HERB experiment (EXP) as a set of control and treatment GEO samples related to an herb/ingredient. The search terms 'berberine' and 'luteolin' were screened on the EXP page of HERB, and the differential gene data were selected for subsequent analysis. The signature genes of CVA were identified in the previous section, and the overlapping genes were then used for KEGG pathway enrichment analysis. WB and qPCR further validated relevant genes in these enriched pathways.

Herb-ingredient and ingredient-target associations were extracted from public databases (Supplementary Fig. S1). To comprehensively compare TCM formulas for CVA, we analyzed five commonly used herbal formulas: SH, HL, SG, SZ, and JW. As shown in Fig. 2a, these formulas are complex systems, containing approximately 523, 440, 767, 836, and 617 ingredients in SH, HL, SG, SZ, and JW, respectively, along with 2,845, 2,679, 3,133, 3,158, and 3,023 ingredient-target pairs. The overlap analysis of ingredients and targets between the five formulas revealed that formulas SZJQT and SGMHT shared the most ingredients (n = 539), and targets (n = 3,067) (Supplementary Fig. S2a and S2b). Notably, these prescriptions share many targets (n = 2,436), with overlap rates ranging from 78% to 90% (Supplementary Fig. S3a).

Figure 2. 

Data collection for five prescriptions. (a) The overlap of herbs and ingredients among the five prescriptions. (b) The overlap of relevant targets between five prescriptions and CVA disease.

To identify the potential therapeutic targets of each CVA formula, we selected 550 downregulated, and 919 upregulated genes with |log₂FoldChange| > 1 and P-adj < 0.05. As shown in Supplementary Fig. S3b–S3f, SZ had the highest number of therapeutic targets (n = 192), followed by SH with 186 targets. Consistent with the high degree of target overlap among the five CVA formulas, these herbal formulas exhibit a high overlap rate in therapeutic targets, ranging from 84% to 94% (Fig. 2b), leading to similar pathways influencing disease treatment (Supplementary Figs S2c and S4). Consequently, these 163 common targets were utilized for further analysis.

Systematic comparison of five prescriptions for CVA on their therapeutic targets

We systematically compared the biological functions of therapeutic targets across the five prescriptions. The 163 common targets are primarily associated with the chemokine-mediated signaling pathway, response to chemokine, and cellular response to chemokine (Supplementary Fig. S5a), which regulate intercellular communication and leukocyte movement during inflammation^[44]. KEGG analysis revealed that the therapeutic effects of these prescriptions on CVA are mainly linked to inflammation and immunity, involving pathways such as viral protein interaction with cytokines and cytokine receptors, the chemokine signaling pathway, and cytokine-cytokine receptor interaction (Supplementary Fig. S5b).

PPI network analysis was also conducted separately for the five prescriptions to identify hub genes (Supplementary Fig. S6a–S6e). We noticed that these five prescriptions shared ten hub genes, namely motif chemokine ligand 10 (CXCL10), C-X-C motif chemokine receptor 3 (CXCR3), C-X-C motif chemokine receptor 1 (CXCR1), C-C motif chemokine receptor 10 (CCR10), C-C motif chemokine receptor 4 (CCR4), C-C motif chemokine receptor 8 (CCR8), G protein subunit gamma 2 (GNG2), G protein subunit beta 1 (GNB1), C-X-C motif chemokine ligand 3 (CXCL3), and neuromedin U receptor 1 (NMUR1) (Supplementary Fig. S5c). Thus, these ten genes can be seen as a typical mechanism of action across different formulas, chemokines and their receptors, and G proteins and neuropeptide receptors (Fig. 3a) for immune responses and inflammatory processes (Fig. 3b). Additionally, the mRNA expression changes of CVA disease models were further validated by these hub genes, with eight genes upregulated and one gene, NMUR1, downregulated (Fig. 3c). NMUR1 downregulation has been reported to be associated with dysfunction in airway smooth muscle contraction^[45].

Figure 3. 

Functional analysis of hub genes for CVA treatment. (a) The protein-protein interaction of the hub genes identified by MCC analysis. (b) A Sankey diagram illustrates the connections between hub-enriched genes and KEGG pathways. (c) Volcano plot showing DEGs in CVA animal models vs controls. (d) Mapping of hub genes to corresponding herbs and prescriptions.

Based on ingredient analysis, these hub genes are primarily ligands for multiple components of Ephedra sinica Stapf. (MA HUANG), Pinellia ternata (Thunb.) Breit. (BAN XIA), and Peucedanum praeruptorum Dunn (QIAN HU) (Fig. 3d). Overall, herbal medicines for CVA tend to share a substantial composition that represents similar mechanisms of action.

Identify synergistic ingredients of CVA prescriptions through a comprehensive combination landscape

We further developed a novel combination-landscape framework, HerbSyner_Finder, to identify synergistic ingredients in herbal medicines for CVA. An overall combination landscape was calculated for the eight herbs and 128 ingredients in the SH formula (Fig. 4a). Specifically, six herbs and two ingredients were identified as significantly associated with CVA (P-value < 0.05). Among them, the ingredient adenosine, and the herb MA HUANG showed the strongest associations with CVA. Notably, the herb-herb pair of MA HUANG and Eriobotrya japonica (Thunb.) Lindl. (PI PA YE) exhibited the closest association with a distance of 1.09, while quercetin and gallic acid were the most closely associated ingredient pairs.

Figure 4. 

Combinational landscape of SH. (a) Overview of the combinational network between herb-herb and ingredient-ingredient pairs, with network distance as edge weight. (b) Louvain community analysis of the proximity interaction network for SH. Ingredients clustered with the CVA disease, with OB ≥ 30% and DL ≥ 0.18, were used to construct an ingredient-disease interaction network. (c) Mapping synergistic ingredients to their respective prescriptions and herbs.

Community analysis was then used to confirm synergistic ingredients in the combined network further. We hypothesize that herbs and ingredients clustered with the disease node are more likely to be active, and exhibit strong synergistic effects. In detail, Modularity Q was used to assess network partitioning quality quantitatively (Supplementary Fig. S7a). For the Louvain algorithm for community analysis, Modularity (Q) was used as an important parameter to evaluate how detected community structures are strong and meaningful. Usually, as long as Modularity (Q) > 0.3, the discrimination is considered robust. Furthermore, to perform optimal community analysis, HerbSyner_Finder could automatically loop through different distance cut thresholds, and return the optimal Q value and the corresponding distance cut threshold. The results of CVA disease show that the Q value increased with decreased distance cutoff for ingredient-ingredient/herb-herb pairs (Supplementary Fig. S7b). Similarly, Q value increased for networks with fewer network edges (Supplementary Fig. S7c). The Q values plateaued at a cutoff of 1.3, indicating an optimal balance for community analysis. When iterating over different resolution parameters (K) on the SH combination network, the core modules with K = 19 achieved the highest Q value of 0.32 (Supplementary Fig. S7b). Among these community modules, 49 ingredients are clustered directly with the CVA node (Fig. 4b) while only three of these ingredients exhibited optimal OB and DL (OB ≥ 30%, and DL ≥ 0.18): luteolin, podophyllotoxin, and (+)-Medicarpin. Importantly, luteolin, present in all three prescriptions, showed the closest interaction with CVA (Distance = 2.19) (Fig. 4c), underscoring its potential therapeutic relevance.

Similarly, we applied HerbSyner_Finder to the other four prescriptions, namely HL, JW, SG, and SZ. Generally speaking, we identified two, six, eight, and six ingredients, and none, five, one, and four herbs, respectively. Kaempferol and quercetin show as significant synergistic ingredients in all five prescriptions (Distance = 1.23, Fig. 5). In addition to kaempferol and quercetin, berberine and luteolin were identified in the SG formulation, which clustered together and exhibited similar interactions with other ingredients and CVA (Distance = 1.24). This finding underscores the importance of using the Louvain community algorithm for combination network analysis.

Figure 5. 

Combinational characteristics via hierarchical clustering and network community analysis. The heatmap displays interaction scores for herb-herb and ingredient-ingredient pairs. The top and left panel show the corresponding herbs and prescriptions associated with these ingredients. The numerical values represent the actual distances, with smaller values indicating closer relationships.

Experimental validation of synergistic ingredients in RAW264.7 macrophages

To validate the prioritized synergistic interactions, LPS-stimulated RAW264.7 cells were treated with various combinations of these ingredient pairs for 24 h. Cellular viability was measured at concentrations ranging from 2.5 to 80 μmol·L⁻¹ using the MTT assay. We found that treatment with 20 μmol·L⁻¹ berberine and 5 μmol·L⁻¹ luteolin for 24 h did not significantly affect cellular viability compared to the control group (Fig. 6a and b). Thus, these non-cytotoxic concentrations were then used for further anti-inflammatory investigations.

Figure 6. 

MTT assay and synergy score of berberine and luteolin in RAW264.7 cells. (a) Chemical structures of berberine and luteolin. (b) Cellular viability with the MTT assay to determine the non-cytotoxic concentrations. (c) Comparison of no inhibition rates among berberine (1 μmol·L⁻¹), luteolin (0.5 μmol·L⁻¹), and Dex (80 μmol·L⁻¹) following LPS stimulation for 24 h. Statistical significance is indicated as *** P < 0.001 vs the control group; and *** P < 0.001 vs the co-administration group. Data are presented as the mean ± SEM (n = 3). (d) Average synergy scores calculated using the HAS, Loewe, ZIP, and Bliss models of berberine and luteolin.

The half-maximal inhibitory concentration (IC₅₀) of berberine and luteolin was determined based on NO production. For synergistic evaluation, we observed that the synergistic effects of 1 μmol·L⁻¹ berberine and 0.5 μmol·L⁻¹ luteolin are equivalent to those of 80 μmol·L⁻¹ Dex. This result suggested that our method is promising for synergistic compounds identification via high-throughput screening (P < 0.001, Fig. 6c). Furthermore, the synergistic score of berberine and luteolin increased with concentration. More importantly, this combination surpassed the additive effects of their individual actions, with average synergy scores of ZIP 16.35, Loewe 22.44, HAS 27.70, and Bliss 16.31 (Fig. 6d). In detail, both compounds demonstrated a concentration-dependent inhibitory effect (Fig. 7a and b). The combination of 1.25 μmol·L⁻¹ luteolin and 0.31 μmol·L⁻¹ berberine achieved optimal synergistic effects, with a ZIP synergy score of 29.84 (Fig. 7c and d). A combination of 0.625 μmol·L⁻¹ berberine and 0.31 μmol·L⁻¹ luteolin produced effects equivalent to those of 5 μmol·L⁻¹ luteolin alone, suggesting the potential to overcome individual drug resistance.

Figure 7. 

Experimental validation of synergistic ingredients in LPS-stimulated RAW 264.7 cells. (a) and (b) Dose-response curves for berberine and luteolin at various concentrations in LPS-stimulated RAW 264.7 cells over 24 h. (c) Dose-response matrix showing the effects of combining different concentrations of berberine and luteolin. (d) 3D plot of combinational landscapes with ZIP synergy scores.

The synergistic effects of kaempferol and quercetin were also studied on RAW264.7 cells. MTT results indicated a safe dosage range for the compounds of 0–4 μmol·L⁻¹ (Supplementary Fig. S8a and S8b). Kaempferol and quercetin exhibited a concentration-dependent inhibitory effect (Supplementary Fig. S8c), with average synergy scores of ZIP 10.25, Loewe 8.01, HSA 8.28, and Bliss 10.72 (Supplementary Fig. S8d–S8f). Because berberine and luteolin showed stronger synergistic effects than kaempferol and quercetin, we selected them for further exploration of their actions and mechanisms.

Experimental validation of synergistic ingredients in ASMCs

Airway hyperresponsiveness (AHR) was also employed to confirm the generalization of synergistic effects of berberine and luteolin. AHR is an important histopathological marker of asthma, and is associated with abnormal proliferation and migration of ASMCs^[46,47]. Studies have found that a LPS concentration of 10 μg·mL⁻¹ effectively stimulates ASMC proliferation^[48].

The effects of berberine and luteolin on LPS-stimulated ASMC proliferation were studied using a cell viability assay. The previously employed concentrations (0–5 μmol·L⁻¹) exhibited no cytotoxicity towards ASMC cells (Fig. 8a). Compared to the positive drug Dex, 1.25 μmol·L⁻¹ berberine and 0.625 μmol·L⁻¹ luteolin can exert a superior synergistic effect (Fig. 8b), achieving an inhibition rate of approximately 40% on cell proliferation (P < 0.001). The average synergy score (ZIP value) was 15.5, indicating a good synergistic effect (Fig. 8c–e).

Figure 8. 

Synergistic effects of berberine and luteolin on LPS-stimulated ASMCs. (a) Assessment of ASMC cell viability using the MTT assay to confirm that all concentrations tested are non-cytotoxic. ** P < 0.01, *** P < 0.001 vs the control group. (b) Comparison of cell proliferation inhibition rate among berberine (1.25 μmol·L⁻¹), luteolin (0.625 μmol·L⁻¹), and Dex (10 μmol·L⁻¹) following LPS stimulation for 24 h. Statistical significance is indicated as *** P < 0.001 vs the control group and *** P < 0.001 vs the co-administration group. Data are presented as the mean ± SEM (n = 3). (c) Dose-response curves for berberine and luteolin at varying concentrations (0–5 μmol·L⁻¹) in LPS-stimulated ASMCs. (d) Dose-response matrix of combined berberine and luteolin in LPS-stimulated ASMCs. (e) 3D plot of synergy landscape with ZIP scores.

Synergistic therapeutic effect of luteolin and berberine on CVA in vivo

To evaluate the effects of berberine and luteolin on CVA, we measured cough frequency and latency in rats. Both berberine and luteolin significantly reduced cough frequency and prolonged cough latency. The combination treatment further enhanced these effects with a synergistic improvement in alleviating cough symptoms (P < 0.001; Fig. 9a and b). Compared to the normal group, H&E staining revealed significant inflammatory cell infiltration around the bronchi and blood vessels in the model group (Fig. 9c and d) (P < 0.001), while the combination treatment group significantly reduced peribronchial inflammation compared to the single-treatment groups (P < 0.05). Masson's staining and PAS staining also demonstrated that the treatment groups effectively alleviated pulmonary fibrosis, and inhibited mucus secretion and goblet cell hyperplasia in the airway epithelium (Fig. 9d).

Figure 9. 

Synergistic therapeutic effects of berberine and luteolin in OVA-induced CVA rats. (a) In vivo experimental procedure. (b) Measurement of cough frequency and cough latency in six groups of rats. (c), (d) Lung tissue pathological sections from rats were stained with H&E, Masson, and PAS, and lung injury was assessed. (e), (f) mRNA expression of CVA-related cytokines (TNF-α, IL-6, IL-1β, MUC5AC, IL-4, IL-5, and COX-2) in rat lung tissues determined by qPCR. Scale bars: 100 μm, ^### P < 0.001 vs the control group; * P < 0.05, ** P < 0.01, *** P < 0.001 vs the OVA group; and * P < 0.05, ** P < 0.01, *** P < 0.001 vs the co-administration group. Data are presented as the mean ± SEM (n = 3).

Finally, qPCR was performed to investigate the combined effects of berberine and luteolin on inflammatory factors in the lungs of CVA rats. The model group exhibited significantly elevated mRNA levels of several inflammatory factors, including interleukin-1β (IL-1β), interleukin-4 (IL-4), interleukin-5 (IL-5), interleukin-6 (IL-6), tumor necrosis factor-alpha (TNF-α), the mucus secretion factor mucin 5AC (MUC5AC), and the inflammatory mediator cyclooxygenase-2 (COX-2) (Fig. 9e and f; P < 0.001). The combination of berberine and luteolin exhibited superior synergistic effects than Dex, alleviating airway inflammation and lung injury in CVA.

Synergistic mechanism of luteolin and berberine

To explore the molecular mechanisms of berberine and luteolin in treating CVA, we analyzed transcriptomic data from the HERB database (EXP IDs: HBEXP000181 and HBEXP000534), and identified 112 genes with CVA-related genes (Fig. 10a). Further KEGG pathway analysis revealed that these two synergistic ingredients may have therapeutic effects for CVA, with a NOD-like receptor, NF-κB, and TNF signaling pathways (Fig. 10b). NOD-like receptors (e.g., NLRP3) could activate the NF-κB pathway and thus trigger the release of pro-inflammatory cytokines and initiate immune-inflammatory responses^[49]. Thus, we hypothesize that the synergistic components alleviate airway inflammation by inhibiting the NLRP3/NF-κB pathway, thereby treating CVA.

Figure 10. 

The mechanism of synergistic treatment of CVA with berberine and luteolin in RAW264.7 cell. (a) Collection and analysis of the transcriptome data of berberine and luteolin from the HERB database. (b) KEGG pathway enrichment analysis of the synergistic components. (c) Measurement of mRNA expression of pro-inflammatory cytokines (TNF-α, IL-6, and IL-1β) in RAW264.7 cells by qPCR. (d) Measurement of mRNA expression of key chemokines and their receptors (CXCL10, CXCL3, and CCR10) in RAW264.7 cells by qPCR. (e), (f) Western blot analysis of the NLRP3 inflammasome and NF-κB signaling pathway. The protein expression levels of NLRP3, p-p65, p65, p-IκBα, and IκBα were detected by Western blot, and quantified using densitometric analysis. ^### P < 0.001 vs the control group; * P < 0.05, ** P < 0.01, *** P < 0.001 vs the LPS group, and * P < 0.05, ** P < 0.01, *** P < 0.001 vs the co-administration group. Data are presented as the mean ± SEM (n = 3)

To confirm our hypothesis, we examined the effects of berberine (1 μmol·L⁻¹) and luteolin (0.5 μmol·L⁻¹) on LPS-induced inflammation in RAW264.7 cells. Both compounds significantly attenuated the LPS-induced upregulation of IL-1β, IL-6, and TNF-α mRNA expression (P < 0.001; Fig. 10c). Notably, the combination of berberine and luteolin elicited a more pronounced reduction in these pro-inflammatory cytokines than either treatment alone (P < 0.01). Then, three key chemokines, CXCL10, CXCL3, and CCR10, through prior hub gene analysis, were further validated. Generally, combination therapy significantly downregulates CXCL10, CXCL3, and CCR10 mRNA levels (Fig. 10d), and thereby suppresses the immune cell recruitment in inflamed tissues. Mechanistically, the combination of berberine and luteolin exhibited stronger suppression of NLRP3 inflammasome components and of phosphorylation of p65 and IκBα in the NF-κB pathway induced by LPS, than single-drug treatments (Fig. 10e and f).

Although many TCM prescriptions show therapeutic effects in clinical practice, the underlying synergies of these complex biological systems remain unclear. In this context, we propose a novel computational model, HerbSyner_Finder, to identify synergistic ingredients from herbal medicines by integrating network proximity and Louvain community detection based on a combined landscape of herbs, ingredients, and diseases within herbal formulas.

Compared to existing frameworks that identify entities associated with specific diseases based on network proximity^[16−18], our study emphasizes the integrated network proximity and community analysis. First, we constructed a multi-layered interaction network as a combinational landscape to capture the complexity of herbal medicines. This integrated network characterized the intricate interplay between herbal medicines and their therapeutic effects. For instance, herb-disease and ingredient-disease interactions facilitate the identification of potential therapeutic targets and the relevance of each herb or ingredient to disease-related pathology. This concept differs from the pharmacological synergistic effect observed in wet experiments. 'Pharmacological synergy score' is calculated from the improvement of bioactive effects of two drug combinations, compared to a single drug treatment. However, to some degree, these two concepts are associated. Here, we hypothesize that the synergistic ingredients will show significant interaction compared with random ingredient pairs, and that these two ingredients interact closely with the disease in protein interactions. As a result, in this study, we aimed to screen for potential synergistic ingredients using computational network models to simulate synergistic effects in wet experiments. However, herbal formulations often yield hundreds of associations, making it challenging to pinpoint the precise molecular-level synergistic mechanisms. Thus, additional in-depth analyses are necessary to prioritize potential synergistic ingredients beyond simple proximity-based assessments. In this study, we employed the Louvain algorithm to identify clusters of tightly connected herbs, ingredients, and disease components within the network, often corresponding to synergistic modules. Louvain is one of the most popular algorithms for community detection from complex networks^[50]. In particular, the Louvain algorithm has been reported to efficiently identify relevant community structures with high connectivity in complex systems such as herbal medicines^[51].

Except for combination herbal medicines, the HerbSyner_Finder is a universal model that can be applied to other combination studies, especially in complex systems with multiple drugs targeting multiple targets. Moreover, while our study focused on CVA as a case study, the framework is broadly applicable to other diseases characterized by complex pathophysiology. HerbSyner_Finder allows us to statistically compare commonly used herbal prescriptions to capture their core pharmacological combination characteristics. Using this strategy, we observed that five CVA formulas share notable similarities in their ingredients and therapeutic targets, despite containing different herbs, suggesting comparable MOAs for CVA. Significantly, we identified ten hub targets across these formulas that have been extensively reported in relation to CVA^[52−59]. To enhance reproducibility and facilitate further research, we have made available our entire methodology, including network construction, Louvain community detection, and synergy score calculations, https://github.com/19900321/HerbSyner_Finder. An interface website has also been developed to support task submission for HerbSyner_Finder calculations at https://herbcomb.com/#/main/finder_synergy.

The HerbSyner_Finder model provides a robust framework for identifying synergistic ingredients in complex biological systems. Using CVA as a case study, we successfully identified the synergistic effects of kaempferol-quercetin, and berberine-luteolin among thousands of potential herbal medicine ingredient pairs. Furthermore, we demonstrated that the combination of luteolin and berberine effectively alleviates inflammation in CVA treatment by modulating the NLRP3/NF-κB signaling pathway. Berberine and luteolin demonstrated diverse therapeutic potential across various diseases, including autoimmune diseases^[60,61], cardiovascular diseases^[62,63], and cancer^[64,65]. Berberine is an isoquinoline alkaloid from Coptidis Rhizoma, and has been reported to have anti-inflammatory and anticancer effects^[66−69]. Luteolin is a flavonoid with strong antioxidant properties^[70−72]. Although similar in their anti-inflammatory, antioxidant, and cytoprotective effects, the synergy between berberine and luteolin has not been reported. The identification of their synergistic effects demonstrated the potential of our framework and provided new insights into the treatment of complex inflammatory diseases.

Given that the TCM prescriptions analyzed are administered orally, we investigate the optimal OB and DL properties. When conducting network analysis without OB and DL filtering, the results showed that the key bioactive compounds were consistently identified, as with OB and DL filtering. However, no OB and DL filtering introduced significant noise, such as non-bioactive organic compounds, and toxic compounds (Supplementary Fig. S9). Therefore, OB and DL thresholds are essential. On the otherhand, this approach may exclude potentially important bioactive compounds with lower OB or DL values, thereby unable to identify synergistic effects in complex herbal mixtures. To address this limitation, data from high-performance liquid chromatography-mass spectrometry could also be involved to assess compound bioavailability in the future.

Despite its promise for identifying synergistic ingredients, the HerbSyner_Finder algorithm is limited by the complexity of herbal medicines. Compound-target interactions are among the foundations of the computational framework for combinational drug discovery^[73]. Our study compiled data from TCM databases through 2022, ensuring that only the most recent data available at that time were included in our analysis. Some TCM databases have since been updated, such as ETCM, TM-MC, and BATMAN-TCM^[74,75], which may introduce additional valuable interactions. Additionally, accurately determining the concentration of each ingredient in each herb remains a challenge. In this study, all herb-ingredient pairs were treated equally on the network, regardless of each herb's proportion in the TCM formula. Future refinements should account for these quantitative aspects to improve the accuracy of synergy predictions. In addition to proximity-based methods, other machine learning-based and target perturbation-based frameworks have shown promise for drug combination studies. Unfortunately, the number of known synergistic ingredient combinations or perturbation data for herbs and ingredients remains limited, making it challenging to apply them to herbal medicines. In the future, as data in this filed accumulates, more advanced models can be incorporated into studies to further improve data quality, and the reliability of results.

In summary, we proposed HerbSyner_Finder, a novel network framework for prioritizing synergistic ingredients by constructing a combinational atlas of herbal medicines to quantify their interactions. This framework would provides deep insights into combinational strategies for disease treatment.