Triterpenes from African Antidiabetic Herbs as Inhibitors of DPP-4 and PTP1B Targets: Molecular Modeling Investigation

1. Introduction

Diabetes mellitus (DM) ranks among the most challenging public health burdens in the world as it affects about 463 million [1] and may increase to 552 million by the year 2030 and 700 million by 2045 as projected by the International Diabetes Federation (IDF) [1,2,3,4]. It is generally characterized by hyperglycemia and impaired glucose homeostasis resulting from defects in insulin secretion, insulin action, or both [5,6]. Type I diabetes mellitus (T1DM) is implicated in about 10% of the diabetic population while Type II diabetes mellitus (T2DM) contributes approximately 90% of diabetes-related cases [7]. While T1DM is often managed by insulin treatments and lifestyle modifications [8], an ideal treatment for T2DM is still elusive as the pathogenesis is highly complicated. Thus, many patients are still unable to meet glycemic goals [9,10]. Therefore, the persistent unaddressed needs of most patients remain an important factor driving the quest for better prevention and treatment options of diabetes mellitus.

Several promising therapeutic agents, which target various players in diabetes are promising oral agents for T2DM patients [11,12]. These include the starch blockers, insulin secretagogues, insulin mimickers and insulin sensitizers [12,13,14,15]. Acarbose, sulfonylureas, meglitinides, thiazolidinediones, biguanides aimed at blunting post-prandial glucose rise, suppressing hepatic glucose output, stimulating insulin release and increasing peripheral glucose utilization respectively are already approved for clinical use [1,16], but not without various adverse effects [1,17]. Beyond alpha-glucosidase, alpha-amylase and other starch blockers widely exploited for antidiabetic drug development [18], dipeptidyl peptidase-4 (DPP-4) enzyme has been recommended as target for antidiabetic agents as it also acts in carbohydrate metabolism by delaying gastric emptying, increasing insulin secretion, and reducing glucagon secretion [12,19]. The DPP-4 inhibitors stimulate the breakdown of endogenous incretin hormones and promote pancreatic glucose-dependent insulin secretion. The enzyme is a serine peptidase which act by rapidly degrading the incretin hormones, such as glucagon-like peptide 1 (GLP 1) that play an important role in blood glucose regulation, causing short life span of the enzyme. Therapeutic strategy that involves inhibiting DPP-4 would help to maintain the endogenous level of GLP 1 and increased β cell mass. Thus DPP-4 inhibitors are capable of enhancing glucose-dependent insulin secretion, slowing down gastric emptying, and reducing postprandial glucagon and of food intake [20,21,22,23]. Metformin is the first-choice medication which is most frequently recommended for treating T2D. This drug act by improving glucose metabolism through both AMP-activated kinase (AMPK) activation and enhanced GLP-1 release [24]. Thus, it improves insulin sensitivity of relevant tissues [25,26]. However, metformin like other antidiabetic drugs is often accompanied by moderate to severe side effects [27]. Another emerging target for enhancing insulin sensitivity in various cells is the Protein tyrosine phosphatase 1B (PTP1B). By improving the sensitivity of the insulin receptor, PTP1B inhibitors have the potential to cure insulin resistance-associated ailments [28,29]. The PTP1B acts by negatively regulating insulin receptor signaling pathway in cells [30]. Agents that inhibit this enzyme can activate phosphorylation of several insulin receptor kinase substrates which in turn stimulate the phosphoinositide-3-kinase as well as subsequent protein kinase B expression. This further promotes translocation of glucose transporter into the cell membrane from intracellular vesicles; and thereby promote uptake of glucose molecules by the cells for energy generation. Due to the major roles played by insulin disturbance in diabetes, insulin signaling is highly recommended for targeting diabetes and explored for developing new antidiabetic agents [10].

African flora possesses a very rich biodiversity of herbs, spices and medicinal which is widely exploited in traditional, alternative and complementary medicine for treating several diseases which include diabetes mellitus (DM) [13,31]. Most traditional healers in Africa who detain the ancestral heritage of these medicinal plants are illiterate with their ethnopharmacological knowledge transmitted verbally from generation to generation which is at risk of disappearing. Thus, the World Health Organization (WHO) has recommended scientific studies to document the folk knowledge and validate the acclaimed therapeutic potential of these plants from the perspective of developing improved medications (WHO, 2013). In this direction, the literature is replete with scientific reports on the anti-diabetic activity of African herbs, spices and medicinal as previously revealed by excellent review articles [31,32,33,34]. Most reported studies mainly focused on plant extracts and/or fractions, with only few having got scientific validation. Several ethnopharmacologically important antidiabetic herbs in Africa have shown potent activity against DPP-4 [31]. The extracts of Antidesma madagascariense Lam, whose leaf and stem bark decoctions are used in folk medicine in Madagascar for the management of DM, showed inhibitory activity against DPP-4 enzyme with an IC50 value of 79.2 mg/ml [35,36]. This inhibition indicated the ability of the extract to improve insulin sensitivity. In the same vein, several bioactive constituents of African indigenous herbs have been screened through several in vitro and in vivo studies for activity against PTP 1B. Ethyl acetate fraction and isolated pimarane diterpenes from Icacina oliviformis, a popular food spice native to the regions of West and Central Africa, have been reported to inhibit PTP1B [31,37,38]. Mohammed and Tajuddeen [31] revealed that the widely reported African antidiabetic plants contain several classes of compounds which include terpenes, alkaloids, flavonoids, and carotenoids. Owing to their wealth of chemical structures, such antidiabetic plants are considered as an inventory of bioactive compounds which might be useful as a basis for drug discovery and development in combating diabetes. Terpene structures, which belong to an important class of natural products have been extensively reported for antidiabetic activity and some of them are already under various stages of pre-clinical and clinical evaluation geared towards developing antidiabetic agents [10,14]. The antidiabetic potential of plant-derived terpenes could be attributed to varying mechanisms of action which include: increasing insulin sensitivity [14,39], inhibiting pancreatic amylase and glucosidase enzymes [40,41], reducing oxidative stress [42], and inhibiting the development of diabetic complications [43,44]. Several studies have revealed that, terpene-rich bioactive extracts of several African medicinal plants have shown antidiabetic activity by targeting key proteins associated with type II diabetes in vitro and in vivo as reviewed by Panigrahy, Bhatt [14]. Several terpene-rich plant extracts and terpene isolates have been widely reported to target DPP-4 [45], PTP1B [46,47] and other targets [14] in diabetes in vitro and in vivo.

