Cancer is considered as one of the worldwide life-threatening and the leading causes of human mortality (Vineis and Wild, 2014; Organization, 2020). According to the latest data released by the International Agency for Research on Cancer (IARC) on 14 December 2020, the annual incidence of cancer in 2020 reached 19.3 million cases and 10 million deaths. Furthermore, evidence based on the World Health Organization (WHO) suggests that there would be 29.5 million new cancer diagnoses and 16.4 million cancer deaths per year by 2040 (Shah et al., 2019; Sung et al., 2021). Accordingly, given the rapid development of oncology researches and the advancement of novel biotechnology approaches, determining different aspects of cancer progression can pave the way for improved cancer prognosis and treatment alternatives (Goyal et al., 2006; Charmsaz et al., 2018; Pucci et al., 2019). Herein, one of the challenges in the field of cancer treatment is the heterogeneity of tumor cells, which may lead to anti-cancer drug resistance or cancer treatment failure. Therefore, a full understanding of tumor heterogeneity can provide a clear picture of cancer progression and lead to the discovery of new cancer therapy options by researchers (y Cajal et al., 2020). Tumor heterogeneity is a condition in which tumor cells differ in various biological aspects such as function, differentiation, tumorigenesis, and sensitivity to anti-cancer therapies (Prager et al., 2019). In addition, depending on the type of heterogeneity, heterogeneous groups of tumor cells can have the same or distinct genomic content (Prager et al., 2019). In addition, heterogeneous populations of tumor cells can have the same or different genome content depending on the type of heterogeneity (Bedard et al., 2013). Hereupon, tumor heterogeneity can be divided into three types: 1) intertumor heterogeneity which is related to the variation of tumor cells among different patients, 2) intersite heterogeneity which is referred to the variation of cells among distinct tumors within a patient such as tumors in the primary site and tumors in the metastatic site, and 3) intratumor heterogeneity which is linked with heterogeneous populations of cells in a single tumor (Piraino et al., 2019). Oncology studies were shown that the cancer stem cells (CSCs) model is one of the models responsible for the generation of heterogeneous populations of cells, especially intratumor heterogeneity type (Prasetyanti and Medema, 2017; Turnquist et al.). Moreover, it can be caused by different factors such as genetics, epigenetics, and various micro-environmental features (Wang et al., 2015). Indeed, CSCs are a subgroup of cancerous tumor cells that display stemness abilities in the same manner as normal stem cells. For instance, they can self-renew to form the same daughter cells and give rise to differentiated multiple lineages of cells which form tumors. Additionally, the quiescence state is one of the distinguishing characteristics of cancer and normal stem cells, and it can play a role in therapeutic resistance and cancer progression (Hung et al., 2019; Lee et al., 2020). Furthermore, CSCs can make the treatment process more challenging because of their resistance to therapeutic approaches such as chemo and radiation therapies. The mentioned therapeutic resistance can be due to a variety of factors and mechanisms, including tumor environment, epigenetic effects, multidrug resistance proteins (MRPs) expression, various signaling pathways, effective mechanisms in DNA damage resistance, and the epithelial-to-mesenchymal transition (EMT) process (Phi et al., 2018). On the other hand, the function of metabolic pathways and processes are crucial in the growth, proliferation and survival of CSCs. In this respect, many investigations at the cellular and molecular level were indicated that unique forms of metabolic processes such as oxidative phosphorylation (OXPHOS), carbohydrate, and lipid metabolisms are observed in CSCs (Chae and Kim, 2018; Yadav et al., 2020). Therefore, the science of metabolomics, as well as the understanding of alterations associated with metabolic processes, could be useful in recognizing CSC behaviors and developing specific therapeutic methods for various types of cancers (Gilany et al., 2018; Rahim et al., 2018; Arjmand, 2019a, 2019b; Goodarzi et al., 2019; Larijani et al., 2019; Tayanloo-Beik et al., 2020), as well as the understanding of alterations associated with metabolic processes, could be useful in recognizing CSC behaviors and developing specific therapeutic methods for various types of cancers (Cuyàs et al., 2017). Additionally, scientists have been pushed to employ targeted approaches for treating cancer due to the problems in CSCs resistance to therapeutic methods. Molecular docking is one of the targeted approaches that play an important role in drug discovery and pharmaceutical researches. This computer-assisted drug design method is based on mathematical algorithms in which the effective biological binding-conformation between the drug and the target molecule can be evaluated. Indeed, the mentioned drug designing is based on the molecular structure that makes it possible to model and predict the molecular interactions as well as evaluate the biochemical processes (Meng et al., 2011; Phillips et al., 2018). Hereupon, in the present study, the cellular and molecular characteristics, signaling pathways, metabolic processes, and drug resistance of CSCs have been reviewed. We have also focused our discussion on molecular docking as a novel therapeutic approach in CSCs targeting.

