We collected the disease-related gene set for HCC from the DisGeNet (http://www.disgenet.org/) and the GeneCards (https://www.genecards.org/) database. Only genes that appeared in both datasets were retained in this study to ensure data reliability. For five TCM syndromes of HCC and other syndromes, their gene sets were collected from the SymMap (http://www.symmap.org/) database, which has already linked TCM symptoms to gene targets through symptom mapping followed by statistical verification (17).We extracted statistically significant relationship data from the database, i.e. screening with a threshold of P-value < 0.05. The datasets of disease/syndrome -related genes/proteins are provided in Supplementary Data S1 .

We collected the disease-related small molecule drugs for HCC, liver disease, and cancer based on the clinical practice guidelines from the European Society for Medical Oncology (ESMO) (18), the European Association for the Study of the Liver (EASL) (19), and the DrugBank database. For five TCM syndromes of HCC, their associated treatments as TCM formulae were collected from the 2019 version of Chinese guidelines for the diagnosis and treatment of HCC (10). The datasets of disease/syndrome-related drugs/formulae are provided in Supplementary Data S2 .

We collected the targets for small molecule drugs of HCC, liver disease, and cancer from the DrugBank database (https://go.drugbank.com/). For TCM formulae, we curated their potential targets through their herb/ingredient constituents based on the curated targets in our previously constructed HERB database (20). We also quantitatively estimated the drug-likeness metric for these herbal ingredients using the RDKit python package (version 2023.03.1) and screened for potential active compounds in the medicinal herbs. The process of determining the quantitative estimation of drug-likeness (QED) threshold is shown in Supplementary Appendix S1 ( Supplementary Figure S1 ). The datasets of drug/formulae targets are provided in Supplementary Data S3 .

We integrated a comprehensive human PPI network from six different sources of experimentally verified data. It includes: 1) HuRI: Known for its binary protein interactions that are typically direct, this bioinformatics resource maps the biophysical human protein interactome using yeast two-hybrid systems, further validated through orthogonal assays, ensuring a high degree of confidence in the interactions (21); 2) BioPlex: With many interactions involving large protein assemblies detected by affinity purification mass spectrometry (22); 3) Interactome3D: This database integrates interaction data from primary pathway repositories, compiling experimentally verified binding events (23); 4) Insider: As a tool that bridges genomic variation information with the structural protein-protein interactome, Insider provides interactions that have been determined experimentally (24); 5) PhosphoSitePlus: Encompassing a dataset of 12,180 high-quality, manually curated kinase-substrate interactions (25); 6) InnateDB: Beyond its 18,000 manually managed interactions, InnateDB also includes 178,000 experimentally verified imported interactions and 3,000 pathways (26).

Then we uniformly mapped all proteins in the network to Entrez gene IDs with official gene symbols based on the R package ‘org.Hs.eg.db’ (version 3.18.0) (27). Next, we calculated the largest connected component of the PPI network using the python package ‘NetworkX’ (version 3.1) (28). As a result, we obtained a large human interactome including 311,499 unique PPI interactions linking 17,345 proteins/genes. We mapped all drug targets in the DrugBank database (29) to this network and found that 94.24% of the drug targets were mappable, indicating that our network is capable of analyzing drugs based on their targets. The PPI network used in this study is provided in Supplementary Data S4 .

Similar diseases exhibit similar patterns in both phenotypes and molecular mechanisms (30, 31). In the disease/TCM syndrome-related network, it was expected that diseases with similar clusters of symptoms will be closely related to adjacent network modules, and vice versa. To verify it, we quantified the distances between the HCC disease (D) and five TCM syndrome (S) related subnetworks by using the separation measure S_DS proposed by Bara´basi et al. (30) (Equation 1). The smaller the value, the closer the two gene modules are in the background network (The principle of the calculation is shown in Figure 2A ).

We employed the disease-drug proximity measure proposed by Bara´basi et al. to quantify the network intersection between disease or TCM syndrome related protein set, and drug or formula related protein set (32). It is a method to calculate distances between two groups of nodes in the network while correcting for degree biases. In this work, the metric calculated the average shortest distance between the drug’s targets and the nearest disease protein on the PPI network, where H represents proteins for disease or TCM syndromes, and T represents drug/formula targets (Equation 2) (The principle of the calculation is shown in Figure 3A ).