Exploring novel phytochemical–target interactions in silico could help to provide useful insights into the mechanism underpinning the biological activities of medicinal plants and their constituent phytochemicals. It also plays an important role in compound screening and drug discovery. Drug-target interactions (DTI) characterize the binding of drug compounds to specific targets. Binding affinity of drugs refers to the strength of DTIs; and thereby greatly influences the drug's efficacy and therapeutic potential. Molecular modeling tools, which provide rapid, economically feasible and environmentally sound techniques for probing DTIs are valuable for screening un-investigated phytochemicals and predicting their interactions with molecular targets using integrated chemical informatics, bioinformatics and systems biology approaches. Employing these tools, the three-dimensional visualization of the electronic and steric molecular properties of the interaction between bioactive compounds and drug targets could be computed. The ultimate goal of such computational modeling in natural product research is to predict, in advance of any targeted laboratory testing, novel biologically active compounds. The modeling methodologies often employ the classical molecular mechanics (MM) and quantum mechanics (QM) computations [48]. Integrated molecular mechanics computations and quantum mechanics computations could provide an accurate and efficient theoretical description of a biomolecular system. Earlier studies have revealed the potential of in silico methods for predicting DPP-4 inhibitors [49,50,51] and PTP1B inhibitors [52,53] from natural sources. Terpenes, a class of compounds with compact molecular structure and less symmetry are known to bind multifarious drug targets to exhibit therapeutic potential in human diseases. Therefore, investigating terpene-target interactions hold great promise in antidiabetic drug discovery [14]. In spite of the reported antidiabetic activity of various African herbs and the therapeutic roles of terpene structure, many terpenes from African herbs remain uninvestigated. The aim of this study was to explore the interactions of terpene structures from African antidiabetic plants with two emerging drug targets in diabetes.

2. Results

2.1. Binding affinity and interactions of terpenes with DPP-4 enzyme

The initial screening based on the Larmackian Genetic Algorithm docking protocol revealed 37 terpene structures with docking scores (-7.7 to -10.0 Kcal/mol) lower or comparable to that of the co-crystalized alogliptin (-7.7 Kcal/mol) (Table S1). Ranking based on the docking scores and screening for favourable interactions with the target enzyme revealed the top 7 compounds (Table 1), most of which belong to the triterpenes sub-class of terpenes.

The topmost 7 terpene-DPP-4 complexes were subjected to post-docking MM-GBSA computations. While MM-GBSA (ΔG_bind) of the reference alogliptin was -37.02Kcal/mol, the MM-GBSA (ΔG_bind) of the top-docked terpene structure ranged from -47.80 kcal/mol (Cucurbitacin B) to -19.72 kcal/mol (1-Deacetylkhivorin) as shown in Table 1. It is worthy of note that, Cucurbitacin B with the lowest binding energy with DPP-4 also exhibited the lowest MM-GBSA free binding energy to the enzyme. While van der Waals, coulombic and non-polar solvation energies made the most notable free energy contributions to Cucurbitacin B binding to DPP-4, van der Waals and non-polar solvation made important contributions to binding of several other terpenes (Table S2). The terpene structures demonstrated high levels of binding affinity, suggesting exceptional potential for sustained interactions with the DPP-4 enzyme. Interaction analysis of the selected terpenes with the binding pocket of DPP-4 revealed that, cucurbitacin B, 6-Oxoisoiguesterin, and 20-Epi-isoiguesterinol had binding poses with DPP-4 similar to that of the co-crystalised alogliptin. These compounds also featured strong and favourable interactions with the critical amino acid residues in the active site of the enzyme in a similar manner to the native compound as shown in Table 2. Table 2 shows that, the hit triterpenes and alogliptin had strong interactions with various residues involved in the catalytic triad, oxyanion cavity, hydrophobic S1 pocket and charged S2 pocket in the binding pocket of DPP-4. The 3-D representations of the interactions is depicted in Figure 1.

2.2. Binding affinity and interactions of terpenes with PTP1B enzyme

Molecular docking, MM-GBSA post docking analysis and molecular interactions were also employed to screen the terpene structures for potential inhibitors of PTP1B enzyme. Although all the docked compounds had docking scores and MM-GBSA binding energy higher than those of the reference isothiazolidinone, the topmost compounds had favourable interactions with the binding pocket the protein. The Molecular docking scores and prime MM-GBSA scores of the top terpenes with PTP1B are shown in Table 3. Favourable interactions were also observed upon visualization of the top terpene-PTP1B complexes as shown in Table 4 and Figure 2.

2.3. Ensemble docking analysis of selected triterpenes against DPP-4

Based on the results from molecular docking and post docking MM-GBSA which showed higher affinity of top-docked terpenes with DPP-4, the top docked terpenes were further subjected to ensemble docking analysis in order to provide more flexibility in the docking experiment. The ensemble docking analysis revealed that, the average docking scores of the hit triterpenes were lower than that of the reference alogliptin.

