Integrated in vitro, microarray, and network pharmacology analysis reveals the multi-target anti-diabetic potential of Vigna unguiculata

Abstract

Introduction: Diabetes mellitus (DM), particularly type 2 DM (T2DM), is a chronic metabolic disorder requiring
novel therapeutic approaches as the available therapies are not meeting the current challenges. This study investigates
the anti-diabetic potential of Vigna unguiculata using a network pharmacology approach, supported by
in vitro and in silico analyses.
Methods: The plant was collected from Khyber Pakhtunkhwa, Pakistan, and subjected to hydroalcoholic extraction
and fractionation. In vitro assays included α-amylase, α-glucosidase, and aldose reductase. Target prediction
using STITCH and SwissTargetPrediction identified 88 common genes linked to T2DM. Protein-protein interaction
(PPI) network analysis highlighted key genes like EGFR, PTGS2, and TLR4 as central nodes in diabetes-related
pathways. Molecular docking was used to study the binding affinities of compounds.
Results: IC50 values were determined using IBM SPSS Statistics 21 software. The data underwent analysis using
one-way ANOVA followed by Dunnett’s multiple comparison test. Significance value was determined at *p < 0.05,
**p < 0.01 and ***p < 0.001. In-vitro assays demonstrated significant α-amylase, α-glucosidase, and aldose reductase
inhibitory activities. Phytochemical screening identified several bioactive compounds. Functional annotation
and KEGG pathway analysis confirmed these genes’ roles in crucial metabolic pathways. Virtual screening
revealed strong binding affinities of compounds like Stigmasterol, Luteoline, and Quercetin with GSK3B, PTGS2,
and TLR4. The Molecular Dynamics (MD) simulation, binding free energy calculations (MM-PBSA and MM-GBSA),
confirmed the results of Virtual screening.
Conclusion: In short, these findings underscore V. unguiculata as a promising source for anti-diabetic agents, supporting
further clinical trials for T2DM management.

Introduction

Type 2 diabetes mellitus (T2DM) is a growing and important health issue worldwide, characterized by continuous high blood sugar levels caused by resistance to insulin and eventual malfunction of pancreatic beta-cells (1). The impact of this prevalent form of diabetes is significant worldwide, accounting for around 90% of cases. This impact is felt deeply in terms of both individual health and healthcare systems (2,3). T2DM’s worldwide occurrence is on the rise due to factors like longer life expectancy, lack of physical activity, and shifts in eating habits, which makes it a significant issue in the epidemiology of non-communicable diseases (4). The mortality rates linked to T2DM are concerning, since it leads to cardiovascular diseases, kidney failure, lower limb amputation, and diabetic retinopathy, as well as other complications (5). The mortality rate linked to T2DM is measured in scientific research using hazard ratios or relative risks when compared to non-diabetic groups. For example, people with T2DM have double the likelihood of dying at any age compared to those without the condition.

Treatment for T2DM mainly focuses on controlling blood sugar levels using lifestyle changes, oral anti-diabetic medications, and insulin treatment (6). The methods mentioned typically do not fully tackle the multifaceted characteristics of the illness. Many traditional treatments focus on only one element of the disease’s intricate pathophysiology, potentially resulting in less than ideal results and failing to stop the advancement of the disease (7). With these restrictions in mind, there is an urgent requirement for a transition towards multi-target drug discovery. This method seeks to create medications that can concurrently impact numerous pathways associated with T2DM, providing a more comprehensive treatment option. These approaches are expected to enhance treatment effectiveness and patient outcomes by tackling the various mechanisms involved in the development and advancement of the disease. Therefore, T2DM continues to be a significant health challenge with substantial morbidity and mortality. The insufficiency of current single-target treatments emphasizes the need for innovative treatment approaches, especially those involving multi-target interventions, to more effectively handle and potentially reverse this intricate and widespread disease.

Conventional treatments primarily focus on managing blood glucose levels and mitigating symptoms, yet they often fail to address the underlying causes and multiple pathways involved in the disease’s progression. Recent advancements in network pharmacology have underscored the importance of targeting multiple biological pathways simultaneously to achieve more effective therapeutic outcomes. Network pharmacology is a cutting-edge approach that integrates systems biology, multi-target drug actions, and bioinformatics to systematically explore the therapeutic potential of drugs, including those derived from medicinal plants (8,9). This approach is especially relevant in the context of T2DM, where the interplay between insulin signaling, inflammation, and oxidative stress contributes to the disease’s complexity. In this study, we explore the potential of V. unguiculata, a plant with a rich history in traditional medicine, as a source of multi-target anti-diabetic agents. V. unguiculata (commonly known as cowpea or lobia) is a leguminous plant belonging to the Fabaceae family (10). Traditionally used in various medicinal practices, it has been suggested to possess anti-diabetic properties, potentially due to its bioactive compounds that may influence glucose metabolism and insulin sensitivity (11,12). However, the exact mechanism is currently unknown.

Following that, this study provides compelling evidence for the anti-diabetic potential of V. unguiculata through an integrated approach that combines traditional knowledge with modern scientific techniques. The study highlights the importance of a multi-target therapeutic strategy, particularly in managing complex diseases like T2DM, where multiple biological processes are dysregulated. The V. unguiculata has been traditionally employed for the ailment of various diseases, including DM. However, the mechanism by which V. unguiculata exhibits anti-diabetic activity has not been explored. Based on the previously reported activities, the bioactive extract of the V. unguiculata has been evaluated against DM using in vitro studies, network pharmacology, molecular docking, MD simulation, and binding free energies with the intention to explore the potential anti-diabetic mechanism. In this study, the key approach was adopted to explore the pathways, which plays crucial role in the pathogenesis of the DM, i.e., by inducing damage to the insulin producing cells and confer insulin resistance such as PTGS2, Glycogen Synthases Kinase 3 Beta (GSK3B) and Toll-like receptor-4 (TLR-4) signaling (13-15). The GSK3B has been focused recently on diabetes management due to its involvement in insulin resistance, insulin synthesis, and protection of the B-cells (Islet of Langerhans cells). Similarly, the role of TLR-4 has also been implicated in the pathogenesis of various diseases, including diabetes; enhanced TLR-4 activity paves the way for inflammation of various tissues, including the pancreas, and cross-talk with multiple signaling pathways such as MAPKs, NF-kB, and JAK/STAT signaling (16,17). Primarily, in this study, V. unguiculata plant bioactive compounds were employed against the multiple signaling pathways involved in Diabetes pathogenesis. By employing comprehensive computational approaches, it was assumed that bioactive compounds from the V. unguiculata interact with the multiple biological signaling pathways concerned with glucose metabolism, insulin resistance, and the inflammatory process, which will eventually improve diabetes symptoms by regulating the blood glucose level (13-15).

Currently, several classes of anti-diabetic drugs are in clinical practice for the management of DM. However, the existing drugs are associated with various unwanted side effects such as weight gain, lactic acidosis, urinary tract infections, and the cases of DM are increasing globally. Thus, there is a need for the exploration of new drugs, which are not only effective but also associated with fewer side effects. Natural products offer a cheap source for new drug development, and they can be used to discover drugs that fight various diseases, including DM (18,19). The use of V. unguiculata plant as an anti-diabetic has been reported traditionally; thus, based on the previous reported studies, it was anticipated that V. unguiculata would be effective against DM.

