Defining the Role of Isoeugenol from Ocimum tenuiflorum against Diabetes Mellitus-Linked Alzheimer’s Disease through Network Pharmacology and Computational Methods

The present study involves the integrated network pharmacology and phytoinformatics-based investigation of phytocompounds from Ocimum tenuiflorum against diabetes mellitus-linked Alzheimer’s disease. It aims to investigate the mechanism of the Ocimum tenuiflorum phytocompounds in the amelioration of diabetes mellitus-linked Alzheimer’s disease through network pharmacology, druglikeness and pharmacokinetics, molecular docking simulations, GO analysis, molecular dynamics simulations, and binding free energy analyses. A total of 14 predicted genes of the 26 orally bioactive compounds were identified. Among these 14 genes, GAPDH and AKT1 were the most significant. The network analysis revealed the AGE-RAGE signaling pathway to be a prominent pathway linked to GAPDH with 50.53% probability. Upon the molecular docking simulation with GAPDH, isoeugenol was found to possess the most significant binding affinity (−6.0 kcal/mol). The molecular dynamics simulation and binding free energy calculation results also predicted that isoeugenol forms a stable protein–ligand complex with GAPDH, where the phytocompound is predicted to chiefly use van der Waal’s binding energy (−159.277 kj/mol). On the basis of these results, it can be concluded that isoeugenol from Ocimum tenuiflorum could be taken for further in vitro and in vivo analysis, targeting GAPDH inhibition for the amelioration of diabetes mellitus-linked Alzheimer’s disease.


Introduction
There have been several insights into the exploration of the biomolecular link between type 2 diabetes mellitus (T2D) and Alzheimer's disease (AD) in the last two decades [1]. However, there is no specific molecular pathway or target that has been identified for the development of AD in diabetic patients. Epidemiological details suggest a strong relationship between high blood glucose levels and the degree of cognitive dysfunction in diabetic patients [2,3]. A few studies have collectively termed these conditions as type 3 diabetes [4,5].
For the treatment of T2D-linked AD, there are no specific therapeutic applications available. However, biguanides, sulfonylureas, glitazones, and DPP-4 inhibitors are known to reduce the symptoms of AD. However, most of these drugs have been reported with adverse effects, such as weight loss, lactic acidosis, ketoacidosis, anemia, bone fractures, and gastrointestinal and cardiovascular complications. Therefore, it is worth testing phytoconstituents with antidiabetic potential, which are less toxic compared to synthetic drugs [6,7].
Ocimum sanctum, commonly known as tulsi, is a fragrant shrub of the basil family Lamiaceae (Tribeocimeae) that is native to the eastern globe tropics, and which is said to have originated in north-central India. Tulsi is supposed to help with disease prevention, general health, wellbeing, and longevity, as well as dealing with the demands of daily life [8]. Specifically, it is also said to have an ameliorating effect on several health maladies [9][10][11]. Several in vitro, animal, and human experiments have proven the pharmacological benefits of tulsi [10][11][12]. Therefore, O. sanctum serves as an important source of phytoconstituents that can be utilized for pharmacotherapeutic applications.
With the advancements of computer-aided drug design technologies, several computational investigations have identified phytoconstituents as an effective replacement for the available therapeutics [13]. Additionally, network pharmacology has been employed to study the biological network of drug candidates, and to construct a poly-target drug molecule to optimize medication efficacy [14][15][16]. Hence, we aim to investigate the antidiabetic as well anti-Alzheimer's potential of O. sanctum phytochemicals through network pharmacology and other computational methods.