Figure 3. (A) Mean scores of docking of top terpenes with structure ensemble of DPP4 (n=4). (B) Interactions of top terpenes with a representative conformer of DPP4.

Figure 3. (A) Mean scores of docking of top terpenes with structure ensemble of DPP4 (n=4). (B) Interactions of top terpenes with a representative conformer of DPP4.

2.4. Molecular dynamics of terpene-DPP-4 interactions

The dynamics simulation was employed to model the interactions of the top-docking terpenes with DPP-4 in a dynamic environment. The thermodynamic parameters that indicate stability of the DPP-4 enzyme back-bone C-α atoms in complex with T1, T4 and T2 were calculated and plotted as shown in Figure 4. The RMSD, a measure of the difference between two sets of coordinates, was computed and depicted in Figure 4a.

The RMSD plot (Figure 4a) revealed that, the DPP-4 backbone, T1-DPP-4, T4-DPP-4, T2-DPP-4 and S1-DPP-4 molecular systems showed a slight rising RMSD trend in initial 0.6 ns, 0.8 ns, 0.1 ns, 0.6 ns and 1.1 ns and then stabilized at mean values of 2.84 Å, 2.85 Å, 3.02 Å, 3.17 Å and 2.50 Å respectively. These values indicated that the terpene–DPP-4 complexes reached a stable conformational trend over 100 ns of the simulation with few conformational transitions as the RMSD values were maintained within 2 Å. Except for the major fluctuations observed with the trajectory of DPP-4-cucurbitacin B complex, all the complexes exhibited relatively stable atomic trajectories displaying minor and consistent fluctuations. Figure 4b shows that the average RoG values for DPP-4 backbone, T1-DPP-4, T4-DPP-4, T2-DPP-4 and S1-DPP-4 biomolecular systems are 26.00 Å, 26.36 Å, 26.25 Å, 26.02 Å and 26.57 Å respectively. The results also revealed that the RoG values were maintained within 2 Å. The result shown in Figure 4c revealed that, the SASA values of DPP-4 backbone, T1-DPP-4, T4-DPP-4, T2-DPP-4 and S1-DPP-4 biomolecular systems are 29052 Å2, 29797 Å, 29669 Å, 29115 Å2 and 30015 Å2 respectively. To provide more insight into the stability of the complexes, the number of H bonds found in the ligand-protein complexes was computed as depicted in Figure 4d. The results revealed a steady trend for the five biomolecular systems which span around 173, 180, 174, 172 and 184 bonds for the DPP-4 backbone, T1-DPP-4, T4-DPP-4, T2-DPP-4 and S1-DPP-4 biomolecular systems respectively. The plot demonstrates large variations in the amino acid residues between 220 and 260 as well as those around 680, with the RMSF reaching notable peaks. Other residues, in addition to these also feature some notable peaks.

2.5. Free energy simulation of DPP-4-terpene complexes

To further validate the docking scores of the ligands to the protein, the trajectories obtained from the MD simulations of the terpene-DPP-4 complexes were subjected to free energy simulation using MM-GBSA. Table 5 shows a number of solvation and molecular mechanics parameters and their energy contributions.

The results show a number of solvation and molecular mechanics parameters and their energy contributions. Among the terpenoids, T1 (ΔGtotal = -11.85 kJ/mol) exhibited the strongest binding affinity as indicated by the lowest free binding energy as compared to other terpenes but with higher free binding energy than the reference agloliptin. The results also revealed that, the main components of the free energy of the top complexes are van der waals energy and electrostatic energy. Amino acids occurring around the binding pocket showed significant contributions to the overall binding affinity to the terpen-DPP-4 complexes as indicated by the MM-PBSA free energy decomposition analysis (Figure 5).

2.6. Frontier molecular orbitals of top triterpenes

To further investigate the behavior of the selected terpene structures as DPP-4 inhibitors, density functional theory computation was performed. Table 6 shows that, the hit compounds had close band gaps which were quite low.

The HOMO and LUMO orbital of the hit compounds are given in Figure 6. The regions of the compounds with red and blue color denote the positive and negative signal/phase of the molecular orbital respectively.

The orbital stands for the likelihood of discovering an electron. The colours red and blue signify the nodal character of the orbital. The HOMO orbital appears to be centralized on the phenyl rings of the hit compounds while the LUMO orbital were predominantly of the side chains of the aromatic ring structures.

2.7. Predicted physicochemical and ADMET properties of hit triterpenes

The drug-likeness characteristics of the top triterpenes were assessed to ascertain the suitability of the compounds as drug candidates. The drug-likeness characteristics of the compunds revealed desirable drug-like properties (Table 7). Using the Lipinski rule of five assessments, it was shown that all the hit compounds are in accordance to Lipinski rule.

3. Discussion

In this study, we used integrated in-silico methods to explore the molecular interactions of terpene structures previously reported from African antidiabetic plants with two emerging drug targets in diabetes. Such methods, which are valuable for sample space minimization and resource maximization, can provide rapid, economically feasible and environmentally sound techniques for probing the biological activities of uninvestigated phytochemicals by identifying their molecular targets; and thereby provide insights into the mechanism underpinning their activity. Besides, they enable more targeted experimental and preclinical evaluations. Herein, structure-based virtual screening was used to screen 107 terpene structures against the 3D-structures of DPP-4 and PTP1B enzymes. Molecular dynamics (MD) simulation, Molecular mechanics with generalized Born and surface area solvation (MM/GBSA) free energy calculation and density functional theory (DFT) were used to clarify the interactions of top the terpenes with the target enzymes. The Structure-based virtual screening helped to identify top terpenes with higher binding affinity with the target enzymes than the reference drugs. Post-docking MMGBSA analysis of the top terpenes revealed good binding affinity scores. It is worthy of note that, Cucurbitacin B (T1) with the strongest binding affinity with DPP-4 as indicated by the lowest docking score also exhibited the lowest post-docking MM-GBSA free binding energy score to the enzyme. Molecular docking, MM-GBSA post docking analysis and molecular interactions also revealed top terpenes as potential inhibitors of PTP1B enzyme.