Methodology

Plant material collection, identification, and extraction

In August 2022, a fresh whole plant (14kg) of V. unguiculata was gathered from Village Rustam, district Mardan, Khyber Pakhtunkhwa, Pakistan. Dr. Mohib Shah from the Department of Botany, AWKUM, authenticated and identified the plant specimen by comparing it with existing specimens at the Department of Botany, AWKUM, and the flora of Pakistan. The Department of Botany, AWKUM, assigned Voucher No. AWKUM.Bot.100.31.3.9 to the plant specimen. After cleaning, fresh plant material was dried in the shade at room temperature. The entire dried plant was crushed, resulting in 838 gm of powdered plant. The powdered plant was macerated using a hydroalcoholic solvent (90 % methanol and 10% distilled water) following the procedure described by Mistriyani et al. (20). The procedure was repeated three times to achieve the highest extraction possible. The solvent was evaporated using a rotary evaporator to concentrate the filtrate. The concentrated extract was then dried in an oven at approximately 40 °C. Ultimately, 186g of oily, viscous extract was acquired and stored in a refrigerator at 2-8 °C.

Fractionation

Distilled water at a warm temperature was added to the hydro-methanolic extract. The mixture was then separated using n-hexane in a separating funnel. The remaining methanol extract was then fractionated using ethyl acetate to yield n-hexane, ethyl acetate, and water fractions. This process was repeated multiple times to obtain the maximum fractions. The fractions were further concentrated using a vacuum rotary evaporator and then dried in an oven at approximately 40 °C.

In-vitro anti-diabetic activity

The in vitro anti-diabetic activities of both the crude hydro-methanolic extract and the subsequent fractions were examined.

α-amylase inhibitory activity

The already established technique was used to conduct the alpha amylase inhibitory assay (21). In a 96-well microplate, 5 μL of the enzyme solution and 15 μl of sample, diluted in phosphate buffer at various concentrations, were added. The reaction was initiated with the addition of 20 μl starch solution and then incubated for 10 minutes at 37 °C, followed by further incubation at 37 °C for 30 minutes. HCl (10 μl of 1M) and iodine reagent (75 μL) were added to each well to stop the reaction. Acarbose at a concentration of 64 μg/mL served as the positive control, and a blank was created using phosphate buffer (pH 6.9). The absorbance at 580 nm was measured, and the percentage inhibitory activity was calculated using the provided equation:

% Inhibition = Abs of Control – Abs of sample /Abs of control x 100

α-glucosidase inhibitory activity

The already established method for measuring alpha-glucosidase inhibitory activity was used (22). Test samples at different concentrations (20 μl), alpha-glucosidase (1 U/ml) (10 μl), and phosphate buffer (100 mM, pH = 6.8) (50 μl) were combined and allowed to preincubate for 15 minutes at 37°C in a 96-well microplate. Following the addition of P-nitrophenyl glucopyranoside (5 mM) as the substrate (20 μl), the mixture was incubated for 20 minutes at 37°C. The reaction was stopped by adding 50 μl of 0.1 M sodium carbonate (Na₂CO₃). The absorbance of the liberated p-nitrophenol was determined at 405 nm using a microplate reader. Acarbose was utilized as a standard in different concentrations. Each experiment was conducted in triplicate, with a parallel setup that did not contain any test sample. The % inhibitory activity was calculated using the provided equation:

% Inhibition = Absorbance of Control – Absorbance of sample /Absorbance of control x 100

Aldose reductase inhibitory assay

The evaluation of aldose reductase inhibitory activity was conducted using a previously established method (23) with some modifications. Upon euthanasia, the eyes of normal Swiss albino mice were promptly extracted, and the lenses were subsequently cleansed with normal saline. The weights of the lenses were then re-determined. Once transparent lenses were obtained, a 10 % homogeneous solution was prepared in 0.1 M phosphate buffer saline (pH 7.4). This solution was then subjected to centrifugation at 5000 × g for 10 minutes. The resulting supernatant was collected and stored on ice. The protein content of the lens homogenate was determined. A sample cuvette containing 0.7 mL of phosphate buffer, 0.1 mL of NADPH (25 × 10−5 M), 0.1 mL of lens supernatant, and 0.1 mL of DL-glyceraldehyde (5 × 10−4 M) was utilized to assess the aldose reductase inhibitory activity of the fractions, with a final volume of 1 mL. Absorbance measurements were taken in comparison to a reference cuvette devoid of DL-glyceraldehyde. The pH of the reaction mixture was adjusted to 6.2, and the enzymatic reaction was initiated upon mixing the substrate with the solution mixture. The absorbance was then monitored at 340 nm for three minutes at 30-second intervals. The extract was dissolved in phosphate buffer saline (PBS) to produce a stock solution containing all the fractions. Cuvettes treated with 0.1 mL of each fraction from various stock solutions, with final concentrations of 10, 25, 50, 100, 200, and 300 µg/mL, were used as reference and standard cuvettes to evaluate aldose reductase inhibitory activity. The process was initiated by adding 0.1 mL of DL-glyceraldehyde. The percentage inhibitions of the aldose reductase assay for the extract were calculated assuming that the lens of a normal mouse exhibited 100 % activity. IC₅₀ values were calculated for each sample.

Phytocompounds Screening

Active compounds of V. unguiculata were identified from existing literature and several databases of biologically active phytochemicals, including Indian Medicinal Plants, Phytochemistry, and Therapeutics (IMPPAT) (24), KNApSAcK (25). The keyword’ Vigna unguiculata’ was used for database searches, while a comprehensive literature review was conducted using PubMed and Google Scholar. The bioavailability (F) and drug-likeness (DL) of all ingredients were assessed to determine their ADME (absorption, distribution, metabolism, and excretion) properties. Compounds were selected based on a DL value of at least 0.18 and an OB value of at least 30%. The OB and DL values for all active compounds were calculated using SwissADME and ADMETlab. Additionally, chemical information such as CID number, structure, and molecular weight was collected from PubChem and ChemSpider.

Retrieval of potential targets of selected compounds

Swiss Target Prediction (26) and the STITCH database (27) were used to predict the potential targets of selected compounds. These compounds were uploaded to the STITCH database with the options set to ‘Homo sapiens’ and using screening criteria as combined scores 0.7 or higher. Further, the SMILES notation of each compound was inputted into the Swiss Target Prediction web application to predict targets by employing the reverse pharmacophore mapping approach and predicted scores with a probability score of 0.7. In network pharmacology, this prediction is very commonly employed for the generation of possible compound-target interactions using available data. The validation of the predicted targets and subsequently comparing them with the known DM target is crucial. Hence, the reliability of the predicted targets was enhanced using publicly available curated databases, i.e., GeneCards and DisGeNET databases and the predicted targets of the STICTH database were compared.

Microarray data analysis

