Table of Contents
Molecular Docking results
Molecular docking was conducted using the PyRx virtual screening tool integrated with AutoDock Vina, a highly efficient molecular docking and virtual screening software. The 3D structures of the target proteins were obtained from the RCSB Protein Data Bank, and the ligands (natural compounds and reference drugs) were sourced from the PubChem and ZINC databases. Before docking, the protein and ligand structures were prepared by adding polar hydrogens, assigning appropriate charges, and removing previously attached ligands, solvents, and water molecules. Grid boxes were defined around the active sites of the target proteins to ensure the correct positioning of the ligands during docking. The docking process involved simulating the binding of each ligand to the target protein’s active sites and calculating the binding affinities, expressed in kcal/mol. These binding affinities were used to predict the interaction strength between the ligands and the target proteins. In molecular docking, more negative binding affinity values (ΔG, in kcal/mol) indicate stronger predicted ligand-target interactions, with values below − 7.0 kcal/mol generally considered strong binders under in-silico conditions. The results were analyzed to identify natural compounds with higher or comparable binding affinities to the reference drugs, indicating their potential as bioactive drug candidates for treating T2DM. Supplementary File 6 (Table S6) encompasses the complete molecular docking record for all genes.
01. GPD2 Docking results
GPD2 was selected for due to its central role in hepatic gluconeogenesis, a process often dysregulated in T2DM. Moracin D and Moracin P bound to active site A with binding affinities of -9.8 and − 9.7 kcal/mol, respectively. These values are higher than the binding affinity of NADH (-9.0 kcal/mol), indicating that Moracin D and Moracin P may serve as better bioactive drug candidates than the drug, NADH. No compounds bound to the active site C that Metformin binds to; therefore, no comparisons could be made for this pocket.
02. IRS1 Docking results
IRS1 was chosen due to its critical involvement in insulin signal transduction. It is a known mediator of insulin action and resistance, making it a particularly relevant target in the context of T2DM. 5-[5-(4-Hydroxy-Benzyl)-4-(4-Methoxy-Benzyl)-1-Methyl-1 H-Imidazol-2-Ylamino]-3-Methyl-Imidazole-2,4-Dione showed a binding affinity of -7.1 kcal/mol. This is competitive with the drug [4-({5-(AMINOCARBONYL)-4-[(3-METHYLPHENYL)AMINO] PYRIMIDIN-2-YL}AMINO)PHENYL]ACETIC ACID (-7.2 kcal/mol), although they bind to different active sites, preventing a direct comparison. Knowing that the compound 5-[1-(4-Hydroxy-Benzyl)-4-(4-Methoxy-Benzyl)-1 H-Imidazol-2-Ylamino]-3-Methyl-Imidazole-2,4-Dione also binds to the other binding pocket than the reference drug, therefore, no comparison can also be made for this compound.
03. PPARG Docking results
Troglitazone, Icosapent, and AMG-131 bound to active site A with binding affinities of -7.9, -6.5, and − 9.1 kcal/mol, respectively. Pyrene, Guggulsterone, and Emodin exhibited stronger binding affinities of -8.6, -8.4, and − 8.1 kcal/mol, respectively, compared to Troglitazone. Ten compounds showed stronger binding affinities than Icosapent. Pyrene emerged as the most potential compound for pocket A; no compounds could present higher values than the reference drug AMG-131.
Pioglitazone and Rosiglitazone bound with affinities of -8.5 and − 6.9 kcal/mol, respectively.
Compounds such as Plantagineoside A (-9.3 kcal/mol), Isosilybin A (-9.1 kcal/mol), Amorphastilbol (-8.9 kcal/mol), Piperitol, Chrysin, and Plantagineoside B (-8.7 kcal/mol), (1R,2 S,5R,6R)-5’-O-Methylpluviatilol, Silibinin, and Bavachinin A (-8.6 kcal/mol) showed stronger binding affinities compared to Pioglitazone. In total, 35 compounds outperformed Rosiglitazone. Plantagineoside A was identified as the most potential compound for pocket B.
04. IAPP Docking results
Due to the 3D structure of Copper only containing an ion (Cu), it was unsuitable for docking. Thus, no comparisons could be made between the only compound found for IAPP (Curcumin (1E,4Z,6E)-5-hydroxy-1,7-bis(4-hydroxy-3-methoxyphenyl)hepta-1,4,6-trien-3-one) with the reference drugs Copper. Despite being unable to compare the compound Curcumin (1E,4Z,6E)-5-hydroxy-1,7-bis(4-hydroxy-3-methoxyphenyl)hepta-1,4,6-trien-3-one with Copper, this compound showed highly promising values of binding affinity (-35.6 kcal/mol) with the protein amylin, making it worth investigating further.
05. GCK Docking results
The drug 3-[(4-fluorophenyl)sulfanyl]-N-(4-methyl-1,3-thiazol-2-yl)-6-[(4-methyl-4 H-1,2,4-triazol-3-yl)sulfanyl]pyridine-2-carboxamide was found to be the most effective with a binding affinity of -8.9 kcal/mol. The compound, Rohitukine showed a higher binding affinity (-7.8 kcal/mol) compared to Beta-D-Glucose (-5.2 kcal/mol), indicating its potential as a more effective drug candidate.
06. ABCC8 Docking results
Glyburide (Glibenclamide) showed a binding affinity equal to Gliquidone (-10.7 kcal/mol), suggesting equal potential. Berberrubine (-8.5 kcal/mol) had a higher binding affinity than Nateglinide (-8.3 kcal/mol), making it a potential drug candidate. Also, Berberine Chloride (-8.2 kcal/mol) exhibited a higher binding affinity than Tolbutamide (-7.5 kcal/mol), indicating its potential efficacy.
07. MAPK8 Docking results
Both compounds Apigenin (-7.1 kcal/mol) and Emodin (-7.0 kcal/mol) displayed higher binding affinities than the reference drug Halicin (-5.8 kcal/mol), therefore, they could be more effective bioactive drug candidates.
08. MTNR1B Docking results
5-Meo-Dmt and Serotonin bound to active sites other than pocket A, thus, no direct comparisons could be made. The drug, Tasimelteon had the strongest binding affinity (-8.2 kcal/mol) among all studied ligands. Compounds N-[2-(5-methoxy-1 H-indol-3-yl)ethyl]cyclopropanecarboxamide (-7.7 kcal/mol), Propionamide (-7.6 kcal/mol), Z02 (-7.4 kcal/mol), and Melatonin (-7.3 kcal/mol) showed high binding affinities after Tasimelteon, respectively.
09. AKT2 Docking results
Chelerythrine (-10.4 kcal/mol) and Alvocidib (-10.2 kcal/mol) showed higher binding affinities than both drugs Capivasertib (-9.7 kcal/mol) and N-[(1 S)-2-amino-1-phenylethyl]-5-(1 H-pyrrolo[2,3-b]pyridin-4-yl)thiophene-2-carboxamide (-9.1 kcal/mol), indicating their strong potential as natural drug candidates. The compound, A-443,654 (-9.3 kcal/mol) exhibited a higher binding affinity than drug N-[(1 S)-2-amino-1-phenylethyl]-5-(1 H-pyrrolo[2,3-b]pyridin-4-yl)thiophene-2-carboxamide, suggesting its potential effectiveness.
10. PTPN1 Docking results
Compounds bound to 11 active sites (A-K). Both studied drugs Trodusquemine (-5.8 kcal/mol) and Ertiprotafib (-7.4 kcal/mol) bound to site A. Eighty-eight compounds exhibited binding affinities of -7.4 or lower, making them competitive with Ertiprotafib. Among these, nine compounds matched the − 7.4 affinity, while seventy-nine compounds had affinities ranging from − 7.5 to -9.0, indicating a stronger potential than Ertiprotafib. Only six compounds (Z02, 325, 263, 216, 394, and 232) showed binding affinities equal to or lower than Trodusquemine. Z02 – Caffeic acid matched the − 5.8 affinity, while the others exceeded this value, making them less comparable. The strongest binding affinities to pocket A were observed with Xambioona (-9.0 kcal/mol), 18Alpha-Glycyrrhetinic Acid (-8.8 kcal/mol), 18Beta-Glycyrrhetic Acid (-8.8 kcal/mol), Taiwaniaflavone (-8.8 kcal/mol), and Oleanderolide (-8.7 kcal/mol), highlighting these as the most potential drug candidates for PTPN1.
11. INSR Docking results
Linsitinib and the compound Altertoxin I exhibited the strongest binding affinity value of -9.3 kcal/mol, indicating that Altertoxin I has equal potential to Linsitinib. Furthermore, the compounds Altertoxin I, Alvocidib, A-443,654, Chelerythrine, Gefitinib, Ellagic Acid, Macrosporin, Alternariol monomethyl ether, Alternariol, Sb-202,190, 6-O-Methylalaternin, and Sp-600,125 demonstrated higher binding affinities than Fostamatinib, suggesting their superior potential compared to this drug.
12. AMPK Docking results
Compounds Chelerythrine, Sp-600,125, Rohitukine, and Sb-202,190 displayed binding affinity values of -10.2, -8.8, -8.4, and − 8.0 kcal/mol, respectively. These values are higher than those for drugs Fostamatinib (-7.8 kcal/mol) and Aspirin (-5.9 kcal/mol), indicating these compounds have higher binding potential.
13. GAA Docking results
The four ligands for the protein GAA bound to four different binding sites. Consequently, no direct comparison of their binding abilities can be made.
14. SLC2A4 Docking results
Nine compounds (C404, C168, C461, C462, C463, C296, C292, C338, and C193) exhibited better binding affinities than the drugs Ascorbic Acid and Glucosamine. Among these, Ursolic Acid and Oleanolic Acid showed the best binding affinities, highlighting their significant potential as alternatives to Ascorbic Acid and Glucosamine.
Table 2 summarizes the top compounds exhibiting potential antidiabetic properties, with superior docking values compared to existing reference drugs. These compounds show significant promise for further evaluation and potential use as alternative medications or natural therapeutics to currently approved or investigational drugs.
The list is later narrowed down from 72 potential ligands to 18 ligands (Table 3) that showed the most remarkable docking scores for each of the studied proteins. The selected 18 ligands were used as a smaller subset of the data for protein-ligand binding site interaction and molecular docking analyses.
Analysis of Docking Ligand-Protein interactions at the binding site
The molecular docking studies revealed specific interactions between the selected subset of 18 ligand molecules and targeted amino acid residues within the binding site of their responsive protein, listed in Table 3.
Supplementary File 7 (Figures S1-S24 and Tables S7-S58) illustrates the docking pose of each protein-ligand interaction complex number to highlight the key interactions with involved residues and show types of interactions and bonds at the binding site, which contribute to the stability and specificity of the ligand binding.
PK/PD and ADME/T analysis results
The PK/PD and ADME/T properties of the 72 selected natural compounds were evaluated and predicted and top ligands were selected for their superior ADME/T characteristics and potential to effectively traverse biological barriers which is essential for oral administration of antidiabetic drugs. The ADME/T results are provided in Supplementary File 8 (Tables S59-S62).
1.
Selection of Antidiabetic Compounds Based on ADME Profiles.
In the selection of potential antidiabetic compounds, a meticulous evaluation was conducted focusing on their Absorption, Distribution, Metabolism, and Excretion (ADME) profiles, with an emphasis on Human Intestinal Absorption (HIA), Caco-2 Permeability, and Blood-Brain Barrier (BBB) permeability. The rationale behind choosing each compound was grounded in their capability to ensure optimal therapeutic efficacy and patient compliance. Below is a detailed justification for the choice of each compound based on these ADME criteria:
Pyrene (C54) and Guggulsterone (C60) exhibit maximum HIA and significant Caco-2 permeability, suggesting excellent oral bioavailability. Their high BBB permeability also suggests potential utility in addressing diabetes-related neurological complications. Despite Pyrene showing favorable ADMET properties, it is flagged as both carcinogenic and mutagenic (Table S62) and was therefore excluded from the lead compound consideration. Melatonin (C103) and Gefitinib (C456) show perfect and near-perfect HIA scores respectively, with considerable BBB permeability, indicating potential benefits beyond glycemic control, particularly in neurological protection. Apigenin (C100), Rotenone (C73), and Curcumin (C96) are characterized by high HIA and Caco-2 scores, ensuring effective absorption and systemic availability, crucial for consistent therapeutic effects. Bavachinin A (C24), Bavachinin (C10), and Quinidine (C61), while displaying lower BBB permeability, maintain high HIA and Caco-2 permeability, focusing their action on peripheral organs involved in glucose metabolism. This selection presumes these compounds are safe, pending further pharmacological and toxicological evaluations to fully ascertain their therapeutic profiles.
2.
Evaluation of Toxicity Profiles for Antidiabetic Compounds.
In parallel with ADME profiling, an assessment of the toxicity profiles was crucial in selecting antidiabetic compounds to ensure patient safety, particularly relevant in chronic conditions like diabetes where long-term medication is common. Here’s the rationale for choosing compounds based on their toxicity characteristics:
2-Tert-Butyl-6-[(3-Tert-Butyl-2-Hydroxy-5-Methylphenyl)Methyl]-4-Methylphenol (C95) and 4-(2-Phenylpropan-2-Yl)Phenol (C45) both fall within Toxicity Class V, indicating they may be harmful if swallowed but possess a relatively low acute toxicity, which is an important consideration for medications intended for chronic use. Their inactive statuses across carcinogenicity, mutagenicity, and nutritional toxicity, coupled with high confidence values, suggest these compounds are associated with minimal long-term adverse effects.
Phenothiazine (C80), although minimally toxic (Class V), shows an excellent safety profile with inactive statuses in all three toxicity indicators. This compound’s high LD50 value and safety markers make it a promising candidate for diabetes treatment, especially considering the potential for long-term therapy. 2-Naphthalen-1-Ylacetic Acid (C26), categorized under Toxicity Class IV, offers a moderate safety margin. Despite its relatively higher acute toxicity, the inactive status in major toxicity categories underscores its potential as a safer alternative for diabetes management, provided that its efficacy and specific diabetic application are validated in clinical settings. The focus on compounds with high safety profiles is aimed at minimizing potential risks associated with long-term drug administration, thereby enhancing patient compliance and improving overall quality of life for diabetic patients. These compounds represent not only effective therapeutic options but also safer alternatives in the development of antidiabetic medications. Further experimental and clinical validation of these compounds could elucidate their efficacy and safety as potential treatments for diabetes.
Network Pharmacology analysis results
The compound-target-pathway network analysis revealed significant interactions between the selected 72 natural compounds and the 14 target proteins associated with T2DM (Fig. 1).
The systematic compound-target-pathway network in this study.
Key findings indicated that several compounds, such as Curcumin, Chelerythrine, and Ursolic Acid, demonstrated strong interactions with critical proteins involved in insulin signaling, such as PPARG, SLC2A4, and PTPN1. These proteins were identified as central hubs in the network due to their high degree of connectivity with multiple compounds.
The identified compounds indicated potential synergistic influences on effective T2DM management. The network visual representation underscored the importance of multi-target interactions of five compounds: Chelerythrine (C105) linked to three target genes AKT2, INSR, and AMPK, and compounds Emodin (C57), Rohitukine (C97), A-443,654 (C107), and Alvocidib (C109), each targeted two genes, simultaneously affecting multiple proteins within the same pathway.
Gene set enrichment analysis
Pathway enrichment analysis further highlighted the involvement of key signaling pathways, including the type II diabetes mellitus pathway (MAPK8, IRS1, ABCC8, INSR, SLC2A4, and GCK), insulin signaling pathway (PTPN1, MAPK8, IRS1, INSR, AKT2, SLC2A4, and GCK), insulin resistance (PTPN1, MAPK8, IRS1, INSR, AKT2, and SLC2A4), adipocytokine signaling pathway (MAPK8, IRS1, AKT2, and SLC2A4), and AMPK signaling pathway (IRS1, INSR, AKT2, PPARG, and SLC2A4). While the pathways with highest combined scores were related to diabetes and obesity, some of the genes also incorporated within other pathways associated with other diseases (Fig. 2). For instance, based on overlap genes, Non-alcoholic fatty liver disease (MAPK8, IRS1, INSR, AKT2, and PPARG), pancreatic and colorectal cancer (MAPK8 and AKT2), Alzheimer disease (MAPK8, IRS1, INSR, and AKT2), and thyroid cancer (PPARG), endometrial and prostate cancer and melanoma (AKT2), tuberculosis and human immunodeficiency virus 1 (HIV) infection (MAPK8 and AKT2) were of other remarkable integrated pathways with these compounds and genes. Supplementary File 9 (Figures S25-S27) offers a detailed Gene Ontology (GO) functional pathway enrichment analysis for selected genes in 2023. In the GO Biological Process category, “Response to Insulin” is the most significantly enriched term, with a p-value of 5.25e-12. The top term for the GO Cellular Component is “Protein Kinase Complex,” having a p-value of 5.27e-05. Lastly, in the GO Molecular Function category, “Insulin-Like Growth Factor Receptor Binding” is highlighted as the most enriched function, with a p-value of 3.03e-05. GO results underscore critical biological functions, components, and interactions associated with the studied gene set to provide valuable insights into their roles in cellular processes in T2DM.
Top enriched terms from the KEGG 2021 Human gene set library37,38,39. The top 10 enriched terms for the input gene set are displayed based on the -log10 (p-value), with the actual p-value shown next to each term. The term at the top has the most significant overlap with the input query gene set.
Molecular dynamics simulation of top-ranked protein-ligand complexes
To assess the dynamic behavior and structural stability of the 17 most promising protein-ligand complexes (Complex 05 excluded) identified via molecular docking (complexes are listed in Table 3), molecular dynamics (MD) simulations were performed for 500 ns under NPT ensemble conditions (310 K and 1 atm). Parameters including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), hydrogen bond formation, and thermodynamic energy profiles were evaluated (Table 4).
1.
Structural Stability and Flexibility (RMSD, Rg, RMSF).
RMSD values ranged from 4.39 Å (Complex 01: Moracin P–GPD2) to 5.33 Å (Complex 18: Oleanolic Acid–SLC2A4), suggesting that all complexes reached equilibrium and maintained overall structural integrity. The average radius of gyration (Rg) across all complexes remained between 30.41 Å and 31.01 Å, confirming that no major conformational collapse occurred. RMSF analysis revealed that most residues remained relatively rigid, with the lowest fluctuation in Complex 06, GCK–Rohitukine (0.67 Å) and highest in Complex 12, AKT2–Alvocidib (2.51 Å), the latter possibly indicating increased flexibility near loop or surface regions.
2.
Hydrogen Bond Dynamics.
Hydrogen bonding, critical for structural stability, varied across the studied complexes. The Complex 02 (Moracin D–GPD2) exhibited the highest average number of hydrogen bonds (439.9), while Complex 04 (Pyrene–PPARG) formed the fewest (249.98). Complexes involving Moracin compounds (Complexes 01 and 02), Berberrubine (Complex 08), and Ursolic Acid (Complex 17) consistently demonstrated high H-bond occupancy which reflects strong ligand-target interaction stability throughout the simulation.
3.
Solvent Accessibility and Buried Surface Area.
SASA analysis provided results on solvation effects and structural compactness. GPD2-bound complexes (Complexes 01 and 02) showed the largest solvent-accessible surface areas (17,273 Ų and 17,643 Ų), consistent with their extensive external contact regions. On the contrary, MAPK8 and PTPN1 complexes (Complexes 09 and 14) exhibited lower SASA values (10,781 Ų and 11,201 Ų), indicating a tighter solvation shell. Buried SASA, which reflects buried interface regions, was highest for Complex 02 (1,076 Ų) which explains its strong complexation potential. These observations were further supported by structural visualizations (Supplementary File 10, Figures S28-S44) illustrating noticeable conformational shifts and compact rearrangements across simulation endpoints.
4.
Energetic Profiles.
Potential energy values were most favorable for Moracin P–GPD2 (Complex 01: − 1.46 kcal/mol) and Moracin D–GPD2 (Complex 02: − 1.15 kcal/mol). Enthalpy and free energy values additionally confirmed these two complexes as energetically stable. Ursolic Acid–SLC2A4 (Complex 17) and Plantagineoside A–PPARG (Complex 03) also ranked among the top performers with free energy values of − 33.94 and − 36.89 kcal/mol, respectively.
5.
Structural and Surface Visualizations.
Visual inspection of initial and final simulation frames highlighted overall structural resilience across all 17 selected complexes, which showcases minimal unfolding and favorable compactness post-equilibration. The relevant visual data for all complexes (Supplementary File 11, Figures S45-S61) support quantitative findings on Rg, SASA, and RMSD stability, visually confirming the adaptive binding nature of the ligands without any major disruptions to the protein scaffold.
