Systems and in vitro pharmacology profiling of diosgenin against breast cancer

Aim: The purpose of this study was to establish a mode of action for diosgenin against breast cancer employing a range of system biology tools and to corroborate its results with experimental facts. Methodology: The diosgenin-regulated domains implicated in breast cancer were enriched in the Kyoto Encyclopedia of Genes and Genomes database to establish diosgenin-protein(s)-pathway(s) associations. Later, molecular docking and the lead complexes were considered for molecular dynamics simulations, MMPBSA, principal component, and dynamics cross-correlation matrix analysis using GROMACS v2021. Furthermore, survival analysis was carried out for the diosgenin-regulated proteins that were anticipated to be involved in breast cancer. For gene expression analyses, the top three targets with the highest binding affinity for diosgenin and tumor expression were examined. Furthermore, the effect of diosgenin on cell proliferation, cytotoxicity, and the partial Warburg effect was tested to validate the computational findings using functional outputs of the lead targets. Results: The protein-protein interaction had 57 edges, an average node degree of 5.43, and a p-value of 3.83e-14. Furthermore, enrichment analysis showed 36 KEGG pathways, 12 cellular components, 27 molecular functions, and 307 biological processes. In network analysis, three hub proteins were notably modulated: IGF1R, MDM2, and SRC, diosgenin with the highest binding affinity with IGF1R (binding energy −8.6 kcal/mol). Furthermore, during the 150 ns molecular dynamics (MD) projection run, diosgenin exhibited robust intermolecular interactions and had the least free binding energy with IGF1R (−35.143 kcal/mol) compared to MDM2 (−34.619 kcal/mol), and SRC (-17.944 kcal/mol). Diosgenin exhibited the highest cytotoxicity against MCF7 cell lines (IC50 12.05 ± 1.33) µg/ml. Furthermore, in H2O2-induced oxidative stress, the inhibitory constant (IC50 7.68 ± 0.51) µg/ml of diosgenin was lowest in MCF7 cell lines. However, the reversal of the Warburg effect by diosgenin seemed to be maximum in non-cancer Vero cell lines (EC50 15.27 ± 0.95) µg/ml compared to the rest. Furthermore, diosgenin inhibited cell proliferation in SKBR3 cell lines more though. Conclusion: The current study demonstrated that diosgenin impacts a series of signaling pathways, involved in the advancement of breast cancer, including FoxO, PI3K-Akt, p53, Ras, and MAPK signaling. Additionally, diosgenin established a persistent diosgenin-protein complex and had a significant binding affinity towards IGF1R, MDM2, and SRC. It is possible that this slowed down cell growth, countered the Warburg phenomenon, and showed the cytotoxicity towards breast cancer cells.


Introduction
Breast cancer spreads through the inner layer of the milk gland or lobules and ducts (Sariego, 2010), and it is one of the second leading causes of mortality for women between the ages of 45 and 55 (Jemal et al., 2009). It includes age, iodine deficiency (Venturi, 2001;Aceves et al., 2005;Stoddard et al., 2008), high hormone levels (Russo and Russo, 2006;Yager and Davidson, 2006), and age-related (Steiner et al., 2008) risk factors. This incident potentially results in the breast being surgically removed entirely or in requiring chemo-, radio-, or hormone-therapy (HeraviKarimovi et al., 2006). However, these practises are preoccupied with multiple side effects that are not specific to the breast tumor. Furthermore, the currently used chemotherapeutic drugs in medical care result in anemia, exhaustion, mouth soreness, vomiting, and diarrhea. This suggests the requirement to discover a novel therapeutic agent against breast cancer.
Diosgenin is a phyto steroid sapogenin obtained from the hydrolysis by strong bases, acids, or enzymes of saponin. It is commercially used as a precursor to synthesis various hormones and steroid products like pregnenolone, and cortisone including progesterone (Marker and Krueger, 1940) for the early manufacture of combined oral contraceptive pills (Djerassi, 1992). Additionally, diosgenin has the potential to inhibit activated pro-inflammatory and pro-survival signaling pathways and promote the death of a variety of cancer cells (Shishodia and Aggarwal, 2006). Additionally, diosgenin inhibits the growth of oestrogen receptor-positive MCF-7 cells by activating caspase three and upregulating the p53 tumor suppressor gene. Further, BCL2 is downregulated in estrogen receptor-negative MDA-MB-231 breast cancer cells (Srinivasan et al., 2009). Although the potential of diosgenin to treat breast cancer has been demonstrated its interaction with proteins involved in the progenesis of breast cancer has yet to be investigated.
As a result, the current study focuses on locating the potential interactions of the diosgenin-regulated proteins and attributable pathways implicated in breast cancer using a range of system biology techniques i.e., gene ontology (GO) analysis, molecular docking, and molecular dynamics (MD) simulations, and reinforcing its findings using diverse functional biomarkers using in vitro experiments in four distinct cell lines.