CSCs are a subset of cancer cells or tumor-initiating cells (TICs) that serve as stem cells and contribute to the original tumor’s phenotypic variety (Lobo et al., 2007). They are found in variable amounts in different tumors. Furthermore, evaluating cell surface markers is the main strategy for detecting CSCs. Normal stem cells and CSCs have many similar characteristics (Jin et al., 2017; Khatami et al., 2019) such as 1) Self-Renewal (Lobo et al., 2007) 2) Differentiation capacity (Mohr et al., 2015) 3) Tumorigenesis (Zhu and Fan, 2018) 4) Capacity of developing resistance to drugs/cytotoxic substances and radiation (Schöning et al., 2017). Despite their similarities, there are some distinctions between cancer and somatic stem cells. The first is the origin of these two types of stem cells: natural somatic stem cells arise during embryonic development and separate from each other. They differentiate and produce a variety of mature cells, while CSCs are differentiated from normal adult stem cells or by multiple mutations in a single cell. The second is the ability to regenerate itself: somatic stem cells regenerate more regular than CSCs, although both types of cells can regenerate themselves. Finally, the organogenesis ability of these two cells is studied: both cells have the ability to organogenesis, but CSCs produce abnormal tissue, whereas somatic stem cells’ organogenesis produces normal tissue (Gjorevski et al., 2014).

Drug Discovery is considered as a multi-process platform in which a specific chemical compound with desired biological activity on the drug target can be selected and eventually enter the drug development as a candidate drug. In this platform, both chemical compounds and biological targets are evaluated from different aspects by using various approaches. Since both drug discovery and development are time-consuming and cost-intensive programs, they pose many challenges for researchers in drug designing and discovering for various diseases such as different types of cancer. Therefore, the use of new technologies can pave the way discovering new drugs with high therapeutic potential and take a big step towards disease treatment. Compound screening assays are one of mentioned methods that can help with hit identification, validation, lead generation, and optimization processes, as well as evaluating the compounds’ effects on the therapeutic target. With the advancement of technology and the integration of computational science with biological and pharmaceutical studies, approaches such as virtual screening are widely applied in drug designing and discovering programs (Reddy et al., 2007; Hughes et al., 2011; Mohs and Greig, 2017; Cui et al., 2020). In this context, virtual screening aims to evaluate and filter a limited number of suitable chemical compounds from large libraries of small molecules by using mathematical calculations. Structure-based virtual screening (SBVS) is one of the virtual screening methods which attempts to model and analyze the efficient biological binding-conformation between a ligand and a target molecule by using the molecular docking technique (Liao et al., 2013). Molecular docking is one of the cutting-edge computational drug designing technologies in which the most effective and stable state form of the ligand-receptor complex can be predicted (Morris and Lim-Wilby, 2008). The determination of the three-dimensional structure of the target and ligand molecules is at the top of the entire process priority list. Therefore, some techniques such as x-ray crystallography, nuclear magnetic resonance (NMR) spectroscopy, cryo-electron microscopy (cryo-EM), and homology modeling are not only useful in determining molecular structure, but also as complementary tools in drug development (Allen and Stokes, 2013; Kershaw et al., 2013; Sturlese et al., 2015; Lohning et al., 2017). Molecular docking includes searching algorithms and scoring function as two fundamental aspects of docking programs. Searching algorithms can be defined as a process that can lead to exploring the predominant and effective matching docking modes of ligand to the molecular target among the myriad configurations. Because a large number of binding modes are actually found between a ligand and a biological target molecule, searching algorithms not only can consider the optimum possible orientations of the ligand with the target but also can be an economical and time-saving solution in the docking process (Meng et al., 2011; Salmaso and Moro, 2018). Molecular dynamics, distance geometry methods, point complementary methods, fragment-based methods, Mote Carlo methods, genetic algorithms, systematic searches, and incremental construction, are some of the examples of search algorithms that can be used in modeling and evaluating the binding form of a ligand molecule to the objective receptor.

