Comprehensive Computational Analysis of Honokiol Targets for Cell Cycle Inhibition and Immunotherapy in Metastatic Breast Cancer Stem Cells

Breast cancer stem cells (BCSCs) play a critical role in chemoresistance, metastasis, and poor prognosis of breast cancer. BCSCs are mostly dormant, and therefore, activating them and modulating the cell cycle are important for successful therapy against BCSCs. The tumor microenvironment (TME) promotes BCSC survival and cancer progression, and targeting the TME can aid in successful immunotherapy. Honokiol (HNK), a bioactive polyphenol isolated from the bark and seed pods of Magnolia spp., is known to exert anticancer effects, such as inducing cell cycle arrest, inhibiting metastasis, and overcoming immunotherapy resistance in breast cancer cells. However, the molecular mechanisms of action of HNK in BCSCs, as well as its effects on the cell cycle, remain unclear. This study aimed to explore the potential targets and molecular mechanisms of HNK on metastatic BCSC (mBCSC)-cell cycle arrest and the impact of the TME. Using bioinformatics analyses, we predicted HNK protein targets from several databases and retrieved the genes differentially expressed in mBCSCs from the GEO database. The intersection between the differentially expressed genes (DEGs) and the HNK-targets was determined using a Venn diagram, and the results were analyzed using a protein-protein interaction network, hub gene selection, gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses, genetic alteration analysis, survival rate, and immune cell infiltration levels. Finally, the interaction between HNK and two HNK-targets regulating the cell cycle was analyzed using molecular docking analysis. The identified potential therapeutic targets of HNK (PTTH) included CCND1, SIRT2, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1, which can potentially inhibit the cell cycle of mBCSCs. Moreover, our results showed that PTTH could modulate the PI3K/Akt/mTOR and HIF1/NFkB/pathways. Overall, these findings highlight the potential of HNK as an immunotherapeutic agent for mBCSCs by modulating the tumor immune environment.


Introduction
Breast cancer was the most prevalent cancer in 2020 (in terms of new cases) and the leading cause of cancer-related deaths among females [1]. According to the World Health Organization, breast cancer has the highest incidence rate in Indonesia, with a mortality rate of 22,692 cases per year [1]. By 2040, the incidence is predicted to reach 89,512 cases [1]. Chemotherapy, along with surgery, radiation, and mastectomy, is the most common treatment [2]. Chemoresistance, or the insensitivity of cancer cells to drug therapy, is a major factor in the failure of chemotherapy against breast cancer.
Breast cancer stem cells (BCSCs) are one of the main factors driving chemoresistance, thereby contributing to poor prognosis and clinical outcomes [3][4][5]. BCSCs can develop into many cell types and repopulate heterogeneous tumors following conventional chemotherapy or radiotherapy [4,6]. BCSCs are mostly dormant and therefore activating dormant cells, and modulating the cell cycle is important for achieving successful BCSCs therapy [7]. Recurrent tumors are highly aggressive, potentially cross-drug resistant, highly metastatic, and have a poor prognosis. A previous study demonstrated that immune cells such as CD8+ lymphocytes induce epithelial to mesenchymal transition of BCSCs [8] Moreover, the tumor microenvironment (TME) promotes BCSC survival and cancer progression [9], and hence it can prevent the success of immunotherapy [10] e use of combination therapy, in which both chemotherapy and natural compounds are used to target metastatic BCSCs (mBCSCs), could be a successful approach to overcome chemoresistance and achieve clinical success in treating breast cancer.
is study aimed to explore the molecular mechanisms underlying HNK-mediated mBCSC-cell cycle arrest, as well as to assess the impact of this compound on the immune environment using bioinformatics studies.