The interactions of the top terpenes were analyzed alongside the reference drugs. The results showed that, the interaction of the co-crystalized compound (alogliptin) with the binding site of DPP-4 involve several types of bonds including Conventional Hydrogen Bonds, Carbon Hydrogen Bonds, Pi-Cation Electrostatic interactions, Pi-Donor Hydrogen Bond, Pi-Sigma Hydrophobic interactions, Pi-Pi Stacked Hydrophobic interactions and Pi-Alkyl Hydrophobic interactions. Key interactions of the co-crystalized compound involve Carbon Hydrogen bond, Pi-Cation electrostatic interactions and Pi-Alkyl hydrophobic interaction with HIS740 as well as Pi-Donor hydrogen bond with SER630. The DPP-4 enzyme, which is catalytically active as dimers in both the soluble and cell-surface form, has 766 amino acids, containing a C-terminal hydrolase domain as well as the N-terminal propeller domain and a C-terminal hydrolase domain. In addition to the catalytic triad comprising Ser630, Asp708 and His740, DPP-4 contains Tyr47 and Ser631 associated with the oxyanion hole; Tyr631, Val656, Trp659, Tyr662, Tyr666, and Val711 functioning in the hydrophobic S1 pocket; and Arg125, Glu205, Glu206, Phe357, Ser209 and Arg358 in the and a charged S2 pocket [54,55]. Interaction analysis of the selected terpene structures with the binding pocket of DPP-4 revealed that, cucurbitacin B, 6-Oxoisoiguesterin (T4), and 20-Epi-isoiguesterinol (T2) had binding poses with DPP-4 similar to that of the co-crystalised alogliptin. Cucurbitacin B, a tetracyclic triterpene containing the cucurbitane nucleus skeleton, namely, 19-(10-9beta)-abeo-10alpha-lanost-5-ene with a several oxygen substitutions in different positions. It was observed that, the oxygen substitutions in the structure played key roles in the interaction of the compound with active site and hydrophobic S1 pocket in the binding site region of DPP-4. There was strong interaction with the catalytic triad residues as important carbonyl oxygen in the compound structure received a strong conventional hydrogen bond from the side chains of SER630 and another from HIS740 side chain. Also, the ether oxygen received a strong conventional hydrogen bond from the side chain of TYR631 in the hydrophobic S1 pocket. In this pocket, the side chain of TYR666 conducted a Pi-Sigma Hydrophobic interaction with a methyl group and Pi-Alky hydrophobic interaction with another methyl group in the compound. The side chain of TYR662 also conducted a Pi-Alkyl hydrophobic interaction with a methyl group in the compound. Cucurbitacin B, was derived from Cogniauxia podolaena [56], a flowering perennial plant that belongs to the family Cucurbitaceae. It grows primarily in the wet tropical biome. The genus Cogniauxia is native to Angola, Cabinda, Cameroon, Central African Republic, Congo, Equatorial Guinea, Gabon, Guinea, Zaïre. The species in this genus are known to possess antidiabetic activity [57]. Four bisnorterpenes derived from Ekebergia capensis [58] viz: 20-Epi-isoiguesterinol, isoiguesterin, 6-oxoisoiguesterin and isoiguesterinol were also observed to show stronger binding affinity than alogliptin as indicated by lower binding energies. The Ekebergia capensis is a flowering plant which belong the family Meliaceae is fondly called Cape ash and its range extends from the South African Eastern Cape to Sudan and Ethiopia. Its antidiabetic potential has been exploited in Eastern Cape KwaZulu-Natal, Limpopo and Mpumalanga [59]. The herbal formulations comprising infusion prepared from the leaves of this plant and taken orally has shown antidiabetic activity streptozotocin-induced diabetic rats [59]. The 6-Oxoisoiguesterin (T4), one of the compounds from this plant interacted with the active site of DPP-4 through a carbon Hydrogen Bond and Pi-Alkyl Hydrophobic interaction with HIS740. As shown in Table 2 and Figure 1. The 20-Epi-isoiguesterinol (T3), another bisnortriterpene, interacted with the catalytic trial through a Pi-Alkyl Hydrophobic interaction with HIS740 as depicted in Table 2 and Figure 1. As extensively reviewed by Hamid, Alqahtani [60], tetracyclic triterpenoids, a class of triterpenoids found in a variety of medicinal plants are notable for treating diabetes and its complications. In a previous study, botulin, a triterpene derived from samples of the root bark of Euclea undulate widely exploited traditionally to treat diabetes and other ailments in Southern Africa, has been reported to show a reduction in level of blood sugar as compared to reference drug glibenclamide [61,62]. This study suggests the molecular mechanism underpinning the activity of such drugs. Also, earlier studies on the structure-activity relationship of terpenes have revealed the few roles of the molecular structure of terpenes and terpenoids in manifesting their antidiabetic activity. For instance, Chen, Lu [63] revealed that, sesquiterpenes are highly effective as antidiabetic agents which may result from their structurally compact configuration, low ramification, and low symmetry. Another author revealed that, the triterpenes with hydroxyl and carboxyl functional groups exhibit high antidiabetic activity which is mediated through various mechanisms [64]. In the review by Panigrahy, Bhatt [14], terpene structures containing the most hydrogen bonds and hydrophobic interactions exhibit enhanced anti-diabetic activity. Other compounds with high affinity for DPP-4 are also reported in this study (Table 1). These include several acyclic triterpenes derived from Glossocalyx brevipes (Siparunaceae) [65]. The native range of Glossocalyx is Southern Nigeria to Western Central Tropical Africa. Glossocalyx brevipes has been documented in the treatment of diabetes and its interconnected diseases [66]. Two important limonoids (1-Deacetylkhivorin and 7-Deacetylkhivorin) derived from Khaya grandifoliola (Meliaceae) [67] also showed low binding energy. Khaya grandifoliola, which is also known as African mahogany belong to the family Meliaceae and found abundantly in the Republic of Benin, Nigeria, the Democratic Republic of the Congo, Ivory Coast, Ghana, Guinea, Sudan, Togo, and Uganda. This plant has been utilized in polyherbal formulations for the treatment of T2DM [68]. Structure-based virtual screening of selected terpenes with PTP1b also revealed strong binding interactions. Although all the docked compounds had docking scores and MM-GBSA binding energy higher than those of the reference isothiazolidinone, the topmost compounds had favourable interactions with the binding pocket the protein. The Beilshmiedic acid derivatives viz: Tsangibeilin B and Cryptobeilic acid C , which were isolated from Beilschmiedia cryptocaryoides (Lauraceae) [69] had the lowest binding score. Several bisnorterpenes (isoiguesterin, 6-Oxoisoiguesterin, isoiguesterinol and 20-Epi-isoiguesterinol) derived from Ekebergia capensis [58] also had low binding scores. In addition, galanolactone, a labdane diterpenoid derived from Aframomum arundinaceum (Zingiberaceae) [70] exhibited low binding energy. Van der Waals, non-polar solvation, and Coulombs energy terms made most significant contributions to the binding affinity of the ligands.