In order to validate the findings of the study, three microarray datasets, including GSE22435, GSE43950, and GSE92724, were obtained from the NCBI GEO database (Online) (28). GEO is a high-throughput gene expression data repository, including microarray and hybridization arrays. In this study, DEGs were analyzed using the Limma R package. Limma is widely known for its efficient tools for reading, normalizing, and analyzing gene expression data for DEGs, and hence, it is the first preferred choice for high-throughput and microarray data. Screening criteria to identify candidate genes were adjusted p-value less than 0.05 and absolute value of log (Fold Change) greater than 1.0 but less than −1.0. In the R program, we used the ggplot2 package to create a volcano plot to differentiate between statistically significant and nonsignificant genes.

Pathway and Functional Enrichment Analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was employed to carry out functional annotation and enrichment analysis (29). Key targets were examined using DAVID to predict functions across three categories: biological process (BP), molecular function (MF), and cellular component (CC). A threshold of an adjusted p-value ≤ 0.05 was used, and the top 10 Gene Ontology (GO) enrichments and the top 10 KEGG pathways with the highest counts were selected for further exploration.

Network Construction

Network analysis was conducted to elucidate the mechanisms through which V. unguiculata affects T2DM. This analysis was facilitated by Cytoscape 3.8.0 (30), a software platform that provides tools for importing, visualizing, and analyzing biomolecular interaction networks. In the constructed network, the nodes are the target genes, and the edges are the interactions between the active compounds and the target genes. To define the network topology, the degree centrality of each compound, target gene, or pathway was calculated using the Network Analyzer of Cytoscape. Targets that showed the maximum level of interconnectedness were termed ‘key targets’.

Protein–Protein Interaction Network Construction

PPIs are essential because they are highly flexible, tunable, and selective processes. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (31) was then used to search for the functional relationship between the above identified key targets, which requires the combined score to be more than 0.9, while the high confidence score was applied to enhance the reliability of the data. A high threshold was selected to ensure that interactions include all the relevant potential interactions to enhance the reliability of the network interaction.

These interactions involved the use of experimental data such as protein homology, or physical type of interaction, and co-expression data, with the intention to make a more biologically relevant network by employing these additional filters. The obtained PPI network was then subjected to the next step of filtration and analysis using the Cytohubba plugin in Cytoscape to identify essential regulatory genes and potential targets.

Virtual screening

The protein targets identified during the Protein-Protein Interaction Network were then further validated through molecular docking, which is a rapid computational technique. Information about the candidate target was obtained from the X-ray crystal structures offered under the RCSB Protein Data Bank (32). Structural optimization of the complex was done by UCSF Chimera (33), while protein targets and probable binding sites were predicted using the CASTp site mapping tool (34). Virtual screening of the reported compounds was performed using the AutoDock Vina in Ubuntu (35,36) in order to assess the binding of the core targets with the active compounds. The favorable docked conformations were considered in terms of binding energies. Docking scores, which reflect the strength and stability of the compound-target interactions, were used as primary criteria for selecting promising compounds and their targets. The results of the AutoDock vina were validated using the PyRx software to cross-check the predictive score of the ligand and protein interaction. This approach ensures that the docking results are robust using multiple docking methods. Visualization of these interactions was achieved using Discovery Studio (37), highlighting the dynamics of compound-protein interactions in a visual format.

MD simulation

After the virtual screening, the molecular dynamic simulation was performed to assess structural as well as dynamic stability of the ligand-protein complexes for 50 ns using GROMACS software 2024.1 (38). The TIP3P water box was used for solvation, and the CHARMM27, i.e., all-atom force field, was applied. The ligand topologies were generated using the SwissParam software. The system was minimized for 50000 steps to remove the steric clashes using the steepest descent model. Following the energy minimization, the system was equilibrated in two steps, i.e., NVT and NPT ensembles by applying the periodic boundary condition in all directions (36). During the NVT and NPT ensembles, the temperature was maintained at 300 K for 100 ps, while the NPT ensemble was maintained at 100 ps by keeping the pressure at 1 bar. The final MD run was performed for 100 ns, and the results were analyzed using various dynamical and structural properties (39).

Binding free energy calculations

To assess the binding free energy for all the complexes, MM-PBSA and MM-GBSA were performed (40). The binding free energy was calculated over the 100 ns simulation using MD trajectories, and the per-residue decomposition analysis was performed for both approaches, i.e., MM-PBSA and MM-GBSA, to determine the individual residues’ contribution towards the binding energy of the complexes. The binding free energy calculation was performed for 1000 frames (3000-4000), keeping the interval of 2. The various parameters that were assessed include EEL (Electrostatic energy), ESURF (Electrostatic surface potential energy), ENPOLAR (Polar solvation energy), GGAS, Van der Waals, and total energy (40).

Statistical analysis

The experiments were performed three times, and the results are presented as mean ± standard error of the mean. The IC₅₀ values were determined using IBM SPSS Statistics 21 software. The data underwent analysis using one-way ANOVA, followed by post hoc Dunnett’s multiple comparison test using GraphPad Prism version 8.0.2. The criteria for the statistical significance were chosen as a p-value less than 0.05. Dunnett’s post hoc was chosen because it adjusts for multiple comparisons inherently and tries to ensure the family-wise error rate control during multiple group comparisons with a single control while requiring no additional corrections. The effect sizes were applied to all ANOVA calculations for the magnitude effects quantifications. The various tests applied include the Epsilon-squared, Eta-squared, Random-effect, and Omega-Squared calculations with confidence intervals of 95 % to give comprehensive insight into the magnitude of effect.

Results

α-amylase inhibitory assay

The inhibitory efficacy of crude hydro-methanolic extract of V. unguiculata and its fractions against α-amylase enzyme was investigated using in vitro activity. Various concentrations, including 8, 4, 2, 1, 0.5, 0.25, 0.125, 0.0625, and 0.03125 mg/mL, were used to test the samples and determine their IC₅₀ values. The crude extract exhibited a significantly lower IC₅₀ value of 0.202 mg/mL compared to the aqueous (0.279 mg/mL), ethyl acetate (0.557 mg/mL), and n-hexane (1.294 mg/mL) fractions. In comparison, Acarbose, used as the standard, demonstrated the highest inhibition of the alpha amylase enzyme with an IC₅₀ value of 0.043 mg/ml, as depicted in Figure 1.

FIGURE 1 -. Representing alpha amylase inhibitory activity of (a) Crude extract, (b) Aqueous fraction, (c) Ethyl acetate fraction, (d) n-Hexane fraction, and (e) standard drug, Acarbose.

α-glucosidase inhibitory assay

The inhibitory potential against α-glucosidase enzyme using in vitro activity was examined on the crude hydro-methanolic extract of V. unguiculata and its fractions. Various concentrations, including 8, 4, 2, 1, 0.5, 0.25, 0.125, and 0.0625 mg/mL, were used to test the samples, and their IC₅₀ values were determined. The crude extract displayed a significantly lower IC₅₀ value (0.111 mg/mL) compared to the ethyl acetate (0.135 mg/mL), aqueous (0.243 mg/mL), and n-hexane (0.629 mg/mL) fractions among all the tested samples. Acarbose, used as the standard, exhibited the highest inhibition of alpha amylase enzyme with an IC₅₀ value of 0.034 mg/mL, as depicted in Figure 2.

