A Comparative Study of Systems Pharmacology and Gene Chip Technology for Predicting Targets of a Traditional Chinese Medicine Formula in Primary Liver Cancer Treatment

Background: The systems pharmacology approach is a target prediction model for traditional Chinese medicine and has been used increasingly in recent years. However, the accuracy of this model to other prediction models is yet to be established. Objective : To compare the systems pharmacology modelwithexperimental gene chip technology by using these models to predict targets of a traditional Chinese medicine formulain the treatment of primary liver cancer. Methods: Systems pharmacology and gene chip target predictions were performed for the traditional Chinese medicine formula ZhenzhuXiaojiTang (ZZXJT). A third square alignment was performed with molecular docking. Results: Identification of systems pharmacology accounted for 17% of targets, whilegene chip-predicted outcomes accounted for 19%.Molecular docking showed that the top ten targets (excludingcommon targets) of the system pharmacology model had better binding free energies than the gene chip model using twocommon targets as a benchmark. For both models, the core drugs predictions were more consistent than the core small molecules predictions. Conclusion:In this study, the identified targets of systems pharmacology weredissimilar to those identified by gene chip technology; whereas the core drug and small molecule predictions were similar.


INTRODUCTION
Traditional Chinese medicine (TCM) has a long history of clinical application in China andforms a complete and independent system of diagnosis and treatment. TCM is efficaciousin the treatment of many diseases, including cancer (Fang et al., 2017). In the treatment of cancer, TCMnot only shows great potential as a means of medical intervention (Yan et al., 2017) but may also be modified according to patient comorbidities. In TCM, sovereign drugs (drugs used to treat the primary disease) may be used together with adjuvant drugs (drugsused to improve the efficacy of the sovereign drugor totreat other symptoms), thus alleviating pain and improving quality of life. TCM usually follows the principle of multiple-herbformulasduring clinical application. In the treatment of a single disease, the core combination of drugs is difficult to change; however, adjuvant drugs may be added to increasethe efficacy of the drug combination or totreat secondary diseases. At Heilongjiang University of Chinese Medicine, our experimental group graduallyformed a unique drug regimen for the treatment of primary liver cancer, by diagnosing and treating primary liver cancer with a large collection of herbs, and then observing the clinical efficacy of these herbs in patients. Using this clinical experience, we finally identified five crucial herbs, includingLigustrum (NZZ), Curcumaerhizoma (EZ), Prunella vulgaris (XKC), Hedyotisdiffusa (BH), and Glycyrrhizae radix (GC). According to theory, TCM formulasconsist of sovereign, adjuvant, assistant, and guide drugs. We applied this principle and assignedNZZ and EZasthe sovereign drugs, XKC and BH asthe adjuvant drugs, and GC as the assistant and guide drug. These five drugs made upthe formula for treating primary liver cancer and we named itZhenzhuXiaojiTang (ZZXJT). We initially administered ZZXJT toa murineH22 hepatocarcinomamodel. TheZZXJT group showedH22 cell degeneration and necrosis and the number of blood vessels was reduced. Additionally, many autophagosomes in H22 cells were observed by transmission electron microscopy. Thesefindings revealed that ZZXJT may induce programmed H22 cell death and inhibit primary liver cancer development (Sun et al., 2017).
TCM uses a multi-component and multi-target approach (Pei et al., 2013), which can be both advantageous and disadvantageous in modern medical research. Advantages includesingle drugsexhibiting multiple mechanisms (Miao et al., 2020;El-Zayat et al., 2021), flexibility of multi-drug treatment (Ren et al., 2020;Zhang et al., 2021), and pairing of drugsto improve drug efficacy (Luo C. H. et al., 2020). However, the disadvantages are consequences of these advantages. A single drug comprisesnumerous compounds which in turn comprise different small molecules that vary in their pharmacological activity; hence, targets are difficult to identify in experimental studies. For example, HoupuDahuang decoction, HoupuSanwu decoction, and Xiaochengqi decoction include the herbsMangnolia officinalis, Rheum palmatum, and Citrus aurantium but the pharmacological effects and indications of the three formulas differ. HoupuDahuang decoction is mainly used to treat cough and exudative pleurisy, while the HoupuSanwu decoction is primarily used to treat paralytic ileus and Xiao Chengqi decoction is mainly used to treat adherent ileus, chronic gastritis, and intestinal paralysis (Kou et al., 2004). It is therefore challenging to regulate results and form a research strategy when core combinations of drugs are studied. Accurate TCM target prediction is crucial due to the rapid development of TCM ingredients in the post-genomic era. The systems pharmacology approach and gene chip technology provide effective target prediction methods for TCM laboratory research, thereby guiding research direction. These prediction models play important roles in drug development and research (Cooke et al., 2009;Boezio et al., 2017;Luo T. T. et al., 2020), as indicated by the increased number of studies in recent years on systems pharmacology particularly (Jiao et al., 2021). However, since the methods and principles of the two prediction models are very different, predicted results heavily influence the research direction.
Systems pharmacologymakes use ofa data platformthat predicts drug-target interactions based on network analyses from previously published research data (Berger et al., 2010). In gene chip technology, a large number of known gene sequences are immobilised and hybridised onto a glass chip, allowing for large-scale prediction (Yanagawa et al., 2005). Molecular docking is used as an auxiliary model for systems pharmacology and gene chip technology. Docking data corroborate the results of the prediction models to check reliability . Although systems pharmacology is cost-effective and saves time compared with the experimental high-throughput screening used in gene chip technology (Yuan et al., 2017), the differences in predictions between the two models have rarely been reported, especially in the research of TCM formulas. This study, therefore, aims toguide the future selection of methods for TCM target prediction.