Based on our results from molecular docking and post docking MM-GBSA which showed higher affinity of top-docked terpenes with DPP-4, the top docked terpenes were further subjected to ensemble docking analysis. Ensemble docking which involves the use of multiple molecular docking simulations was employed to improve the accuracy and reliability of predicting the top terpenes binding mode and affinity with DPP-4 enzyme. To achieve this, a DPP-4 ensemble was constructed by subjecting the experimental structure to a 100 ns MD simulation which generated three clusters. A representative structure was selected for each cluster and docked against the hit triterpenes (T1, T4, T2) along with the reference compound (alogliptin) using AutoDock Vina software [71,72].The result of the ensemble docking are shown in Figure 3. The study showed that, docking of top triterpenes with the active site of DPP-4 conformational representatives was consistent with the docking result obtained from docking with the single DPP-4 experimental structure. It is also worthy of note that, interaction analysis of the docked conformation with representative clusters revealed that the amino acid interaction were preserved and reproducible (Figure 5).

The top terpene-DPP-4 complexes were subjected to MD simulation alongside the reference alogliptin and unliganded enzyme to capture a more comprehensive picture of biomolecular systems considering the inherent flexibility and dynamics of proteins. The MD simulation is often integrated with docking simulation and other molecular modeling techniques to further validate virtual screening through static docking and assess the stability of drug-receptor complexes. Computation of various thermodynamic parameters from the MD simulation trajectories revealed the structural stability and conformation flexibility of the ligand-DPP-4 complexes alongside the unliganded protein as depicted in Figure 4. The RMSD values which indicate the overall structural stability of the DPP-4 enzyme back-bone C-α atoms in complex with alogliptin, T1, T4 and T2 are shown in Figure 4a. Except for the major fluctuations observed with the trajectory of DPP-4-cucurbitacin B complex, all the complexes exhibited relatively stable atomic trajectories displaying minor and consistent fluctuations. During MD simulation, low levels of RMSD which are characterized by such consistent fluctuations indicate equilibration and stability of the system. However, higher fluctuations indicate low stability. However, highly deviated structures after equilibration such as the observed RMSD of DPP-4-cucurbitacin B complex may imply major conformational transitions by DPP-4 to attain stable conformation of the enzyme with cucurbitacin B. In the context of molecular modeling and structural biology, RMSD is often used to compare the structures of two molecules, such as a crystal structure and a molecular dynamics simulation trajectory. To calculate RMSD, one first aligns the two sets of coordinates, typically by superimposing them using a least-squares fitting algorithm. Then, the distance between each corresponding pair of atoms in the two sets of coordinates is calculated, and the average of the squared distances is taken. Finally, the square root of the average squared distance is computed to obtain the RMSD value. The RMSD is a widely used metric in molecular modeling and structural biology, as it provides a quantitative measure of the similarity between two structures. A low RMSD value indicates a high degree of similarity, while a high RMSD value indicates a greater difference between the two structures. The radius of gyration calculation helps to assess how spread out the mass of the complexes is from its rotation axis. It is often employed in MD simulation to estimate the compactness of the bound structures. The RoG values of the apoprotein and the protein-ligand complexes were maintained within 2 Å indicating stability of the molecular system. Also the SASA, a measure of the total surface area of a molecule that is available for interaction with other molecules, such as enzymes or receptors indicate stability of the complexes (Figure 4c). It is calculated by measuring the area of the molecule that is exposed to solvent, such as water. The SASA is important because it can help predict the binding affinity of a molecule to its target. Generally, molecules with a larger SASA are more likely to bind strongly to their targets because they have more surface area available for interaction. Since the RoG values indicate compactness of the complex structures and those of SASA reflect the extent of the solvent accessible surface tendency of the protein structure. Taken together, the terpene-DPP-4 complexes with reduced RoG and SASA values reflect a low tendency of unfolding deviation from their initial structure. Hydrogen bond number existing in the molecular systems was calculated in order to provide more insight into the stability of the protein-ligand complexes. The study revealed that, DPP-4-ligand complexes exhibit a high tendency to maintain compact and well-folded biomolecular structure characterized by stable intermolecular bonds. The high number of H-bonds observed in the ligand-protein complexes (Figure 4d) further indicate well bounded ligand and compact structures. In this study, the values of RMSF for DPP-4 enzyme residues towards the triterpenes were calculated and shown in Figure 4e. The plot demonstrates large variations in the amino acid residues between 220 and 260 as well as those around 680, with the RMSF reaching notable peaks. Other residues, in addition to these also feature some notable peaks. This suggests that they have a strong affinity for interacting with terpenes and that the structures may be able to bind strongly to the protein at these points. Computation of the RMSF from the trajectory files helps to evaluate the crucial contributions that a protein's amino acid residues make towards achieving stable conformations for protein-ligand complexes. The RMSF analysis often reveals how a portion of the protein deviates from its mean structure upon interactions with a ligand. The RMSF plot displays the degree of fluctuation seen for each protein residue, providing insight into the degree of adaptability they have in a changing environment. As a result, motifs or individual amino acid residues that exhibit higher RMSF and flexibility have a stronger propensity of interactions with the ligand.

Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) is often employed in molecular modeling investigations to assess the binding free energy of a compounds with a receptor. It is often considered as a post-processing technique for analyzing trajectories obtained from an MD simulation run. Given the limitations of molecular docking simulation, which is mostly reliant on molecular mechanics calculations, the trajectories obtained from the MD simulations of the terpene-DPP-4 complexes were subjected to free energy simulation using MM-GBSA, which combines conventional MM with implicit solvation models. The main components of the free energy of the top complexes are van der waals energy and electrostatic energy. Amino acids occurring around the binding pocket showed significant contributions to the overall binding affinity to the terpen-DPP-4 complexes as indicated by the MM-PBSA free energy decomposition analysis (Figure 5).