FIGURE 2 -. Representing α-glucosidase inhibitory activity of (a) Crude extract, (b) Ethyl acetate fraction, (c) Aqueous fraction, (d) n-Hexane fraction, and (e) standard drug, Acarbose.

Aldose reductase inhibitory assay

V. unguiculata crude hydro-methanolic extract and its fractions were examined for their ability to inhibit aldose reductase enzyme through an in vitro assay. Various concentrations, including 300, 200, 100, 50, 25, and 10 µg/mL, were tested, and the IC₅₀ values were determined. The crude extract displayed the lowest IC₅₀ value (15.823 µg/mL) among all the samples tested, followed by ethyl acetate (25.806 µg /mL), aqueous (52.376 µg /mL), and n-hexane (106.857 µg /mL) fractions. Quercetin, used as the standard, exhibited the highest inhibition of alpha amylase enzyme with an IC₅₀ value of 8.545 µg /ml, as indicated in Figure 3.

FIGURE 3 -. Representing aldose reductase inhibitory activity of (a) Crude extract, (b) Ethyl acetate fraction, (c) Aqueous fraction, (d) n-hexane fraction, and (e) standard drug, Quercetin.

Network Pharmacology

Screening of potential compounds

The screening of potential compounds from V. unguiculata identified several candidates with drug-like potential, as summarized in Table 1. The compounds belong to various classes, including flavonoids, phenolic acids, isoflavones, triterpenoids, phytosterols, and alkaloids, each exhibiting different molecular weights, drug likeness (DL), and oral bioavailability (OB). Among the flavonoids, Quercetin, Kaempferol, Catechin, Luteolin, and Apigenin were identified, with molecular weights ranging from 270.24 to 302.23 g/mol. These compounds showed moderate drug likeness scores, with Catechin having the highest DL of 0.64, and all of them exhibited an oral bioavailability of 0.55, indicating their potential for oral administration. Protocatechuic acid, a phenolic acid with a molecular weight of 154.12 g/mol, displayed a lower drug likeness score of 0.23 but a slightly higher oral bioavailability of 0.56, suggesting that while its drug-like properties may be limited, it could still be bioavailable if administered orally. Isoflavones such as Genistein, Daidzein, and Glycitein were also included in the screening. These compounds have molecular weights between 254.24 and 284.26 g/mol and exhibit moderate drug likeness, with Genistein having a DL of 0.44. Like the flavonoids, they all have an oral bioavailability of 0.55. In the category of triterpenoids, Betulinic acid stands out with a significantly higher molecular weight of 456.7 g/mol, a drug likeness score of 0.25, and a notably higher oral bioavailability of 0.85, suggesting its strong potential as an orally administered drug candidate. The phytosterols identified include Beta-sitosterol, Campesterol, and Stigmasterol, with molecular weights around 400 to 414 g/mol. These compounds showed varying drug likeness scores, with Beta-sitosterol achieving the highest DL of 0.78, and all three exhibiting an oral bioavailability of 0.55. Finally, Vicine, an alkaloid with a molecular weight of 304.26 g/mol, demonstrated the highest drug likeness score of 0.84 among the compounds screened, coupled with an oral bioavailability of 0.55, highlighting it as a particularly promising candidate for further drug development. Overall, the screening results indicate that several compounds from V. unguiculata possess favorable drug-like properties and oral bioavailability, making them viable candidates for further investigation in drug development efforts. Similarly, the comprehensive pharmacokinetic studies of the selected compounds, including Stigmasterol, Quercetin, Luteoline, and Acarbose, are shown in the Supplementary Table 1.

Table 1 -. Summarizes the list of compounds obtained from V. unguiculata with drug-like potential

Compound	Compound Class	Molecular Weight (MW)	Drug Likeness (DL)	Oral Bioavailability (OB)
Quercetin	Flavonoid	302.23	0.52	0.55
Kaempferol	Flavonoid	286.24	0.5	0.55
Catechin	Flavonoid	290.27	0.64	0.55
Luteolin	Flavonoid	286.24	0.38	0.55
Apigenin	Flavonoid	270.24	0.39	0.55
Protocatechuic acid	Phenolic Acid	154.12	0.23	0.56
Genistein	Isoflavone	270.24	0.44	0.55
Daidzein	Isoflavone	254.24	0.29	0.55
Glycitein	Isoflavone	284.26	0.37	0.55
Betulinic acid	Triterpenoid	456.7	0.25	0.85
Beta-sitosterol	Phytosterol	414.7	0.78	0.55
Campesterol	Phytosterol	400.7	0.59	0.55
Stigmasterol	Phytosterol	412.7	0.62	0.55
Vicine	Alkaloid	304.26	0.84	0.55
Acarbose	Oligosaccharide	645.608		0.17

Microarray Data Analysis

Microarray gene expression data were analyzed using datasets obtained from the GEO database. Additionally, the GSE22435 dataset contained data from 7 affected individuals and 10 healthy controls (Fig. 4A and 4D). Differential expression analysis was conducted using the LIMMA package, resulting in the identification of a significant number of differentially expressed genes (DEGs) across these datasets. Specifically, the analysis of GSE22435 revealed 671 upregulated and 198 downregulated genes, leading to a total of 870 DEGs. In the GSE43950 dataset, 2,133 DEGs were identified (Fig. 4B and 4E), and the GSE92724 dataset yielded 1,160 DEGs (Fig. 4C and 4F). All identified DEGs met the stringent criteria of an adjusted p-value less than 0.05 and an absolute log fold change (|log(FC)|) of greater than or equal to 1. These findings underscore the substantial alterations in gene expression between affected individuals and controls, providing a robust basis for further biological interpretation and validation.

FIGURE 4 -. Differential Gene Expression Analysis for GEO Datasets GSE22435, GSE43950, and GSE92724. (A), (B), and (C) show UMAP plots, boxplots of gene expression, and adjusted p-value distributions for each dataset. (D), (E), and (F) present volcano plots highlighting upregulated (red) and downregulated (blue) genes.

Target Identification

A comprehensive target identification process was carried out, resulting in the identification of 1,310 targets for 14 active compounds using the STITCH and SwissTarget Prediction databases. These targets were then cross-referenced with the DEGs identified through the microarray data analysis. This comparison revealed 88 common genes that are implicated in both the disease context and the plant-derived compounds (Fig. 5A). The compound-target network analysis demonstrated that these 88 genes were targeted by multiple compounds, indicating their potential as druggable targets.

Network Analysis

Following the identification of the 88 common genes, a protein-protein interaction (PPI) network was constructed (Fig. 5B). The network consisted of 88 nodes and 144 edges, illustrating the complex interactions among these genes. Further analysis identified the top 10 genes with the highest degree of connectivity within the network (Fig. 5C), which included EGFR (degree score of 43), SRC (41), PTGS2 (39), MMP9 (33), BCL2 (33), TLR4 (29), GSK3B (26), PTPRC (24), KIT (22), and BCL2L1 (20). Due to their high connectivity, these genes may play central roles in the disease mechanism and serve as critical targets for therapeutic interventions.