Prediction of Target Genes
The genes associated with type 2 diabetes and AD were obtained using DisGeNET (https://www.disgenet.org/search) (Accessed on 5 January 2022) and the GeneCards database (https://www.genecards.org/) (Accessed on 3 January 2022). The target genes were searched for using keywords "T2D" and "Alzheimer's", using "Homo sapiens" as a filter for both.

Molecular Docking Simulation
The 3D crystal structure of the target protein GAPDH (PDB ID: 1U8F) was downloaded from the RCSB Protein Data Bank (https://www.rcsb.org/) (Accessed on 10 January 2022). The molecular docking protocol was validated by redocking the co-crystallized ligand with the target protein and evaluating the RMSD between the bound and re-docked ligand. The molecular docking simulation was carried out according to the previous study conducted by   [21]. For the molecular docking simulation, AutoDock Tools 1.5.6 was used to prepare the protein structures. AutoDock is a suite of automated docking tools. It is designed to predict how small molecules, such as substrates or drug candidates, bind to a receptor of a known 3D structure [22]. To purify the protein structures, water and heteroatoms were removed. Polar hydrogens, on the other hand, were added to stabilize the same. The energy of the protein structures was reduced by using Kollmann-united charges and Gasteiger charges. After energy minimization, all atoms were assigned an AutoDock 4 atom type prior to obtaining the prepared protein structure in the PDBQT format for molecular docking simulations. The PDBQT format saves the protein's atomic coordinates, partial charges, and AutoDock 4 atom type [23]. Further, OpenBabel 2.3.1 was used to transform the 3D SDF structures into PDBQT format during ligand preparation [24]. For ligand preparation, the structures were imported into AutoDock Tools 1.5.6. To reduce the energy of the ligands, Kollman-united and Gasteiger charges were added. As indicated by AutoDock Tools 1.5.6, the number of torsions was left at the default [25].
Binding site prediction for the target proteins was completed using a literature analysis. The binding residues present in the binding pocket were identified using the literature available for GADPH [26]. Using AutoDock Tools 1.5.6, the binding pockets of the respective proteins were set in the grid boxes. The grid box measuring 11.20 Å × 11.20 Å × 11.20 Å was coordinated at x = 12.128023, y = 24.467386, and z = 29.015455. Virtual screening of the phytocompounds was completed with a command line-based software known as AutoDock Vina 1.1.2. For the perturbation and local optimization of ligands into the target site, it employs the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, which analyses the scoring function of each ligand conformation [27]. Because of the number of torsions generated during ligand production, ligands were assumed to be flexible throughout the docking simulation, while the proteins were believed to be rigid. However, 10 degrees of freedom were permitted for ligand molecules. The first binding poses with zero rootmean-square deviation (RMSD) of atomic positions are deemed to be highly valid out of ten binding poses generated [28]. They also possess the strongest binding affinity of all the positions, implying that the binding is stronger. Biovia Discovery Studio Visualizer 2021, an open-source visualizing GUI software, was used to visualize the results. The extent of ligand interaction was determined using binding affinity, total number of bonds, and respective hydrogen bonds [29].

Molecular Dynamics Simulation
The docked conformations of GAPDH protein and respective ligands (isoeugenol and myrene) with the lowest negative binding affinity were chosen for the molecular dynamics simulation, which was carried out according to the previous study conducted by   [30]. Myrene was chosen as the negative control, since the molecule was predicted to have the highest positive binding affinity (inferior binding) with GAPDH. For the simulation, the GROMACS-2018.1 biomolecular software suite was utilized [31]. All of the protein-ligand complexes were assigned using the CHARMM36 force field, and the ligand topology was acquired using the CGenFF server [32]. The GROMACS pdb2gmx module was used to add hydrogen atoms to the heavy atoms present. The steepest descent technique was then utilized to achieve 5000 vacuum minimization steps. The protein-ligand complexes were placed in a box with a 10 Å perimeter. The solvent was included in the water model TIP3P. The entire system was neutralized by adding the proper amount of Na + and Cl − counter ions. A total of 3 simulations were run for 100 ns, including a GAPDH protein backbone atom (5051 residues), GAPDH protein-isoeugenol complex (5065 residues), and GAPDH protein-myrene complex (5062 residues). The energy of the generated systems was decreased using the steepest descent and conjugate gradient techniques. The NVT ensemble was then followed by an NPT ensemble, which was followed by a brief (1000 ps) equilibration (1000 ps). All simulations took 100 ns at 310 K temperature and 1 bar pressure. A trajectory analysis consisting of root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), and solvent accessible surface area (SASA), and ligand hydrogen bond parameters was performed using the XMGRACE software, with the results displayed in graphical form [33,34].

Binding Free Energy Calculations
The determination of a protein-ligand complex's binding free energy is another application of molecular dynamics simulations and thermodynamics for determining the extent of ligand binding with protein. The binding free energy calculations for all the protein-ligand complexes were evaluated according to the previous study by   [30]. The molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA) technique was used in this work. It is another application of molecular dynamics simulations and thermodynamics for determining the extent of ligand binding with protein [35,36]. The g_mmpbsa program with MmPbSaStat.py script, which utilizes the GROMACS 2018.1 trajectories as input, was used to determine the binding free energy for each ligand-protein combination. In the g_mmpbsa program, three components are used to calculate the binding free energy: molecular mechanical energy, polar and apolar solvation energies, and molecular mechanical energy. The calculation is performed using MD trajectories of the last 50 ns, which were considered to compute ∆G with dt 1000 frames. It is evaluated using molecular mechanical energy and polar and apolar solvation energies. The equation (I) and (II) that are used to calculate the free binding energy are given below [21,30]. respectively, ∆G: standard free energy, ∆E MM : average molecular mechanics potential energy in vacuum, G Solvation : solvation energy, ∆E: total energy of bonded as well as nonbonded interactions, ∆S: change in entropy of the system upon ligand binding, and T: temperature in Kelvin.

Target Prediction and Active Compound Library
A total of 61 compounds of Ocimum tenuiflorum were retrieved from the PubChem database. On the basis of the ADMET (Table 3) screening, an active compound library was built, and a total of 26 compounds passed the screening criteria. Further, the target genes were fished for using the DisGeNET database, and a total of 6531 were considered, where 3134 belong to T2D and 3397 belong to AD, respectively. However, only 1381 genes were found to be common in both the diseases (Figure 1). in Kelvin.

Target Prediction and Active Compound Library
A total of 61 compounds of Ocimum tenuiflorum database. On the basis of the ADMET (Table 1) screeni built, and a total of 26 compounds passed the screenin were fished for using the DisGeNET database, and a t 3134 belong to T2D and 3397 belong to AD, respective found to be common in both the diseases (Figure 1).

Construction of PPI Network
The interaction network was constructed for the target proteins. Out of 1381 common genes, the network was constructed for 1287 targets which were specific to Homo sapiens. The constructed network contained 1287 nodes and 5685 edges with the 8.83 as the average node degree and 0.564 as the coefficient. The pictorial representation of the diseases interaction network is given in Figure 2.   The interaction network was constructed for the target proteins. Out of 1381 common genes, the network was constructed for 1287 targets which were specific to Homo sapiens. The constructed network contained 1287 nodes and 5685 edges with the 8.83 as the average node degree and 0.564 as the coefficient. The pictorial representation of the diseases interaction network is given in Figure 2.

Network Analysis
The disease gene interaction network obtained from STRING was imported to Cytoscape 3.8.2 for further analysis. The interaction network that was imported to Cytoscape contained 1025 nodes and 5685 edges. By using the CytoCluster plugin, the gene cluster was obtained to find the highly significant genes. A total of 157 clusters were obtained, out of which 19 had significant p values of less than 0.05, as depicted in Figure 3. Further, using the CytoNCA plugin, the top 15 genes were selected (Table 2), and again the selected genes were screened using the CytoHubba plugin to obtain the critical genes. After all three analyses, the top 14 genes were selected for the construction of a C-T network. The compounds were selected on the basis of the above ADMET screening.

Network Analysis
The disease gene interaction network obtained from STRING was imported to Cytoscape 3.8.2 for further analysis. The interaction network that was imported to Cytoscape contained 1025 nodes and 5685 edges. By using the CytoCluster plugin, the gene cluster was obtained to find the highly significant genes. A total of 157 clusters were obtained, out of which 19 had significant p values of less than 0.05, as depicted in Figure 3. Further, using the CytoNCA plugin, the top 15 genes were selected (Table 4), and again the selected genes were screened using the CytoHubba plugin to obtain the critical genes. After all three analyses, the top 14 genes were selected for the construction of a C-T network. The compounds were selected on the basis of the above ADMET screening.    After the discovery of core target and active compounds, the C-T network was constructed, as shown in Figure 4. The interaction network contains 35 nodes and 196 edges, which were analyzed using Cytoscape. The analysis shows that the phytocompounds of Ocimum tenuiflorum are connected to the target genes. Among these 14 genes, GAPDH and AKT1 were the most significant, according to the analysis.  After the discovery of core target and active compounds, the C-T network was constructed, as shown in Figure 4. The interaction network contains 35 nodes and 196 edges, which were analyzed using Cytoscape. The analysis shows that the phytocompounds of Ocimum tenuiflorum are connected to the target genes. Among these 14 genes, GAPDH and AKT1 were the most significant, according to the analysis.

Gene Ontology and Pathway Enrichment Analysis
The top 14 core genes were then analyzed using ClueGo and CluePedia for the GO functional annotation and KEGG pathway analysis. On the basis of the analysis result, a total of 148 pathways were obtained, which are divided into 17 functional groups that represents the pathways with GO terms as depicted in the pie chart ( Figure 5). Figure 6 shows the graphical representation of the cellular location of the potential marker using the interaction score from the GO term. The outcome of GO and pathway enrichment reveals the AGE-RAGE signaling pathway to be the most significant target pathway with

Gene Ontology and Pathway Enrichment Analysis
The top 14 core genes were then analyzed using ClueGo and CluePedia for the GO functional annotation and KEGG pathway analysis. On the basis of the analysis result, a total of 148 pathways were obtained, which are divided into 17 functional groups that represents the pathways with GO terms as depicted in the pie chart ( Figure 5). Figure 6 shows the graphical representation of the cellular location of the potential marker using the interaction score from the GO term. The outcome of GO and pathway enrichment reveals the AGE-RAGE signaling pathway to be the most significant target pathway with 50.53% probability in comparison with the other pathways analyzed ( Figure 5).

Molecular Docking Study
Following the network analysis, the most significant gene glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was considered for the docking study. The structural details of the selected protein molecule (PDB ID: 1U8F) have been depicted in Supplementary Figure S1A. The results from the molecular docking protocol reveal that the protocol was suitable for the screening of the phytocompounds. The RMSD between the bound and re-docked ligand was found to be 2.10 Å (Supplementary Figure S1B). We then selected the active compounds for the docking study on the basis of the C-T network and ADMET screening, and a total of 11 compounds were considered for the study. The virtual screening result is given in Table 3. The screening result showed that the GAPDH-bound isoeugenol complex showed a better binding affinity of −6.0 kcal/mol with the highest non-bonded interaction than other compounds. However, myrcene was found to have the lowest binding affinity (−3.6 kcal/mol) and was considered to be a negative control for further analyses. Both isoeugenol and myrcene were bound to the same binding site, where the co-crystallized NAD + was bound. Isoeugenol formed a total of eight non-bonding interactions, with three of them being hydrogen bonds. It bound to ARG 13 (3.20 Å), GLY 12 (3.41 Å), and SER 98 (2.03 Å) with hydrogen bonds. It also interacted with the protein using hydrophobic pi-alkyl bonds with ARG 13 (4.27 Å and 5.10 Å), ILE 14 (5.24 Å and 4.98 Å), and an alkyl bond with ALA 183 (3.90 Å). However, myrcene was found to have only two non-bonding interactions, both being alkyl bonds. It bound only with

Molecular Docking Study
Following the network analysis, the most significant gene glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was considered for the docking study. The structural details of the selected protein molecule (PDB ID: 1U8F) have been depicted in Supplementary Figure S1A. The results from the molecular docking protocol reveal that the protocol was suitable for the screening of the phytocompounds. The RMSD between the bound and re-docked ligand was found to be 2.10 Å (Supplementary Figure S1B). We then selected the active compounds for the docking study on the basis of the C-T network and ADMET screening, and a total of 11 compounds were considered for the study. The virtual screening result is given in Table 3. The screening result showed that the GAPDH-bound isoeugenol complex showed a better binding affinity of −6.0 kcal/mol with the highest non-bonded interaction than other compounds. However, myrcene was found to have the lowest binding affinity (−3.6 kcal/mol) and was considered to be a negative control for further analyses. Both isoeugenol and myrcene were bound to the same binding site, where the co-crystallized NAD + was bound. Isoeugenol formed a total of eight non-bonding interactions, with three of them being hydrogen bonds. It bound to ARG 13 (3.20 Å), GLY 12 (3.41 Å), and SER 98 (2.03 Å) with hydrogen bonds. It also interacted with the protein using hydrophobic pi-alkyl bonds with ARG 13 (4.27 Å and 5.10 Å), ILE 14 (5.24 Å and 4.98 Å), and an alkyl bond with ALA 183 (3.90 Å). However, myrcene was found

Molecular Docking Study
Following the network analysis, the most significant gene glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was considered for the docking study. The structural details of the selected protein molecule (PDB ID: 1U8F) have been depicted in Supplementary Figure S1A. The results from the molecular docking protocol reveal that the protocol was suitable for the screening of the phytocompounds. The RMSD between the bound and re-docked ligand was found to be 2.10 Å (Supplementary Figure S1B). We then selected the active compounds for the docking study on the basis of the C-T network and ADMET screening, and a total of 11 compounds were considered for the study. The virtual screening result is given in Table 5. The screening result showed that the GAPDH-bound isoeugenol complex showed a better binding affinity of −6.0 kcal/mol with the highest non-bonded interaction than other compounds. However, myrcene was found to have the lowest binding affinity (−3.6 kcal/mol) and was considered to be a negative control for further analyses. Both isoeugenol and myrcene were bound to the same binding site, where the co-crystallized NAD + was bound. Isoeugenol formed a total of eight non-bonding interactions, with three of them being hydrogen bonds. It bound to ARG 13 (3.20 Å), GLY 12 (3.41 Å), and SER 98 (2.03 Å) with hydrogen bonds. It also interacted with the protein using hydrophobic pi-alkyl bonds with ARG 13 (4.27 Å and 5.10 Å), ILE 14 (5.24 Å and 4.98 Å), and an alkyl bond with ALA 183 (3.90 Å). However, myrcene was found to have only two non-bonding interactions, both being alkyl bonds. It bound only with ILE 14 (5.14 Å and 4.51 Å) (Figure 7). A protein-ligand interaction fingerprint (PLIF) was conducted to interpret the docking results (Supplementary Figure S2).

Molecular Dynamics Simulation
A molecular dynamics simulation was used to validate the docking investigation and determine the degree of stability of the docked complex along with the target protein.
Graphs for RMSD, RMSF, Rg, SASA, and H-bond are represented in Figure 8 to analyze the trajectories.
The root-mean-square deviation (RMSD) graph represents the protein-ligand complex's stability throughout the course of a 100 ns simulation. By examining the plot, it can be said that the complex and target protein show a similar pattern of stability. The complex is within the range of 0.2-0.5 nm and the protein backbone was found within 0.2-0.35, and both were found to be stable after 80 ns. However, the protein-myrcene complex was found to have an RMSD value of 0.6 nm, and the plot was not concordant with that of the protein backbone atoms. In the RMSF analysis, both the isoeugenol complex and the protein backbone atoms were on par, with an almost similar fluctuation pattern. In the protein-myrcene complex, a higher number of fluctuations was observed in the loop regions (200 residues). Fluctuations were also found between 50 and 100 residues in 200-250 regions. Further, Rg and SASA plots were analyzed to show the structural compactness of the structure formed. The Rg plot analysis shows that the protein and the complex were found to be ranging between 2.0 and 2.1 nm, whereas the SASA value was found to be on par, with a similar pattern. The same trend of incompatible plot was predicted for protein-myrcene complex in Rg analysis, which was observed in the case of RMSD and RMSF. It was equilibrated at 2.2-2.3 nm. In the SASA plots, the value of protein-myrcene complex was found to be ranging between 130 and 135 nm 2 , which made the plot incompatible with that of the protein backbone atoms. Finally, the H-bond was assessed to determine the structural re-agreement, and it can be seen that the complex may have undergone structural alterations. In ligand hydrogen bond formation, isoeugenol formed more hydrogen bonds (six) in comparison with myrcene (four). Table 5. Binding affinity and non-bonding interactions of compounds with the target proteins.

Binding Free Energy Calculations
To determine the energy generated during the MD simulation, the free binding energy calculation was performed for both protein-ligand complexes. From the analysis, it can be seen that, aside from polar solvation energy and electrostatic energy, other energies contribute negatively to the complex development in the case of both the complexes. From the results, it is clear that the GAPDH-isoeugenol complex is more stable than the GAPDH-myrcene complex. Apart from the experimental values, the standard deviations also indicate the deviation of the GAPDH-myrcene complex from stabilization. Table 6 gives the details of the energy formed, along with their values.

Discussion
Over the past few decades, there have been a substantial number of studies conducted on the relation between type 2 diabetes (T2D) and AD (AD), as both T2D and AD are age-related conditions [37][38][39]. According to a study that has been recently conducted [40], it is known that there is about a 50-150% chances of an increase in cases of AD with people suffering from T2D compared with the normal population. T2D affects, directly or indirectly, various metabolic, inflammatory, signaling, and hyperinsulinemia functions, which can also be a related risk factor leading to AD. One of the most common features of T2D is impairments in insulin actions and signaling, which leads to hyperglycemia and hyperinsulinemia [40]. The most common factor linked to T2D and AD is the role of insulin resistance [3,40].
Ocimum tenuiflorum, which is commonly called holy basil or tulsi, belongs to the Lamiaceae family. Over the past few decades, many researchers have studied the pharmacological properties of tulsi, and it has been used in Ayurveda up until now. In accordance with [41][42][43], tulsi has exhibited a beneficial effect on both diabetes and AD. With this in consideration, we used phytocompounds from tulsi to study the effect of potential anti-diabetic and anti-Alzheimer's drugs.
To understand the relationship between T2D and AD, the alternative herbal medicine network pharmacology approach was used. The core target and the active compounds were identified to understand the disease network and their pathway from the network perspective. The PPI network for the common genes from both T2D and AD was constructed using the STRING database, which was later imported to Cytoscape for analysis. The main focus of this objective was to find the relationship through network construction and to find the most significant genes and their pathway using GO terms.
The constructed PPI contained a total of 19 significant clusters that had p values of less than 0.05 (p < 0.05). Further, using CytoNCA and CytoCluster, the most significant GAPDH gene was considered for a docking and dynamic study to understand the interaction of active compounds on GAPDH as target. In an in vivo study, GAPDH was reported to play a significant role in the development of diabetic retinopathy by elevating retinal mitochondrial superoxide levels. It is also reported to elevate the AGE-RAGE pathway, protein kinase C, and the hexosamine pathways, which are considered to be the hallmarks of T2D [44][45][46]. However, in the presence of αβ-42, the nitrosylation of GAPDH can increase tau acetylation, causing a disruption in the microtubule association process, and the amount of nitrosylated GAPDH is elevated in AD brains. These studies indicate that GAPDH is linked with both T2D and AD and could be used as a potential target. Thus, the inhibition of GAPDH would be an essential step to reduce the pathogenesis of diabetes-linked AD [47].
We screened the core 14 genes for the pathway analysis with GO term using ClueGo and Cluepedia plugins. The core 14 genes are correlated with multiple biological processes and pathways, which includes the AGE-RAGE signaling pathway, the interleukin-4 and interleukin-13 signaling pathways, the photodynamic therapy-induced NF-kB survival signaling pathway, the negative regulation of the lipid catabolic process, and pancreatic cancer, which are the top interacting pathways. Our study shows that the AGE-RAGE signaling pathway may play a key role in drug development for both T2D and AD.
In both hyperglycemic and calcification settings, AGE/RAGE signaling promotes both cellular and systemic responses to increase bone matrix proteins via the PKC, p38 MAPK, fetuin-A, TGF-, NFB, and ERK1/2 signaling pathways. Through the activation of Nox-1 and the decreased expression of SOD-1, AGE/RAGE signaling has been demonstrated to increase oxidative stress and promote diabetes-related vascular calcification. Increased oxidative stress caused by AGE/RAGE signaling in diabetes-related vascular calcification was also linked to the phenotypic transformation of VSMCs to osteoblast-like cells in AGEinduced calcification. Researchers discovered that pharmacological drugs and antioxidants reduced calcium deposition in the diabetic vascular calcification caused by AGEs [48]. However, GAPDH is linked with the activation of the AGE-RAGE pathway, which could be the initiating point of diabetes-linked AD pathogenesis [44,45].
The most significant gene, GAPDH, was selected as a target of an interaction study after C-T network construction. Before going for the molecular docking simulation, AD-MET profiling of the retrieved compounds was performed. Compounds that obey the desired ADMET properties were selected for further analysis. The Lipinski rule of five was considered as an important parameter for this analysis, as it deals with the druglikeness parameters. Although most of the compounds accepted the rule, other pharmacokinetic rules, such as oral bioavailability, the blood-brain barrier, drug half-life, and drug-induced liver injury, were considered for the screening of the compounds. This was conducted because natural products are often cited as an exception to Lipinski's rules. Phytochemicals are considered to be the product of nature's 'evolution', which enhances the functionality of the compounds in specific metabolic pathways [49]. When it comes to making physiologically active molecules with high molecular weights and a large number of rotatable bonds, nature has learned to maintain low hydrophobicity and intermolecular H-bond donating potential. Natural products are also more likely than pure synthetic compounds to resemble biosynthetic intermediates or endogenous metabolites, allowing them to benefit from active transport systems. Surprisingly, the natural product leads in the Lipinski universe all delivered an oral medication with a 50% success rate [50].
Only selected compounds from the ADMET screening were docked with the target protein. From the docking study, we see that isoeugenol had a better binding affinity of 6.0 kcal/mol with the highest non-bonded interaction and hydrogen bonds of seven and two. Where the hydrogen bond formed is within the active site of the structure. The binding of the isoeugenol to the hydrophobic site, which is supposed to be occupied by nicotinamide adenine dinucleotide phosphate (NADPH), suggests that isoeugenol inhibits the enzyme activity. Moreover, other known details of the hydrophobic site also mention the same [51]. Among the docked compounds, myrcene was predicted with the lowest (more positive) binding affinity (−3.6 kcal/mol). It had only two non-bonding interactions with zero hydrogen bonds. Therefore, it was further considered as a negative control for MD simulation. Although geranyl acetate was predicted with −5.9 kcal/mol of binding affinity, which is close to that of the isoeugenol's, the compound was not considered to be a computational hit. This was because of the comparatively fewer non-bonding interactions (six) and hydrogen bonds (one) than isoeugenol. In the PLIF analysis, binding residues interacting with isoeugenol and myrcene were used, for which residues interacting with NAD + were used as a reference (Supplementary Figure S2). However, ALA 183 was not considered for PLIF, since it does not interact with NAD + . A study conducted on GAPDH inhibition by 3-bromopyruvate and its derivatives showed that molecular docking can be performed for this protein [52]. However, the study failed to provide sufficient proof for the stability of the experimental compounds through MD simulation and binding free energy calculations.
Both isoeugenol and myrcene were found to obey Lipinski's rule of five. Both the compounds passed the oral bioavailability parameter. Isoeugenol was predicted to have negative drug-induced liver injury (DILI), yet myrcene was found to have moderate DILI. Therefore, isoeugenol could be a better drug candidate than myrcene in terms of cytotoxicity. Additionally, the CYP450 metabolism of isoeugenol generates a corresponding epoxide, which could act as a potential inhibitor of GAPDH. This is due to GAPDH's increased proclivity for reacting with epoxide as a result of its contact with the phosphate moiety, which causes the enzyme to change its conformation [53].
Further, to know the stability of the docked complex, an MD simulation was performed. From the MD result, we can see that the complex formed is stable, and a structural re-arrangement may have taken place with the docked structure. In the MD result analysis, isoeugenol-bound GAPDH was found to be more stable in comparison with myrcenebound GAPDH. The graphical analysis of the MD trajectories, including RMSD, RMSF, Rg, SASA, and ligand hydrogen bonds, depict the stability of isoeugenol inside the binding pocket over myrcene. The binding free binding energy was calculated, as it is known that the greater the negative value, the higher the stability. The calculated result had the free binding energy of −143.458 kj/mol; thus, we can conclude that the docked structure is stable. However, the GAPDH-isoeugenol complex is chiefly formed by the van der Waal's energy, which is −159.277 kj/mol. However, the binding free energy values of the GAPDHmyrcene complex were found to be inferior in comparison with the GAPDH-isoeugenol complex. Although both the complexes were formed majorly using van der Waal's energy, the GAPDH-isoeugenol complex was predicted to have more stability. Results obtained from both molecular dynamics simulation and binding free energy calculations support the molecular docking simulation results. The present investigation is the first ever to be conducted on GAPDH inhibition using plant-based compounds through network pharmacology and computational studies. Until now, there have been no studies on GAPDH inhibition through an in-silico approach using phytochemicals. With these findings, we report the pharmacological significance of isoeugenol from Ocimum tenuiflorum against diabetes-linked AD targeting GAPDH.

Conclusions
A phytochemical intervention into T2D has already been proven to ameliorate the lifestyle-based disorder. However, diabetes-linked AD and other diabetic complications have not been focused on so far. Therefore, in this study, we report the virtual screening of isoeugenol from Ocimum tenuiflorum as an inhibitor of GAPDH, a common target for both T2D and AD. The study systematically identifies GAPDH as a promising target through network pharmacology tools. However, the phytocompounds have also been screened through a druglikeness and pharmacokinetic approach. Moreover, the compound target network also identifies isoeugenol as the best inhibitor of GAPDH. This hypothesis is supported by the computational tools, including a molecular docking simulation, a molecular dynamics simulation, and binding free energy calculations, where the binding interactions and stability of the isoeugenol with GAPDH have been proven to be satisfactory. With these results, we can conclude that isoeugenol could be the lead potential inhibitor of GAPDH. In future, isoeugenol could be used for in vitro and in vivo studies against GAPDH, targeting T2D-linked AD.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules27082398/s1, Figure S1: (A) Structural details of protein GAPDH (PDB ID: 1U8F) bound to its co-crystallized NAD + cofactor. N-terminal of the protein is highlighted in blue, whereas the red terminal is highlighted in red. (B) Result from the validation of the docking protocol: Crystallized and bound NAD + (green) complexed with re-docked NAD + (red) with an RMSD of 2.10 Å; Figure