Density functional theory is a quantum mechanical approach that calculates the electronic density of a system, which can then be used to calculate other properties such as total energy, electronic charge distribution, and molecular geometry. The essential chemical characteristics of the terpene structures were analysed using frontier molecular orbital analysis of the completely optimised structures as shown in Table 6. According to FMO theory, one of the key factors that influence the biological activities of drug-like compounds is the energy level of their HOMO and LUMO orbitals [73]. The HOMO energy indicates the compound's ability to donate electrons, whereas the LUMO energy quantifies the compound's capacity to accept electrons. All the triterpene structures displayed higher HOMO energy than the reference alogliptin, which is consistent with the molecular docking, ensemble docking and post-docking MM-GBSA. In simple term, a compound with small energy gap is more polarizable and often associated with high reactivity [74]. These suggest that the hit compounds have high tendency of participating in molecular reactions. In a previous study, it was observed that clinically approved therapeutic showed good correlation with the experimental data having electron affinity value ranging between −1.5 and 2.0 eV [75]. The fact that all our hit compounds EA value lies within this range suggest a good druglikeness property. In addition, it was observed that the highest electrophilicity index was observed in cucurbitacin B, oxo-isoguesterin and Epi-isoguestrin. The electrophilicity index depicts information about electron system structure, stability, bonding, reactivity, and dynamics at ground and excited states [76]. Ganesan, Raja [76] further explained that, compounds with a high electrophilicity index are more likely to interact with biomolecules, suggesting a preferential reactivity of these compounds. Overall, all the hit compounds show promising reactivity indices, which emphasizes significant potential as lead candidates.

The drug-likeness characteristics of the top triterpenes were assessed to ascertain the suitability of the compounds as drug candidates. These compounds possess desirable properties (Table 7), which are important to the ultimate nutraceuticals and drug development. Using the Lipinski rule of five assessments, it was shown that all the hit compounds are in accordance to Lipinski rule. Based on Lipinski rule, orally active drugs have no more than one violation of the set criteria viz: less than 5 hydrogen bond donors, less than 10 hydrogen bond acceptors, molecular mass below 500 daltons and octanol-water partition coefficient less than 5 [77,78]. The fact that the compounds did not exceed one of the Lipinski violations suggest they can be readily absorbed and in therapeutic doses. The lead-likeness property was further corroborated by the Veber, Egan and Muegge rule. Further emphasizing the druglikeness property of the drugs is the result from pain (Pan-assay interference compounds) alert prediction. While the previous rules focused on the oral bioavailability of the drug, pain alert helps to identify the tendency of the compounds to participate in non-specific binding to multiple targets. The PAINS is one of the major results for false positive in drug development. Based on these prediction, it was observed that most of the top compounds do not contain scaffolds capable of non-specific binding as such the high binding affinity obtained is less likely due to false positive. Furthermore, most of hits were not predicted as p-glycoprotein substrate. The Pgp is responsible for binding foreign particles and facilitating their rapid elimination from the body. Since most of the compounds are not pgp substrate, the implication is the compounds are likely to have long circulating time thereby achieving maximum therapeucity before clearance from the body.

4. Materials and Methods

4.1. Preparation of protein structures