FIGURE 5 -. Network and Target Analysis of Common Genes. (A) Venn diagram showing the overlap between compound-related targets (3,186) and disease-related targets (327), identifying 88 common genes. (B) Protein-protein interaction (PPI) network of the 88 common genes, illustrating the complex interactions among these targets. (C) Sub-network highlighting the top 10 genes based on degree of connectivity

Functional Annotation

Later, the functional annotation of these 88 common genes was performed using the DAVID software, which revealed their involvement in multiple diabetes-related pathways. The analysis identified significant clusters of genes associated with key biological processes (Fig. 6), including metabolic pathways such as prostanoid and fatty acid biosynthesis. Specifically, genes were enriched in processes like long-chain fatty acid metabolic processes, organic acid biosynthetic processes, and eicosanoid metabolic processes. Additionally, another cluster highlighted the role of these genes in retinoid metabolic processes and response to retinoic acid, emphasizing their involvement in important signaling pathways. Moreover, a smaller subset of genes was associated with oxidoreductase activity, particularly acting on CH-OH groups of donors and the incorporation of two atoms of oxygen, indicating their role in oxidation-reduction processes.

In terms of molecular function, the common genes were linked to activities such as ligand-activated transcription factor activity and nuclear receptor activity, which are crucial for regulating gene expression, particularly in response to retinoic acid. Other identified functions include protein autophosphorylation, essential for signal transduction, and protein tyrosine kinase activity, which plays a vital role in various cellular signaling pathways. The GO analysis also revealed that these genes are integral components of critical cellular structures, such as the Schaffer collateral-CA1-CA1 synapse, highlighting their role in synaptic transmission and plasticity in the brain. Additionally, other genes were located in the tertiary granule lumen and the intrinsic component of the presynaptic membrane, emphasizing their involvement in neurotransmitter release and synaptic function.

Further KEGG pathway analysis was also performed (Fig. 7). One of the most prominent pathways identified is the Arachidonic acid metabolism, which is crucial for the synthesis of bioactive lipid mediators that play significant roles in inflammation and cellular signaling. Another critical pathway is the HIF-1 signaling pathway, which is essential for cellular responses to hypoxia and is often implicated in cancer progression. The PPAR signaling pathway was also depicted, which plays an important role in regulating fatty acids and glucose homeostasis and is relevant to metabolic diseases. It also marks down pathways that relate to particular diseases, for example, non-small cell lung cancer or Bladder cancer; this means that these genes are involved directly in causing these types of cancer or in their advancement. Moreover, the information about the discovered VEGF signaling pathway, which plays an important role in the formation of new blood vessels and is widely used in cancer treatment. Other pathways under study include NF-kappa B signaling, via which most immune responses are controlled, Th17 cell differentiation during autoimmune diseases, etc. The occurrence of the Adherens junction pathway was also observed for cell-cell adhesion, underlining the importance of these genes in supporting the cell’s structure.

In general, the KEGG pathway analysis shows that the common genes are found in various signaling pathways and biological processes that are known to be important for maintaining cellular homeostasis and regulating metabolic and disease processes. This type of analysis is helpful in understanding the possible functions of these genes from a human disease perspective and might inform future directions on treatments and cures. Apparently, the network representation strengthens the claims made with respect to the connectivity of these genes in the designated GO classes and proves the multifaceted nature of these genes in terms of their functions in various biological processes, molecular functions, and cellular components.

FIGURE 6 -. Gene Ontology Analysis of 88 common genes representing their involvement in multiple biological processes, cellular components, and molecular functions.

Compound-target pathway network

The compound-target pathway network analysis (Fig. 8) revealed that GSK3B, PTGS2, and TLR4 are the primary target proteins, displaying the highest connectivity within the network and interacting with the most bioactive compounds. The high degree of connectivity values means that these proteins are core to the network and highlight them as potential targets for treatment. Furthermore, the investigation also pointed to the fact that many compounds, including those from V. unguiculata, interact with these proteins, making them valuable targets for therapeutic agents. Notably, GSK3B, PTGS2, and TLR4 are known to be associated with type 2 diabetes; GSK3B was reported to be a key protein regulating insulin signaling and glucose homeostasis, PTGS2 participates in inflammation aggravating insulin resistance, and TLR4 was reported to be involved in immune response and chronic inflammation in the development of T2D. These observations indicate that modulating the expression of GSK3B, PTGS2, and TLR4 using the active compounds isolated from V. unguiculata may be a potential therapeutic approach for type 2 diabetes.

Molecular docking analysis

The docking results provide insights into the binding affinities of various compounds with three target proteins: GSK3B, TLR4, and PTGS2, as shown in Table 2. These interactions are critical for evaluating the efficacy of these compounds as inhibitors or modulators of these proteins. The Virtual screening results showed Stigmasterol had the highest binding affinity with the 8AV1 (GSK3B), Luteolin exhibited the highest binding affinity with the 4PH9 (PTGS2), while Quercetin showed the highest binding energy with the 3FXI (TLR4). The ligand-protein complexes were further validated using MD simulation and binding free energy calculations. During the virtual screening of the compounds with the 4PH9, the center of the grid box includes x = 20.685, y = 31.547, and z = 27.584, while the size of the grid box includes x = 80, y = 102, and z = 72. Similarly, the virtual screening of the compounds with the 8AV1, the center of the grid box was kept as x = −3.818, y = 0.936, and z = −2.868, while the size of the grid box was kept as x = 60, y = 64, and z = 72. Furthermore, the virtual screening of the compounds with the 3FXI, the center of the grid box was kept as x = 4.536, y = −17.363, and z = 20.067, while the size of the grid box was kept as x = 60, y = 40, and z = 100. The molecular docking of Luteoline with the 4PH9 showed three hydrogen bonds (ASNA39, GLUA466, and GLYA136), two pi-alkyl bonds (LEUA153 and PROA154), two carbon-hydrogen bonds (CYSA47 and ARGA470), and one unfavorable acceptor-acceptor bond (TYRA131). The molecular docking of the Stigmasterol with the 8AV1 showed seven pi-alkyl bonds (ILEA37, TYRA109, LEUA163, ALAA58, VALA45, LYSA60, and LEUA153), and one carbon-hydrogen bonds (ASPA175). The molecular docking of Quercetin with 3FXI showed one hydrogen bond (LYSA204), two pi-anion bonds (ASPA155 and ASPA183), one carbon-hydrogen bond (HISA153), and one unfavorable acceptor-acceptor bond (SERA183). The molecular docking of Acarbose was conducted against the 8AV1, 3FXI, and 8AV1, and results showed that the binding energy of Acarbose with the 8AV1 is −6.9 kcal/mol, with 3FXI is −6.7 kcal/mol, and with the 4PH9 is −8.9 kcal/mol. The Acarbose showed the highest binding energy against the inflammatory target, i.e., PTGS2; however, against the GSK3B and TLR-4, the Acarbose showed lower binding energy as compared to the other compounds Table 3. The compounds and protein interactions are shown in three-dimensional and two-dimensional views using Discovery Studio visualizer_16 (Fig. 9).