After the algorithm searching, it is time for the scoring function to step into the docking arena to find the good pose between the ligand and the target molecule. Scoring function refers to a process in which putative docking modes are ranked by evaluating their binding affinity and lowest binding energy to achieve top-ranked poses between a ligand and a target molecule. Force field function, Empirical scoring functions, knowledge based scoring functions, knowledge-based potentials, machine learning based scoring functions, comparative assessment of scoring functions, physics-based methods, and descriptor-based scoring functions are some of the examples of scoring function classifications in molecular docking (Madhavilatha and Babu, 2019; Sethi et al., 2019).

Regarding scoring function, the study of Wang et al. (2003) is one of the best examples of meticulous studies of popular scoring functions in molecular docking. In this study, the authors compared 11 popular scoring functions, including four scoring functions of the LigFit module in Cerius2 (LigScore, PLP, PMF, and LUDI), four scoring functions of the CScore module in SYBYL (FScore, G-Score, D-Score, and ChemScore), the scoring function of the AutoDock program, and two stand-alone scoring functions (DrugScore and X-Score) by performing them on 100 protein−ligand complexes to scrutinize their performance and sift the most effective and efficient methods among them. In this regard, after generating a set of docked conformations for each ligand by Autodac software, each 11 scoring function was tested and implemented on the maintained set and significant results were obtained. In this study, the authors used root-mean-square deviation ≤2.0 Å as a criterion for examining those scoring functions. Based on the analysis of the mentioned criterion, it was concluded that six scoring functions, including PLP, F-Score, LigScore, DrugScore, LUDI, and X-Score achieved high success rates (about 66–76%). In addition, the study implied that the combination of some of those scoring functions and generating the consensus scoring scheme can also increase the success rate (about 80%) (Figures 1, 2). In addition to success rates, the authors examined the correlation between 100 complexes’ binding scores and experimentally determined protein-ligand binding affinities. As a result of this experiment, X-Score, PLP, DrugScore, and G-Score could represent correlation coefficients of more than 0.50, which demonstrate superiority over other scoring functions. Since the best scoring function should perform excellently in both docking and scoring, the three scoring functions, including X-Score, DrugScore, and PLP, can be considered the top scoring functions in molecular docking, according to the study by Wang et al. (Wang et al., 2003).

