We searched for the active ingredients of thefive herbs, EZ, NZZ, BH, XKC, and GConthe traditional Chinesemedicine systems pharmacology database andanalysis platform (TCMSP) (Ru et al., 2014). We screenedthese active ingredients by specifyingtwo ADME-related properties as the screening criteria, namely, oral bioavailability (OB) ≥ 30% and drug-likeness (DL) ≥ 0.18. Active ingredients without targets were removed and we integratedthe identified active ingredient targets. After determining the target protein information, the collected target information was unified in UniProt protein database ¹ , to normalise it (UniProt Consortium, 2021). We then constructed the “drug active ingredient-target” network.

Using “liver cancer” as the keyword, potential liver cancer-related targets were mined from the OMIM database ² and GeneCards database ³ (downloaded December 2020). Liver cancer-related genesidentified by the GeneCards database were filtered andonly results with relevance score ≥15 wereretained. Data collected from the two databases weremerged and denoted as Dataset 1.

To further clarify the mechanism and relevance of the interaction between ZZXJT and hepatocellular carcinoma (HCC), intersecting targets between ZZXJT (drug), and HCC (disease) were obtained by drawing Venn diagrams ⁴ . These intersecting targets were submitted to the STRING database ⁵ (Szklarczyk et al., 2019) and a protein-protein interaction (PPI) network was constructed by setting “Organism” to “Homo sapiens” and confidence level to 0.4. This network was imported to Cytoscape software (version 3.70) and core potential proteins were analysed by adjustingvisualisation parameters according to the “combined degree” values of individual nodes. The functions of the target proteins involved in the biological process were described.

We selected Metascape database ⁶ (Tripathi et al., 2015; Zhou et al., 2019) as the gene-list analysis portal because it is updated regularly and is data comprehensive. The relevant ZZXJTand HCC targets were entered into the Metascape database and screening criteria were set as follows: statistical difference at p < 0.01 and mode of analysisas custom analysis. The following analyses were then carried out: Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO). GO has three subdivisions that were conducted, namely, molecularfunction (MF), biologicalprocess (BP), and cellularcomponent (CC). KEGG pathway analysis yielded the top 20 results and BP, CC, and MF yielded the top 10 results. Results were displayed in bubble plots, created using ImageGP ⁷ .

The cells were routinely cultured in 25 cm² culture flask in DMEM medium (Boster, #PYG0103) supplemented with 5% fetal bovine serum and 1% antibiotics (penicillin/streptomycin). When the cells placed under the culture flask about 70–80%, the original culture medium was discarded, and the culture flask was rinse with PBS gently. Then trypsinization was added in the flask. When the intercellular space of cells enlarged, fresh culture medium was put in the flask to finish trypsin digestion. After centrifugation, the culture medium was extracted carefully on the cell upper layer, and then the cells were made a single cell suspension by the fresh culture medium. The cells were cultured in a constant temperature (37°C) with a volume fraction 5% CO₂.

Six 25 cm² cell culture flasks of SMMC-7721 cells were cultured together for 24 h. Three flasks were used as the control group (sample numbersA12854, A12855, and A12856) and the remainder as the drug group, in which 20% ZZXJT-containing serum was added to each flask (sample numbersA12857, A12858, and A12859). RNA was extracted from samples using the TRIzol reagent (Invitrogen, United States), according to the manufacturer’s instructions. We measured the A260/A280 ratio using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, United States). The required reagents were then left at room temperature for 30 minand the sample and RNA ladder were placed on ice. A 550 µl RNA 6000 Nano gel matrix (Agilent Technologies, United States) was placed in a centrifuge tube and centrifuged at 1500 g for 10 min at room temperature. We then added 1 µl of dye to 65 µl of the gel, shook it well, and centrifugedit at 13000 g for 10 min at room temperature to make the gel-dye mix. Next, we added 9 µl of the gel-dye mix to the GN-genechipClariom™ S Array, human gene chip (catalogue number 902927, Affymetrix, United States) by pressingthe gel-dye mix into the gene chip capillary with a piston. 5 µl of RNA 6000 Nano maker (Agilent Technologies, United States)was added to the sample well and ladder well. Finally, we added 1 µl of denatured ladder into the Agilent 2,100 instrument (Agilent Technologies, United States) to assess the RNA integrity number (RIN) and the 28S/18S ratio. Data were analysed using the Agilent 2,100 Expert software. Quality control standards were as follows: A260/A280 ratio of 1.7–2.2, RIN≥7.0, and 28S/18S > 0.7.

The first and second strands of the cDNA were synthesised. Labelling cRNA was synthesised by in vitro transcription. The synthesised cRNAs were then purified and quantified. Single-stranded cDNA was synthesised in the second cycle. cRNA was hydrolyzed by RNase H and single-stranded cDNA remained. After the second cycle of single-stranded cDNA was purified, its concentration was measured. The purified ss-cDNA was transformed into dUTP residual fragments and broken DNA strands, which were covalently linked to biotin, Affymetrix proprietary DNA labelling reagent, to complete cDNA fragmentation and labelling, using. The gene chip was removed and hybridisation and washing were performed.