FIGURE 7 -. KEGG Pathway analysis demonstrating the involvement of selected genes in diabetes relevant pathways.

FIGURE 8 -. Compound-Target-Pathway network. The green octagon represents the compound, the red circle represents the hub genes, and the blue square represents the pathways targeted by hub genes.

Table 2 -. Docking results of various compounds with 8AV1 (GSK3B), 3FXI (TLR4), and 4PH9 (PTGS2) proteins. The table summarizes the binding affinities (in kcal/mol)

Compound Name	Binding Affinity (kcal/mol)	Binding Affinity (kcal/mol)	Binding Affinity (kcal/mol)
Stigmasterol	−8.854	−8.679	−7.244
Campesterol	−8.833	−8.062	−7.084
Beta-sitosterol	−8.739	−8.052	−7.149
Daidzein	−8.198	−7.088	−6.667
Quercetin	−8.146	−9.122	−7.237
Apigenin	−8	−8.941	−7.066
Luteolin	−7.982	−9.297	−6.885
Kaempferol	−7.782	−8.619	−6.732
Genistein	−7.685	−8.81	−6.865
Glycitein	−7.645	−7.297	−6.864
Catechin	−7.637	−8.098	−7.159
Vicine	−6.17	−7.04	−6.479
Acarbose	−6.9	−6.7	−8.9

MD simulation analysis

During the entire 100 ns, it was found that all ligand-protein complexes remained stable during the simulation, and no significant deviation in the structural conformation of the ligand-protein complex was noticed. The ligand-protein complex stability was measured using various key parameters such as RMSD, RMSF, RoG, etc., and it was noticed that within the first few ns, the system reached equilibrium and thereafter remained stable. The RMSD values were calculated to assess the convergence of the system to its equilibrium conformation. During the RMSD calculation, initially, the system showed deviation for a few ns due to thermodynamic shock, but thereafter, the RMSD showed no significant fluctuations for all the complexes, as shown in Figure 10. The RMSD analysis showed that the average RMSD of the 3FXI_Quercetin was 0.25828, the average RMSD of the 4PH9_Luteoline was 0.45871, while the average RMSD of the 8AV1_Stigmasterol was 0.28030 nm. The RMSF analysis was conducted to study the behavior, i.e., flexibility of the individual amino acid of the protein following binding with the ligand. The RMSF analysis showed that complexes remain stable during the simulation, and no significant fluctuations were noticed, as shown in Figure 11. The RMSF analysis showed that the average RMSF of the 3FXI_Quercetin was 0.14277, the average RMSF of the 4PH9_Luteoline was 0.10841, while the average RMSF of the 8AV1_Stigmasterol was 0.11547 nm. The compactness of the protein structure following ligand binding was observed using the radius of gyration (RoG). The RoG analysis showed that the complexes remained stable during the entire simulation, which further validated the results of the RMSD and RMSF. Minor fluctuations have been observed in the complexes, which correlate with the small-scale motions and conformational changes associated with the ligand binding with proteins Figure 12. The RoG analysis showed that the average of the 3FXI_Quercetin complex was 3.1447, the average RMSF of the 4PH9_Luteoline was 2.4579, while the average RoG of the 8AV1_Stigmasterol was 2.1969 nm. The changes in the ligand-protein complex surface following exposure to the solvent during the 100 ns simulation were studied using SASA. Initially, the SASA showed minor fluctuations, and thereafter the complexes achieved stability as shown in Figure 13. The hydrogen bonds analysis was conducted during the entire simulation to study the stability of ligand-protein complexes. The number of hydrogen bonds fluctuated for all complexes during the course of the simulation; however, the maximum number of hydrogen bonds was retained during the simulation, as shown in Figure 14.

FIGURE 9 -. Molecular docking analysis of key compounds with target proteins TLR4 (3FX1), PTGS2 (4PH9), and GSK3B (8AV1). Based on the virtual screening, stigmasterol showed binding affinity with the GSK3B, luteoline exhibited the highest binding energy with the PTGS2, and quercetin showed the highest binding energy with the TLR4. The interaction of the ligands with the target proteins was presented in two-dimensional and three-dimensional views using Discovery Studio Visualizer.

Binding Free Energy parameters

The MM-PBSA and MM-GBSA analyses were performed for all the complexes to determine their stability in terms of binding energy. During the binding free energy calculations of the 8AV1_Stigmastrol using MM-PBSA revealed that the van der Waals energy, EEL, ENPOLAR, and GGAS remain negative and exhibited favorable binding energy. However, EPB energy and GSOLV energy of the system remain positive and indicate unfavorable binding free energy. However, the total energy of the system remains negative and indicates favorable binding energy, as shown in Table 3 and Supplementary Figure 1. Similarly, the MM-PBSA results were validated using the MM-GBSA analysis.

FIGURE 10 -. The RMSD parameter following the MD simulation using the compound with the highest binding affinity. The stigmasterol showed the highest binding affinity with the 8AV1, Luteolin with 4PH9, and Quercetin with the 3FXI. The RMSD analysis showed that all complexes showed an acceptable value of RMSD, and no significant changes were observed.

Figure 11 -. The RMSF analysis was performed for all the complexes following the simulation. The results showed no major fluctuation in the RMSF value, and the results indicate that all complexes remain stable.

During the binding free energy calculations of the 8AV1_Stigmastrol using MM-GBSA revealed that the van der Waals energy, EEL, ESURF, and GGAS remain negative and exhibited favorable binding energy. However, the EGB energy and GSOLV energy of the system remain positive and indicate unfavorable binding free energy. However, the total energy of the system is negative and indicates favorable binding energy, as shown in Table 4 and Supplementary Figure 2.

FIGURE 12 -. The RoG of all the complexes was determined following the simulation to assess the stability of the complexes. The RoG results showed that 8AV1_Stigmasterol, and 4PH9_Luteolin complexes are lowered as compared to the 3FXI_Quercetin complex.

FIGURE 13 -. The 8AV1_Stigamsterol result revealed the lowest SASA as compared to the 4PH9_Luteolin and 3FXI_Quercetin complexes.

FIGURE 14 -. The 3FXI_Quercetin showed the highest number of hydrogen, i.e., five, as compared to the 8AV1_Stigmasterol and 4PH9_Luteolin complex.

Table 3 -. The MM-PBSA result of the 8AV1_Stigamsterol complex

Frames	VDWAALS	EEL	EPB	ENPOLAR	GGAS	GSOLV	TOTAL
Average	−39.97	−8.1	29.08	−4.52	−48.07	24.56	−23.51
SD	3.34	5.49	4.09	0.21	6.2	4.03	4.8
SEM	0.15	0.25	0.18	0.01	0.28	0.18	0.21

Table 4 -. The MM-GBSA result of the 8AV1_Stigamsterol complex

Frames	VDWAALS	EEL	EGB	ESURF	GGAS	GSOLV	TOTAL
Average	−39.97	−8.1	21.8	−5.26	−48.07	16.54	−31.52
SD	3.34	5.49	3.5	0.45	6.2	3.36	4.2
SEM	0.15	0.25	0.16	0.02	0.28	0.15	0.19

