Identification of targetsThe keyword “glioblastoma” was utilized to retrieve data from multiple databases. A comprehensive search was conducted based on the data obtained from Swiss Target Prediction, and 196 genes have been identified to be associated with naringin. DisGeNET identified glioma-related targets (9039 genes). One hundred twenty-nine intersection targets were obtained from the two gene sets for further analysis. An online tool, Venny 2.1.0 (Fig. 1), was utilized to display the intersection of naringin targets and GBM-related differentially expressed genes. Categorical analysis of matching proteins between naringin and GBM indicated the top three categories to be affiliated as protease (21.0%), enzymes (17.0%), and lyase (11.0%) to be the top three classes. Figure 2 represents the pie chart of the predicted target of naringin based on the protein family.Fig. 1Identification of targets related to Naringin and GBM: Venn diagram of naringin targets and GBM targets.Fig. 2Pie chart of the predicted target of naringin based on protein family.Analysis of protein–protein interaction, protein selection for molecular docking and network constructionA database called “STRING” was employed to analyze the relationship between 115 potential therapeutic targets. The Cytoscape version 3.9.0 software was employed to represent the protein–protein interaction network visually. The visualization indicates that the PPI network comprised 570 edges and 129 nodes, with an average node degree 8.84. The nodes are intended to represent the individual proteins, while the edges are meant to depict the relationships that exist between the proteins. Figure 3 represents the PPI networking of common genes from disease and naringin. An increase in the degree value indicates that the protein plays a more substantial function in the network. There were 10 hub genes found when the aforementioned criteria were used to filter the data (Fig. 4): ALB, CASP3, ESR1, PARP1, HSP90AB1, CDK2, CDK1, ABL1, CSP8, and ACE. When it comes to the PPI networks, these genes are considered to be core targets since they possess superior network topological metrics (CC, MF, and BP). A variety of biological regulating activities, including human signal transcription and protein phosphorylation, are regulated by these proteins, which are biological enzymes and cytokines. There were five nodal targets that were shown to be the most essential: ALB, CASP3, ESR1, PARP1, and HSP90AB1. These nodal targets are highly connected with other prospective pharmacological targets and have a vital role in the treatment of GBM. For the purpose of conducting a cluster analysis of protein groups, the MCODE plugin was utilized, and the result was the acquisition of six cluster modules (Fig. 5). There were sixteen significant protein targets that were grouped together in the top functional modules. These targets included MMP7, MMP12, TCF12, ADAM17, MMP3, ALB, MME, FGF2, TYR, CAB39, CXCL10, FOLH1, PON1, ACHE, SLC5A2, and MMP13.It was shown that MMP7 and BACE1 are the core nodes in each functional module, which suggests that these major targets have a considerable influence on the effects that naringin has against GBM.Fig. 3PPI (protein–protein interaction) networking of common genes from disease and naringin.Fig. 4Top10 proteins or genes by degree analysis from common targets.Fig. 5Cluster analysis of the PPI network.GO (gene ontology) enrichment analysisIn order to investigate the molecular mechanisms by which naringin acts against GBM, we conducted a GO enrichment analysis on 101 possible therapeutic targets of naringin that are involved in the treatment of GBM. This analysis was performed using the DAVID database. The evaluation was categorized into three distinct groups: BP (Biological Process), MF (Molecular Function), and CC (Cellular Component). We acquired 137 base pairs (BP), 10 MF, and 3 CC terms. The bar graph (Fig. 6) illustrates the top 5 Gene Ontology (GO) items. It provides the gene enrichment and corresponding p-value. A greater abundance of therapeutic genes within that GO term shows a stronger association between the detected GO term and the treatment of GBM compared to other GO terms. The findings showed that primary enriched biological process categories were the response to naringin, regulation of apoptosis, response to an organic substance, response to chemicals, and response to oxygen-containing substances. Endopeptidase activity, hydrolase activity, peptidase activity, catalytic activity, and activity acting on proteins were the primary enriched MFs. The terms with the highest enrichment in CC analysis were the extracellular exosome, extracellular space, and caspase complex.Fig. 6Gene ontology analysis of the top five biological processes, molecular functions, and cellular component categories associated with naringin and GBM.KEGG pathway analysisThe enrichment studies of the KEGG pathway were used to identify the metabolic pathways that could be possible therapeutic targets for treating GB. Utilizing the DAVID database. Table 1 highlights the signalling pathways associated with naringin and GBM, arranged from least to greatest by P value. The findings demonstrated the enrichment of the major targets, which included numerous signalling pathways associated with cancer, such as the TNF, IL-17, and p53 pathways. Supplementary sheet 1 illustrates the network connectivity between the forecasted pathways of glioblastoma and the target genes obtained from the KEGG pathway. With Cytoscape v3.10.0, we created a visual network to examine the connections between active chemicals, possible therapeutic targets, and important pathways in more detail (Supplementary sheet 2). ABL1, CSP8, ACE, CDK2, CDK1, ABL1, PARP1, HSP90AB1, and CASP3 were among the core GBM targets with which the visual network shown a strong correlation, making it the most crucial drug. Consequently, we have chosen these specific targets to be examined in more detail through molecular docking investigations.Table 1 Identified KEGG pathways modulated by naringin against GBM.Molecular binding analysisMolecular binding is a fundamental technique that is essential in drug discovery processes. It involves predicting the interactions between a protein (receptor) and a small molecule (ligand) to understand binding energetics, molecular interactions, and induced conformational changes. If the configuration of the ligand and the receptor is stable, then the potential for action will increase while the binding energy will decrease simultaneously. An interaction between the ligand and the receptor that occurs spontaneously is indicated by a binding energy that is less than 0 kcal/mol, but a binding energy that is less than − 5 kcal/mol implies a snug binding. A selection of the top 10 proteins (ALB, CASP3, ESR1, PARP1, HSP90AB1, CDK2, CDK1, ABL1, CASP8, ACE) were made for conducting molecular docking studies. These proteins were then compared with the drugs the FDA approved. Evidence demonstrates that naringin had superior binding activity compared to the FDA-approved medicines talazoparib and imatinib. Table 2 depicts the docking score of naringin with its corresponding protein and the binding interaction of FDA-authorised drugs with proteins. Docking interactions of naringin with PARP-1 (4ZZZ) and ABL1 (2ABL) proteins and with their standard drugs (talazoparib and imatinib) were illustrated in Supplementary sheet 3 as well as in Fig. 7. Among the PARP-1 inhibitors approved by the FDA, Naringin demonstrated superior PARP binding activity to Talazoparib. The docking score for Naringin was − 12.90 kcal/mol, indicating that it exhibited excellent binding interactions with PARP-1 (4ZZZ). Amino acids such as ILE872,TRP861, LEU877, TYR889, ILE879, ALA880, ILE895, TYR896, LEU984, PHE897, MET890,ALA898, TYR907, CYS908, and, TYR689 showed hydrophobic interaction. It was shown that other amino acids, such as SER904, TYR907, ARG878, ASP770, ASP776, and SER864, had a role in the inhibition of PARP-1 (4ZZZ), which indicates that they interact with the ligand through hydrogen bonds and form π–π bonds with TYR 907. It was observed that the amino acid complexes of the enzyme PARP1 inhibitors (4ZZZ) exhibited polar interaction with ligands. These complexes included ASN767, SER864, HIS862, ASN868, HIE909, and SER904. Talazoparib, an FDA-approved PARP inhibitor, had a − 4.45 kcal/mol binding score, which is lower than naringin’s affinity for the 4ZZZ protein. Only one amino acid, the ASP766 amino acid of the PARP1 inhibitor protein, was involved in the formation of an H-bond interaction with this molecule. Moreover, hydrophobic interaction (TYR907, ILE895, TYR896, PHE897, VAL762, MET890, TYR889, and LEU877), polar interaction (GLN759, HIS862, SER864, SER904), pi–pi stacking (HIS862, TYR889) was observed. Naringin exhibited superior binding activity compared to FDA-approved Tyrosine kinase protein (2ABL) inhibitors (imatinib) in the instance of Protein ABL1 (2ABL) with a score of binding energy of − 8.4 kcal/mol. In comparison, imatinib showed a binding score of − 5.12 kcal/mol. Naringin having established H-bonding interaction with (GLU142, THR136, ILE135, ASN115, GLU174, ASN83) amino acid residue of 2ABL protein. We observed that lipophilic interactions (LEU141, PRO137, ILE135, TYR112, PRO150, TYR147, LEU88, PHE85, and PRO82), polar interactions (SER140, THR136, ASN133, SER132, ASN115, SER173, SER175, SER176, ASN83) and no π–π stacking involved with amino acid complex of 2ABL protein. The molecular docking interactions of imatinib with 2ABL proteins were illustrated in Supplementary sheet 3. Among the top 10 proteins screened, the top two leads with the highest binding activity than FDA-approved drugs were selected for further studies (Molecular dynamic simulation studies).Table 2 Docking score of naringin and FDA-approved drugs with respective protein.Fig. 72D &3D interaction of naringin (A,E) and talazoparib (B,F) with protein 4ZZZ and naringin (C,G) and imatinib (D,H) with 2ABL.Molecular dynamics simulation analysisComplex of naringin with 4ZZZThe RMSD (root mean square deviation) of the complex exhibited a peak variation of around 11.1 Å. It demonstrated deviations below 0.7 Å, showing the consistent stability of the complex for the whole duration of the experiment. After a stabilization time of 50 ns, the RMSD remained constant and showed variations of less than 0.3 Å. The root mean square fluctuation (RMSF) of the complex atoms exhibited a peak fluctuation of approximately 10.2 Å, which was the same as the residues that showed the highest variation in the naringin–4ZZZ complex (Fig. 8; Movie 1—Supplementary sheet 6). These findings suggest that the binding between naringin and PARP-1 forms a durable complex, and that this interaction occurs specifically at the active site of the protein. The Radius of Gyration (RoG) for the complex exhibited fluctuations ranging from around 20.7 Å to 21.5 Å, with slight deviations confirming the protein’s compactness throughout the entire process. The surface area of the solvent-accessible region (SASA) exhibited variability between approximately 168 nm2 and 186 nm2. Initially, during the first 25 ns of the molecular dynamics (MD) simulation, the SASA was unstable. However, subsequently, there was a progressive increase in the SASA, which might be attributed to the formation and deformation of bonds throughout the MD run. During the MD run, a maximum of 10 hydrogen bonds were observed, while a minimum of 8 bonds remained constant. However, the production and deformation of hydrogen bonds occurred intermittently throughout the run. The energy decomposition analysis per residue revealed that ARG218 and TYR247 had the lowest energy contributions of − 2.13 and − 2.54 kcal/mol, respectively. It was discovered that for the approximately 95 ns of the MD run, the ligand’s (Naringin) total energy contribution was − 34.71 kcal/mol (Fig. 9).Fig. 8Molecular dynamic simulation result of PARP1-Talazoparib and PARP1-Naringin complex where: (a) Root mean square deviation (RMSD) of the backbone (black) and complex for PARP1-Talazoparib (red), and (green) backbone and complex (Blue) for PAPRP1-Naringin, (b) (root mean square fluctuations) RMSF of c-alpha atoms (black) and complex (red) for PARP1-Talazoparib and c-alpha atoms (green) and complex (blue) for 4ZZZ-Naringin, (c) Radius of gyration of complex for PARP1-Talazoparib (Black) and PARP1-Naringin (Red), (d) Number of hydrogen bonds formed for PARP1-Talazoparib (Black) and 4ZZZ-Naringin(Red), (e) SASA of complex for PARP1-Talazoparib (Black) and PARP1-Naringin (Red).Fig. 9The decomposition energy of Naringin (blue) and Talazoparib (Orange) with PARP1 where (a) Total energy decomposition vs time and (b) Total energy decomposition vs residue.Complex of talazoparib with 4ZZZThe root mean square deviation (RMSD) of the complex exhibited a peak variation of around 9.4 Å. However, the RMSD of the complex demonstrated deviations below 0.9 Å, suggesting that the complex remained stable throughout the whole run. After a stabilisation time of 100 ns, the RMSD was steady and showed variations of less than 0.3 Å. The RMSF of the complex atoms exhibited a maximal fluctuation of approximately 12.4 Å (Fig. 8 and Movie 2—Supplementary sheet 6). The Radius of Gyration (RoG) for the complex exhibited variations between approximately 20.6 Å and 21.4 Å, with slight deviations suggesting that the protein maintained a compact structure throughout the experiment. The surface area of the SASA exhibited variability between about 169 nm2 and 190 nm2. Initially, during the first 25 ns of the molecular dynamics (MD) simulation, the SASA was unstable. However, it gradually increased subsequently. This might be attributed to the formation and deformation of bonds occurring during the MD run. Throughout the MD run, the strongest hydrogen bonds out of the three was consistently observed. However, the formation and breaking of hydrogen bonds were also observed during the run. The energy breakdown analysis per residue revealed that TYR236 and LYS243 had the lowest energy contributions of − 1.63 and − 5.51 kcal/mol, respectively. The total energy contribution of the ligand was determined to be − 19.12 kcal/mol for the MD run that lasted approximately 100 ns, as depicted in Fig. 9.Complex of naringin with tyrosine kinase protein (2ABL)Over the course of the 200 ns MD run, the complex’s RMSD fluctuated at a maximum of ~ 5.2 Å, falling between ~ 2.5 Å and 5.2 Å, suggesting that the complex is stable. Following an approximately 120-ns MD run, the complex’s RMSD was shown to stabilise, and it subsequently had RMSD deviations of less than 2.0 Å. A maximum fluctuation of approximately 7.4 °C was seen in the RMSF of the complex atoms (Fig. 10; Movie 3—Supplementary sheet 6). There was variation in the complex’s Radius of Gyration (RoG), which ranged from about 16.5 °C to about 18 °C. For the first twenty ns, the SASA displayed instability, fluctuating between approximately 90 nm2 and 108 nm2. After around 80 ns of simulation, there was a progressive increase in SASA with little fluctuations. A total of six H bonds were seen during the whole run. Initially, for approximately 75 ns, the bonds were not stable, but they became stable after approximately 125 ns of molecular dynamics (MD) simulation. Consistently, at least four hydrogen bonds were seen over the duration of approximately 100 ns. The analysis of energy decomposition, per residue demonstrated that PRO77 and PRO64 had the lowest energy contributions of − 1.04 and − 0.57 kcal/mol, respectively. Ligand’s total energy contribution was − 31.79 kcal/mol, after an MD run of approximately 40 ns, as depicted in Fig. 11.Fig. 10Molecular dynamic simulation results of Tyrosine-protein kinase-imatinib and Tyrosine-protein kinase—Naringin complex were, (a) Root mean square deviation (RMSD of complex (red) and backbone (black) for Tyrosine-protein kinase-imatinib, and complex (blue) and backbone (green) for Tyrosine-protein kinase-Naringin, (b) RMSF of carbon-alpha atoms (black) and complex (Red) for Tyrosine-protein kinase-imatinib and c-alpha atoms (green) and complex (blue) for Tyrosine-protein kinase—Naringin, (c) RoG of complex for Tyrosine-protein kinase-Imatinib (Black) and Tyrosine-protein kinase -Naringin (Red), (d) Number of H bonds formed for Tyrosine-protein kinase-Imatinib (Black) and Tyrosine-protein kinase-Naringin (Red), (e) SASA of complex for Tyrosine-protein kinase-Imatinib (Black) and Tyrosine-protein kinase-Naringin (Red).Fig. 11The energy decomposition of Naringin (blue) and Imatinib (Orange) with tyrosine kinase protein where (a) Total energy decomposition vs time and (b) Total energy decomposition vs residue.Complex of imatinib with tyrosine kinase protein (2ABL)During the MD run of 200 ns, RMSD of the complex exhibited a maximum fluctuation of around 9.4 Å, which fell within the range of approximately 4.2 Å to 9.48 Å. This fluctuation serves as an indication of the complex stability. The RMSD value of the complex reached a stable state after approximately 100 ns of MD simulation. Subsequently, the RMSD exhibited variations of less than 3.0 Å. The RMSF of the complex atoms exhibited a maximal fluctuation of approximately 7.0 Å (see Fig. 10; Watch Movie 4—Supplementary sheet 6). The complex and backbone exhibited a fluctuation in the Radius of Gyration (RoG) between approximately 16.5 Å and 18 Å. After 20 ns of unstable behaviour, the SASA fluctuated between ~ 94 nm2 and ~ 110 nm2, and then began to gradually increase with very slight fluctuations after ~ 125 ns of simulation. During the run time a maximum of seven H bonds were observed; the bonds were unstable for around 75 ns at first and became stable after about 125 ns of MD run. Consistently, three H bonds were observed for a duration of approximately 100 ns. PRO09 and THR124 were found to have the lowest energy contributions of − 2.21 and − 0.65 kcal/mol, respectively, according to the total energy decomposition per residue analysis. The ligand’s overall energy contribution was determined to be − 39.95 kcal/mol during a ~ 70 ns MD simulation. According to the results of the total energy decomposition for each residue, PRO09 and THR124 are more likely to interact with one another. They have the lowest energy contribution, which is − 2.21 and − 0.65 kcal/mol, respectively (Fig. 11).Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) analysis of naringin and talazoparib with PARP1 (4ZZZ)Based on the MMPBSA analysis, which was presented in Supplementary sheet 4, it was found that the naringin-PARP1 complex exhibited the lowest values of VDWAALS, EPB, and GSOLV, with respective values of − 60.37 ± 0.42, 101.20 ± 1.02, and 95.20 ± 1.01 kcal/mol. Similarly, the conventional Talazoparib provided a total binding energy of − 31.80 ± 0.45 kcal/mol. This value is lesser than the naringin–4ZZZ complex, which reveals that the complex is highly stable and has a higher affinity for binding with protein.MM-PBSA analysis of naringin and imatinib with ABL-1 (2ABL)The MMPBSA analysis, which was presented in Supplementary sheet 5, demonstrated that the Naringin-Tyrosine kinase protein complex exhibited the lowest values of VDWAALS, EPB, and GSOLV, with respective values of − 31.18 ± 0.62, 46.80 ± 1.63, and 43 ± 1.60 kcal/mol. In the same vein, the complex exhibited the least amount of electrostatic and non-polar contributions from solute–solvent interactions to the solvation energy and total gas phase molecular mechanics energy. The energy decomposition for the complex was − 3.81 ± 0.06, and − 59.34 ± 1.87 kcal/mol, respectively. FDA-approved drug Imatinib exhibited a total relative binding energy of − 27.05 ± 0.58 kcal/mol, which is lower than that of the naringin-tyrosine kinase protein complex.Principal component analysisThe principal component analysis (PCA) study helps in examining the important coordinated movements that occur during ligand binding. Only the first few eigenvectors are known to define the overall motion of the protein. In this study, the eigenvectors are obtained by diagonalizing the matrix. This study used the top 40 and 200 eigenvectors to estimate coordinated movements. Figure 12a,b show the eigenvalues generated from the diagonalization of the covariance matrix of atomic fluctuations versus the appropriate eigenvector in decreasing order. Out of the 40 eigenvectors, it can be safely inferred that both of the studied complexes had lesser movements and formed a stable association with naringin than imatinib and talazoparib. This could only be possible if the ligand association has altered protein structure and dynamics. Using PCA to construct 2D projection plots is another approach for obtaining the dynamics of protein–ligand complexes. Figure 12b,d depict a two-dimensional projection of the trajectories in phase space for the first two main components, eigenvector1 and eigenvector2, for both complexes. In general, the complexes occupying a compact phase space along with a stable cluster are considered more stable, whereas those complexes taking greater space with unstable clusters are considered less stable. As they took up minimal area in the phase space, all of the complexes were found to be very stable. The above PCA and other MD simulation findings were also in accordance with the 2D PCA reports.Fig. 12Principal component analysis. The plot of the eigenvalues vs first 40 eigenvectors (a,b), the first two eigenvectors showing the protein motion in phase space for naringin (black) and Imatinib (red), (c,d) naringin (black) and Talazoparib (red).In vitro anticancer effect of NaringinThe cytotoxic effects of Naringin on C6 rat glioma cells and U-87MG human glioblastoma cells were assessed using the MTT assay. Figure 13 demonstrated a dose-dependent reduction in cell viability for both cell lines upon treatment with naringin. The EC50 value, representing the concentration of naringin required to inhibit cell viability by 50%, was determined to be 34.1 µg/ml. This suggests that C6 cells exhibit moderate sensitivity to Naringin treatment. A more pronounced cytotoxic effect was observed in U-87MG cells, with an EC50 value of 10 µg/ml. This indicates that U-87MG cells are significantly more sensitive to Naringin compared to C6 cells.Fig. 13In vitro cytotoxicity assay of naringin against (A) C6 cells and (B) U87MG cells.