Molecular docking is a substantial method for estimating the interaction between macromolecules such as protein and small molecules such as ligand. On the other hand, molecular docking is also capable of analyzing the molecular kinds of behavioral variability for those molecules which are located at the binding site of a targeted protein. Also, molecular docking is a computational approach and some docking programs are required to carry out its many functions. Some of the most considerable docking programs are Gold, Fred, and Flex. Moreover, they are useful in the prevision of protein and ligand`s binding conjunction (Kumar et al., 2013; Pagadala et al., 2017). In this respect, molecular docking is applied in different CSCs-related pathways including metabolic pathways and signaling pathways. Shedding light on metabolic pathways, the activity modification which can be applied by cancer cells, result in the production of metabolic precursors which leads to cancer cells anabolic and energetic requirement fulfilling. Furthermore, different metabolic pathways take a part in tumor progression and malignant tumor alterations. Accordingly, metabolic reprogramming is considered as one of the cancer insignias (Jagust et al., 2019b).

Herein, there will be a few examples in the context of some molecular docking usages in metabolic pathways such as 1) mitoketoscins application in targeting metabolic tumor promoters (OXCT1 and ACAT1) in both ketone re-utilization and mitochondrial function (Ozsvari et al., 2017). 2) Survivin protein interaction with andrographolide, which can lead to having an influence on human breast CSCs apoptosis (Wanandi et al., 2020). 3) S-phase kinase-associated protein-2 (Skp2) inhibition process by compound #25 which can result in CSCs survival suppression (Chan et al., 2013).

The other molecular docking-affected pathway in CSCs is the signaling pathway. This pathway is useful for targeted CSCs therapies expansion (Koury et al., 2017), embryonic evolvement, maintaining CSCs, etc. (Karamboulas and Ailles, 2013). Some molecular inhibitor agents of Wnt, Notch, Hh, and some other signaling pathways are implying as one of the important effects of molecular docking process on signaling pathways (Yang Y et al., 2020).

Modulating some target proteins is a striking aspect of molecular docking which has done by natural products. Natural products are able to be considered on the ground of multi-targeting drugs. As such, alkaloids are one of the natural products that have the strength to act as an anticancer lead molecule in the molecular docking process of CSCs. In this regard, Jaitak et al. provided an in-depth analysis of multitargeting drugs as an effective strategy to fight against CSCs and prevent disease recurrence. In this study, the authors examined the effect of some alkaloids that have anticancer potentiality by focusing on the Hh pathway in cancer stem cells. After selection and preparation of target ligands and proteins, Grid parameter selection and validation, implementation of glide docking module of Schrö6; dinger Maestro 9.6 suite, and determination of ADME profile for the studied alkanoid ligands, significant results were discovered. For instance, according to the findings of this study, emetine, and cortistatin, were able to target CSCs maintenance feature by binding to sonic Hh, smoothened (Smo) and, gli protein. Therefore, these two drugs can be applied as multi-targeting drugs in a combination with cancer chemotherapy compounds. Moreover, solamargine alkaloid could also have a good effect on gli protein and sonic hedgehog due to its pharmacophores. Furthermore, both solasonine and tylophorine modulated the Hh pathway and exert anticancer effects on CSCs by affecting only gli proteins. However, unlike other drugs, solamargine and solasonine need to improve the properties of ADME features (Jaitak, 2016).

In addition to alkaloids, the targeting of overexpressed CD44 surface marker in triple-negative breast cancer (TNBC) tissues can have an anticancer effect on CSCs. Regarding targeting CD44 surface markers, Yang et al. determined the positive role of drug carriers including Gambogic acid (GA)-loaded, zirconium-89 (89Zr)-labeled, chitosan (CS)-decorated multifunctional liposomes (MLPs) on TNBC CSCs by designing two in vitro and in vivo experiments. In this study, researchers examined 3D mammospheres and TNBC tissues of 32 women who were diagnosed with TNBC and found that the CD44 surface marker was overexpressed in the disease. Therefore, in this study, 89Zr@CS-MLPs were constructed and predicted how the drug carriers interact with CD44 surface markers in TNBC by applying molecular docking and dynamics simulations methods. The results obtained from the in vitro (examination on tumor cell lines) and in vivo (examination on mice) experiments were implied that 89Zr@CS-MLPs has a great potentiality for TNBC-targeted therapy as a drug carrier. Additionally, since Zr has a long half-life, it can also be used as an ideal radiolabel for positron emission tomography (PET) imaging of cancer. Moreover, 89Zr@CS-GA-MLPs have a high ability to target CSCs in vivo (Yang R et al., 2020).

In 2021, Hongwiangchan et al. demonstrated that hydroquinone 5-O-cinnamoyl ester of renieramycin M (CIN-RM) can be recognized as a fundamental approach in targeting lung CSCs which has been confirmed by molecular docking computational analysis. The effect of CIN-RM is based on the reduction of CSCs markers and upstream inhibition of the AKT pathway. As a result of inhibition of Akt, the expression level of transcription factors involved in self-renewal, such as c-Myc, Nanog, Oct4, and Sox2, are decreased. It is, therefore, CSCs are suppressed and tumor growth can be inhibited. To a lesser extent, CIN-RM also induces its inhibitory effect on the mTOR pathway, but the inhibitory effect of CIN-RM on the protein kinase C (PKC) signal pathway was not significant. Another important result obtained in this study is that CIN-RM even has an effect on inactivating the AKT pathway related to c-Myc regulation of lung non-stem cancer cells, which can be a promising therapeutic approach in cancer treatment (Figure 3) (Hongwiangchan et al., 2021).

A comprehensive analyzing of CSCs, including signaling pathways and metabolic activities, as well as recognizing their distinctions with normal cells, is an essential technique in cancer targeted treatment (Bjerkvig et al., 2005; Shyh-Chang and Ng, 2017). Major obstacles to molecular binding as a targeted therapeutic technique include receptor flexibility, ligand flexibility, and drug resistance, all of which have contributed to cancer therapy failure (Meng et al., 2011). Computational methods or computer tools, which are a form of artificial intelligence (AI), have recently proven to be a useful approach in a variety of domains, including structure prediction, molecular bond modeling, and junction prediction (Menke et al., 2021). The mentioned methods are divided into two groups: classical and machine learning (ML), of which ML is widely used in molecular binding. Computational methods play a key role in molecular binding for drug design and discovery, and results can be analyzed cheaper and often faster than other conventional methods (Torres et al., 2019). Eventually, with the advancement of science and the identification of various computational methods, computing software and hardware are still being updated, and researchers are looking for the most accurate way to target CSCs for cancer treatment (Phillips et al., 2018).