The 3-dimensional structure of human DPP-4 (PDB ID: 3G0B) which was co-crystallized with T22800 (alogliptin) and that of the human PTP1B (PDB ID: 2CM7), which was co-crystalized with isothiazolidinone, were retrieved from the Protein Data Bank (https://www.rcsb.org). From these protein structures, the co-crystalized compound and water molecules were deleted. Each hydrogen atom was integrated into protein the structures using Autodock Vina 4.2 program.

4.2. Ligand preparation

An in-house library of 107 terpenes previously reported from African medicinal plants was created after a thorough search of the literature. The structures of the terpenes from the in-house library and the protein native ligands were obtained in the structural data format (SDF) from the PubChem database at www.pubchem.ncbi.nlm.nih.gov. Discovery studio was used to convert SDF structures of ligands and the standard drugs into PDB chemical format. The ChemDraw version 19 was used to prepare unavailable structures, which were then converted to mol2 chemical format and then to PDB chemical format. As previously demonstrated, the non-polar hydrogen molecules were combined with the carbon atoms, and the polar Gasteiger-type hydrogen charges were assigned to the atoms in the chemical structure [79]. Internal degrees of freedom and torsions were both set to 0. The structures were then transformed using AutoDock Tools into the dockable PDBQT format.

4.3. Molecular docking calculations

The structures of the terpenes were imported through OpenBabel [80] into AutoDock Vina incorporated in PyRx 0.8 [71]. Energy minimization of the ligands was performed using Universal force field (UFF) using conjugate gradient descent as the optimization algorithm. The structures of the terpenes were then docked with the binding pockets of both DPP-4 and PTP1B. The active site regions of these enzymes were defined and selected using grid boxes with parameters depicted in table 1. The virtual screening was performed keeping all docking parameters and exhaustiveness as default. After the scoring, the docked poses were retrieved and the ligand-protein interactions visualized using Discovery Studio Visualizer version 16.

Table 8. Grid box parameters for defining active site regions of DPP-4 and PTP1B.

Dimensions	DPP-4 (Å)	PTP1B (Å)
center_x	34.88	14.79
center_y	29.17	-3.40
center_z	16.00	1.69
Size x	21.92	16.50
Size y	17.77	19.72
Size z	22.83	27.64

Table 8. Grid box parameters for defining active site regions of DPP-4 and PTP1B.

Dimensions	DPP-4 (Å)	PTP1B (Å)
center_x	34.88	14.79
center_y	29.17	-3.40
center_z	16.00	1.69
Size x	21.92	16.50
Size y	17.77	19.72
Size z	22.83	27.64

4.4. Prime MM-GBSA post docking calculations

Post docking free energy simulation through MM-GBSA was applied to the top ligand-protein complexes in order validate the docking score. To achieve this, the MMGBSA panel in Maestro was employed to compute the binding affinity of the top terpenes with the target enzymes as demonstrated in a previous study [81]. The MMGSA helped to estimate the difference in free binding energy between the terpene structures and the enzymes in both the non-bound state and complexed form after the process of energy minimization. The OPLS3forcefield parameters were used for the MM-GBSA, while VSGB was employed to model the continuum solvent. The computations were done with all other options set as default. The binding free energy computations were undertaken based on the following equations

$$∆G_{bind}=∆E+∆G_{solv}+∆G_{SA}$$

(i)

$$∆E=E_{complex}-E_{Enzyme}-E_{Terpene}$$

(ii)

$$∆G_{solv}=∆G_{solv (complex)}-∆G_{solv \left(Enzyme\right)}-∆G_{solv (Terpene)}$$

(iii)

$$∆G_{SA}=∆G_{SA (complex)}-∆G_{SA \left(Enzyme\right)}-∆G_{SA (Terpene)}$$

(iv)

Where

E is the minimized energies for molecular systems

$∆G_{SA}$ represents the non-polar contribution to the energy of solvation resulting from the surface area, while $G_{SA}$ is the surface energy of the systems.

4.5. Ensemble docking analysis

Three different coordinates of DPP-4 resulting from clusterization of the MD simulation trajectories were docked against the selected top three terpenes (T1, T4, T2) along with the reference compound (alogliptin) using AutoDock Vina software [71,72]. The selected terpene structures were docked following the docking protocol. The docking complexes were analyzed using Protein-Ligand Interaction Profiler (PLIP) web server and visualized using PyMOL 2.4 software [82].

4.6. Molecular dynamics simulation

Initially, the human DPP-4 (PDBID: 3G0B) protein was subjected to a 100 ns full atomistic MD simulation in order to generate an ensemble of the protein for performing the ensemble docking analysis. Then, the apo enzyme, and the DPP-4-terpene complexes obtained from docking calculations were subjected to MD simulation. These were performed using GROMACS 2019.2 and GROMOS96 43a1 selected as forcefield on the WebGRO [83,84,85]. The preliminary topology files of the terpene structures were prepared using PRODRG webserver (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg) [86]. The apo DPP-4 and the DPP-4-terpene complexes were solvated within a cubic box of the transferable intermolecular potential using a four-point (TIP4P) water model. The periodic boundary conditions were applied and a physiological condition of 0.154 M concentration set by neutralized NaCl ions. The minimization of the apo DPP-4 system and the complex systems was carried out in 10000 steps employing the steepest descent algorithm and applying the NVT ensemble for 0.3 nanosecond. This was followed by 0.3 ns of equilibration under NPT condition. The system temperature was set up and maintained at 310 K with the velocity rescale. The system pressure was also set and maintained at 1 atm employing the Parrinello-Rahmanbarostat. The integrator used for the computation was the Leap-frog integrator with a time step of 2 femtosecond. During the 100 ns full atomistic MD simulation production run, 0.1 ns snapshot was saved with a total of 1000 frame for each biomolecular system. Various thermodynamic parameters viz: RMSD, RMSF, SASA, RoG, and hydrogen bond number were computed from the trajectory files obtained from the dynamic simulation runs using VMD TK console scripts [87].

4.7. Clustering analysis

The trajectory files obtained from the initial MD simulation of the DPP-4 protein was clustered to generate the ensemble of the protein which were used for the ensemble docking simulation. Also, the trajectory files obtained from the simulation runs of DPP-4-terpene complex systems were clustered. The clustering analysis was performed employing the elbow method using TTClust V 4.9.0. Following cluster generation, a representative structure of each cluster of the complexes was chosen for interaction analysis with the aid of the Protein-Ligand Interaction profiler (PLIP) [88].

4.8. Binding free energy calculation using MM-GBSA

The Gmx_MMPBSA algorithm was used to calculate the binding free energy for each complex using Molecular Mechanics Generalized Born Surface Area (MM-GBSA). The salt concentration and the solvation method were 0.154 and 5, respectively. The internal and external dielectric constants were set to 1.0 and 78.5, respectively with other options as default. The decomposition of free energy was calculated to find which amino acids within 10 Å contribute most to the binding [89,90]. The MM-GBSA approach utilized is depicted in equation (v):

$$∆\mathrm{G}=<\mathrm{Gcomplex}-\mathrm{Greceptor}-\mathrm{Gligand}>$$

(v)

The ‘< >’ indicates the means of the free energies of complex, protein, and terpene structure over the frames used for the calculations. Several frames (200) were used to calculate the free energy. Various energy terms were computed using Equations (vi) to (x) as follows:

$$∆\mathrm{Gbind}=∆\mathrm{H}-\mathrm{T}∆\mathrm{S}$$

(vi)

$$∆\mathrm{H}=∆\mathrm{Egas}+∆\mathrm{Esol}$$

(vii)

$$∆\mathrm{Egas}=∆\mathrm{Eele}+∆\mathrm{EvdW}$$

(viii)

$$∆\mathrm{Esolv}=\mathrm{EGB}+\mathrm{ESA}$$

(ix)

$$\mathrm{ESA}=\mathsf{\gamma }.\mathrm{SASA}$$

(x)

Where:

∆H represents the enthalpy which is computed from gas-phase energy (Egas) and solvation-free energy (Esol). The T∆S term represents the entropy contributions to the total binding affinity which was not included because of the aim of comparing the relative binding free energies. Egas is comprised the electrostatic and van der Waals energy terms; Eele, EvdW, respectively. Esol was computed from the polar solvation energy (EGB) and nonpolar solvation energy (ESA) which was estimated from the solvent accessible surface area [91,92].

4.9. Density functional theory

On the basis of the electron density associated with the ligands, DFT asserts that ground state energy and other molecular attributes are determined only by this density. As a basis set for a DFT calculation using the B3LYP functional and 6-31G** as a single point energy calculation, the Jaguar panel of Maestro was employed in this work. A variety of molecular reactivity indicators, including electrophilicity, hardness, softness, electron affinity, ionization potential, and molecular electrostatic potential, were assessed. The following mathematical equations were used to perform the computations as discussed previously [93].

𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑛 𝑎𝑓𝑓𝑖𝑛i𝑡𝑦 (𝐸𝐴) ≈ − (LUMO)

(xi)

𝐼𝑜𝑛𝑖𝑧𝑎𝑡𝑖𝑜𝑛 𝑝𝑜𝑡𝑒𝑛𝑡𝑖𝑎𝑙 (𝐼𝑃) ≈ − (HOMO)

(xii)

$$𝐻𝑎𝑟𝑑𝑛𝑒𝑠𝑠 (𝜂) \approx (\frac{IP-EA}{2})$$

(xiii)

$$𝑆𝑜𝑓𝑡𝑛𝑒𝑠𝑠 (𝜎) \approx (\frac{1}{2n})$$

(xiv)

$$\mathrm{Ell}𝑒𝑐𝑡𝑟𝑜\mathrm{ga}𝑡𝑖𝑣𝑖𝑡𝑦 (𝜘) \approx (\frac{IP+EA}{2})$$

(xv)

Chemical potential (µ) = - (I + A)/2

(xvi)

Electrophilicity (ω) = µ²/2η

(xvii)

4.10. Physicochemical admetSAR analysis

The hit triterpene structures were passed through predictive physicochemical analysis. The drug-like properties were obtained from the SwissADME web werver platform (http://www.swissadme.ch/index.php) [94]. The ADMETox prediction was performed on the Protox-II webserver (https://tox-new.charite.de/protox_II/index.php) [95]. To perform these analyses, The SDF file format as well well as the canonical SMILES of the of the selected terpenes were retrieved from the PubChem Database and imported into the various webservers for druglikeness and ADMET computations using default parameters.

5. Conclusion

Drug-target investigations geared towards providing underlying mechanism of herbs, spices and medicinal plants are highly recommended for validating their acclaimed therapeutic potential. Catalogue of West African indigenous plants are extensively reported for protective and therapeutic potential against diabetes. However, most of the constituent phytochemicals that may account for the therapeutic potential are largely uninvestigated. We employed computational modeling in this study to assess the interactions of selected plant-derived terpenes with two emerging drug targets in diabetes. In this study, several triterpenes were reported to exhibit strong binding tendency with DPP-4 and PTP1B enzymes. Three hit triterpenes viz: Cucurbitacin B, 20-Epi-isoiguesterinol and 6-Oxoisoiguesterin are suggested as promising inhibitors of DPP-4. These compounds maintained strong binding affinity and stable interactions tendency with DPP4 as indicated by ensemble docking, post-docking MMGBSA, MD simulation, MD-based MMGBSA and density functional theory. The binding interactions of these triterpenes were sustained by interactions with several amino acid residues in the catalytic triad, oxyanion cavity, hydrophobic S1 pocket and charged S2 pocket in the active site region of DPP-4. In addition, the hit compounds exhibit good drug-likeness and moderate ADMET properties. Therefore, these triterpene structures are suggested for experimental validation towards their development as potent inhibitors of DPP-4 enzyme in diabetes.