The obtained gene chip microarray data were combined with Dataset 1, and a Venn diagram was drawn to determine the intersectinggene chip-predicted targets and HCC targets. The intersecting targets were submitted to the STRING database and the PPI network model was constructed by setting“Organism” to “Homo sapiens” and confidence level to 0.4. This network was imported to Cytoscape software (version 3.70) and core potential proteins were analysed by adjusting visualisation parameters according to the “Combined degree” values of individual nodes. The functions of the target proteins involved in the biological process were described.

The gene chip-predicted targets of ZZXJT and relevant targets of HCC were entered into the Metascape database. Screening criteria were set as follows: statistical difference as p < 0.01 and the analysis mode as custom analysis. KEGG pathway analysis and GO analysis at MF, BP, and CC levels, among others, were then performed.

The total RNA of drug group and control group was extracted with Trizol (Invitrogen, United States). Use Prime Script first Stand cDNA Synthesis Kit to reverse transcribe RNA into cDNA according to the instructions. Adopt SYBR pre-mixed ExTaq kit (Takara, Dalian, China), PCR amplification was conducted at 40 cycles of denaturation at 95°C for 30s, after pre-denaturation treatment at 95°C for 5 min, annealing treatment 58°C for 30 s and extension at 72°C for 30 s. The temperature at the end was 4°C. GAPDH was used as the internal reference gene to detect the relative expression levels of AKT, FOXO1, FOXO3, and FOXO4. The primers of these genes were displayed in Table 1. Each sample was repeated 3 times to ensure the accuracy of the data. The relative expression levels of these genes were calculated in accordance with 2^−∆∆ct.

According to the TCM theory of sovereignand adjuvant drugs, the active ingredients of sovereigndrugs in Table 2 and the active ingredients shared by Ligustrum (NZZ) andCurcumaerhizoma (EZ) were integrated. The mol2 file corresponding to these active ingredients was downloaded from the TCMSP database for later use.

We downloaded the protein crystal structures corresponding to the common target genes of the two models from the RCSB PDB database ⁸ , using the search criteria “Homo sapiens” and “protein”Thetarget gene proteins were then ranked by “Combined degree” values of the two sets of models. Crystal structures of the top tenwere selected (excluding the common target proteins) and the PDB numbers were recorded.

AutoDock Vina (version 1.5.6) was the software used for molecular docking (Trott and Olson, 2010). According to the AutoDock report, 78% of molecular target docking results had aroot mean square deviation (RMSD) < 2. RMSD <2 was considered feasible, therefore, molecular docking predictions conducted using AutoDock could be considered accurate. The collected small molecules and target proteins of ZZXJT were prepared by removing water molecules, adding hydrogen atoms, and setting semiflexible docking and blind docking methods as parameters. The remaining parameters were set to default values. Given that this study already covered experimental microarray prediction, molecular docking was used as a third-party reference model. The possibility of interfering conditions such as hydrogen bonds was therefore ignored andonly maximum binding energy was taken as a reference. The free energy results of the system pharmacology model, the gene chip model, and the common target gene proteins were then statistically aligned using SPSS 26.0 and the independent samples t-test method was used.

A total of 126 chemical components were initially extracted from ZZXJT and 103 active components were identified after ADME screening and normalisation of the UniProt protein database, including poriferosterol, glycyrol, and hederagenin. As shown in Table 2, these small molecules, numbered A1 (BH, GC, XKC, and NZZ), A2 (BH, XKC), A3 (BH, XKC, and NZZ), B1 (GC, XKC, and NZZ), and C1 (XKC,NZZ) are common to many drugs. After merging and removing duplicates, the number of ingredient targets was 235. Figure 1 shows the “drug active ingredient-target” network constructed usingCytoscape software.

HCC-associated targets from the OMIM database and GeneCards database were merged, resulting in a total of 1,067 targets (after duplicates were removed).

The PPI network from the systems pharmacology modelhas a total of 145 nodes (Figure 2A), which represent predicted targets and 3,167 edges which represent protein-protein interaction relationships (Figure 2B). Darker nodes and larger circles represent greater degree of connectivity. Visualisation revealed that five targets, namely, phosphorylated protein kinase (AKT1), tumour suppressor gene (TP53), interleukin (IL)—6, mitogen activated protein kinase 3 (MAPK3), and vascular endothelial growth factor A (VEGFA) had the highest connectivity scores. This indicates that these five targets were most strongly associated with this model and can be considered as the core five targets of ZZXJTin HCC in the systems pharmacology model.

The enrichment analysis of ZZXJT in relation to HCC targets was conducted using the Metascape data platform, which resulted in 462 KEGG pathways, 755 GO molecular functions, 5670 GO biological processes, and 426 GO cellular components. Only the top 20 results from each group were considered. In the KEGG pathway analysis (p < 0.01), pathways related to HCC were PI3K-AKT signalling pathway, apoptosis, and transcriptional misregulation in cancer. Figures 2C,D show the distribution and significance of the proteins involved in each pathway. The exported results for BP, CC, and MF are presented in Figure 2E. BP is mainly involved in the response to toxic substances, response to inorganic substances, cellular response to lipids, response to lipopolysaccharide, and apoptotic signalling pathways. CC is mainly involvedin membrane rafts, vesicle lumens, receptor complexes, protein kinase complexes, and RNA polymerase II transcription factor complexes. MF is mainly involved in transcription factor binding, protein kinase binding, nuclear receptor activity, protein kinase activity, and protein homodimerization activity.