Computational pharmacology
The detected diosgenin-regulated proteins were enriched to determine the altered pathways. Later, diosgenin was docked with proteins, and the proteins with the maximum diosgenin binding affinity were chosen for MD simulation. Furthermore, the diosgenin-regulated proteins involved in breast cancer were examined for survival and gene expression in normal and tumor cell lines, including the three hub proteins determined by molecular docking.
diosgenin-modulated targets were compared with DisGeNETrecorded targets to trace the reported targets involved in oestrogen receptor-positive breast cancer using PivotTable (Microsoft excel 2007; https://www.microsoft.com/en-in/ microsoft-365/excel).

GO and cluster analysis
The diosgenin-modulated targets involved in oestrogen receptor-positive breast cancer were queried in STRING (Szklarczyk et al., 2021; https://string-db.org/) ver 11.5 for "Homo sapiens" to trace three GO terms i.e., cellular components, molecular function, and biological processes. In addition, probable regulations of multiple pathways were also traced concerning the KEGG database (https://www. genome.jp/kegg/pathway.html) with whole genome statistical background. Also, the regulated proteins were concerned with subcellular location (COMPARTMENT), protein domains and features (InterPro), protein domain (Pfam), and tissue expression (TISSUES). Later, the protein-protein interaction (PPI) was also assessed for the cluster analysis via k means clustering to identify three lead distinct clusters.

Network construction and analysis
The network between diosgenin, its targets (involved in oestrogen receptor-positive breast cancer), and the regulated pathways were constructed using Cytoscape ver 3.5.1 (Shannon et al., 2003; https://cytoscape.org/). The constructed network was recognized as directed and inspected by translating node size and color to low values corresponding to small sizes and low values corresponding to bright colors toward edge count. In addition, the edge size and color were mapped to edge betweenness, with low values corresponding to small sizes and low values equating to bright colors.

Molecular docking
The regulation of insulin-like growth factor 1 receptor (IGF1R), E3 ubiquitin-protein ligase Mdm2 (MDM2), and proto-oncogene tyrosine-protein kinase Src (SRC) were majorly triggered in the network interaction among diosgenin-target(s)-pathway(s). As a result, these three proteins were considered for molecular docking studies.
Ligand-protein docking: Diosgenin was docked against the aforementioned targets using AutoDock Vina, which was run through the POAP pipeline (Trott and Olson 2010;Samdani and Vetrivel 2018;Patil et al., 2022). Nine different poses of ligand were obtained after docking. Docking results were analysed based on the binding affinity, and number of interactions as explained previously (Dwivedi et al., 2021;Badraoui et al., 2022). Further, the diosgenin pose with the lowest binding energy was chosen to visualize the ligand-protein interactions and perform MD simulation (Samdani and Vetrivel 2018).

MD simulation
To examine the structural stabilities and intermolecular interactions of diosgenin with IGF1R, MDM2, and SRC, an all-atom MD simulation in an explicit solvent was performed for 150 ns. We used the GROMACS software package, ver. 2021.3 (https://www.gromacs.org/) and Amber ff99SB-ildn force field to run MD simulations (Berendsen et al., 2005;Saeed et al., 2020;Dwivedi et al., 2022). The topological parameters of the ligands and the entire complex were calculated using the AmberTool's xleap module (https://ambermd.org/AmberTools.php) and the partial charges of the small molecules were calculated using an antechamber with a "bcc" charge model. The built systems were solvated in a rectangular box with 10.0 Å boundary conditions from the protein's borders in all directions using the TIP3P water model. The required amounts of counter ions were introduced to the prepared systems to neutralize the charges. The steepest descent and conjugate gradient energy reduction methods were used to discover the least energy conformations of the nearglobal state. To equilibrate the systems for 1 ns, "Canonical (NVT) and isobaric (NPT) ensembles" were used. A modified Berendsen thermostat approach was utilized to keep the volume and temperature consistent during NVT equilibration (300 K). During NPT equilibration, a Parrinello-Rahman barostat was used to keep the pressure constant at 1 bar. Furthermore, the Particle Mesh Ewald approximation with a cut-off value of 1 nm was used to calculate the long-range electrostatic interactions, van der Waals, and Coulomb interactions. A similar LINear Constraint Solver method was used to constrain bond length. Every complex went through a 150 ns production run with coordinates recorded every 2 fs. Other software programs, in addition to the built-in gromacs tools, were utilized to perform specific analyses on the acquired trajectories as required.
Analysis of free binding energy of complexes utilizing molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) In MD simulations and thermodynamic calculations, the relative binding energy of a ligand-protein complex was employed to investigate the binding free energies. The relative binding free energy and its contribution to individual residues were calculated using the MM-PBSA method and the "g_ mmpbsa" tool (Kumari et al., 2014;Dwivedi et al., 2022). The parameters for binding free energy calculations were taken from our previous study . The binding free energy (ΔG) was calculated using 50 representative snapshots taken throughout the stable trajectory observed between 100 and 150 ns. The change in entropy (ΔS) was calculated using the Schlitter formula and finally, accurate binding free energy was calculated using the formula, ΔG = ΔH-TΔS, where, ΔG = Gibbs free energy, ΔH = enthalpy change, T = temperature (Kelvin), and ΔS = entropy change.
Principal component (PCA) and dynamic crosscorrelation matrix (DCCM) analysis PCA was performed over the stable MD trajectory to examine various forms of molecular motion (Amadei et al., 1993;Amadei et al., 1996;Bhandare and Ramaswamy 2018). To accomplish this, the "least square fit" to the reference structure was used to account for the molecular translational and rotational motion. The collection of eigenvectors obtained by diagonalizing a covariance matrix was produced by a linear transform of cartesian coordinate space to reflect the direction of the molecular motion. The energy contribution of each eigenvector to the motion was presented by the eigenvalue associated with the respective eigenvector. The "timedependent motions" that the components carry out in an atomic vibrational mode were demonstrated by projecting the trajectory onto a particular eigenvector. The atomic vibrational components' contribution to this form of coordinated motion was shown by the projection's time average (Amadei et al., 1993;Van Aalten et al., 1995;Amadei et al., 1996;Bhandare and Ramaswamy 2018).
DCCM evaluates the magnitude of each pairwise crosscorrelation coefficient to determine whether or not atomic pair motion is correlated i.e., positive or negative (Khanal et al., 2021). We examined each DCCM component in this section, where Cij = 1 denotes the same period and phase (positive correlation), Cij = 0 indicates a lack of correlation, and Cij = −1 indicates a negative correlation between the fluctuations of i and j (Khanal et al., 2021).

Gene expression analysis in tumor, normal and metastatic tissues
Herein, we evaluated the top three genes based on the logrank test from survival analysis, and in silico molecular docking was evaluated for the gene expression in normal, tumor, and metastatic tissues using RNA sequence data in the platform of tumor, normal and metastatic samples. These data were analyzed using Kruskal Wallis and Dunn test.

Experimental pharmacology
Through different in vitro pharmacology profiling, we investigated the influence of diosgenin on breast cancer cell Frontiers in Pharmacology frontiersin.org 04 lines. To begin, brine shrimp lethality (BSL) bioassay was used to assess its cytotoxicity and the effect was compared with doxorubicin. Later, the effect of diosgenin on breast cancer cell lines (MCF7, MDA-MB-231, SKBR3, and T47D) was compared to normal epithelial cell lines (Vero) by emphasizing on cell viability, proliferation, and Warburg effect.

BSL bioassay
The brine shrimp lethality bioassay was performed as explained by McLaughlin et al. (1998) with minor modifications. Here, Artemia salina Leach. eggs from Seamonk international Artemia cyst 003 were used for the assay. Briefly, 10-12 brine shrimps were incubated within the different concentrations of diosgenin and doxorubicin (prepared in seawater) for 24 h. Controls were used without the test agents. After 24 h, the survived shrimps were counted and the % cytotoxicity was calculated as % cytotoxicity Total shrimps added − live shrimps Total shrimps × 100 The LC 50 was calculated using a linear regression curve.

In vitro MTT assay tumor and non-tumor cell lines
The cytotoxic activity of diosgenin and doxorubicin on tumor and normal cell lines was performed using an MTT assay (Mosmann, 1983) with minor modifications. Briefly, cell lines were plated onto 96-well flat-bottom plates, maintaining the cell density (20,000 cells/well), and were allowed to proliferate (24 h). After that, the cells were treated with different concentrations of diosgenin and doxorubicin maintaining the final volume of 250 µL after adding DMEM media (supplemented with 3% FBS) and incubated (37°C, 48 h, 5% CO 2 ). Next, 20 μL of MTT reagent was added and incubated (37°C in 5% CO 2 , 4 h). After incubation, the wells were washed (PBS, 3X) to discard the MTT. Then, formazan crystals were dissolved in DMSO (99.5% v/v, 100 μL) by gentle shaking. The absorbance was then recorded (570 nm) using an ELISA plate reader. The cell viability was calculated as % viability Absorbance of control − Absorbance of sample Absorbance of control × 100 In vitro scratch assay In vitro scratch assay was performed as explained by Bolla et al (2019) with minor modifications. Briefly, a stock solution (200 μg/ml) of diosgenin and doxorubicin was prepared and sterilized by filtering using a sterile membrane filter (0.22 µm). Later the solution was diluted using geometric series up to 1.56 μg/ml. This series of concentrations were chosen based on the number of experiments (trial and error). All the cells were seeded (2 × 10 5 cells/well) in 12-well tissue culture plates to obtain the confluence of 80-90% after 24 h of culture. After 24 h of post-seeding, the cell monolayers were scraped to create scratches of (300 µm). The detached cells and debris were washed with phosphate buffer. The media containing the samples were added to each well. Suitable controls were used by adding the minimal media and the scratch coverage was recorded at 0, 12, 24, 48, and 72 h after sample addition. The percentage scratch coverage was calculated using the following formula % scratch coverage Scratch length at 0 min − t min ( ) Scratch length at 0 min × 100

Effect of diosgenin on Warburg effect
The effect of diosgenin on the Warburg phenomena was evaluated by evaluating glucose uptake via the above-mentioned cancer cell lines vs. normal cell lines (Vero). Initially, the cells were grown in six well plates and incubated (37°C, 48 h) in a CO 2 incubator. After the formation of the confluent monolayer, the culture was renewed (DMEM consisting of 0.2% BSA) and again incubated (37°C, 18 h) in the CO 2 incubator. After incubation, the media was discarded and washed with KRP buffer. The cells were then treated with diosgenin and metformin in the presence of insulin followed by the addition of glucose (1 M) and incubated (30 min). The remaining amount of glucose was quantified from the supernatant. The percentage glucose uptake was calculated as the difference between the initial and final glucose content in the incubated medium (Gupta et al., 2009).

Statistical analysis
For the enrichment analysis, the whole-genome statistical background was used. The inhibitory constant (IC 50 ) and effective concentration (EC 50 ) were calculated using linear regression in GraphPad Prism (https://www.graphpad.com/) ver 5.0.

Computational pharmacology
Diosgenin ADMET profile and its targets Diosgenin has molecular weight of 414.62 g/mol with 3 H-bond acceptors and 1 H-bond donor. Diosgenin was predicted to possess the high human intestinal absorption (>30). In addition, it also showed the plasma protein binding, fraction unbound in plasma, and its volume of distribution 97.743%, 1.872%, and 1.695 L/kg, respectively. It was predicted as no inhibitory action towards CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 and showed Frontiers in Pharmacology frontiersin.org 05 high clearance rate of 23.332 ml/min/kg. In addition, it was not predicted for any side effects based on the ADVERPred server projection.
Diosgenin was predicted to regulate 105 targets in SwissTargetPrediction and four in DIGEP-Pred and BindingDB. The predicted targets in BindingDB were in common with SwissTargetPrediction. Similarly, a total of 510 different targets were recorded in DisGeNET for the neoplastic process of oestrogen receptor-positive breast cancer targets (UMLS CUI: C2938924). In general, the Diosgenin targeted 21 different proteins involved in breast cancer compared to recorded targets in DisGeNET (UMLS CUI: C2938924); Supplementary Figure S1.

Gene set enrichment analysis
The interaction of diosgenin-targeted 21 proteins had a total of 57 edges, 5.43 average node degree, 0.701 average local clustering coefficient, 17 edges (expected), and 3.83e-14 enrichment p-value ( Figure 2).

FIGURE 2
Protein-protein interaction of the diosgenin-triggered protein. Node color; colored nodes: query proteins and first shell of interactors, white nodes: second shell of interactors, Node content; empty nodes: proteins of unknown 3D structure, filled nodes: some 3D structure is known or predicted, Known Interactions; from curated databases, experimentally determined, Predicted Interactions; gene neighborhood, gene fusions, gene co-occurrence and Others; text mining, coexpression, protein homology.

FIGURE 3
Network interaction of diosgenin-regulated proteins and respective pathways. Pink color presents a low edge count and green color presents a high edge count. A large node size presents a higher edge count and a smaller edge count presents a lower edge count.
Frontiers in Pharmacology frontiersin.org 07 GO analysis GO analysis identified 12 GO terms for cellular components in which receptor complex; GO:0043235 had the lowest false discovery rate i.e., 0.00025 regulated seven genes i.e., PDGFRB, IGF1R, APP, RET, GRM1, FGFR2, and LYN against 381 genes at 1.23 strength. Also, a total of 21 genes were triggered in multiple cellular components in which tyrosine-protein kinase Lyn (LYN) was majorly triggered in 11 cellular components except for cytoplasmic vesicles; GO:0031410 (Supplementary File S1; Supplementary Sheet S1). In addition, the Pearson p-value for cellular components was 0.795 strength vs. false discovery rate, 0.0003 strength vs. observed gene count, 0.8093 false discovery rate vs. observed gene count, 2.684e-004 observed gene count vs. strength, and 0.809 observed gene count vs. false discovery rate. In addition, the minimum observed Pearson r was −0.084 and the maximum was 1.000 (Supplementary Figure S3). The proteins modulated by the diosgenin in the different cellular compartments are presented in Figure 5.
Likewise, 27 different molecular functions GO terms were identified via the PPI in which protein tyrosine kinase activity; GO:0004713 had the lowest false discovery rate i.e., 3.53E-05 to trigger six genes i.e., PDGFRB, IGF1R, RET, SRC, FGFR2, and LYN against 137 background genes at 1.61 strength. Herein, a total of 22 genes were triggered in 27 molecular functions in which protooncogene tyrosine-protein kinase Src (SRC) was involved in 18 molecular functions except hormone binding; GO:0042562, transmembrane receptor protein tyrosine kinase activity; GO: 0004714, transition metal ion binding; GO:0046914, steroid hydroxylase activity; GO:0008395, oxygen binding; GO:0019825, nuclear receptor activity; GO:0004879, peptide hormone binding; GO:0017046, platelet-derived growth factor receptor binding; GO: 0005161, and steroid-binding; GO:0005496 (Supplementary File S1; Supplementary Sheet S3). Herein, the correlation p values were 0.130 strength vs. false discovery rate, 2.0510e-009 strength vs. observed gene count, and 0.753 false discovery rate vs. observed gene count. Within the interaction of strength, false discovery rate, and observed gene count, −0.877 and 1.000 were minimum and maximum r values for molecular function (Supplementary Figure S3).
Similarly, a total of 307 biological processes was traced due to PPI in which cellular response to oxygen-containing compound; GO:1901701 had the lowest false discovery rate i.e., 2.72E-09 in modulating 14 genes i.e., NR3C1, CDK4, MDM2, PDGFRB, STAT3, IGF1R, APP, RET, PTPN1, SRC, AR, CRHR1, FGFR2, and LYN against 1,055 background proteins at 1.09 strength. In addition, a total of 21 genes were triggered for 307 biological processes via the 21 proteins associated with breast cancer in which SRC was involved and linked with 213 biological processes. In addition, the order of triggered genes in multiple pathways was as CYP17A1 (17) The detailed biological processes of participation of the above genes are supplemented in Supplementary File S1; Supplementary Sheet S4. The Pearson correlation p-value was found to be 0.034 strength vs. false discovery rate,

FIGURE 4
The associated protein-pathways interaction with respective false discovery rate and strength.
Frontiers in Pharmacology frontiersin.org 08 1.0622e-063 strength vs. observed gene count, and 8.1767e-017 false discovery rate vs. observed gene count. In addition, the minimum Pearson r values for biological processes were -0.779 and 1.000 respectively (Supplementary Figure S3).

Compartments enrichment analysis
The PPI reflected the modulation of proteins in 16 different compartments in which the intrinsic component of plasma membrane; GOCC:0031226 had the lowest false discovery rate i.e., 0.00027 in regulating nine genes i.e., PDGFRB, STAT3, IGF1R, APP, RET, GRM1, CRHR1, FGFR2, and LYN against 841 background genes at 1 strength. Herein, a total 21 were triggered in 16 compartments in which LYN was traced in all; the order of genes regulation in multiple enriched compartments was

InterPro enrichment analysis
Protein-protein enrichment analysis traced seven InterPro in which tyrosine-protein kinase, catalytic domain; IPR020635 was Frontiers in Pharmacology frontiersin.org 09 traced with minimum false discovery rate i.e., 8.42E-06 to regulate six genes i.e., PDGFRB, IGF1R, RET, SRC, FGFR2, and LYN against 79 background genes at 1.85 strength. Herein, in the enrichment of InterPro analysis, a total of nine genes were triggered in order as MDM2 (1)

Tissues enrichment analysis
In cluster 1, the interaction of seven nodes corresponded to 11 edges with 3.14 average node degree, 0.8 average local clustering coefficients, 1 expected edge, and 1.07e-09 PPI enrichment p-value. Herein, a total of nine molecular functions were traced in which steroid hydroxylase activity (GO:0008395) and oxygen binding (GO:0019825) scored the lowest false discovery rate i.e., 0.0013 via the regulation of three genes i.e., CYP3A4, CYP17A1, and CYP19A1 against 36 background genes at 2.37 strength. Similarly, Steroid binding (GO:0005496) showed 0.0082 false discovery rate and regulated three genes i.e., NR3C1, CYP3A4, and AR against 104 background genes at a strength of 1.91. Similarly, two KEGG pathways i.e., steroid hormone biosynthesis (hsa00140) and ovarian steroidogenesis (hsa04913) were traced. Herein, steroid hormone biosynthesis was identified with a false discovery rate of 0.00048 to regulate three genes i.e., CYP3A4, CYP17A1, and CYP19A1 against 59 background genes at 2.15 strength. In addition, ovarian steroidogenesis regulated two genes i.e., CYP17A1 and CYP19A1 against 50 genes at 2.05 strength and 0.0327 false discovery rate. Also, two different biological processes i.e., androgen metabolic process (GO:0008209) and organic cyclic compound biosynthetic process (GO:1901362) were identified to regulate three genes i.e., CYP3A4, CYP17A1, and CYP19A1 against 27 background genes at 2.47 strength and 0.002 false discovery rate and six genes (NR3C1, FASN, CYP3A4, CYP17A1, AR, and CYP19A1I) against 1,211 background genes at 1.14 strength and 0.0033 respectively (Supplementary File S2; Supplementary Sheet S1).
In cluster 3 the interaction of the nodes had 1 edge count, average node degree, and average local clustering coefficient with 0.0398 PPI enrichment p-values. Here, only two KEGG pathways and eight biological processes were traced. In KEGG pathways, the p53 signaling pathway; hsa04115, and microRNAs in cancer; hsa05206 were associated with two genes i.e., MDM2 and MDM4 against 72 and 160 background genes, 2.43 and 2.09 strength and 0.0064 and 0.0155 false discovery rate respectively.

Diosgenin-targets-protein network analysis
The combined interaction between the diosgenin, its targets, and regulated pathways traced IGF1R, MDM2, SRC, CDK4, and PDGFRB as the top five lead hub proteins. In addition, pathways in cancer; hsa05200, EGFR tyrosine kinase inhibitor resistance; hsa01521 prostate cancer; hsa05215, viral carcinogenesis; hsa05203, and PI3K-Akt signaling pathway; hsa04151 were traced as the top five lead hub pathways modulated within diosgenin-targets-pathways interactions.

FIGURE 6
Edgebetweenness of network interaction and node interaction count. The highest edge betweenness was 51 (1 interaction) and the lowest was 1 (125 interactions).

Stability of diosgenin-IGF1R complex
The diosgenin-IGF1R complex showed stable dynamics up to 150 ns after a 20 ns equilibration phase. Initial backbone and complex RMSD values climbed steadily increased from 1.0 Å to 3.1 Å and~1.4 Å to~3.7 Å, respectively from 0 to 20 ns After 20 ns, it was discovered that the backbone and complex RMSD (~2.5 Å and 3.0 Å, respectively) were stabilized with lesser fluctuations ( Figure 8A). The loopforming residues Leu1064 to Pro1077 showed comparatively greater fluctuations (7.0 Å). On the other hand, residues Leu975, Val983, Ala1001, Glu1050, Asp1123, and Ile1130 that interacted with diosgenin during docking studies didn't exhibit fluctuation because they were involved in stable non-bonded interactions. Additionally, it was discovered that residues Gly1122 to Tyr1131 forming the loop region were involved in ligand binding and showed the least RMS fluctuation (~2.0 Å) ( Figure 8B, Supplementary Movie S1). There was a formation of a compact globular shape, which was supported by a gradual drop in Rg value from 20.5 Å to 19.6 Å and was further found to be stable at 20.0 Å ( Figure 8C). The initial and final surface  Figure 8D). The complex formed 3 H-bonds of which two were consistent during the simulation ( Figure 8E). It was discovered that diosgenin and IGF1R had relative binding energy of −35.143 ± 3.03 kcal/mol. Supplementary Table S2 summarizes the free energy contribution of diosgenin with IGF1R, MDM2, and SRC. The per residue contribution energy revealed that 16 residues from the binding pocket i.e. Glu974, Leu975, Val983, Tyr984, Glu985, Ala1001, Glu1020, Val1033, Leu1051, Met1052, Asp1056, Met1112, Asp1123, Met1126, Ile1130, and Tyr1131 significantly contributed to the formation of a stable complex. These residues also scored the least per residue decomposition/contribution energy which ranged from −2.0 kJ/mol to −6.8 kJ/mol whereas the positive contribution energy of 4.0 kJ/mol was achieved by residues Lys1003, Lys1058, and Arg1109 ( Figure 8F).

Stability of diosgenin-MDM2 complex
Throughout a 150 ns production run, the Diosgenin-MDM2 complex exhibited stable dynamics. The initial and final backbone RMSD values were 0.62 and 1.21, respectively, and an average was 1.40 Å. Similar to this, the initial and final RMSD of the docked complex MDM2 and diosgenin were 0.92 and 2.13 respectively with an average of~2.30 Å ( Figure 9A). The loop forming Met17 to Ser22 residues at the N-terminus had substantially greater fluctuations (~4.2 Å). Further, residues interacting with diosgenin during docking studies (Leu57, Tyr67, Phe91, Val93, and Ile99) did not exhibit variations as they were involved in stable non-bonded interactions ( Figure 9B and Supplementary Movie S2). By monitoring a stable Rg, a more compact and stable complex was formed. The complex's initial and final Rg values were determined to be 13.2 Å and 1.31 Å. Similarly, the initial and final surface area occupied by IGF1R and diosgenin docked complex was 63.90 nm 2 and 61.57 nm 2 . The complex had an average surface area of 62.89 nm 2 ( Figure 9D). Two H-bonds were established by the complex of which one was consistent during MD simulation ( Figure 9E). The relative binding energy between diosgenin and IGF1R was discovered to be −34.619 ± 2.81 kcal/mol. Further, seven residues (Leu57, Gly58, Ile61, Met62, Tyr67, Val93, and Ile99) from the binding pocket scored the lowest per residue decomposition/ contribution energy, ranging from −3.65 kJ/mol to −8.87 kJ/ mol. These residues contributed significantly to forming the stable complex according to the per residue decomposition/ contribution energy. The residue Glu69 had positive contribution energy of 2.68 kJ/mol ( Figure 9F).

Stability of diosgenin-SRC complex
Another diosgenin-SRC complex also showed stable dynamics with the RMSD value of < 2Å after the equilibration period of 50 ns( Figure 10A). In addition, the RMSD of the complex showed a sharp increase in its values after 40 ns to 8 Å ( Figure 10A). The careful observation of the entire complex trajectory and representative snapshots extracted from the region reveal the structural transition of diosgenin from its primary binding pocket to the neighboring alternate pocket which formed a stable complex during the entire simulation period. Thus, for the first time, we report the existence of an alternate binding site other than the primary binding site on SRC. From the simulation movie (Supplementary Movie S3) it was observed that the complex undergoes major conformational changes in the Frontiers in Pharmacology frontiersin.org 13 secondary structure of the protein including the primary binding pocket, leading to the increased conformational flexibility of SRC, hence the diosgenin switches its original position from the primary binding pocket to the reported alternate binding pocket (formed by the residues Tyr152 to Leu164). The residual fluctuations observed show fewer C-alpha fluctuations (<2.2Å), interestingly the residues from the primary binding pocket as well as new reported alternate binding pocket show the least RMSF values~0.5Å ( Figure 10B), mainly due to the stable non-bonded interactions shown by them. The complex diosgenin-SRC gains a compact globular shape during the simulation, thus we propose the stable complex formation was favored after the equilibration period of MD simulation ( Figure 10C). The solvent-accessible surface area represents an exposure of the hydrophobic residues to the solvent, here in this complex SASA exhibited a similar trend as that of Rg values. SASA value decreased gradually till the equilibration state and further, it stabilized till the simulation end ( Figure 10D), suggesting the proper folding of the hydrophobic core including both the binding sites. Two H-bonds were established by the complex of which one was consistent during MD simulation ( Figure 10E). The relative binding energy between diosgenin to SRC was discovered to be −17.994 ± 5.67 kcal/mol. In addition, three residues i.e., Ile156, Thr157, and Leu164 from the binding pocket scored the lowest energy contribution per residue. These residues contributed significantly to forming the stable complex according to the per residue decomposition/ contribution energy. The residue Ile156 had the least contribution energy of −8.9 kJ/mol ( Figure 10F).

Principal component and dynamic crosscorrelation matrix
PCA is a statistical technique being used to study the dynamics of bimolecular complexes as it limits the 3 N (N = number of atoms in the protein) degrees of freedom describing functionally crucial motions of the protein. It was observed that in all the complexes maximum dynamics during the simulation have been captured by the first 10 eigenvectors, of which the first two contributed significantly to the collective motions exerted by all the simulated complexes ( Figures 14A-C). Hence, we examined the collective motion sampled by the first two principal components (PCs), and 2D projections for PC1 and PC2 were plotted ( Figures 11D-F). The complexes of diosgenin with MDM2 and SRC express the compact clusters in the conformational spaces those range from -1.5 to 1.5. In the MD trajectory of complexes Diosgenin-SRC and Diosgenin-MDM2, PC1 and PC2 (top two modes) showed the uniform distribution across the configurational space while the remaining complex of Diosgenin-IGF1R showed a large diversity in the conformational space and was widely clustered in the range of -4.5 to 4.5 ( Figure 11D). MD trajectory sampled three states of the protein as seen by the three individual clusters in the scatterplot of PC1 v/s PC2. Herein, we propose that the Diosgenin-IGF1R complex has undergone significant conformational changes in the secondary structure during the simulation those favored in forming a stable complex. However, other complexes namely Diosgenin-SRC and Diosgenin-MDM2 were well stabilized and undergone comparatively lesser conformational changes in the secondary structure hence exhibited compact cluster in the conformations space. Further, the convergence of sampling was also analyzed by calculating the   Frontiers in Pharmacology frontiersin.org 15 close to 1 indicated bad sampling of the trajectory pointing to all our simulation trajectories being well converged and properly sampled in the free energy landscape.
The concerted motion exerted by the three complexes during the simulation was examined using DCCM (dynamic crosscorrelation matrix). The calculation of the correlation matrix is utilized to depict the dynamical information of proteins in two dimensions. To observe the correlation in the dynamics of the binding site correlation matrix over the stable trajectory of all the complexes were plotted (Figures 11F,G). The diagonal orangered line indicates the self-correlation of the individual residues with themselves. The orange-red region in the correlation map signifies the concerted movement of the residues in the same direction whereas the dark blue region represents anti-correlated fluctuations. The N-and C-terminal region of the complex Diosgenin-MDM2 represents strong anticooperative movement with each other. Comparatively, this complex exhibited a maximum amplitude of negative correlation across all the residues. The maximum region in the complex Diosgenin-IGF1R showed positive correlation mainly at the binding pocket region (residues Glu974, Leu975, Val983, Tyr984, Glu985, Ala1001, Glu1020, Val1033, Leu1051, Met1052, Asp1056, Met1112, Asp1123, Met1126, Ile1130, and Tyr1131). This signifies the closure movement observed at the binding pocket region facilitating stable complex formation (Supplementary Movie S1). Similarly, in the Diosgenin-SRC complex, the binding pocket residues (Arg158, Glu181, Thr182, Cys188, Lys203-Lys206) showed moderately positive cooperative motion compared to other regions in the SRC structures. Interestingly, the residues from the alternate binding region show the cooperative motion results in stable non-bonded contact with the newly reported alternate binding pocket.
In general, the conformational flexibility of all the complexes i.e., Diosgenin-IGF1R, Diosgenin-MDM2, and Diosgenin-SRC varies greatly as observed in the DCCM plot and the collective dynamics nature observed in the PCA and DCCM plot favors the stable complex formation during the simulation.

In vitro pharmacology BSL bioassay
Exposure to the different concentrations of diosgenin and doxorubicin showed concentration-dependent brine shrimp lethality. In addition, the LC 50 was found to be 19.15 and 0.71 μg/ml respectively reflecting the doxorubicin to be 27.06 times more potent than diosgenin (Supplementary Figure S4).

In vitro MTT assay
Here, the MTT assay showed the IC 50 of the diosgenin to be significantly higher than that of the doxorubicin in MCF7 (p < 0.001), MDA-MB-231 (p < 0.05), SKBR3 (p < 0.05) and Vero (p < 0.001) compared to the doxorubicin. In vitro MTT assay on MCF7 cell lines reflected the doxorubicin (IC 50 3.21 ± 0.29 μg/ ml) to be 3.7 times more potent than the diosgenin (IC 50 12.05 ± 1.33 μg/ml). In addition, over the MDA-MB-231 cell lines, the IC 50 of diosgenin was found to be (45.54 ± 23.41) µg/ml compared to doxorubicin (6.30 ± 2.67) µg/ml which suggests the doxorubicin to be 7.2 times more potent than diosgenin. In addition, the diosgenin and doxorubicin had the IC 50 (15.11 ± 5.32) µg/ml and (5.47 ± 1.09) µg/ml respectively over SKBR3 reflecting the doxorubicin to be 2.76 times more potent than that of diosgenin. In addition, the IC 50 of diosgenin and doxorubicin was found to be (17.78 ± 7.86) Frontiers in Pharmacology frontiersin.org 16 and (8.18 ± 8.18) µg/ml respectively over the T47D cell lines to point doxorubicin to be 2.17 times more potent. Further, diosgenin and doxorubicin had the IC 50 of (38.59 ± 4.03) and (7.05 ± 0.69) µg/ml respectively in which the doxorubicin was found to be 7 times more cytotoxic than Vero cell lines; Supplementary Figure S5; Supplementary Table S4.

In vitro scratch assay
After the 72-h exposure to different cells, it was observed that diosgenin (200 μg/ml) had the highest effect on SKBR3 cell lines (18.22 ± 1.237%) to prevent cell migration compared to the rest. In addition, we observed a significant difference in percentage scratch closure within MCF vs. SKBR3 (p < 0.01), MDA-MB-231 vs. SKBR3 (p < 0.05), SKBR3 vs. T47D (p < 0.01) and SKBR3 vs. Vero (p < 0.001) cell lines (Supplementary Figure S6). In addition, the concentration and time-dependent effect of the diosgenin on scratch closure are presented in Supplementary Figure S7.

Effect of diosgenin on Warburg effect
It was observed that diosgenin had efficacy in dealing glucose uptake in tumor cells. The EC 50 of the diosgenin was observed to be significantly higher in glucose uptake in MCF7 (p < 0.001), SKBR3 (p < 0.05), T47D (p < 0.001), and Vero (p < 0.01) cell lines compared to metformin in the presence of insulin. It was observed that diosgenin had a comparatively higher EC 50 (26.19 ± 2.77) µg/ml compared to metformin (3.10 ± 0.99) µg/ml which showed the diosgenin to possess a comparatively lower glucose uptake efficacy than metformin in

FIGURE 12
Differential gene expression analysis of PDGFRB, FBFR2, and STAT3 in tumor, normal and metastatic tissues. These genes were identified to possess the least log-rank test.

FIGURE 13
Differential gene expression analysis of IGF1R, MDM2, and RSC in tumor, normal and metastatic tissues. Diosgenin was predicted to possess the highest binding affinity with these targets.
Breast cancer is one of the main causes of death among women. Chiefly, the cells lining the milk-forming duct of the mammary glands are the origination of breast cancer (Herbein and Kumar, 2014) which can be further subdivided based on the presence or absence of the hormone receptors i.e., estrogen and progesterone subtypes and human epidermal growth factor receptor-2 (HER2). In addition, the estrogen receptor pathway triggers hormone receptor-positive breast cancer (Rani et al., 2019). Similarly, in HER2-positive breast cancer, HER2 triggers RAS/RAF/MAPK and PI3K/AKT signaling pathways that stimulate cell growth, survival, and differentiation (Dittrich et al., 2014). We discovered 21 diosgenin-regulated proteins that are involved with oestrogen receptor-positive breast cancer (UMLS CUI: C2938924) targets in this investigation. In addition, three pathways i.e., PI3K-Akt (modulated five genes i.e., CDK4, MDM2, PDGFRB, IGF1R, and FGFR2), Ras (modulated three genes i.e., PDGFRB, IGF1R, and FGFR2), and MAPK (modulated three genes i.e., PDGFRB, IGF1R, and FGFR2) signaling pathways associated to HER2 positive breast cancer were modulated. Since diosgenin modulated three pathways closely associated with HER2-linked pathogenesis, it can be speculated that diosgenin could act via the manipulation of HER2.
It has previously been proposed that Akt activation influences endocrine resistance in metastatic breast cancer. In addition, Akt activation in the downstream pathway of HER2 could resist the endocrine therapy of breast cancer (Tokunaga et al., 2006). Furthermore, Ras proteins activate the cytoplasm and extracellular signaling networks via receptor tyrosine kinase and are involved in cell proliferation, survival, growth, metabolism, motility, and apoptosis, and their hyperactivation promotes the growth and progression of breast cancer. In addition, Ras's intracellular localization, activation, and signaling have been used in developing therapeutic candidates against breast cancer via the enzymes involved in posttranslational modification of Ras e.g., farnesyltransferase and geranylgeranyltransferase 1 (Moon, 2021). Further, MAPK links the extracellular mitogenic signals to cell proliferation which may be concerned with or act independently towards estrogenmediated events in breast cancer cells (Yue et al., 2002). In the present study, we identified the modulation of the abovemodulated pathways i.e., PI3K-Akt signaling pathway (false discovery rate: 0.0011, strength: 1.12), Ras signaling pathway (false discovery rate: 0.0263, strength: 1.09), and MAPK signaling pathway (false discovery rate: 0.0446, strength: 0.99). This suggests the probability of involvement of these pathways linked to HER2 function which could be modulated by diosgenin in breast cancer.
The role of EGFR dysregulation or mutation in cancer etiology, particularly breast cancer, has been proposed previously. However, resistance towards EGFR tyrosine kinase inhibitors may occur due to secondary mutations (T790M), activation of secondary pathways (AXL, c-Met, HGF), aberrant downstream pathways (K-RAS mutations, loss of PTEN), deregulation of the EGFR tyrosine kinase-mediated apoptosis, histological transformation, and ATP binding cassette transporter effusion, etc (Huang and Fu 2015). The FoxO signaling system interacts with the PI3K-Akt signaling pathway and relates to cancer progression, particularly breast cancer advancement (Farhan et al., 2017). In addition, FoxO negatively regulates activated EGFR signaling which was demonstrated via the in vitro cell line culture method and in vivo models (Sangodkar et al., 2012). Here, in the present study, we identified the regulation of the FoxO signaling pathway via the modulation of four genes i.e. MDM2, STAT3, IGF1R, and GRM1 which could support the functioning of the PI3K-Akt and EGFR signal against breast cancer pathogenesis.
Frontiers in Pharmacology frontiersin.org 18 Furthermore, the Rap1 signal has been linked to tumor cell proliferation, invasion, and metastasis through regulating integrin-or cadherin-mediated cell function, cytoskeletal alterations, protease (metal metalloprotease) production, and cell adhesion (Zhang et al., 2017). In addition, Rap1 has been traced to attenuate metastasis and EGFR-triggered carcinoma (Huang et al., 2012). Since breast cancer progression is closely linked to EGFR tyrosine kinase signal; it can be speculated that diosgenin-mediated Rap1 signal could attenuate tumor invasion and metastasis; was observed to be modulated via the modulation of three genes i.e., PDGFRB, IGF1R, and FGFR2.
The p53 acts as a transcription factor for p21 via cyclin-CDK interactions which is important for the transition of the G2 phase to the mitosis phase (Alam et al., 2015). In addition, p21 protects cells from apoptosis, regulates the cell cycle, causes apoptosis, and decreases cell proliferation in tumor cells . Mutant p53 has been pointed to as the guardian of the cancer cells (Mantovani et al., 2019) and is also associated with worsening breast cancer affecting overall survival (Gasco et al., 2002). In the present study, we identified the regulation of the three genes i.e., CDK4, MDM2, and MDM4 which could activate the p53 signal against breast cancer. Furthermore, excessive plasma prolactin levels have been associated with an increased risk of breast cancer in both premenopausal and postmenopausal women; is more prominent in estrogen or progestogen receptor cancer type (Tworoger et al., 2004;Tworoger et al., 2006). Here, in the present study, we identified the diosgenin to regulate the prolactin signaling pathway via the regulation of three genes i.e., STAT3, CYP17A1, and SRC which could avoid estrogen or progestogen receptor-mediated breast cancer progression. In addition, chemokine signals are not limited to tissue differentiation, hematopoiesis, inflammation, and immune regulation but also process tumor development by triggering angiogenesis, metastasis, drug resistance, and immunity of breast cancer (Liu et al., 2020). In the present study, we identified the diosgenin to trigger the chemokine signaling pathway, and regulated three genes i.e., STAT3, SRC, and LYN at 1.18 strength and 0.0168 false discovery rate). In addition, in enrichment analysis i.e., tissues, Pfam, InterPro, and compartments enrichment analysis we observed the multiple proteins that are concerned with the breast cancer prognosis.
Also, in the present study, we identified three proteins i.e., IGF1R, MDM2, and SRC in diosgenin-protein(s)-pathway(s) interaction. Hence, these were further considered for post-network analysis i.e., molecular docking and molecular dynamics simulation. Apart from handling the transcription, IGF1R can trigger the growth and metastasis of malignant melanoma cells through the PI3K-Akt signaling pathway (Ekyalongo and Yee, 2017). Diosgenin was anticipated to interact with IGF1R in the current investigation, potentially preventing breast cancer metastasis and tumor invasion, which could be PI3K-Akt driven, as previously discussed.
The IGF1R is a transcription factor that binds to DNA and influence transcription. Both ERK1/2 and AKT are involved in the transcriptional control of the IGF1R gene. MicroRNA-139-5p modulates the growth and metastasis of malignant melanoma cells via the PI3K/AKT signaling pathway by binding to IGF1R (binding energy -8.6 kcal/mol). Previously, MDM2 amplification has been reported to relate to estrogen receptor status and its presence has been indicated in human breast cancer cell (Quesnel et al., 1994); was observed to be manipulated with diosgenin in the third cluster of the PPI binding (binding energy −8.5 kcal/ mol). Likewise, another modulated protein i.e., SRC has been reported to increase its expression in human breast cancer by 4-30 fold which was evidenced via both immunohistochemistry and immunoblotting (Verbeek et al., 1996); was also predicted to be modulated by diosgenin interaction (binding energy -7.4 kcal/ mol). In addition, since 17 genes i.e. APP, AR, CDK4, CYP19A1, FGFR2, GRM1, IGF1R, LYN, MDM2, MDM4, NR3C1, PDGFRB, PRCP, PTPN1, RET, SRC, and STAT3 were observed to have a significant role in disease prognosis, it can be speculated that the above-modulated proteins are primarily concerned with breast cancer management with diosgenin treatment.
The docking study revealed diosgenin to interact with active site residues of three potential targets involved in breast cancer via IGF1R, MDM2, and SRC. Diosgenin formed stable intermolecular interactions throughout 150 ns MD simulation revealing them as the best lead. Among the interactions of diosgenin with IGF1R, diosgenin interactions with Leu975, Val983, Met1112, Met1126, and Ile1130 were consistent in both docking and MD simulation. Multiple studies have demonstrated these residues involve the pocket as a primary binding site (also validated by PrankWeb server; https:// prankweb.cz/) for inhibition of IGF1R (Munshi et al., 2002;Li et al., 2009;Guo et al., 2015). This indicates, that diosgenin as a potent lead hit against IGF1R. Similarly, diosgenin scored the lowest binding energy of -8.5 kcal/mol and binding free energy of -34.619 kcal/mol and formed stable interactions Leu57, Gly58, Ile61, Met62, Tyr67, Val93, and Ile99 throughout the 150 ns MD production run. Both docking and MD simulation revealed diosgenin as a potent lead hit for targeting MDM2. Interestingly, a study by Li et al. (2021) identified AG-690/ 37072075 and AO-022/43452814 as potent anticancer lead hits against MDM2. These molecules were predicted to interact with the residues "Leu57, Gly58, Ile61, Met62, Tyr67, Val93, and Ile99" and were also found to inhibit p53-MDM2 interaction in wild-type p53 cells. Hence, we believe that diosgenin may interfere with the p53-MDM2 interaction (MDM2 inhibits the transcriptional activity of p53 by attaching to its transactivation domain), and further experimental studies are required to validate our findings. Further, diosgenin scored -7.4 kcal/mol binding energy against SRC and formed interaction with Ala168, Arg172, Phe194, and Leu200. The structural shift of diosgenin from its initial binding pocket to the neighboring alternative pocket was revealed through MD simulation, and the complex remained stable throughout the entire simulation. As a result, we are the first to disclose the discovery of a secondary binding site in the SRC. The simulation movie (Supplementary Movie S3) shows that the complex undergoes significant conformational changes in the protein's secondary structure, including the primary binding pocket, increasing the conformational flexibility of SRC and causing diosgenin to switch positions from the primary binding pocket to the alternate binding pocket that has been reported. The PCA and DCCM revealed that diosgenin with IGF1R, MDM2, and SRC exhibited significant differences in conformational flexibility and also support the stable complex formation during the simulation.
During the differential gene expression analysis of diosgeninregulated genes, we observed that PDGFRB, FBFR2, and STAT3 possess maximum gene expression in tumor and metastatic vs. normal tissues. In addition, molecular docking and simulation also identified IGF1R, MDM2, and SRC could be the prime diosgenin-modulated targets against breast cancer. Hence, to confirm this we further assessed the functional role of each target using cell line studies. The PDGF family (PDFGRB) has been reported to use PDGF ligands released by cancer stromal cells from breast cancer cells to drive cell proliferation (Farooqi and FIGURE 14 Mechanism of action of diosgenin against oestrogen receptor-positive breast cancer. Diosgenin acts on the two cell surface proteins IGF1R and PDGFRB and inhibits invasion, metastasis, angiogenesis, and cell proliferation. In addition, it acts on the four cytoplasmic proteins i.e., MDM2, STAT3, FGFR2, and SRC. Diosgenin may inhibit the cell to escape p53 surveillance by binding with MDM2, inhibit cell proliferation and promote apoptosis by binding to STAT3, interfere with NF-kB signals check the cancer cell self-renewal, and inhibit the metastatic spread by inhibiting the SRC action. Frontiers in Pharmacology frontiersin.org 20 Siddik, 2015). Similarly, FBFR2 has been reported to promoting cell self-renewal by interacting with NF-kB signals (Kim et al., 2013). Likewise in previous studies, STAT3 inhibition has been reported to decrease cell proliferation and apoptosis promotion in various cancers including breast cancer (Kanai et al., 2003;Pancotti et al., 2012;Chen et al., 2013). In this regard, because these targets were directly triggered by diosgenin and were also significantly higher in the tumor than in the normal, diosgenin may limit cell proliferation and promote apoptosis by interfering with the self-renewal process. This hypothesis was further supported by the MTT and scratch assays. Previously, an increase in glucose uptake has been reported in the cancer cell which supports ATP production and acts as a fuel (Adekola et al., 2012). This glucose uptake can be reduced by blocking insulin's impact on cancer cells. In comparison to metformin, the current investigation found a considerably greater effective concentration to limit glucose uptake. As a result, diosgenin may inhibit glucose uptake in cancer cells. In the current investigation, however, the glucose absorption assay was done in the presence of insulin. As a consequence, more research is needed to determine its role in glucose uptake in the absence of insulin. Similarly, MDM2 helps cancer cells to escape p53 surveillance and avoid cellular apoptosis (Momand et al., 1998). In the present study, since, diosgenin had a significant binding affinity with MDM2, it could probably promote the apoptosis that needs to be further confirmed. In addition, SRC has been indicated for cancer metastasis and also helps in cancer progression and development which is an indicator of cell proliferation (Wheeler et al., 2009). In the current investigation, we discovered that diosgenin inhibited cell proliferation in a variety of cell lines. This may be due to the blockage of the SRC generated by diosgenin by binding to it. This could have prevented cell multiplication, as revealed by the scratch assay. Based on the information presented above, it is

FIGURE 15
Probable checkpoints affected by the diosgenin in the pathogenesis of breast cancer (Homo sapiens (human); hsa05224) targets/ pathways. In the KEGG pathway analysis, it was diosgenin-regulated targets involved in breast cancer were also identified to trigger EGFR, PI3K-Akt, MAPK, and p53 signaling pathways.
Frontiers in Pharmacology frontiersin.org 21 reasonable to believe that diosgenin may act against breast cancer by reducing glucose utilization, invasion, and metastasis via the surface protein IGF1R. Furthermore, it may limit cancer cell proliferation by acting on another surface protein, PDGFRB. In addition, the diosgenin may also act over the MDM2 and prevent the cancer cell to escape the p53 surveillance in the cytoplasm. Also, diosgenin may primarily act over two cytoplasmic targets FGFR2 and SRC to prevent the interaction with NF-kB (prevents cell self-renewal) and metastatic spread respectively ( Figure 14) and may also modulate other pathways within the breast cancer pathogenesis (Figure 15) which was evidenced during KEGG pathway analysis.
Although diosgenin inhibited cell growth, was cytotoxic, and had an effect on the various cellular compartments of the tumor cell, the inhibitory constant efficacy was found to be greater than that of doxorubicin. In addition, the previous drug discovery process utilized the concept of the "lock and key" approach in which a designed drug is specific to a single protein e.g. doxorubicin complexes with DNA by intercalation and inhibits topoisomerase II. However, this approach has often failed many times and molecules were potent cytotoxic to normal cells which was also observed in the present study. However, a single compound preferably from a natural source tends to act through a polypharmacology approach in which a single molecule can target multiple proteins based on the concept "master key can unlock multiple locks" to target multiple proteins and pathways (Chandran et al., 2017). If it happens, the amount of concentration required is a bit high; however, approaches should be made through the targeted drug delivery to breast cancer which is yet to be studied. This issue could be remedied by increasing diosgenin's cellular permeability and promoting anti-cancer action against breast tumors via a novel drug delivery mechanism, as established in previous studies (Dhamecha et al., 2015;Jagwani et al., 2020;Dhamecha et al., 2021;Ramasamy et al., 2021). In addition, the effect of diosgenin on the hub genes (SRC, MDM2, and IGF1R) expression is based on the computational models which need to be further evaluated using real-time or reverse transcriptase polymerase chain reaction even though the present study pointed their functional effect via the glucose uptake, cell proliferation, and apoptosis which is the perspective of the present findings.

Conclusion
The present study utilized a series of system biology tools to trace the potential action of diosgenin against breast cancer. Herein, we identified the probable action of the diosgenin against breast cancer via FoxO, PI3K-Akt, p53, Ras, and MAPK signaling pathways. In addition, we traced the selectivity of the diosgenin to manipulate the action of three hub genes i.e., IGF1R, MDM2, and SRC. Our molecular modeling study reveals that the stable complex formation is primarily facilitated by the cooperative closure motion exerted by the primary binding pocket residues in diosgenin-MDM2, and diosgenin-IGF1R complex while, in Diosgenin-SRC similar compact closure dynamics are also observed at the alternate binding pocket. Conformational flexibility and convergence of the trajectories during MD have been investigated using PCA. Diosgenin showed stable non-bonded interactions forming stable binary complexes with all three screened targets namely SRC, MDM2, and IGF1R. Thus, we proposed that these stable interactions of Diosgenin would trigger the successful inhibition of SRC, MDM2, and IGF1R; the newly identified targets in breast cancer. The results of our computer modeling investigations agree well with the results of the experimental cell line tests. As a result, we hope that this work will open the way for the development of novel therapeutic techniques and/or medication candidates for breast cancer.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions
PK: Conceptualized, experimented, drafted, edited, reviewed, and supervised the entire work. VP: Equal contribution in experimentation, drafting, editing, and reviewing the draft. VVB, PPP, PSRD, KB: Assisted in experimentation and editing. BMP, DRH, and SR: Reviewed, edited, and approved the entire work.

Funding
The work was partially supported by the Indian Council of Medical Research-National Institute of Traditional Medicine (ICMR-NITM), Belagavi, Karnataka, India (Intramural funds).