Introduction

Breast carcinoma is a major global health concern that affects millions of women globally and is the primary cause of fatalities associated with cancer. Breast cancer accounts for 2–3% of all cancer diagnoses and 14% of annual cancer deaths worldwide. It is also the most frequently identified cancer and the main factor in women’s cancer-related mortality (Bray et al., 2018). It is estimated that around 83,000 incidents of breast cancer are recorded annually in Pakistan. Approximately 40,000 women suffer from the disease each year. It does seem that the number of deaths from breast cancer is declining, though. According to the American Cancer Society predictions, 530 men and 43,700 women would have lost their lives to breast cancer in 2023 (Giaquinto et al., 2022). In spite of structural diversity, anticancer drugs have adverse effects. These adverse effects are mainly due to the detrimental effects that spread to normal cells. Over 80% of the global population, particularly in underdeveloped nations like Pakistan, depends on traditional medicines to fulfill all their medical needs. As one of the most frequent and deadly type of cancer in the world, breast cancer requires ongoing research to find new treatment approaches, including preventative strategies. Because of their various pharmacological characteristics, such as their anti-inflammatory, antioxidant, and anti-cancer effects, research on natural compounds, especially medicinal plant flavonoids (MPFs), has become more important in this regard (Nurgali et al., 2022; Ullah et al., 2020). This paper uses an in vitro investigation against MCF-7 and HepG2 cell lines, in addition to pharmacological networking, to validate MPFs as possible agents for the management of breast cancer. Approximately 60% of the currently used anticancer drugs are obtained from natural sources. Secondary metabolites exhibit potential to treat cancer and are still used widely as a source of therapeutic and preventive anticancer drugs (Rayan et al., 2017). During the last 10 years, a lot of research has been done to discover novel treatments that can lessen the adverse effects of the existing drugs (Dagogo-Jack and Shaw, 2018). Flavonoids (polyphenolic compounds) are produced by plants as secondary bioactive metabolites. They are responsible for pigment, taste, and pharmacological activities of plants (Scarano et al., 2018). Flavonoids plants possess antioxidant potential and shield plants from abiotic stress. Hence, various studies have been conducted to determine their beneficial role in treating a variety of acute and chronic human ailments (Vrhovsek et al., 2004). Pharmacological networking entails a thorough examination of the complex relationships that exist between flavonoids and important biological targets related to the advancement of breast cancer. The goal of this approach is to identify the intricate molecular processes that underlie flavonoids’ anticancer properties. This may yield important information for the logical development of designed treatments (Naponelli et al., 2015). Concurrently, MCF-7 and HepG2 cell lines research conducted in vitro provides essential platforms for assessing the effectiveness of herbal flavonoids in cancer prevention. A thorough evaluation of the broad-spectrum action of these drugs is provided by MCF-7, which represents breast cancer cells, and HepG2, which represents hepatocellular carcinoma cells. This study aims to find particular herbal flavonoids with strong anticancer potential by utilizing pharmacological networking and in vitro experiments. As phytochemicals are increasingly being considered as additional or alternative cancer treatments, the ultimate goal is to support the development of focused and effective therapeutic approaches for breast cancer (Martinez-Perez et al., 2014). To test our hypothesis, ethnobotanical data from Southern Punjab were collected and five plants from family Amaranthaceae were selected for further quantitative validation. Phytochemicals of the selected plants were retrieved from published resources while their molecular targets and anticancer mechanisms of action linked to the flavonoid detection were studied using network pharmacology. Seven extracted flavonoids are polar in nature. The most suitable solvent for the extraction of flavonoids as per literature were searched. For further verification, plant ethyl acetate extract was tested in vitro against MCF-7 and HepG2 cell lines. The extracts showing the best results were subjected to liquid chromatography–mass spectrometry (LCMS) studies to verify the presence of retrieved flavonoids from in silico studies. Collectively, this research validates the discovery of bioactive compounds, specifically flavonoids, for the treatment of breast cancer and provides scientific validation for the traditional usage of plants in cancer treatment.

Materials and Methods

Ethnobotanical data assemblage and analysis

Description of the study area