Then, 100 random pairs of protein targets with the same size and same degree distribution were sampled from the whole network for the permutation test. Based on this, the z-scores and P-values for the proximity measures between drug-HCC disease, formula-HCC disease, and formula-TCM syndrome were calculated respectively (Equation 3). The criteria for significance is set to be z-score<−1.5 and P-value< 0.05 (33).

where d ¯ and σ are the mean and standard deviation of the proximity measures based on 100 random samples.

The relevant computational codes for network separation and proximity measures can be accessed via the following GitHub repository: [Barabasi-Lab/COVID-19] (https://github.com/Barabasi-Lab/COVID-19).

For gene sets related to either HCC diseases, TCM syndromes, and formulae, we analyzed their potential functions through the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. This analysis assesses the biological pathways of a particular gene set and was performed using the clusterProfiler R package with the following parameters: pvalueCutoff = 0.05, pAdjustMethod = “BH”, qvalueCutoff = 0.2, and the background gene set was the human gene set (34).

The R package used for KEGG enrichment analysis is clusterProfiler, version 4.10.1, available on GitHub at: [YuLab-SMU/clusterProfiler](https://github.com/YuLab-SMU/clusterProfiler).

We also prioritized the active herbal ingredients that constitute TCM herbs. Functionally related nodes tend to be close to each other in the network, which can be calculated based on the guilt-by-association method. In this study, we constructed a comprehensive heterogeneous network of drug-target-PPI, using a random walk with restart (RWR) algorithm (35), and took disease genes as seed nodes, to explore key ingredients in TCM formulae. We used ingredient targets or disease genes as mediators to connecting with the proteins constructing the PPI network. The TCM syndrome related genes were used as seed nodes for the random walk algorithms, and the probability of restart was set to 0.7. Networks were visualized using Gephi 0.9.3 (https://gephi.org/).

The RandomWalkRestartMH R package, version 1.22.0, was utilized for key molecule mining based on RWR methodology and can be found at: [alberto-valdeolivas/RandomWalkRestartMH] (https://github.com/alberto-valdeolivas/RandomWalkRestartMH).

In our previous work (20), we comprehensively curated the differentially expressed genes (DEGs) related to 211 ingredients after the cell line or animal models were perturbed by these ingredients. In this work, we firstly identified key ingredients for each formula based on network analysis, and then we connectively mapped these key ingredients to drugs based on the DEGs from both HERB and CMap databases (36), with the latter provided pharmacotranscriptomics dataset for modern drugs. As a result, the ingredients with unknown functions can be functional inferred based on the drugs that were connectively mapped to them with scores > 90. A score of 90 indicates that only 10% of reference perturbations showed stronger connectivity to the query, thus is recommended as the default threshold by CMap. In addition, we also conducted data mining on the key small molecules from the website of ClinicalTrials.gov to prove the druggability of these natural compounds.

According to the 2019 edition of the “Guidelines for the Diagnosis and Treatment of Hepatocellular Carcinoma”, issued by the National Health Commission of the Chinese government (10, 11), five distinct types of TCM syndromes were published and treated using five different TCM formulae. These syndromes have been categorized as common TCM classifications for HCC. The composition of these five formulae is shown in Table 1 , with each formula corresponding to one of the five syndromes. The formulae involve a total of 44 herbs, and most of them appeared only once. And the information about their ingredients and targets is provided in the supplementary information ( Supplementary Appendix S1 , Supplementary Table S1 , Supplementary Figure S2 ).

TCM syndromes are subtypes of a particular disease with a cluster of symptoms that can be stably observed in long-term clinical practice. In the HCC disease, 5 TCM syndromes were curated by national experts. We thus compared the molecular similarities among the HCC disease and its 5 different TCM syndromes respectively, based on their associated genes in the large human PPI network ( Figure 2A ). We firstly retrieved 2,835 syndrome-associated genes for five TCM syndrome types ( Supplementary Appendix S1 , Supplementary Table S2 ). Then we set out to further quantify the distances of the gene sets associated with HCC disease and five TCM syndromes by network-based separation analysis. The separation measure is formed as S_DS , with smaller S_DS indicating a tighter network-based overlap between two gene sets. It is shown that the mean separation scores were fairly small (0.10, green box in Figure 2B ) when comparing the genes related to the HCC disease to genes for 5 TCM syndromes, indicating that the TCM syndromes of HCC are similar to the HCC disease in both phenotypes and molecular mechanisms (30). Furthermore, we also benchmarked this analysis by introducing random controls. To this end, we selected other 195 unrelated TCM syndrome types from SymMap database ( Supplementary Data S2 ). We also computed the separation measures between genes of HCC and genes for these 195 other syndromes. And their mean separation scores are significantly larger (0.24, orange box in Figure 2B ) than that for the comparison of HCC syndromes-related genes, indicating that diseases or TCM syndromes with similar clusters of symptoms will be closely related to adjacent network modules. To further demonstrate the robustness of the collected syndrome-related genes, we conducted the same calculations using different P-value thresholds. The results obtained were similar, to some extent confirming the reliability of the syndrome-related genes identified through SymMap mining in the context of traditional Chinese medicine ( Supplementary Appendix S1 , Supplementary Figure S3 ).

In addition to the topological approach, we further compared HCC and its related TCM syndromes based on the functional analysis of genes. We performed functional enrichment analysis on the HCC disease/HCC syndromes gene sets based on the KEGG annotation, and found 35 classes of enriched metabolic pathways for these genes ( Figure 2C ). And we showed that the associated pathway for the HCC disease overlapped with that for the five TCM syndrome types of HCC. Among them, syndrome 1 (syndrome of liver depression and spleen deficiency) demonstrated the highest degree of overlap with HCC. By integrating clinical data, some researchers have found that this syndrome is one of the most common syndromes of liver disease (37, 38). Meanwhile, the integrated application of intestinal microbiota and metabolomics research further provides experimental basis for the efficacy of TCM formulae in the treatment of hepatocellular carcinoma of the syndrome type (39). In addition, each syndrome overlapped to the HCC disease in a different set of pathways, indicating that the differentiation of a disease into several TCM syndrome subtypes captures the different aspects of the disease at the molecular level. The perspective promoted us to further untangle the molecular basis of different TCM formulae that were used to treat different HCC syndromes.

Different from modern drugs, TCM formula is consisting of multiple herbs, with each herb constituted by multiple ingredients. The mixture nature of TCM formula brought huge difficulties in studying its mechanism of action. Fortunately, with the introduction of systems biology as well as our previous study of the HERB database (20), TCM formulae can be connected to gene targets in recent years based on the global viewpoints of network medicine (40). We thus gathered potential gene targets for the 5 TCM formulae of HCC, and then compared them with the genes for HCC and its TCM syndromes in the large human PPI network. We used the network-based proximity measure to compute the closest distance between a gene set related to drugs/formulae and that for disease/TCM syndromes ( Figure 3A ), as the metric is shown to be suitable for quantifying treatment effect (32). Inspired by several studies using this metric in drug discovery and drug repurposing, we expected that effective drugs would demonstrate more close gene modules to the disease than irrelevant drugs (33, 41–43).

Similarly, we benchmarked the analyzing results for TCM formulae by introducing other drugs that are approved for treating several disease types, including HCC, liver disease, and cancers. We performed the network-based proximity analysis for several pairwise comparisons of gene sets. For example, TCM formulae vs. HCC, anti-HCC drugs vs. HCC, other anti-cancer drugs vs. HCC, other liver disease drugs vs. HCC, and drugs unrelated to either liver disease or cancer vs. HCC. For each comparison, we conducted permutation test to calculate the relative z-score values and P-values, and defined the statistically significant pairs when it meets a stringent criteria: z-score <−1.5 and P-value < 0.05 (33). As a result, we found that four first- and second-line anti-HCC drugs, named sorafenib (DB00398, z-score = -4.27), cabozantinib (DB08875, z-score = -3.28), regorafenib (DB08896, z-score = -3.79), and lenvatinib (DB09078, z-score = -4.18), all showed good proximity to HCC, which is consistent with their clinical indications (2). In contrast, other types of drugs demonstrated a fairly decreased ratio of significant proximal gene modules to HCC. As expected, the proximity analyzing results for five TCM formulae are all statistically significant ( Figure 3B ). Based on these findings, we expanded our scope of analysis to include direct comparisons with other database resources. This benchmarking process allowed us to further validate the effectiveness of our methodology. Different formula-target data from two reputable high-quality databases, including HIT 2.0 (44) and ITCM (45), provided a solid foundation for our network-based proximity analysis. We then proceeded to conduct a series of pairwise comparisons of the genomes, employing the same stringent statistical significance criteria previously mentioned. Our comparative analysis revealed that the results from different data sources were largely consistent overall. In terms of topological analysis ( Supplementary Appendix S1 , Supplementary Figure S4 , Supplementary Data S5 ), the targets of TCM formulas, as extracted from these databases, exhibited significant proximity to HCC genes, underscoring the inherent capacity of TCM to engage with the intricate pathophysiology of HCC. This topological congruence across different data sources not only validates the consistency of our results but also reaffirms the systemic approach inherent to TCM. The holistic nature of TCM is reflected in its ability to address HCC through the modulation of multiple biological pathways concurrently. This polypharmacological approach is particularly pertinent in the context of complex diseases like HCC, where a multi-target intervention strategy may offer superior therapeutic outcomes compared to single-target treatments. This result demonstrated the efficiency of TCM formulae in treating HCC, which aligns with the national guidelines for the treatment of HCC in China.

Next, we examined the proximity metrics between each TCM syndrome and its corresponding TCM formula ( Figure 3C ). It is showed that all five recommended syndrome-formula pairs (cells on the main diagonal in Figure 3C ) exhibiting significant proximity, further proved their clinical efficacy. Besides, among them, the syndrome 5 is exclusively associated with formula 5, implying that the dysregulated state of syndrome 5 is distinct from other syndromes in the biological network. A clinical study utilized peripheral blood mononuclear cells from HCC patients as test specimens and identified differences between the liver-kidney yin deficiency syndrome group and the non-liver-kidney yin deficiency syndrome group through bioinformatics analysis. The result further supports our findings (46). That is to say, the molecular analysis of the TCM clinical experiences across millennia can shed light on better precision medicine in the future.

We then performed KEGG enrichment analysis on gene sets related to HCC, its 5 TCM syndromes, drugs, and formulae ( Supplementary Appendix S1 , Supplementary Figure S5 ). We found that the gene targets for the TCM formulae for HCC shared several common KEGG pathways with that for anti-HCC drugs, including cancer: specific types, signal transduction, cancer: overview, endocrine system ( Figure 4A ). Compared with modern drugs, TCM formulae targeted larger number of pathways, not only includes cancer-related pathway but also embraces other pathways, such as viral infectious disease. This result is interpretable because HCC frequently originated from viral infections, including hepatitis B virus (HBV) and hepatitis C virus (HCV) (47). In contrast, small molecule drugs are more focused on cancer-specific pathways and signal transduction processes. These results are consistent with the target-based drug discovery mode of small molecule drugs (48) and phenotypic drug discovery for TCM formula (17). This consensus is further reinforced by the convergence of findings from diverse sources, which not only validates the robustness of our results but also underscores the potential of integrating systems biology with TCM ( Supplementary Appendix S1 , Supplementary Figure S5 ).

Furthermore, we analyzed the related pathways for HCC, and its TCM syndromes, formulae, and drugs in detail, based on 42 curated pathways mined from the PubMed database ( Figure 4B ). We then performed KEGG enrichment analysis for gene sets related to HCC, its TCM syndromes, drugs or TCM formulae, respectively. We also performed KEGG enrichment analysis for gene sets concurrently related to the disease-drug or the syndrome-formula pairs. The results demonstrated similar patterns as that in Figure 4A . Small molecule drugs showed fewer but more specific related pathways, while TCM formulae seem to act on more broad set of regulated pathways. We also noticed that different TCM formulae showed distinct major regulated pathways. For example, formula 1 regulates similar pathways with small molecule drugs, formula 2 mainly targets MAPK and RAS signaling pathways, formula 3 specifically targets purine metabolism and TGF-beta signaling pathways, formula 4 exclusively targets tyrosine metabolism, and formula 5 only targets ribosome-related processes. The results further provide functional evidence for understanding the usage of different treatment material for each detailed differentiated subtype of disease in a traditional TCM way.

Furthermore, we intended to dissect active small molecules in TCM formulae, in order to elucidate the key material basis from the mixture and provide support for modern drug discovery. To this end, we assembled a new heterogeneous network covering herbs, ingredients, and protein interactions by integrating the PPI network to the HERB target curations ( Figure 1D ). Then we navigated the network to explore its global topology by using the random walk with restart (RWR) algorithm (35). In a systematic effort to validate and demonstrate the advantages of our approach over traditional network pharmacology computational methods, we compared our drug screening algorithm with these methods. This systematic comparison was conducted to evaluate the efficacy and accuracy of the RWR algorithm in identifying the most relevant ingredients within the TCM formulae. The detailed results of the systematic comparison can be found in the supporting information ( Supplementary Appendix S1 , Supplementary Figure S6 ).

We used genes of the HCC-related TCM syndromes as the seed nodes to perform the RWR algorithm, and then measured the association degree of each ingredient in the heterogeneous network. Finally, the functions of top predicted ingredients were explored based on literature as well as connectivity mapping analysis. The 46 key ingredients were provided in Table 2 . The networks of these ingredients with the TCM formulae and syndromes were shown in Figure 5A .

Several ingredients were prioritized repeatedly in five formulae, such as citric acid, quercetin, coumarin, and curcumin. Among them, citrate is an important metabolic regulator, which is involved in inflammation and cancer (49, 50). It is shown that citrate uptake broadly affects tumor cell metabolism through citrate-dependent metabolic pathways (51), because citrate promotes lipid excess biosynthesis and disrupts lipid metabolism in cancer cells, resulting in cancer cell senescence and growth inhibition (52). Curcumin and quercetin not only have inhibitory effects on HCC (53, 54), but also have the potential to interfere with HCC comorbidities. Curcumin can also inhibit HCV entry of all major genotypes into primary human hepatocytes, thereby impairing viral binding and fusion (55). Quercetin also exhibited hepatoprotective effects on acute hepatitis by inhibiting the TRAF6/JNK pathway (56). To substantiate the reliability of our results about the small molecules screened, we also have conducted an extensive review of the most recent experimental literature, encompassing a variety of in vitro cellular assays and in vivo animal studies ( Supplementary Appendix S1 , Supplementary Table S3 ). This comprehensive analysis has provided a wealth of evidence supporting the efficacy of the compounds we have identified. These results demonstrated the power of our network medicine analysis for decoding promising material basis for TCM.

Some ingredients have evidence that they were absorbed into the blood or used for quality inspection based on SymMap annotation (17). Therefore, these ingredients are more likely to be the pharmacological substances of the formulae. So, we functionally annotated these ingredients based on their pharmacotranscriptomics datasets curated in the HERB database ( Figure 5B ). The results showed more evidence for their potential pharmacological effects. For example, nobiletin from Aurantii Fructus is a MEK inhibitor, which has been shown to significantly inhibit HCC in vitro and in vivo (57), by regulating cell proliferation, differentiation, angiogenesis, and survival (58). Note that the first-line drug, sorafenib, exerted its anti-HCC effects by targeting the RAF/MEK/ERK pathway (59). Besides, baicalein, a flavonoid from Scutellariae barbatae herba, may act as VEGFR inhibitor, FGFR inhibitor, and JAK inhibitor, which were significantly associated with anti-HCC effects (60). Note that VEGFR is the primary target of most FDA-approved HCC drugs such as sorafenib, lenvatinib, and regorafenib. We also utilized RDKit to visualize the structural similarities between nobiletin, baicalein, and sorafenib ( Supplementary Appendix S1 , Supplementary Figure S7 ). The result revealed that baicalein exhibits a higher degree of structural similarity to the core scaffold of sorafenib. This finding suggests that baicalein and sorafenib may share molecular fragments that are crucial for their pharmacological actions, potentially indicating a common mode of interaction with biological targets. Although nobiletin showed a lower overall similarity to sorafenib, the visualization provided by the similarity maps was invaluable. It shed light on the specific molecular features that contribute to the observed similarity, offering a roadmap for the discovery of new structural scaffolds with potential therapeutic relevance. Therefore, baicalein may exert its anti-angiogenic effect as a new VEGFR inhibitor in HCC. In conclusion, these active herbal ingredients exert diverse pharmacological activities in treating different TCM syndromes of HCC, and ultimately jointly intervene in the progression of HCC.