Table 5 -. The MM-PBSA result of the 4PH9_Luteolin complex

Frames	VDWAALS	EEL	EPB	ENPOLAR	GGAS	GSOLV	TOTAL
Average	−30.05	−20.24	43.36	−3.45	−50.29	39.91	−10.38
SD	2.24	3.66	4.34	0.07	3.99	4.32	2.78
SEM	0.1	0.16	0.19	0	0.18	0.19	0.12

Table 6 -. The MM-GBSA result of the 4PH9_Luteolin complex

Frames	VDWAALS	EEL	EGB	ESURF	GGAS	GSOLV	TOTAL
Average	−30.05	−20.24	35.12	−4.74	-50.29	30.38	−19.91
SD	2.24	3.66	3.43	0.17	3.99	3.34	2.03
SEM	0.1	0.16	0.15	0.01	0.18	0.15	0.09

Table 7 -. The MM-PBSA result of the 3FXI_Quercetin complex

Frames	VDWAALS	EEL	EPB	ENPOLAR	GGAS	GSOLV	TOTAL
Average	−21.08	−46.49	48.11	−2.8	−67.58	45.31	−22.27
SD	2.85	5.84	4.45	0.06	4.95	4.45	2.8
SEM	0.13	0.26	0.2	0	0.22	0.2	0.13

The MM-PBSA and MM-GBSA analyses were performed for all the complexes to determine their stability in terms of binding energy. During the binding free energy calculations of the 4PH9_Luteoline using MM-PBSA revealed that the van der Waals energy, EEL, ENPOLAR, and GGAS remain negative and exhibited favorable binding energy. However, EPB energy and GSOLV energy of the system remains positive and indicate unfavorable binding free energy. However, the total energy of the system remains negative and indicates favorable binding energy, as shown in Table 5 and Supplementary Figure 3. Similarly, the MM-PBSA results were validated using the MM-GBSA analysis.

During the binding free energy calculations of the 4PH9_Luteoline using MM-GBSA revealed that the van der Waals energy, EEL, ESURF, and GGAS remain negative and exhibited favorable binding energy. However, the EGB energy and GSOLV energy of the system remains positive and indicate unfavorable binding free energy. However, the total energy of the system is negative and indicates favorable binding energy, as shown in Table 6 and Supplementary Figure 4.

The MM-PBSA and MM-GBSA analyses were performed for all the complexes to determine their stability in terms of binding energy. During the binding free energy calculations of the 3FXI_Quercetin using MM-PBSA revealed that the van der Waals energy, EEL, ENPOLAR, and GGAS remain negative and exhibited favorable binding energy. However, EPB energy and GSOLV energy of the system remains positive and indicate unfavorable binding free energy. However, the total energy of the system remains negative and indicates favorable binding energy, as shown in Table 7 and Supplementary Figure 5. Similarly, the MM-PBSA results were validated using the MM-GBSA analysis.

During the binding free energy calculations of the 3FXI_Quercetin using MM-GBSA revealed that the van der Waals energy, EEL, ESURF, and GGAS remain negative and exhibited favorable binding energy. However, the EGB energy and GSOLV energy of the system remains positive and indicate unfavorable binding free energy. However, the total energy of the system is negative and indicates favorable binding energy, as shown in Table 8 and Supplementary Figure 6. The binding free energy calculations revealed that 8AV1_Stigmaterol showed the highest total binding free energies (MM-PBSA = -23.51, MM-GBSA = −31.52), followed by the 3FXI_Quercetin (MM-PBSA = −22.27, MM-GBSA = −30.56), and the lowest binding free energies were shown by the 4PH9_Luteoline (MM-PBSA = −10.38, MM-GBSA = −19.91).

The contribution of the individual residues was determined using per-residue decomposition analysis for all the complexes. The per-residue decomposition analysis showed that VAL45 and ILE37 were the highest energy contributing amino acids in the 8AV1_Stigmasterol complex using MM-PBSA and MM-GBSA, respectively, Supplementary Figure 7. Meanwhile, in 4PH9_Luteoline complex decomposition analysis, it was revealed that CYS36 and PRO154 were the highest contributing amino acids using MM-PBSA and MM-GBSA parameters, respectively, Supplementary Figure 8. The per-residue decomposition analysis of the 3FXI_Quercetin showed that GLU128 was the highest contributing residue using MM-PBSA and MM-GBSA analysis, as shown in the Supplementary Figure 9.

Discussion

The use of natural remedies as additional treatments alongside conventional therapies for the treatment of diabetes and its associated complications is increasing worldwide, with many plants in different countries known to have effects in lowering blood sugar (41). Indian literature from the past has documented the anti-diabetic properties of more than 800 plants, while ethnopharmacological research has shown that over 1200 plants can be used for their anti-diabetic effects (42). The primary reason for postprandial hyperglycemia is the action of α-amylase and α-glucosidase enzymes, which break down carbohydrates. The process of breaking down carbohydrates starts with α-amylase, which breaks the 1,4-glycosidic bonds in polysaccharides, converting them into disaccharides. Subsequently, α-glucosidase transforms the disaccharides into monosaccharides, resulting in elevated blood glucose levels after meals (43). Delaying the digestion of carbohydrates and reducing postprandial glucose levels, inhibitors of alpha-amylase and alpha-glucosidase aid in controlling high plasma glucose levels. Investigators typically seek out natural substances that can control elevated blood sugar levels by reducing the activity of alpha-glucosidase and alpha-amylase, leading to decreased adverse effects such as diarrhea, nausea, and abdominal pain (44).

Table 8 -. The MM-GBSA result of the 3FXI_Quercetin complex

Frames	VDWAALS	EEL	EGB	ESURF	GGAS	GSOLV	TOTAL
Average	−21.08	−46.49	40.85	−3.84	−67.58	37.02	−30.56
SD	2.85	5.84	4.6	0.12	4.95	4.61	2.3
SEM	0.13	0.26	0.21	0.01	0.22	0.21	0.1

During an in vitro study, researchers analyzed cowpea protein hydrolysates (PH) and ultrafiltered peptide fractions (UFPF) for their ability to inhibit enzymes like α-amylase, α-glucosidase, and dipeptidyl peptidase IV. The findings indicated that both PH and UFPF demonstrated the most potent inhibitory activity against the enzymes under investigation (45). A different research study found that cowpea seed protein hydrolysates demonstrated greater effectiveness in inhibiting α-amylase (46). The efficacy of sorghum-cowpea composite biscuits in inhibiting starch hydrolysis enzymes like α-amylase and α-glucosidase was examined. This combination could potentially contribute to the reduction of postprandial hyperglycemia by targeting specific enzymes (47). Aldose reductase can be found in almost all mammalian cells, with significantly higher levels in organs such as the cornea, retina, lens, myelin sheath, kidney, and sciatic nerves. These organs are prone to diabetic complications (48). The accumulation of polyol in lens fibers, due to the heightened activity of the polyol pathway, results in the influx of water and the generation of osmotic stress. Consequently, this leads to tissue damage through osmotic swelling, changes in membrane permeability, and oxidative stress. The polyol pathway has garnered considerable focus in the clinical treatment of secondary diabetes complications because of these outcomes (49). The breakdown of glucose through the polyol pathway by enzymes leads to the production of sorbitol, which is a harmful substance for tissues as it increases osmotic pressure and can cause tissue swelling. By inhibiting the activity of aldose reductase, which is mainly present in the lens’s epithelial cells, it is possible to prevent the development of cataracts (23). The literature does not contain any data on the in vitro anti-diabetic properties, such as the inhibition of α-amylase, α-glucosidase, and aldose reductase activities, of the entire V. unguiculata plant. Therefore, our current study aims to assess the ability of the whole V. unguiculata plant to inhibit alpha amylase, alpha glucosidase, and aldose reductase enzymes. Among all the samples tested, the crude extract exhibited the most significant inhibitory activity against alpha amylase, alpha glucosidase, and aldose reductase with IC₅₀ values of 0.202 mg/mL, 0.111 mg/mL and 15.823 µg/mL, respectively.