This multiregional study was carried out in different parts of the Southern Punjab, Pakistan. It was chosen as the study area because of its unique geographical location (Figure 1). The ethnobotanical investigation was performed in the following different regions of Southern Punjab: Layyah (area: 6,291 km²), Khanewal (area: 4,349 km²), Muzaffargarh (area: 8,435 km²), Vehari (area: 20 km²), Multan (area: 3,721 km²), Dera Ghazi Khan (area: 70 km²), Lodhran (area: 1,790 km²), Bahawalnagar (area: 8,878 km²), Rajanpur (area: 12,318 km²), Bahawalpur (area: 246 km²), and Rahim Yar Khan (area: 92.71 km²). In research regions, 990 persons of various ages were interviewed between 2019 and 2022. The criteria of selection of plant species was based upon the information we gathered from the respondents. First, all the plants were collected and quantitative analysis (use value [UV], informant consent factor [ICF], fidelity level [FL], and relative frequency citation (RFC) were applied. Plant species showing greater use value were selected for further analysis. Multiple visits were planned to cover all seasonal variations of the study areas.

Figure 1. Geographical location of Southern Punjab, Pakistan.

Ethnobotanical data collection and preservation

Ethical approval was sought from the Departmental Ethical Committee, LCWU, Lahore, Pakistan, while approval from local government was also sought for ethnobotanical data collection. The face-to-face interview method was used to acquire ethnobotanical data. This procedure was implemented due to low literacy rate in the study area. As a result, the same methodology was used to collect data for therapeutic herbs as previously reported (Ahmed et al., 2007). The standard technique was employed for the preservation of plant material. Proper pressing keeps the necessary plant parts visible for identification and keeps them from curling throughout the drying process. Herbarium specimens that have been carefully pressed are more functional and attractive to the eye. The procedure involves arranging the plant specimens on a pressing frame that has straps tightened around it after they have been arranged in folded newsprint sheets partitioned by cardboard sheets (Balick, 1996).

Data analysis

Medicinal plants collected during this study were examined using a variety of quantitative characteristics (UV, ICF, FL, and RFC) and sorted in alphabetical order. The collected ethnobotanical data include names of plants, family, life form, part used, and diseases cured. First, plants were dried by pressing between the layers of newspaper to remove its moisture at room temperature for 1–3 weeks, and then mounted on herbarium sheets. The Flora of Pakistan (http://www.efloras.org/index.aspx) was employed to confirm the nomenclature, while for identifying the correct botanical name, the International Plant Name Index (IPNI) (www.ipni.org) was utilized.

Quantitative analyses of ethnobotanical data

Informant consent factor, fidelity level, use value, and relative frequency citation

Quantitative analysis of ethnobotanical data (ICF, FL, UV, and RFC; Bennett and Prance, 2000; Tabuti et al., 2003) was computed as cited in the literature. ICF was examined to demonstrate the uniformity of the information for various disease categories. The calculation was as follows:

where “Nur” and “Nt” represent the number of total use reports for a particular ailment category, and accordingly the number of taxa used for a specific ailment category. This consensus result varies from 0 to 1.

Fidelity level predicts the preference of one species over another for curing a particular disease. Higher value suggests that a plant is more effective at treating a certain condition whereas lower FL values indicate that a plant is less effective at treating diseases. It was derived using the following equation:

where “Ip” is the number of informants providing information about the use of a certain species for a particular disease category, and “Iu” represents the number of respondents reporting using the plant for any condition category.

The use value of reported species was computed using the following formula

where “U” is the total number of use reports per species, and “n” represents the total number of informants surveyed for a given plant. UVs are high if there are more usage reports of a plant, indicating that the plant is significant, whereas they are close to zero if there are few use reports.

The RFC was calculated using the following equation:

This relation indicates the local prevalence of each species and is determined by dividing “FC,” the number of informants reporting the usage of species by the total number of informers contributing to the survey (N) without considering the use-categories.

Selection of plants for analysis

Plants with higher UV values were chosen for further examination, while diseases with a higher ICF index were chosen for additional research.

Network pharmacology studies

Mining of chemical constituents of selected plants

We obtained information regarding the phytochemicals of the selected plants from published resources. The plant’s name was entered into several search engines, including KNApSAcK (http://www.knapsackfamily.com/), PubMed (https://pubmed.ncbi.nlm.nih.gov/), GoogleScholar (https://scholar.google.com/), Tropicos (https://www.tropicos.org/), etc., to find active components. The PubChem database (available at: https://pubchem.ncbi.nlm.nih.gov) was explored to find the molecular formula, molecular weight, simplified molecular-input line-entry system (SMILES), etc. Using these canonical SMILES, the pharmacokinetic properties of complete active chemical ingredients were determined.

Drug-likeness (DL) and chemical absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties

The online admetSAR tool (http://lmmd.ecust.edu.cn/admetsar2) (Yang et al., 2019) was utilized to evaluate the ADMET parameters for each component (Li, 2001). The oral bioavailability (OB) evaluates the oral bioavailability of pharmacological compounds, while DL reveals whether a medication and a component are sufficiently comparable to be considered potential pharmaceuticals. Then, the active chemicals that fulfill drug-likeness (DL ≥ 0.18) and oral bioavailability (OB > 30) requirements were utilized for further screening (Li, 2001).

Possible cancer target screening

SwissTargetPrediction (http://www.swisstargetprediction.ch/) is a computational tool used to guess the protein targets of bioactive compounds, such as small molecules and natural products. It was used to determine the possible targets of certain active components by providing it with their canonical SMILES strings from PubChem. The species under consideration was Homosapiens, and the expected outcomes were obtained. The possible targets were selected from GeneCards (https://www.genecards.org/) and database of gene disease association (DisGeNET) (https://www.disgenet.org/) using “BreastCancer” as the search keyword (Zhang et al., 2020). Targets from both datasets were consolidated, and redundant gene entries were excluded to enhance the integrity of the subsequent analysis. Using UniProtKB, the target’s standard name was ascertained, with Homosapiens being the selected organism. A Venn diagram was created to show the mutual objectives of stability and drug for prevalent diseases and chemicals that were identified as potential study targets.

Gene expression and enrichment analysis

Gene Ontology (GO) function and pathway enrichment for breast cancer targets were examined using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) (Huang et al., 2009). The output data with P < 5×10^-2 were selected for further analysis, with smaller P values indicating greater enrichment. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of DAVID was used to identify biological pathways in the gene targets of breast cancer, and pathways with P ≤ 0.05 were analyzed further.

Formation of protein association network

The breast cancer gene targets were analyzed for protein–protein interactions using the online search tool for the retrieval of interacting genes/proteins (STRING) program (version 11.5) (https://string-db.org/). The interactions were then visualized and analyzed using Cytoscape (version 3.9.1) (Bhagya, 2023).

Analysis of protein–protein interaction (PPI) network

Cytoscape (version 3.9.1) was employed to generate a network of bioactive chemicals, proteins, and pathways, and duplicate data were deleted. In this network, nodes symbolize bioactive compounds and target genes, while edges symbolize connections between the compounds and targets. The significance of a component, target, or pathway in the network was measured using the degree, a topological metric that was calculated using Cytoscape’s analysis tool. The CytoHubba plugin for Cytoscape was used to identify target genes in the network.

Molecular docking

The Protein Data Bank (PDB) provided the PDB formatted protein crystal structures of the top six putative breast cancer core targets. The water and other small molecules from the protein crystal structure complexes were removed using PyMOLv 2.4 (Schrödinger, L.L.C) (Schrödinger and Delano, 2020) and saved in the PDB format. The three-dimensional (3D) configurations of selected phytomolecules were attained from the National Center for Biotechnology Information (NCBI) PubChem (Kim et al., 2021), an online database in Spatial Data File (SDF) format. Furthermore, energy minimization was performed using the Yet Another Scientific Artificial Reality Application (YASARA) software (Krieger et al., 2009), subsequent docking calculations, adding Gasteiger charges, using AutodockTools1.5.6 (Morris et al., 2009), the structures were translated to PDBQT format. Following the acquisition of protein PDB files, AutoDock Tools 1.5.6 was used to process them, combining nonpolar hydrogens and adding Gasteiger charges before converting them to the PDBQT format. AutoDock Vina software (Scripps Research Institute) (Eberhardt et al., 2021) was then used to simulate molecular docking between the target proteins and potential bioactive compounds to determine their binding affinity. The binding affinity was used as an evaluation criterion, with lower values indicating better docking.

Simulation studies

MD simulation

Molecular dynamics (MD) simulations were carried out to evaluate the binding affinities of the top-performing compounds following docking. The GROningen MAchine for Chemical Simulations (GROMACS) 2020 software was used for these simulations, which lasted for 200 ns. The protein topology was generated using the CHARMM36 force field (Huang and Mckerell 2013). While using the CHARMM general force field (CGenFF) server (https://cgenff.umaryland.edu), for MD simulation, the necessary ligand topology and parameters were produced. All the systems were solvated using the TIP3P water model and then neutralized with the necessary amounts of Cl⁻ and Na⁺. Subsequently, the steepest descent minimization technique was employed to minimize the energy of every system, with a maximum of 50,000 iterations and less than 10.0 kJ/mol of force. A LINear Constraint Solver (LINCS) holonomic constraints, a 2 femtosecond (fs) timestep, the number of atoms, volume, and temperature (NVT) ensemble with a leap frog integrator, and position constraints were applied to the receptor and ligand of both systems for 100 ps during heating (300 K). During the number of atoms, pressure, and temperature (NPT) equilibration phase, the NPT ensemble was employed for 100 ps at a temperature of 300 K with a 2 fs timestep. The structure’s coordinates were saved every 10 picosecond (ps) during an MD production run lasting 200 ns at a timestep of 2 fs, which came after all systems were optimized and energy minimized. The trajectories were used for different dynamics evaluations after a 200 ns MD simulation, including root mean square deviation (RMSD) of ligands relative to the backbone of proteins. The amount of H-bonds between the ligand and proteins was estimated over a 200-ns period. The Coulombic short range (Coul-SR) and Lennard Jones short range (LJ-SR) ligand–protein interaction energies were also calculated.

Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA)

A method for determining end-state free energy using GROMACS MD trajectory data is MM-PBSA. The user can choose from several parameters when using gmx MMPBSA, comprising estimates of the binding free energy using different solvation models (PB, GB, or 3D-RISM), computational alanine scans, entropy adjustments, and stability calculations (Valdés-Tresanco et al., 2021).

In vitro studies

Anticancer activity

The 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl tetrazolium bromide (MTT) assay was used for analysis, and the plant extracts’ cell viability was recorded. HepG2 and MCF-7 human cancer cell lines were employed. The cell lines were grown in 10% fetal bovine serum (FBS) RPMI-1640 medium. The cells were kept in a humidified environment with 5% CO₂ at 37°C.

MTT assay

The assay targets the mitochondrion’s nicotinamide adenine dinucleotide (NAD) + hydrogen (H) (NADH) dehydrogenase pathway in living cells. During this colorimetric approach, yellow tetrazoliren crystals breakdown to generate purple formazan. Since MTT formazan is insoluble in water, it crystallizes in cells into pro-pie needle-shaped crystals. MCF7 and HepG2 cells were used. Primary human hepatocytes (HH) are considered as the most reliable model for xenobiotic metabolism and cytotoxicity studies. However, a number of issues, such as insufficient availability of fresh human liver, difficulties in isolation procedures, brief life span, individuals’ heterogeneity, and expensiveness of the method, can limit its use in in vitro screening processes (Madan et al., 2003). Hepatocellular carcinoma cell line HepG2, the immortalized liver-derived cell line, is the best alternative proposed thus far, because it has unlimited availability and phenotypic stability. They have close genotypic resemblance to the normal liver cells with distinct differentiation (Sassa et al., 1987). Thus, HepG2 cell line provides the ideal screening approach for cytotoxicity potential of new chemical entities at the lead generation phase.

MCF-7 and HepG2 cells were grown in RPM11640 media supplemented with 10% FBS, 100 mg/mL of streptomycin, and 100 units/mL of penicillin. Cell line subculturing was kept in an incubator at a temperature of 37°C and 5% CO. After repeated doses of 24 and 48 h, cultures were observed (Alley et al., 1998). For 24 h, cell culture was carried out in a 96-well plate at 37°C with varying concentrations of foetal bovine serum (FBS). Cells in the negative control received simple medium treatment. After removing the supernatant from each well, cells were twice cleaned with phosphate buffer saline (PBS). The MTT solution was then added to the cells. After 4 h of incubation, formazan crystals started to grow on the cells. The resulting formazan crystals were dissolved in 100 µL of dimethyl sulfoxide (DMSO), and a microplate reader (Bio-RAD680, USA) was used to detect the absorbance intensity at 570 nm. The experiment was set up in randomized complete block design (RCBD), and the proportion of viable cells was calculated in comparison to the cells that were left untreated (Shahraki et al., 2016).

LCMS analysis

LCMS of selected plant species

The LCMS analysis was executed for the separation of bioactive compounds from the ethyl acetate and n-hexane extracts of Salsola spp. and Digera muricata. The mobile phase consisted of two solvents: (a) deionized water with 0.04% acetic acid, (b) acetonitrile. For the purpose of ion identification, a full total ion chromatogram (TIC) positive mode scan was acquired. All the samples were eluted within 24.3 min. The TIC scan summary for all extracts revealed a similar pattern of peaks eluted. Retention time is the amount of time a solute remains in the column or in the polar mobile phase and non-polar stationary phase. Because of their physiochemical characteristics, different components in the solution have different flow velocities (Farag et al., 2007).

Results

Ethnobotanical data results

Sociodemographic profile of respondents

Different field visits were conducted to collect voucher specimens. Following careful drying and identification, all of the voucher specimens were adhered to the herbarium sheets and then stored in Herbarium, LCWU. Ethnobotanical information about plants for all seasonal variations from 2019 to 2022 was covered. There were 990 individuals surveyed, and the sample size was determined in the same manner as stated previously (Kadam and Bhalero, 2010). Nearly all of respondents were natives (Figure 2), and were aged 30–40 years (27%) and 40–50 years (23%). Most of the respondents were males (55%) while rest were females (45%); only 2% respondents were illiterate.

Figure 2. Demographic profile of respondents.

Life form of plant species and plant parts used for curing diseases

This study reveals that 51% of the flora comprised herbs, followed by shrubs (28%), trees (15%), and weeds (6%) as shown in Figure 3. Herbs may predominate because bioactive compounds help them to adapt to their surroundings (Howard et al., 2012). Herbs are followed by shrubs, whose dominance may be attributed to their multifaceted environmental benefits, as shrubs improve air quality by absorbing dust and contaminants. Shrubs are dependable and easy to cultivate if grown in a suitable climate and good soil, and they prevent erosion, which reduces the amount of storm water runoff and harmful contaminants in waterways (Howard et al., 2012).

Figure 3. Life form of plant species.

Leaves, stems, roots, flowers, fruits, and other plant components are used to treat various ailments. However, according to the current study, leaves are the most prevalent plant portion utilized to treat various diseases, accounting for around 32% of the research area. Because of the ease with which leaves could be obtained and handled, they were widely used to cure a variety of diseases. The leaves are followed by the entire plant (26%), roots (18%), fruit (14%), flowers (9%), and stem (3%); also, most leaves were used to manufacture herbal remedies (Figure 4). Leaves are increasingly being used in the production of herbal remedies not only in Southern Punjab, but across the province (Santosh Kumar et al., 2015).

Figure 4. Plant parts traditionally used for treating diseases.

Informant consent factor

The ICF values ranged from 0.08 to 0.96 for 11 conditions, such as fever, flu, analgesics (pain), toothache, gastrointestinal, urinary, respiratory, muscle, skin, digestive, heart diseases, and tumors (Table 1). Tumors had the highest ICF of 0.96, indicating that it is more prevalent at the research site, indicating that people were more responsive about the practice of medicinal plants to treat this disease. Analgesics (0.08) and urinary diseases (0.05) had the lowest ICF. Lower ICF values suggest that people in this area know little about using plants to treat these diseases. Typically, types of plants native to a region dictate the efficiency of ICF for treating diseases (Rajakumar and Shiwanna, 2009). It was observed that natives were using plants (Abutilon indicum (L.) Sweet, A. javanica, Alhagi maurorum Medik, Artemisia vulgaris L., Calotropis procera (Aiton) W.T. Aiton , D. muricata, Fagonia indica Burm f., H. curassavicum, H. strigosum, S. baryosma, Salsola aphylla, S. kali, Tamarix a phylla (L.) H. Karst, and Thuja occidentalis L.) for the treatment of tumors.

Use value, relative frequency of citation, and fidelity level

In the present study, UV ranged from 0.37 to 0.090 (Supplementary Table S1). The maximum UV was noted for C. album (0.090). Many regions of Pakistan also utilized UV-rich plants. Future herbal drug development would emphasize on plant species with greater UV levels, improve the sustainability and preservation of plant resources, and follow pharmacological and phytochemical screening. While the minimum UV was reported for Flueggea leucopyrus Willd (0.37) and Sophora millus (Royle) Baker (0.37), the majority of respondents knew little to nothing about these plant species or their ethnobotanical applications. Less information about a given species in the study area was indicated by lower UV levels, despite the fact that it was not possible to correlate quantitative data within the region prior to the first quantitative ethnobotanical documentation in this area, notably in Southern Punjab. The value of RFC ranged from 0.1 to 0.86 (Supplementary Table S1). The highest value of RFC was discovered in C. album (0.86). The most prevalent plants at that location were those with the highest RFC, and the majority of individuals feel that they had therapeutic potential. While the lowest value of RFC was present in Capparis decidua Edgew. In the present study, FL extended from 10.0% to 75.8% (Supplementary Table S1). The higher the FL value, the more would be the usage of plant. Maximum FL was present in S. baryosma (75.8%), while minimum FL was observed in C. decidua (10.0%). Respondents utilized these herbs to treat illnesses and for other uses. Informants dealing with particular illnesses revealed importance of higher FL value (Islam et al., 2014).

Network pharmacology analysis

Active constituents selection of collected plants

After locating, filtering, and removing duplicates, approximately 494 active compounds of S. Kali, S. baryosma, D. muricata, C. album, and A. javanica were obtained through a published literature and two databases, ADMET analysis and drug likeliness and OB (DL ≥ 0.18 and OB ≥ 30). Through the results of ADMET analysis of these 494 compounds, 15 active compounds (Figure 5) were chosen as effective components (DL ≥ 0.18 and OB ≥ 30). The molecular structure of these peculiar compounds was established using PubChem.

Figure 5. Screened compounds with oral bioavailability, drug likeliness, and chemical structures.

Identifying potential targets

The SwissTargetPrediction database was utilized to build the 700 potential targets from seven chemical components. The GeneCard (Stelzer et al., 2016) and DisGeNet (https://www.genecards.org/) databases were explored to find prospective breast cancer targets, yielding 15,447 and 184 potential breast cancer targets, respectively. Prospective mapping of these 700 active protein targets in breast cancer revealed 16 common targets, which were classified as prospective breast cancer targets (Figure 6).

Figure 6. Venn diagram showing common targets.

Protein association network

A PPI network was acquired by STRING to illustrate relationships between 16 common targets of breast cancer, as shown in Figure 7.

Figure 7. Protein–protein interaction (PPI) network of common genes.

Network Compound Target

Cytoscape was used to construct a network to investigate the relationship between potential targets and active compounds. Seven active ingredients and 16 putative target genes were utilized to create a compound network target (Figure 8).

Figure 8. Disease compound network target.

The core orange node of the network symbolizes the Amaranthaceae family whereas purple nodes represent plant components found in the Amaranthaceae family. The remaining yellow nodes, on the other hand, are the possible breast cancer targets. A molecular docking investigation was performed on each target. Research of the Target–Compound network found that while the identical target may intermingle with several active compounds, one vigorous component could have an impact on many targets. This reveals that the plants in the Amaranthaceae family have numerous targets and multi-component effects in the treatment of breast cancer.

Exploration of protein–protein intersections network

Protein-protein interactions are crucial because they reveal the relationship between targets and are adaptable, flexible, and selectable. The CytoHubba plugin, which has roughly 12 topological analysis methodologies, was used to identify target genes, such as androgen receptor (AR), ESR1, EGFR, and Cytochrome P450 family 1 subfamily A member 1 (CYP1A1), with the highest degrees. These genes are represented by green nodes in the network (Figure 9).

Figure 9. Target PPI network analysis. Topological analysis found four green targets as key nodes (green nodes are target genes with higher degrees, and orange nodes are other potential targets).

Gene ontology analysis and KEGG pathway

Gene Ontology annotations and KEGG pathway analysis on 16 anti-breast cancer targets revealed the molecular mechanism of selected plants in the management of breast cancer. The GO analysis resulted in the identification of 61 biological processes (BP), including apoptotic process, estrogen biosynthetic process, steroid biosynthetic process, mammary gland alveolus development, intracellular steroid hormone receptor signaling pathway, etc. Thirteen cellular components (CC) (Figure 10), including macromolecular complex, cytoplasm, nucleus, intracellular membrane-bounded organelle, etc., and 35 molecular functions, including enzyme bind KEGG study (Figure 11) suggested 22 anti-breast cancer pathways.

Figure 10. Gene ontology term biological process (GOTERM_BP) and gene ontology term cellular component (GOTERM_CC) pathways showing richment factor.

Figure 11. Gene ontology term molecular function (GOTERM_MF) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway showing richment factor.

Molecular docking

It was used to find possible targets for drugs that could reduce the risk of developing breast cancer. The docking study accurately anticipated the components and the target gene’s great binding affinity. Linking the target genes with the outcomes of KEGG analysis in the breast-cancer pathway led to the identification of four target genes for docking: AR (PDBID:2AMA) with the resolution of 1.90, EGFR (PDBID:1IVO) with the resolution of 3.30, ESR1 (PDBID:1UOM) with the resolution of 2.28, and CYP1A1 (PDBID:6UDL) with the resolution of 2. Molecular docking was performed using AutoDock tools, and RMSD was determined using docking pose distance calculation (DockRMSD) program. The docking results of these four genes are displayed in Table 2.

The method of molecular docking was used to identify potential targets for substances that could decrease the probability of breast cancer. The strong binding affinity between the component and the target gene was correctly predicted by the docking analysis. Four target genes for docking were found by comparing the target genes with the results of KEGG analysis in the breast-cancer pathway: AR (PDBID:2AMA) with the resolution of 1.90 Å, EGFR (PDBID:1IVO) with the resolution of 3.30 Å, ESR1 (PDBID:1UOM) with the resolution of 2.28 Å, and CYP1A1 (PDBID:6UDL) with the resolution of 2.85 Å. AutoDock Tools was used for molecular docking, and RMSD was calculated from DOCKRMSD (Bell and Zhang, 2019). Compounds were docked with their respective targets as shown in the results of CytoHubba (Table 3).

The top six phytochemicals predicted by CytoHubba were docked against genes AR, EGFR, ESR1, and CYP1A1. The affinity of docked phytochemicals with their respective genes are presented in Table 3. By using the BIOVIA Discovery Studio program (version 4.5), protein ligand complex interactions were determined, which primarily evaluated various interactions (hydrogen-, hydrophobic-, pi-alkali-bond, and so on). Docking structures are shown in Figure 12. Docking structure of apigenin with AR along with interacting residues are shown in Figures 12A–E.

Figure 12. (A) Docking structure of apigenin with AR along with interacting residues. (B) Docking structure of luteolin with AR along with interacting residues. (C) Docking structure of chrysin with CYP1A1 along with interacting residues. (D) Docking structure of 5,7,2,3 tetrahydroxyflavone with ESR1 along with interacting residues. (E) Docking structure of ferulic acid with CYP1A1 along with interacting residues.

As Autodock vina perform rigid docking, where the receptor (usually a protein) is kept static while the ligand (the molecule being docked) is allowed to move. This approach does not account for the flexibility of receptor, which can lead to an incomplete representation of molecular interactions. These protein–ligand interaction includes ionic interactions, hydrogen bonds, and van der Waals interactions. Although quantum mechanics (QM) offers the most accurate estimation of these interactions, QM methods are often too computationally expensive for docking. To accelerate the process, simpler potential energy functions, typically related to force fields or statistical potentials, are used. However, despite improvements in force fields and scoring functions, they still lack detailed polarization effects and accurate proton affinity estimation. Strong binding contacts were observed between major active compounds and targets, as validated by molecular docking. The top three combinations are thought to be important in treating breast cancer because they could be effective on the main targets. Our results are supported by our docking study, and we can suggest these complexes as a treatment approach for breast cancer. Moreover, docking results were confirmed using MD simulations.

Simulation studies

Molecular dynamics simulation

Molecular dynamics simulations are used to comprehend binding dynamics between ligands and proteins. One of the most useful and often utilized computer programs for exploring biological macromolecules is MD simulation. This is very useful for comprehending the dynamic behavior of proteins at various time scales, such as fast internal motions, slow structural alterations, and protein folding processes. On the basis of the highest scores from docking results, MD simulations were performed on the following three complexes: apigenin-AR, apigenin-ESR1, and luteolin-EGFR.

GROMACS was used to simulate the MD of apigenin-AR, apigenin-ESR1, and luteolin-EGFR complexes with all atoms for 200 ns. RMSD, hydrogen bonding, and interaction energies, such as Columbic short-range (Coul-SR) and Lennard–Jones short-range (LJ-SR) were examined for the trajectories. The RMSD values demonstrated that the ligands had strong interactions with the proteins. RMSD values for all three ligands ranged between 0.1 and 0.38 (Figures 13B, 14B, and 15B). During the simulation period of 200 ns, these data suggest that ligands shifted somewhat from their initial positions. Similarly, ligands exhibited a high number of hydrogen bonds with protein structures throughout the same time. This also indicates that the ligands bind to their respective proteins more effectively. H-bonds stabilize protein–ligand complexes. H-bonds interaction exploration was carried out on MD trajectories to conclude the total number of H-bonds generated between protein–ligand complexes to comprehend the binding affinity of ligands to proteins. The apigenin-AR, apigenin-ESR1, and luteolin-EGFR complexes showed H-bonds between 0 and 4.1, 0 and 4, and 0 and 4, respectively. The outcomes revealed that throughout the simulations, the total number of H-bonds produced by all protein–ligand complexes remained stable, as shown in Figures 13C, 14C, and 15C. In addition to this analysis, the Coul-SR and LJ-SR interaction energies for MD trajectories were also estimated. Apigenin with AR (Figure 13A), apigenin with ESR1 (Figure 1A), and luteolin with EGFR (Figure 15A) were shown to have average Coul-SR interaction energies of –30.6433, –62.3507, and –95.9883 KJ/mol, respectively. The LJ-SR interaction energies were found to be –146.702, –120.451, and –127.218 KJ/mol. These findings demonstrated that the interaction between luteolin and EGFR is relatively stronger than that of the other two complexes. Complexes are shown in Figures 13D, 14D, and 15D.

Figure 13. Complex of apigenin with AR. (A) Time (ns) vs. interaction energies (KJ/mol) plot for the MD simulation of docking complex involving receptor AR and ligand apigenin. (B) Time (ns) vs. root mean square deviation plot for the MD simulation of docking complex involving receptor AR and ligand Apigenin. (C) Time (ns) vs. the number of hydrogen bonds plot for the MD simulation of docking complex involving receptor AR and ligand apigenin. (D) Movement of ligand during MD simulation. The position of ligand at 10 ns, 50 ns, 100 ns, 150 ns, and 200 ns is shown in green, cyan, magenta, yellow, and blue color, respectively.

Figure 14. Complex of apigenin with ESR1. (A) Time (ns) vs. interaction energies (KJ/mol) plot for the MD simulation of docking complex involving receptor ESR1 and ligand apigenin. (B) Time (ns) vs. root mean square deviation plot for the MD simulation of docking complex involving receptor ESR1 and ligand apigenin. (C) Time (ns) vs. the number of hydrogen bonds plot for the MD simulation of docking complex involving receptor ESR1 and ligand apigenin. (D) Movement of ligand during MD simulation. The position of ligand at 10 ns, 50 ns, 100 ns, 150 ns, and 200 ns is shown in green, cyan, magenta, yellow, and blue color, respectively.

Figure 15. Complex of luteolin with EGFR. (a) Time (ns) vs. interaction energies (KJ/mol) plot for the MD simulation of docking complex involving receptor EGFR and ligand luteolin. (b) Time (ns) vs. root mean square deviation plot for the MD simulation of docking complex involving receptor EGFR and ligand luteolin. (c) Time (ns) vs. the number of hydrogen bonds plot for the MD simulation of docking complex involving receptor EGFR and ligand luteolin. (d) Movement of ligand during MD simulation. The position of ligand at 10 ns, 50 ns, 100 ns, 150 ns, and 200 ns is shown in green, cyan, magenta, yellow, and blue color, respectively.

MM-PBSA calculations

Using this method, the binding free energies of the complexes were determined, and thermodynamic parameters ΔG, ΔH, and –TΔS were calculated. Here, ΔG represents the Gibbs free energy change upon ligand binding, ΔH denotes the enthalpy of binding, and –TΔS corresponds to the entropic contribution, specifically the conformational entropy associated with ligand binding. The calculated value is the effective free energy when the entropic factor is eliminated, and this is typically adequate for comparing the relative binding free energies of related ligands. The apigenin_AR showed ΔG = –14.43 kcal/mol, apigenin_ESR1 showed ΔG = –13.85 kcal/mol, and luteolin_EGFR showed ΔG = –8.13 kcal/mol. Stability of the complex depends upon the value of ΔG. More negative the value of ΔG, more stable is the complex. Hence, apigenin_AR and apigenin_ESR1 are more stable, compared to the third complex luteolin_EGFR. TΔS values of all complexes were positive as shown in Figure 16. For apigenin_AR, ΔH = –25.65 and -TΔS = 11.22. apigenin_ESR1 showed ΔH = –28.12 with -TΔS = 14.27. luteolin_EGFR showed ΔH = –21.49, while this complex showed –TΔS = 13.36.

Figure 16. Binding energies of apigenin_AR, luteolin_EGFR, and apigenin_ESR1 complexes.

Anticancer activity

The anticancer activity of S. kali, S. baryosma, D. muricata, C. album, and A. javanica extracts was evaluated through MTT assay. Two breast cancer cell lines were used, that is, HepG2 and MCF-7. N-hexane and ethyl acetate extracts were used for the calculation of these selected plants (Table 4). The percentage of viable cells was negatively correlated with the concentration of plant extract. There was a decline in the cell viability of cancer cells with an increase in concentration of plant extracts. Half-maximal inhibitory concentration (IC₅₀) of ethyl acetate extracts showed that D. muricata L. showed more cytotoxic activity against MCF-7 in ethyl acetate extract as compared to n-hexane extract, that is, 30 µg/mL. A study showed the IC₅₀ value of 50 µg/mL against prostate cancer cell line using leaf extract of D. muricata. Ethyl acetate extracts of S. kali, S. baryosma, C. album, and A. javanica showed the IC₅₀ values of 34.3 µg/mL, 40.1 µg/mL, 37 µg/mL, and <30 µg/ml, respectively, against MCF-7, while the IC₅₀ values of these plants in N-hexane extract was 48.1 µg/mL, 54.4 µg/mL, 42 µg/mL, and 36 µg/mL. S. kali also proved to be more cytotoxic against MCF-7 cell line. Any substance with anticancer properties may do so by either destroying cancer cells or altering the genetic make-up of healthy cells (Asare et al., 2015). While in case of HepG2 cell line in ethyl acetate plants, extracts showed the IC₅₀ values of 42 µg/mL (C. album), 39 µg/mL (A. javanica), 46 µg/mL (D. muricata), 51 µg/mL (S. baryosma), and 54 µg/mL (S. kali).

LC-MS study

Liquid chromatography–mass spectrometry was performed for the isolation of compounds from ethyl acetate and n-hexane extracts of best anticancer plants among the five selected plants. All of the compounds’ mass spectra were determined by measuring their retention time and mass-to-charge ratio (m/z). D. muricata and S. Kali were subjected to LC-MS analysis for confirmation of flavonoids among the two extracts, i.e., ethyl acetate and N-hexane. Consistent with other findings, the LC-MS results showed that the extracts include a variety of phytocompounds, including fatty acids, alkaloids, terpenes, and flavonoids (Al-Dalahmeh et al., 2022). It was discovered that flavonoids account for a greater reaction to different diseases and damage than any other chemical. Traditional folk remedies for gastrointestinal, vascular, and respiratory conditions use flavonoids. Ethyl acetate extract of D. muricata (Figure 17) contains various compounds, such as phenol,2,6-dibromo-4-nitro, benzene,4-chloro-2-iodo-1-methyl, cyclopropane carboxyanilide,3’,4’-dichloro,1,1-dichloro-2,2,2-trifluoroethyl chlorodifluoromethyl ether, l-methionine, N-methoxycarbonyl-, dodecyl ester, didocosyl succinate, l-methionine, n-heptafluorobutyryl-, octadecyl ester, tetratriacontyl pentafluoropropionate, tetrapentaconytyl benzene, 1-hexacontanethiol, lupeol, and luteolinetc. While its n-hexane extract showed the presence of barban, benzamide, N-(-4-flurophenyl)-4-methyl, cyclopropane, carboxyanilide, 3,4-dicloro, L-tyrosine, N-(trifluoroacetyl)-, butyl ester, trifluoroacetate (ester), octadecanoic acid, and 1,2-ethanediyl ester (Figure 18) In case of S. Kali, its ethyl acetate extract (Figure 19) showed various compounds, such as beta-cadinene, 4-aminobenzoic acid, N,N-bis(pentafluoropropionyl)-, tert.-butyldimethylsilyl ester, ethyl 5,7-dichloro-4-hydroxyquinoline-3-carboxylate, L-methionine, acetamide, N-(4-bromophenyl)-2-bromo, fumaricacid, D-alanine, fumaricacid, silane, L-tyrosine, N,O-bis(2,6-diflurobenzyl)-, methyl ester, beta-alanine, N-(2,6-difluorobenzoyl)-octyl ester, L-valine, n-heptafluorobutyryl-, octadecyl ester, ergotaman-3’,6’,18-trione, 12-hydroxy-2’-(1methylethyl)5’-(2-methylpropyl)-(5’alpha), Tetratriacontylheptafluoro butyrate, aluminium palmitate or hen pentacontylbenzene, L-methionine, L-leucine, n-heptafluorobutyryl, octadceyl ester, etc. Its n-hexane extract showed sarcosine, N-(2-thienylcarbonyl)-, pentodecylester, 3-phenoxybenzyl alcohol picolinyl oxy dimethyl silyl ester, 2-acetamino-5-bromobenzoic acid, trichlormethiazide, hexan amide, N-ethyl-N-(3-methhylphenyl)-6-bromo, triphenylphosphine selenide, tricosanoic acid, 1-dotetracontanethiol, hen tetracontyl cyclohexane, luteolin, etc. (Figure 20).

Figure 17. Chromatogram of ethyl acetate extract of Digera muricata.

Figure 18. Chromatogram of n-hexane extract of Digera muricata.

Figure 19. Chromatogram of ethyl acetate extract of Salsola kali.

Figure 20. Chromatogram of n-hexane extract of Salsola kali.

Discussion

Breast carcinoma is a frequent cancer in women and a primary cause of death globally (Darooei et al., 2017). Medicinal plants exhibit a wide range of pharmacological actions due to the presence of bioactive metabolites (Fagbemi et al., 2003). Ethnopharmacological research for cancer treatment is continuously performed globally. In this study, ethnobotanical data were acquired through interviews from selected sites in Southern Punjab (Pakistan). In all, 990 respondents were chosen to gather data about the customary uses of plants used in Southern Punjab. The objective of this study was to calculate the anticancer potential of various plant species by using various quantitative indices. Five plant species (S. kali, S. baryosma, D. muricata, C. album, and A. javanica) with high ICF, UV, and RFC as potential candidates for further analysis were identified using network pharmacology. The outcomes indicated that the active components of the selected plant species had the ability to target key proteins involved in breast cancer signaling pathways, but further research is needed to confirm their effects on functional alterations and treatment outcomes. The findings indicate that the active components of the chosen plants can target key proteins linked to a variety of biological processes and signaling pathways involved in breast cancer. However, the functional alterations in these factors and the profound effects of the treatment still need to be clarified in additional animal studies and human clinical trials.

The GeneCard and DisGeNet databases identified an overall of 15,447 and 184 possible targets for breast cancer, respectively. The SwissTargetPrediction database gathered an additional 700 potential targets from seven chemical compounds. By mapping the data from these databases and phytochemical targets, we were able to identify 16 common breast cancer genes. These genes were used to create a network of phytochemicals and gene targets using Cytoscape. Among these, the genes with the highest degree of values in CytoHubba were AR, ESR1, EGFR, and CYP1A. KEGG enrichment analysis was conducted and several signaling pathways were identified that play a part in the growth and advancement of the disease. Gene ontology analysis was also performed, which identified 61 biological processes related to breast cancer, including processes involving apoptosis, steroids, and estrogen synthesis and mammary gland development. Based on their PPI network and KEGG enrichment analysis, we selected the following four key target proteins that were closely related to breast cancer: AR, ESR1, EGFR, and CYP1A1. These proteins are linked to breast cancer in various ways. For example, AR expression is associated with increased survival in ER+/AR+ breast cancer, and ESR1, which comes in two forms, called ERα and ERβ, is most closely correlated with carcinogenesis, while lack of ERβ expression is linked to tumor progression (Caiazza et al., 2015). EGFR is overexpressed in around 50% of triple-negative breast sarcoma and inflammatory breast growth (Masuda et al., 2021), and CYP1A1 is suspected to be convoluted in breast cancer due to its part in the metabolism of certain compounds that may increase the danger of oxidative stress and cancer.

Several active ingredients that target these proteins, including apigenin, luteolin, ferulic acid, chrysin, and N-trans-feruloyl-4-O-methyldopamine were identified. These phytochemicals have demonstrated therapeutic effects and anticancer properties, particularly in breast cancer. Molecular docking studies were performed and found that the binding affinities of the docked complexes ranged from –6.3 to –10.7 kcal/mol, indicating stable binding.

We also conducted MD simulations on three top-scoring docked complexes for 200 ns. In this study, MD simulations were also conducted to comprehend the behavior of the final compounds. MD simulations are a powerful tool in drug innovation and proposal because they offer detailed insights into the behavior of molecules and their interactions at anatomic level (Liu et al., 2018). MD simulations are particularly useful for studying the dynamics of proteins, which play a central role in many biological processes and are often targeted by drugs. In drug discovery, MD simulations are used to identify potential drug candidates and to optimize their properties, such as binding affinity and specificity. They can also be used to study the mechanisms of drug action and to predict the potential adverse effects of drugs. For example, MD simulations can help researchers understand how a drug molecule interacts with its target protein, which can provide valuable information about how the drug is likely to behave in the body. In addition to their use in drug discovery, MD simulations can also be used to study the arrangement and functioning of proteins to calculate the stability of protein–ligand complexes, and to design new proteins with desired properties (Ashraf et al., 2022). Overall, MD simulations are an important instrument for understanding the behavior of molecules and their interactions, which is critical for the improvement of new drugs and therapies.

In our study, all three complexes, such as apigenin_AR, luteolin_EGFR, and apigenin_ESR1 showed considerably lower LJ-SR and Coul-SR interaction energies. This indicated the strong binding of ligands with the receptors. Furthermore, the RMSD values for all three ligands with reference to the backbone were in the range of 0.1–0.36 Å, indicating the small change in the ligands’ position, compared to their original docked position. Using the MM-PBSA method, the complexes’ binding free energies were determined, ranging from –8.13 kcal/mol to –13.85 kcal/mol.

Based on the MD simulation data, we can say that shortlisted ligands could be good therapeutic candidates to be used in the future treatment of breast cancer. Even though many compounds isolated from plants are tested thoroughly for their anticancer properties, it is becoming more widely accepted that the beneficial effects of plants are caused by an intricate interaction of a complex combination of compounds present in the whole plant (additive/synergistic, and/or antagonistic), rather than constituents in single agents alone (Karna et al., 2011). The anticancer activity of ethyl acetate and n-hexane extracts of S. kali, S. baryosma, D. muricata, C. album, and A. javanica was evaluated through MTT assay. Results showed that D. muricata and S. kali are more cytotoxic among selected plants. Ethyl acetate extract of D. muricata observed cytotoxic (30 µg/mL) against breast cancer MCF7 cell line. In order to resist apoptosis and escape controlled growth, cancer has developed a number of strategies. Therefore, utilizing entire cell extracts, which have multiple constituents with various potential intracellular targets, may be superior to employing a single isolated plant ingredient.

In LCMS, the extracted ion chromatogram, which was created by extracting each peak from TIC, revealed that each peak’s intensity varied at the same location. Each compound’s MS spectra were acquired by the individual data acquisition (IDA) technique. Since temperature affects the thickness of the mobile phase, it is crucial for all chromatography modes. Adropin liquid viscosity is a result of rising temperature. As a result, temperature has an impact on both chromatogram retention and solute partitioning (between mobile and stationary phases) (Milićević et al., 2010). The retention time increases with increasing ratio (k). As a result, retention of each metabolite, as well as total flow rate, is determined by column temperature. The ability of the solute to dissolve in the melted phase and adsorbed in the dense phase determines the solute’s dispersal between two phases (mobile and stationary). The proportion of solute absorption in the mobile phase to the ratio in the stationary phase is known as the distribution ratio. Luteolin and other flavonoids identified by LC-MS analysis are common flavonoids. As common flavonoids are found in a variety of medicinal plants, traditional Chinese medicine has utilized luteolin-rich plants to cure variety of illnesses (Su et al., 2003). Based on the results of LC-MS analysis data, the selected plants are greatly endorsed as potential anticancer therapeutics.

Conclusions

The outcomes of this study has paved the approach for new breast cancer targets, multi-target chemical regimes, and representations of treatment efficacy. Experimental validation of ethyl acetate extracts of all screened plants (Salsola kali L., Salsola baryosma (Schult.) Dandy, Digera muricata (L.) Mart, Chenopodiumalbum L., and Aerva javanica (Burm. f.) Juss. ex Schult) has demonstrated a significant growth inhibitory effect on MCF-7 and HEPG2 cancer cell lines. The network analysis results suggest that these plants contain flavonoids with numerous targets that may act on a variety of pathways connected to breast cancer. Furthermore, the genes AR, EGFR, ESR1, and CYP19A1 have been recognized as prospective and real therapeutic targets for the inhibition and decrement of breast cancer.