Screening of Active Ingredients and Related Targets ofZZXJT
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 1 , to normalise it (UniProt Consortium, 2021). We then constructed the "drug active ingredient-target" network.

Screening for Liver Cancer Targets
Using "liver cancer" as the keyword, potential liver cancer-related targets were mined from the OMIM database 2 and GeneCards database 3 (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.

Construction and Analysis of "Drug Active Ingredient-Target" Network
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 4 . These intersecting targets were submitted to the STRING database 5 (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.

Gene Enrichment Analysis
We selected Metascape database 6 (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 7 .

Cell Culture
The cells were routinely cultured in 25 cm 2 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 2 .

RNA Extraction and RNA Quality Control
Six 25 cm 2 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.

Gene Chip Preparation and Hybridization
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. Singlestranded 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.

Construction and Analysis of the "Drug Active Ingredient-Target" Network
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.

Gene Enrichment Analysis
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.

qRT-PCR Assay
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 .

Identification Treatment of Two Models
The resulting sets of data from Section 2.1.3 and Section 2.2.3 were combined and a Venn diagram was drawn to determine intersection of predicted targets. The intersecting data were imported into the STRING database and the PPI network model was constructed by setting "Organism" to "Homo sapiens" and confidence level to 0.4. The resulting network was imported into Cytoscapesoftware (version 3.70) for analysis. The free proteins were identified by screening core potential proteins by adjusting the visualisation parameters according to the "Combined degree" values of individual nodes. We then performed KEGG pathway analysis using the Metascape database.

Main Active Ingredients of ZZXJT
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.

Selection of Docking Targets
We downloaded the protein crystal structures corresponding to the common target genes of the two models from the RCSB PDB database 8 , 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.

Molecular Docking Software and Protocol
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.

Identified Active Ingredients and Targets of ZZXJT
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.

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

Genes
Sequences

Construction of the ZZXJT-Liver Cancer Target Interaction PPI Network
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.

Enrichment Analysis of Target Pathways and Functions
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.

Chip RNA Quality Control Assessment and Results
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.

Construction of the ZZXJT-HCC Target PPI Network
A total of 118 nodes with 870 edges representing proteinprotein 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.

Enrichment Analysis of Target Pathways and Functions
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.

qRT-PCR to Detect Changes in the Level of Genes in Cells
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).

Identity Treatments for Two Models
Using Cytoscape software, a total of 25 nodes with 147 edge lines were obtained. Darker nodes and larger circles represent a greater degree of connectivity. Visualisation results showed that the top five targets, namely, CCND1, vascular endothelial growth factor A (VEGFA), MAPK1, EGF, and FOS, had the highest connectivity values. This indicated that these five targets were the most highly correlated in this model and can therefore be considered as the core five targets for ZZXJT action inHCC in this model. The intersecting targets were further analysed using the Metascape data platform, resulting in eight KEGG pathways. The distribution and significance of the proteins involved in each pathway were analysed. BP, CC, and MF analyses were not performed because of the small number of common targets and the lack of statistical significance of less than 20 pathways ( Figure 6).

Designationof Small Molecules for Docking
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.

Collection of Docking Targets
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).

Results of Docking
The heatmap (Figure 7) displays the docking resultspresented as free binding energy values (kcal·mol −1 ) 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 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).

DISCUSSION AND CONCLUSION
In this study, the targets of HCC were predictedusing a gene chip model and systems pharmacology model. The results of this paper have been analysed and discussed from the perspective of the contrasts and similarities between the targets predicted by thegene chip and systemspharmacology models, respectively.
It is evident from the results that the system pharmacology approach neglects the relationship between the activity of a drug's small molecules and the dosage of the drug. Numerous  relationshipsbetween small molecules andtargets may be predicted by the systems pharmacology approach, but this method may lead to errors since it simply uses superimposition based on past research results. This error may be augmented in formulation research, which focuses on the relationship between drug dosage, and drug compatibility. The results of this study showed that thecore targetspredicted by systems pharmacology and gene chip accounted for only 17   CCND1  2W9F  VEGFA  1VPF  EGF  SL0T  MAPK1  6G54  CCNA2  1OIY  STAT1  3WWT  PRKCA  2GZV  CAT  1QQW  MCL1  6OQC  RB1  7CZG  CASP8  4ZBW  FOS  1FOS  EDN1  1T7H  AHR  5Y7Y  CHEK1  2HOG  MET  3EFJ  TOP2A  5NNE  SERPINE1  3PB1  CASP9  3YGS  CAV1  7LUD  ERBB3  3KEX  SREBF1  1AM9  FASN  2JFD HIF1A 4H6J  and 19% of targets, respectively. The similaritiesintargets identified by the two models differed. The top 10 core targets identified by the system pharmacology model only accounted for 2 common targets. Counting the top 10 targets without including the common targets would increase this to 13 targets. The common targets in the top 10 core targets predicted by the gene chip model accounted for 5 targets.
Counting the top 10 targets without including the common targets would increase this to 22 targets. Itcan therefore be concluded that the two models were poorly similar based on the target aspect. In terms of KEGG pathway analysis, the top identicalpathwayof the two models was pathways in cancer and the top fiveincluded the PI3K-AKT signalling pathways. While there were common KEGG pathways, except for the inclusion of pathways in cancer, the two modelsdid not share the same pathways, and so the use of KEGG pathway analysis was limited. In conclusion, the gene chip modelis more accurate in target predictions than the systems pharmacology model (Heber and Sick, 2006). Thismay be a direct result of gene chips having been developed as an assay for genetic diseases and as an alternative to clinical disease testing (Huang et al., 2004;Shi et al., 2021). In addition, gene chips usually employ normalisation, which directly expands the target genes that may be involved. In comparison, the predicted results of the systems pharmacology model in this study were poor. Molecular docking was introduced in this comparative study as a reference prediction model. The binding free energy was used as a benchmark and the top 10 targets from the two models (without co-targets) were individually aligned. Results showed that the docking scores of system pharmacology and co-targets were not statistically significant (p > 0.05), while those of genechipsand co-targets were statistically significant (p < 0.05). The mean values were larger than the mean values of the common targets. From the docking data of 10 small molecules docked with two model core targets, it can be derived that, because the numerical value is smaller, and the binding is more stable (Hsin et al., 2013). Based on the docking results, system pharmacology was better than the binding of the gene chip. However, experimental results are insufficient to fully reject the prediction that molecular docking is more accurate, possibly due to too little docking data, not only in actual drug effect results but also not in full agreement with binding energy values (Ye et al., 2021). For example, in the vicinity of the target tissue, the lowest drug concentration isa factor that affects the regulation of gene expression (Shin et al., 2020), or when the same drug molecule acts on targets of different cells, and results should differ (Akita and Sliwkowski, 2003). These situations are currently difficult to simulate using a computer model. Although drug absorption was previously modelled through the specification of OB and DL values, such predictions by systems pharmacologymay result in extremely low confidence for a sophisticated and complex network of target structures, especially after multiple error factors are incorporated. Therefore, systems pharmacology corroborated by molecular docking is risky as results may direct researchers incorrectly. Nevertheless, based on the docking data, the predictions of the core small molecules from the two models were consistent, mainly focusing on five small molecules, namely, A1, A2, B1, C1, and NZZ4. Among these, the four consensus small molecules are part of the drug XKC while the remaining four small molecules and A2 are part of NZZ. Based on the concept of"sovereign and adjuvant" in TCM, it is evident that NZZ and XKC play a major role in ZZXJT. Thissupported the TCM theory that adjuvant drugs may also besovereigndrugs. It was inferred that system pharmacology and molecular docking models were of positive significance to identify the main drugs and molecular components in the formula. This was highly consistent with the gene chip.
Systems pharmacology predicts that the targets of the core small molecules are the same, therefore the main treatments should be the same. This prediction in terms of TCM, may wrongly direct researchas important drug-target relationships are circumvented. The integrity of formula research is affected by the simple superposition of a single molecule, which may lead to the low credibility of predicted results. This assumption is made based on systems pharmacology for ZZXJT target prediction of a single formula (ZZXJT), hencemore data is required to firmly establish the reliability of using systems pharmacology predictions in TCM formula research. This can be addressed in future studies.
Since predictions by the systems pharmacology approachare poorly similarto the gene chip technology predictions, we can therefore conclude thatsystems pharmacology has low model confidence and is less reliable. However, the consistency between the core drug predictions versus the core smallmolecule predictions was greater in the systems pharmacology model.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories, which can be found in the article and Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by Heilongjiang University of Chinese Medicine Institutional Animal Care and Use Committee (IACUC).

AUTHOR CONTRIBUTIONS
YaS and SL contributed conception of the study. YaS provided financial support. SL and YuS collected the dataset. All authors wrote sections of the manuscript, contributed the manuscript revision, and approved the submitted version.

FUNDING
This study was supported by the National Nature Science Foundation of China (Grant No. 81704054) and Postdoctoral Scientific Research Developmental Fund (Grant No. LBH-Q19184).