RNA quality control information of the normal and drug groups was examined using Nanodrop 2000 and Agilent 2,100 (Table 3). Raw data for this experiment were acquired using the GeneChip Scanner 3,000 (Affymetrix, United States) (Figure 3). The signal intensity profiles of the 6 sample sets were calculated. It is assumed that the plots of the signal intensity profiles can demonstrate the signal intensity profiles of all chip probes. As shown in Figure 4A, the abscissa represents the probe signal intensity interval, and the ordinate represents the number of probe sets within the signal intensity interval. The better the sample coincidence of the signal intensity distribution curves, the greater the reliability of the gene chip model. Data was taken for normalisation, resulting in a total of 1,543 differential genes (| fold change | ≥ 2.0 and FDR <0.05), of which 194 were upregulated genes and 1,349 were downregulated genes. The quality control requirements for RNA were metfor subsequent experiments to be conducted and the results of the gene chip analysis were highly reliable and reproducible.

A total of 118 nodes with 870 edges representing protein-protein interaction relationships were obtained using Cytoscape software. Darker nodes and larger circles represent a greater degree of connectivity. After visualisation, the greatest connectivity values were observed for the following five targets: vascular endothelial growth factor A (VEGFA), EGF, MAPK1, CCND1, and CYCS. These five targets were most strongly associatedwith this model and can therefore be considered as the core five targets of ZZXJT activityin HCC, as shown in Figures 4B,C.

The enrichment analysis of ZZXJT and HCC-related targets was conducted using the Metascape platform, which resulted in 461 KEGG pathways, 686 GO molecular functions, 5114 GO biological processes, and 484 GOcellular components. Only the top 20 resultsdisplayedby KEGG were considered. In KEGG pathway analysis (p < 0.01), pathways related to cancer were the PI3K-AKT pathway, FOXO pathway, p53 pathway, and viral carcinogenesis. Figures 4D,E show the distribution and significance of the proteins involved in each pathway. The export results for BP, CC, and MF are presented in Figure 4F. Each group shows only the top 10 results. BP is mainly involved in blood vessel development, cellular response to growth factor stimulation, glandular development, response to inorganic substances, and positive regulation of the cell cycle. CCmainly involves the DNA recombinase mediator complex, protein kinase complex, membrane raft, mitochondrial envelope, and cell base. MF mainly involves DNA secondary structure binding, phosphotransferase activity, alcohol group as acceptor, protein heterodimerization activity, single-stranded DNA binding, and kinase binding.

In KEGG pathway analysis, pathways related to cancer were the PI3K-AKT pathway and FOXO pathway. As shown in Figure 5, AKT mRNA expression downregulated compared with the control group (p < 0.05). FOXO pathway related genes (FOXO1, FOXO3 and FOXO4) mRNA expression gradually upregulated (p < 0.05).

A total of 10 small molecules were involved in docking, among which the smallmolecules of the sovereign drug (small molecules common to the drugs were not counted), and the other drugs were numbered EZ1, NZZ1, NZZ2, NZZ3, and NZZ4. The five multidrug common small molecules were numbered A1, A2, A3, B1, and C1.

The two models had 25 common targets. One targetprotein structure was removed as it was not searched. The final total of common targets was 24, corresponding to the PDB numbers shown in Table 4. There were three common targets among the top10 target genes of the systems pharmacology model according to “Combined degree”, therefore the targets orderedfrom 1 to 13(PDB numbersshown in Table 5). The top 10 target genes predicted by the gene chip model according to “Combined degree”, included 13 common targets, 2 proteins without crystal structures, and 1 without Homo sapiens protein crystal structure, therefore the number of target genes was ordered from 1 to 26(PDB numbers shown in Table 6).

The heatmap (Figure 7) displays the docking resultspresented as free binding energy values (kcal·mol⁻¹) ranging from-10 to 0. The 10 small molecules were individually docked with the target protein crystals, and 440 docking results were obtained, comprising 240 for the common target gene proteins and 100 each for the systems pharmacology model and the gene chip model, respectively. VMD software was used to visualise the docked conformation stack plots according to the obtained data to visualise the three optimal strips in each group. The action profiles and visualisations of candidate compounds with each target gene protein are shown in Figure 7 the common target pairsare shown in Figures 7A,a–c, the systems pharmacology target pairs are shown in Figures 7B,d–f, and the gene chip docking targets are shown in Figures 7C,g–i. The three sets of binding free energy results were integrated for statistical comparison, and it was found that there was no statistical difference (p > 0.05) in binding free energy for the common target gene proteins (mean = −7.05)and the systems pharmacology model (mean = −7.04). There was a significant statistical difference (p < 0.05) in the binding free energy between the common target gene proteins (mean = −7.05) and the genechip model (mean = −6.74).