Data Collection and Differentially Expressed Genes'
(DEGs) Identification. Proteins that interact with HNK were searched using STITCH (https://stitch.embl.de), [23]. Swisstargetprediction (https://www.swisstargetprediction. ch), [24] canSAR Black (https://cansarblack.icr.ac.uk/) [25], and SEA (https://sea.bkslab.org/) [26]. e retrieved proteins were considered as HNK-mediated proteins (HMPs) and were included in the subsequent analyses. e microarray data of metastatic breast cancer stem cells were collected from the GEO database (https://www.ncbi.nlm. nih.gov/geo) using keywords such as metastatic breast cancer stem cells, and Homo sapiens. e inclusion criteria were: use of patient samples or patient-derived xenografts; focus on metastatic breast cancer; characterization of breast cancer stem cells; and clear description of the identity of samples in the GSE datasets. e exclusion criteria included: use of breast cancer cell lines; no emphasis on metastatic breast cancer; no characterization of breast cancer stem cells; and ambiguity around the identity of samples in the GSE datasets. One GSE Dataset (GSE151191) was selected among the 62 datasets for this study (Supplementary Figure 1). GEO2R, a web-based interactive program (https://www. ncbi.nlm.nih.gov/geo/geo2r) that compares two groups of samples under the same conditions, was used to identify the DEGs between primary and metastatic tumors, on the basis of the following criteria for significance: P < 0.05 and log | Fold Change| > 1. Using [27], the overlapping proteins between the HMP and those encoded by the DEGs were identified and further analyzed using a protein-protein interaction (PPI) network.

Construction of the PPI Network.
e PPI was constructed and displayed using STRING-DB v11.0 and Cytoscape software, respectively [28,29]. Proteins included in the top-10 rank according to the Maximal Clique Centrality (MCC) score determined by the Cyto-Hubba plugin were considered hub genes [30]. e hub genes were subjected to GO and KEGG enrichment analyses using the tools [31] and WebGestalt [32]. Statistical significance was set at P < 0.05.

Genetic Alterations Analysis.
e proteins nominated by GO and KEGG enrichment analyses and an in-depth literature study on the hub genes were used to determine the potential therapeutic targets of HNK (PTTH). In this study, genes encoding PTTH such as CCND1, SIRT2, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1 were screened for genetic changes in all breast cancer studies available in the cBio-portal database ((https://www. cbioportal.org) [33]. e studies with the highest number of genetic changes were selected for analysis further connectivity.

Survival Rate and Immune Cell Infiltration Level.
e online database Gene Expression Profiling Interactive Analysis (GEPIA, https://gepia.cancer-pku.cn) was utilized to analyze the contribution of PTTH to the overall survival (OS) [34]. Tumor Immune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/) was used to analyze the correlation between PTTH expression levels and immune cell infiltration level [35].

Validation of the mRNA and Protein Expression Levels of PTTH.
e mRNA and protein expression levels of PTTH were determined using TNMPlot and Human Protein Atlas Evidence-Based Complementary and Alternative Medicine (HPA). Differentially expressed genes and mRNA levels in tumor, normal, and metastatic tissues were analyzed using TNMplot (https://www.tnmplot.com/) [36], in which the database utilized data from GEO or RNA-seq libraries from e Cancer Genome Atlas (TCGA), erapeutically Applicable Research to Generate Effective Treatments (TAR-GET), and e Genotype-Tissue Expression (GTEx). e protein expression levels were analyzed using the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) [37], an online database that contains a wide range of transcriptomic and proteomic data from various tissues and cells.

Molecular Docking.
To predict the binding properties of HNK to AURKB and RAC-1 through molecular docking, computational prediction was conducted on a Windows 10 operating system, Intel Core (TM) i5-10th Gen with 8 GB of RAM. MOE 2010 (licensed from Faculty of Pharmacy UGM) was used for docking simulation, RMSD-docking score calculation, and visualization interaction. e PDB IDs of the proteins AURKB and RAC-1 (3ZCW and 3TH5, respectively) were searched for in https://rcsb.org. e HNK structure was obtained from PubChem, subjected to conformational search, and minimized in the MOE using Energy Minimize Menu. For the docking simulation setting, London dG was used for both Rescoring 1 and Rescoring 2. Triangle Matcher was used for score function and placement setting, and Forcefield was used to refine the docking results from 30 retained settings. e results of this method will determine which conformation has the lowest binding interaction between the ligand and its receptor.

DEG and HMP Identification.
e DEGs are considered to be the molecular drivers and/or molecular biomarkers of Evidence-Based Complementary and Alternative Medicine 3 various phenotypes [38]. e identification of DEGs was carried out to determine the genes that act as biomolecular markers of metastatic breast cancer stem cells (mBCSCs). In total, 6,970 DEGs in the GSE 151191 dataset were found to be up/downregulated in metastatic breast cancer stem cells, according to the adjusted P value of <0.05, and a |logFC| ≥ 1.0 (Supplementary Table 1). Subsequently, proteins that interact directly and/or indirectly with HNK, referred to as HMPs, were identified. A total of 128 HMPs were retrieved from Swisstargetprediction, STITCH, canSAR Black, and SEA (Supplementary Table 2). Finally, 30 overlapping genes (OGs), including 18 upregulated and 12 downregulated genes, were identified to be both HMP and DEGs (Figure 1(b); Table 1).

GO and KEGG Pathway Enrichment Analysis.
e functions of the OGs were further investigated using GO and KEGG pathway enrichment analyses. e biological processes in which these OGs were implicated are summarized in Figure 2(a). Among the identified biological processes, cell communication, metabolic process, cellular component organization, multicellular organismal process, developmental process, response to stimulus, and biological regulation are strongly linked to cancer progression. According to the enrichment analysis of cellular components, the OGs were abundant in the nucleus, cytosol, membrane-enclosed lumen, and protein-containing complex (Figure 2(a)). Finally, the OGs were enriched in the molecular function protein binding (Figure 2(a)). KEGG pathway enrichment analysis demonstrated that the OGs were particularly enriched in the PI3K-Akt signaling pathway, pathways associated with cancer, and the regulation of the cell cycle (Supplementary Table 3).

Genetics Alteration Analysis.
Eight genes (SIRT2, CCND1, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1) that play an essential role in the growth and development of mBCSC were selected from hub genes and referred to as potential therapeutic targets of honokiol (PTTH). Genetic variation within these genes was analyzed using cBioportal. e breast cancer study with the highest number of genetic changes was selected for further analysis (Figure 2(b)). Oncoprint was used to determine the percentage of PTTH gene alterations in patients with mBC. Genetic alterations in PTTH ranged from 1.1% to 35% in the 180 patient samples analyzed (Figure 2(c)), with amplification being the most common gene alteration. e genes that were most frequently mutated were CCND1 (35%), HSP90AB1 (11%), and VEGFA (8%). Mutual exclusivity analysis showed that VEGFA mutations significantly cooccurred with HSP90AB1 mutations (Table 3). Copy number alterations (CNAs) are particularly common in cancer and play a significant role in its development and progression. CNA status can be homozygously deleted (shallow deletion), heterozygously deleted (deep deletion), diploid, gained (amplification event with relatively few copies), or amplified (amplification event with many copies). CNAs analysis revealed that SIRT2 mRNA expression was lower in shallow deletion cases and higher in amplification cases than in diploids (normal/without change) (Figure 2(d)). HDAC1 and HSP90AB1 mRNA expression was lower in cases with shallow deletions and higher in cases with gain. HSP90AA1 mRNA expression was lower in patients with gain than in those with diploid gain. CASP9 mRNA expression was lower in the gain than in the in shallow deletion cases, but not significantly different from that in diploid cases. All CNAs other than those mentioned were not differently expressed. Finally, the cBioportal pathway analysis showed that the cell cycle pathway is the main pathway that is disrupted by PTTH genetic alterations. Among the genes involved in the regulation of cell pathway, CCND1, encoding cyclin D1, was identified as PTTH ( Figure 2(e)).

Survival Rate and Immune Cell Infiltration Level.
To assess the clinical value of PTTH genes' expression levels, we examined whether they are associated with the OS or prognosis of patients with breast cancer. Low expression levels of CASP9 and HSP90AB1 were significantly associated with poor OS (P < 0.05) (Figure 3(a)). To understand the role of the immune microenvironment in the development and prognosis of patients with BRCA mutations, we analyzed the correlation between the expression levels of PTTH and immunocyte infiltration. e expression level of PTTH was either positively or negatively related to the infiltration level of different immune cells, indicating that PTTH modulated the immunologic microenvironment by influencing immune cell infiltration. e expression levels of CCND1, VEGFA, AURKB, HDAC1, HSP90AA1, and HSP90AB1 were positively correlated with immune cell infiltration levels, whereas the expression levels of SIRT2 and CASP9 were negatively correlated with the tumor purity of BRCA (Figure 3(b)). Additionally, the B cells' infiltration level was positively correlated with the expression levels of HSP90AB1 and AURKB, and negatively correlated with the expression level of CCND1. Moreover, we observed a positive correlation between CD8+ T cells' infiltration levels and SIRT2, HDAC1, CASP9, and HSP90AA1 expression levels, and between CD4+ T cells' infiltration levels and SIRT2, HDAC1, and CASP9 expression levels. However, CD4+ T cells' infiltration levels were negatively correlated with HSP90AA1 expression levels. Infiltration levels of macrophages were positively correlated with CCND1, SIRT2, CASP9, and HSP90AA1 expression levels, whereas they were negatively correlated with AURKB expression levels. Neutrophil infiltration levels were positively correlated with SIRT2, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1. Dendritic cell infiltration levels were negatively correlated with CCND1 expression levels and positively correlated with SIRT2, AURKB, and HDAC1 expression levels. Other non-mentioned data are not statistically significant.

Validation of the mRNA Expression Level and the Protein Expression Level of PTTH.
e expression levels of CCND1, AURKB, HDAC1, VEGFA, HSP90AA1, and HSP90AB1 were increased in the BC tissue, and even higher in mBC (Figure 4(a)). CASP9 expression levels were not significantly different between normal, tumor, and metastatic breast cancer tissues. Interestingly, SIRT2 expression was decreased in breast cancer tissues, but increased in metastatic tissues compared to normal tissues. ese results were supported by immunohistochemical data from HPA, that showed that CCND1, AURKB, and HDAC1 were overexpressed in the nucleus, while SIRT2, VEGFA, HSP90AA1, and HSP90AB1 were overexpressed in the cytoplasmic/membranous region (Figure 4(b)). Finally, CASP9 was not differentially expressed between the normal tissue and the tumor tissue.

Molecular
Docking. Molecular docking analysis revealed that AURKB and RAC-1 could bind to their respective native ligands and to HNK ( Figure 5). e affinity of the interaction between these proteins and HNK was similar to that of their natural ligands. e interaction between AURKB and its native ligand ADP was stronger than that between AURKB and HNK according to the docking score (−12.89 and−8.68, respectively, Table 4). Furthermore, ADP interacted with several amino acids of AURKB, such as     Evidence-Based Complementary and Alternative Medicine Gly108, Gly110, Lys111, and r112, while HNK interacted with only one amino acid (Pro27). Likewise, the binding interaction between Rac-1 and its native ligand (phosphoaminophosphonic acid-guanylate ester/GNP; docking score −21.38) was stronger than that between Rac-1 and HNK (docking score −18.17).
is was a result of the number of amino acids that the compounds interacted with and the distance between the interacting amino acids and the compound. For instance, the distance between r17 and GNP was much closer (1.82Å) than that with HNK (3.35). However, despite the lower affinity of the interaction, HNK could potentially compete with the native ligands to inhibit the function of these proteins.

Discussion
is study identified eight PTTHs, including CCND1, SIRT2, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1. e expression levels of these genes were strongly associated with immune infiltration levels. e tumor microenvironment influences angiogenesis and the immune response, and has long been recognized as a primary determinant of long-term tumor progression [39][40][41].
In addition, it can greatly impact the effectiveness of immunotherapy, highlighting the need of its further understanding [42]. Cyclin D1, encoded by CCND1, is considered an oncogene that promotes cell proliferation, growth, angiogenesis, and resistance to chemotherapy and radiotherapy [43,44]. In this study, we revealed that CCND1 expression was positively correlated with BRCA purity and macrophage infiltration levels, and negatively correlated with B cell and dendritic cell infiltration levels. Many studies have shown that tumor-associated macrophages play an important role in the proliferation, invasion, angiogenesis, and metastasis of human breast carcinoma, and that increased macrophage tumor infiltration confers metastatic potential and is associated with poor prognosis in breast cancer [42]. Our findings are in line with those of Pestell et al., who demonstrated that cyclin D1 expression was increased in human cancer stroma, and promoted tumor inflammation, angiogenesis, and stem cell expansion in advanced breast cancer [41]. Interestingly, a previous study reported that HNK could inhibit cyclin D1 expression [14].
SIRT2, an NAD-dependent histone deacetylase, has been suggested to be a promising therapeutic target in cancer  deep deletion, 2: shallow deletion, 3: diploid; 4: gain; 5: amplification. Statistical analyses were done by one-way ANOVA using Tukey's multiple comparison test. e symbol * or * * or * * * or * * * * symbolizes P < 0.05 or P < 0.01 or P < 0.001 or P < 0.001, respectively. (e) Pathways related to genetic alterations predicted by cBioportal. e results showed that genetic alterations of the PTTH disrupted the pathways regulating the cell cycle.  treatment [47]. SIRT2 is thought to affect carcinogenesis in a context-dependent manner, affecting epigenetic pathways implicated in cancer initiation, development, and progression [48][49][50]. SIRT2 expression level was negatively correlated with BRCA purity, but positively correlated with CD8+ T cell, CD4+ T cell, macrophages, neutrophils, and dendritic cells' infiltration levels. ese results confirm those of previous studies, in which SIRT2 expression level and CD8+ T cell infiltration level were positively correlated in breast cancer patients [51]. In addition, systemic SIRT2 has been suggested to promote tumor development by suppressing NK cells [42]. Interestingly, SIRT2 expression was significantly lower in breast cancer than in normal breast tissue, suggesting that SIRT2 may act as a tumor suppressor during the initiation of tumorigenesis. Moreover, a previous study reported that high SIRT2 expression in advanced tumor tissues is associated with poor prognosis, suggesting that SIRT2 may function as an oncogene [50].
VEGFA is a cytokine that promotes vascular development and the formation of new blood vessels from preexisting vascular networks during embryogenesis [52][53][54]. In addition, VEGFA can also be released by cancer and stromal cells [55]. In several murine and human cancer models, it has been demonstrated that VEGFA stimulates the tumor-initiating epithelial-mesenchymal transition and metastasis, and that VEGFA expression levels are positively correlated with BRAC purity and neutrophil infiltration levels [56][57][58][59][60][61][62][63]. In line with these results, another study reported that patients with mBC had higher levels of circulating VEGFA than patients without metastases [64]. VEGF can stimulate neutrophil migration through the activation of VEGFR1, [65] and can prevent dendritic cells from maturing, resulting in cytotoxic T cells' inactivation [66]. Tregs, tumor-associated macrophages, and myeloid-derived suppressor cells are all highly induced by VEGF, resulting in an immunosuppressive TME [67]. Furthermore, VEGF increases the expression of PD-1 on CD8+ CTLs and Tregs in a VEGFR2-dependent manner, [64] as well as the expression of Fas ligand, interleukin (IL)-10, and prostaglandin E3, leading to cytotoxic T cells' depletion [65]. Hence, VEGF-A can be used as a biomarker for immune-targeting therapy in breast cancer patients [66].
Histone deacetylase 1 (HDAC1) is overexpressed in breast cancer cells and human breast cancer tissues and can trigger the proliferation and migration of these cells via activation of Snail/IL-8 signals [41]. HDAC1 suppression has been reported to reduce the invasion of breast cancer cells by inhibiting matrix metalloproteinase-9, [67] and to reduce PD-L1 and HLA-DR expression and Treg frequency in triple negative breast cancer [68]. In our study, there was a positive correlation between HDAC1 expression levels and tumor purity, and CD8+ T cells', CD4+ T cells', neutrophils', and dendritic cells' infilitration levels.
HSP90α, encoded by HSP90AA1, is the stress-inducible isoform of HSP90. Previous studies have shown that high expression levels of HSP90 (HSP90α and HSP90β) increase the likelihood of recurrence and distant metastases in triple negative and ER+/HER2-breast cancer, and are associated with higher mortality [69]. Overexpression of HSP90 in human breast cancer cells has been linked to enhanced cell proliferation [42] and metastasis, [70] as well as to short OS and aggressive clinicopathological characteristics, such as high clinical stage, large tumors, and lymph node involvement [71]. Lin et al. reported that elevated HSP90AB1 expression was linked to a better overall survival of ER-and Basal-like breast cancer patients [55]. However, we found that high HSP90AB1 expression was associated with poor prognosis in BRCA patients. Additionally, the expression of HSP90AA1 and HSP90AB1 was positively correlated with tumor purity; and the expression of HSP90AA1 was positively related to CD8+ Tcells, macrophages, and neutrophils, but negatively correlated with CD4+ T cells. Further studies for exploring the infiltration of CD8+, macrophages, neutrophils, CD4+, and the effects of HNK on HSP90AA1 and HSP90AB1 are warranted.
CASP9 encodes caspase-9, an initiator of the intrinsic apoptosis pathway [79]. When the apoptosome, a multimolecular complex comprising cytochrome c and the apoptotic peptidase activating factor 1 (Apaf-1), is formed, it cleaves pro-caspase-9, forming caspase-9, triggering the caspase activation cascade by activating executor caspases, including caspase 3 and caspase 7 to cleave other cellular targets [80]. e aurora kinase family includes Aurora kinase B (AURKB), a mitotic serine/threonine protein kinase, and aurora kinase A (AURKA), which is a member of the Chromosomal Passenger Complex (CPC). e CPC plays a role in cell cycle progression and is a prognostic marker of    Evidence-Based Complementary and Alternative Medicine breast cancer arising from BRCA2 mutation [81]. Deregulation of AURKB is observed in several tumors, and its overexpression is frequently linked to tumor cell invasion, metastasis, and drug resistance [82]. Hence, AURKB has emerged as an attractive drug target for the development of small-molecule inhibitors [83]. RAC1 and AURKB were selected for in-depth molecular docking analysis because they were connected to the PTTH-  Figure 6: Proposed mechanism of honokiol (HNK) in mBCSCs (created with https://BioRender.com). 14 Evidence-Based Complementary and Alternative Medicine mBCSCs-cell cycle axis.
e Rho-GTPase family, including Rho, Rac1, and Cdc42, regulates the cytoskeleton [84] and thus modulates cell motility, migration, and invasion [85]. Rac-1/Cdc42 activation can induce cell growth by activating the PAK1/cyclin D1 pathway, or cell death by activating the PAK1/Akt/BAD pathway [86]. Rho has been suggested to be a potential therapeutic target, since Rho and VEGFA crosstalk leads to cancer progression and metastasis [87].
is study revealed the potential targets and molecular mechanisms of HNK on the cell cycle of mBCSCs (Figure 6). It is known that HNK exerts anticancer effects by suppressing angiogenesis, migration, invasion, and proliferation in a variety of cancer cell lines and tumor models [12]. HNK inhibits the cell cycle via the PI3K/Akt/mTOR pathway by upregulating PTEN and P21, and suppressing p-Akt, cyclin D/CDK4, c-Myc, Rac1, and AURKB [13,14,22,71]. Angiogenesis is inhibited through the HIF1/NFkB pathway, which is activated under hypoxic conditions and blocks the release of VEGF. Immune infiltration analysis showed that HNK is correlated with VEGFA inhibition, suggesting HNK can effectively block VEGFR2. HNK was also found to reduce HIF-induced VEGFR/VEGF activation and inhibit matrix metalloproteinases activity and cell migration [88]. In addition, HNK can induce apoptosis through the upregulation of BAD, caspase-9, caspase-3, and caspase-8 [89]. In the tumor microenvironment, oncogenic drivers such as β-catenin, STAT3, PI3K/PTEN/AKT/mTOR, p53, NF-kB, and RAS/RAF/MAPK are activated to suppress the production of chemokines, reduce the recruitment of dendritic cells, macrophages, T cells, and NK cells to tumor sites, and to suppress the immune system of these immunocytes [90]. Furthermore, tumor-intrinsic signaling can cause tumor cells to express PD-L1, resulting in T cell dysfunction in the tumor microenvironment. is study highlights the potential of HNK as an immunotherapeutic agent for mBCSCs by modulating the tumor immune environment. However, the results of this study were obtained through bioinformatics studies; therefore, further in vitro, in vivo, and clinical trials are needed to validate the findings.

Conclusions
is study identified eight PTTHs consisting of CCND1, SIRT2, AURKB, VEGFA, HDAC1, CASP9, HSP90AA1, and HSP90AB1, which can inhibit the mBCSC-cell cycle axis. In addition, PTTHs may regulate the PI3K/Akt/ mTOR and HIF1/NFkB pathways. is study is speeding up the development of HNK as anti-mBCSCs by targeting certain genes. However, this study have several limitations; for example, the targets of HNK are predicted from database. Other additional machine learning algorithm will provide more validated candidates of HNK targets. Another limitation of this study is that we used a bioinformatics approach; therefore, more needs to be explored further for validation and clarification in laboratory experiments.
Data Availability e data generated or analyzed during this study are included within the article and its supplementary files.

Conflicts of Interest
e authors declare that they have no conflicts of interest.