Moreover, the current study also used a network pharmacology approach to come up with a complex network connecting the compounds of V. unguiculata and the potential protein targets. Unlike most previous studies that focus on single-target drugs, our approach considers multiple targets and their interactions, including those with T2DM. For instance, previous works have defined involvement of proteins such as GSK3B in insulin signaling and glucose metabolism, but this work affords elaborate information to how multiple compounds influence GSK3B, PTGS2, and TLR4 at the same time. This multiple-target approach is in harmony with modern practices in network pharmacology, where the impact of drug candidates on multiple interconnected biological targets is researched. The GO and KEGG study provides insight into important signaling, including PTGS2, GSK3B, and TLR4, which are influenced by the V. unguiculata constituents and evidently support the potential multi-target anti-diabetic activity. However, the multi-target analysis conducted in this study increases the likelihood of off-target effects and paves the way for undesired biological responses. Furthermore, the non-specific interactions of the studied protein with multiple signaling protein may alter biological responses and produce adverse effects. This study focused on the computational predictions and in vitro analysis; however, to validate the safety and specificity of the compounds, in vivo and in vitro studies will be required. GSK3B is implicated in insulin signaling and glucose homeostasis. Several studies have revealed that GSK3B knockdown promotes insulin signaling and corrects glucose intolerance in diabetic animal models. Inhibition of GSK3B is important for regulating hyperglycemia and insulin resistance since it is targeted by these inhibitors (50). Similarly, PTGS2, commonly known as COX-2, is an enzyme that plays a significant role in inflammation (51). Inflammatory processes play a crucial role in the evolution and progression of insulin resistance and consequent type 2 diabetes. Select members of the COX-2 enzyme have been studied as an anti-inflammatory, which may have beneficial effects on other metabolic complications related to diabetes.

In contrast, TLR4 contributes to the immune system, and its overactivation is suggested to induce the inflammation that drives insulin resistance (52). It has been proposed that using specific inhibitors or antagonists of TLR4 may present a useful strategy for the modulation of chronic inflammation and the improvement of insulin resistance in type 2 diabetes. Multi-target therapy in this case entails the identification of a number of pathways or proteins that are related to the disease and treating all of them at once. This is thought to be a better approach to treating various intricate conditions like type 2 diabetes by combating various causes.

In the future, the findings of this study can provide a good foundation for subsequent research into the medical value of compounds of extracts from V. unguiculata. In vivo studies can also be undertaken to validate the findings made in silico and further examine the PK/PD characteristics of these compounds. We also supposed that the other target proteins may be required for further research, and they may discover the new multi-target drug for the better treatment of T2DM, especially GSK3B, PTGS2, and TLR4. Further insights on how these compounds collectively function might be helpful in designing improved combination therapies, which could be beneficial in patients. Hence, this work not only opens more evidence about network pharmacology as a strategy for searching the targets and new drugs’ designs but also shows the application of plant extracts as a treatment for chronic diseases.

Conclusion

In the present study, the V. unguiculata was evaluated for anti-diabetic potential using in vitro, network pharmacology, and computational approaches. The results of the in vitro analysis using V. unguiculata phytochemicals showed significant activity against the diabetes targets such as α-amylase, α-glucosidase, and aldose reductase enzymes concerned with glucose metabolism. Furthermore, the network pharmacology studies showed 88 shared target genes associated with T2DM, and the possible drug targets, i.e., GSK3B, PTGS2, and TLR4 were deemed as hub genes in the PPI network. The in vitro and network pharmacology studies were followed by the computational studies using virtual screening, MD simulation, and binding free energies calculations. The virtual screening showed the phytochemicals, i.e., Quercetin, Luteolin, and Stigmasterol, with the highest binding energy against the diabetes targets. The MD simulation studies were conducted following the virtual screening to assess the dynamic stability, and the results showed that the complexes remain stable. The binding free energy studies using MM-PBSA, and MM-GBSA exhibited negative and favorable binding energies, hence indicating the stability of the complexes. In conclusion, the study showed that V. unguiculata phytochemicals exhibited significant activity against the diabetes targets using in vitro, network pharmacology, and comprehensive computational approaches; however, to employ the V. unguiculata phytochemicals clinically additional studies such as comprehensive in-vivo approaches will be required.

Limitations of the study

The limitations of this study include the lack of in vivo validation despite the fact that the present study includes in vitro analysis, network pharmacology approaches, and comprehensive computational approaches. In order to validate the results of the study, the in vivo studies, including pharmacokinetics, pharmacodynamics, and toxicities associated with the top hits, should be conducted. Furthermore, the study mainly focused on identifying the molecular targets and binding affinities; the stability of the compounds in the in vivo system, such as metabolism, protein binding, and solubility, was not studied. Similarly, the computational studies offer key insight into the ligand’s interaction with the target protein; however, the computational studies are influenced by various factors such as the type of force field used and the input of the structure provided. Furthermore, the study also lacks in vitro and in vivo studies to assess the off-target effect and to explore the safety profile of the compounds; thus, to establish the specificity, selectivity, and safety of the studied compounds, comprehensive in vitro and in vivo studies will be required.

Authors contributors

HS: Data curation, Writing – original draft; Fatima Noor: Data curation, Writing – original draft; SMS: Data analysis, Writing – original draft; JSAO: Data curation, Writing – original draft; AUK: Data curation, Writing – original draft, MD: Data curation, Writing – final draft; HK: Designing, supervision, Writing – final draft

Acknowledgment

The authors are thankful to HEC Pakistan for financial support under project No: (Ref No. 20-16097/NRPU/R&D/HEC/2021 2021). The authors extend their appreciation to the Princess Nourah Bint Abdulrahman University Researchers Supporting Project number (PNURSP2025R13), Princess Nourah Bint Abdulrahman University, Saudi Arabia.

Data availability statement

The data that support the findings of the present study are available from the corresponding authors upon reasonable request.

Other information

This article includes supplementary material.

Corresponding author:

Haroon Khan

email: haroonkhan@awkum.edu.pk; haroonkhan@korea.ac.kr

Disclosures

Conflict of interest: The authors of this article have a conflict of interest

Financial support: This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors