Finding Natural Compounds That Help Activate the GLP-1 Receptor Using Computer Screening and Analysis

The molecular modeling software Maestro by Schrödinger (Version 2024-4 software release) [53] (RDIA Grant 12990-iau-2023-iau-R-3-1-HW: P.O. 6947 License key: 03cb87b8-723c-4fec-9b8c-8a58137d7a76) was used for computational studies. The desktop workstation was equipped with Intel^® Core™ i7-10700F Processor, Linux Ubuntu 22.10 operating system and an RTX 5000 graphics card.

The natural products (NPs) database consisting of 695133 NPs was downloaded from the Coconut website “https://coconut.naturalproducts.net/↗ (accessed on 20 June 2024)”. The obtained structures were filtered using Schrödinger’s Canva based on their physicochemical descriptors including MW, logP, number of rings, total heavy atoms count, total charge, number of H-donors, number of H-acceptors, number of rotatable bonds, total charge, and QED score. Final database filtration comprised no violation of rule of 5 [54] and a natural product-like (NPL) score <2 [55].

Shape screening was performed using the shape screening tool on Schrödinger Maestro. The filtered database was first energy minimized in Maestro using the OPLS3 force field and a default value of 0.30 Å for rmsd for non-hydrogen atoms and then subjected to shape screening using the structure of the crystal ligand as the reference structure. The typed pharmacophore technique was used for volume scoring. This technique evaluates the compounds based on their pharmacophore features. The shape similarity (SS) score was used to rank the compounds and a threshold of 0.3 was utilized to further filter the database for virtual screening [56].

The Cryo-EM structure of the GLP-1 receptor in complex with G protein, GLP-1 peptide, and a positive allosteric modulator (PDB ID: 6VCB) were extracted from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) [57,58].

The protein structure was prepared for docking using the Protein Preparation Workflow on Maestro. The preparation and minimization process were carried out at a pH of 7.4, and with adjusting ionization states. Polar hydrogens were added, and non-essential water molecules were deleted from the structures. Finally, the receptors were optimized using OPLS3 force field. In the final stage, the optimization and minimization on the ligand–protein complexes were carried out with the OPLS3 force field and the default value for rmsd of 0.30 Å for non-hydrogen atoms was used [59]. The receptor grids were then created at the center of the bound ligand utilizing a 1.00 van der Waals (vdW) radius. The scaling factor for vdW radius was used, along with a cutoff of 0.25 for partial charges. The binding sites were contained inside a grid box of 20 ų using default parameters and without any constraints.

The final filtered ligands were prepared using the built-in LigPrep tool in the virtual screening workflow available in Maestro. The ligands’ three-dimensional structures were created by adding missing hydrogens and generating the most likely ionization states at a pH of 7 ± 2 using Epik. Ligands’ geometry was then optimized with the OPLS3 force field with tautomer generation, desalting and producing 32 isomers per ligand at most [59]. The produced conformations represented the initial input structures for the virtual screening workflow.

The validation of the molecular docking method was performed to evaluate the accuracy of Maestro Glide to predict docking poses for the studied protein [60,61]. The crystal ligand was docked into GLP-1 receptor using the same criteria used for library screening. The docked pose with the lowest binding energy was aligned with the crystal structure conformation using Maestro’s structure superimposition feature. The root mean square deviation (RMSD) of the alignment was then calculated.

The virtual screening workflow involved three steps: high throughput virtual screening (HTVS), standard precision docking (SP), and extra precision docking (XP). For each docking protocol, flexible docking was selected with performing post-docking minimization, generating three poses per compound and holding 10% of the best scoring compounds. No more filters or constraints were applied in the docking process. The specific parameters used were a van der Waals (vdw) radius scaling factor of 0.80 and a partial charge cut-off of 0.15. The XP score was utilized to rank ligands and determine hit molecules. Molecular mechanics–generalized Born surface area (MM-GBSA) binding free energy calculation was applied from the virtual screening workflow to estimate the binding affinity for the top hit compounds using the Prime MM-GBSA module.

The docking output from the SP and XP docking was subjected to docked poses filtration based on previous experimental findings in the crystal structure using the pose filter tool in Maestro. The filtration criteria utilized were selected to must be all fulfilled and included contact with Leu 142 within 5 Å, π-π contact with Tyr 145 within 5 Å, H bond with Lys 202 within 2.5 Å, maximum H bond distance of 2.5 Å, minimum H bond-donor angle of 90 and H bond-acceptor angle of 60, and contact with Phe 12, Val 16, and Leu 20.

Molecular dynamic (MD) simulation studies were conducted using the Desmond Module on Schrödinger’s Maestro platform as previously described. In brief, the ligand–protein complex in its optimal docked configuration was first minimized using the Protein Preparation Wizard and prepared for MD simulation using the system builder application of Desmond. The simulation environment generated contained a water-based solvent system: the TIP3P water model. An orthorhombic simulation box with a 10 Å buffer parameter from the protein surface was generated; the entire system was neutralized by calculating and adding the required number of counter ions and 0.15 M NaCl in order to attain isosmotic conditions. MD simulation was carried out at a temperature and atmospheric pressure of 300 K and 1.013 bar, respectively. The simulation was run for a total of 100 nanoseconds (1001 frames were saved in order to compile the trajectory). Analysis was run and results presented using the simulation interaction diagram tool of Desmond.

We used the pkCSM web server (http://biosig.unimelb.edu.au/pkcsm/prediction↗ accessed on 11 November 2024) [62], in predicting the ADMET absorption, distribution, metabolism, excretion, and toxicity descriptors together with drug-likeness properties for the finally selected potential inhibitors. A total of eight molecular descriptors were generated in calculating the ADMET attributes in the potential KHK hits. In addition, Swiss ADME www.swissadme.ch/↗ accessed 11 November 2024 [63] was implemented for the calculation of the physicochemical properties, medicinal chemistry aspects, and drug-likeness attributes.

The virtual screening of the final database involved three stages including high throughput screening (HTVS), standard precision (SP), and extra precision (XP) docking. The protocol used a flexible docking method which incorporated post-docking energy minimization and was set to keep 10% of the best scoring compounds in each phase while retaining scoring hierarchy. The final output after the XP docking stage included 136 natural products. These compounds were further subjected to a second round of XP docking using the same criteria to produce a final hits list of 10 natural products. The binding energies (MM-GBSA) were calculated for the best 10 hits and were correlated with their XP score (Figure 3). The docking protocol was initially validated by calculating the RMSD of the crystal and docked poses of the co-crystallized ligand that was found to be 1.77 Å, a figure that establishes a good accuracy for glide to predict binding to the GLP-1 allosteric site.

The XP scores ranged from −9.6 to −7.6, with a maximum and minimum score difference of two units. Compared to the XP score (−3.3) of the crystal ligands, the identified hits showed better affinity to the GLP-1 allosteric site based on their docking scores. The calculated binding energies (MM-BBSA DG Bind) ranged from −50.8 to −32.4 and showed acceptable correlation with the XP scores especially for hits 1, 5, 6, 7, 8 that have SP scores below −8.0 and MM-BBSA DG Bind below −40.0 (Table 1). Given the fact that the hits are not a congeneric series, the ranking based on MM-GBSA DG Bind are not expected to perfectly correspond with the ranking based on docking scores and experimental binding affinity. MM-GBSA is more reliable in prioritizing compounds for experimental testing compared to alternative docking calculations [68]. It represents the calculated total free energies of the protein–ligand complex which is the sum of the van der Waals, electrostatic, General Born solvation, and surface area energies [69].

The crystal ligand binds in a pocket formed by residues within the interface between TM1 and TM2 and forms van der Waals (vdW) contacts with Leu 142 and π-π stacking with Tyr 145 [46]. In a mutagenesis study by Bueno et al. [46], Leu 142 and Tyr 145 single mutants reduced the potentiation effect of the allosteric modulator on GLP-1 binding to its receptor in cAMP assays. The Lys 202 mutant also showed a slight reduction in the potentiation effect. Triple mutation (Leu 142, Tyr 145, Lys 202) completely abolished the positive modulation activity and indicated π–π stacking is the primary interaction between the ligand and Tyr 145 [46].

In addition to the ligand–receptor interactions, ligand–GLP-1 contacts were also found to be detrimental. Establishing contacts with Phe 12, Val 16, and Leu 20 of the GLP-1 peptide is of great importance. Such interaction with Phe 12 of the GLP-1 peptide stabilizes its vdW contact with Leu 384 and Leu 388 of the receptor. The vdW interactions between the peptide and ligand are crucial for affinity as the affinity of the mutant peptide GLP-1 Phe 12, or Val 16, or Leu 20 can only be weakly potentiated by the allosteric modulator when compared to native GLP-1 [46].

The analysis of the binding pattern and interactions with the key residues revealed that the hits were able to bind the key residues Leu 142, Tyr 145, and Lys 202, combined or selected ones. The analysis also uncovered the importance of other important residues in the allosteric pocket that would increase the binding affinity of potential ligands and revealed additional important binding interactions with residues including Ser 206, Glu 138, and Asp 198 (Table 2). Interestingly, the hit compounds were able to establish H bonds with the key residues, an interaction that was absent from the resolved crystal ligand to 3.3 Å.

Comparative docking of the top ten hits with the crystal ligand showed that they became deeper in the binding site and closer to the GLP-1 peptide. This resulted in better interactions with key residues in both the GLP-1 receptor and GLP-1 peptide and an overall better XP score (Figure 4).

Hit 1 docking pose showed interesting, centered orientation of its fused ring that enabled vdW interactions with the receptor key residues; Leu 142, and Tyr 145 and the peptide key residues; Phe 12, Val 16, and to a lesser extent Leu 20. Other important noted electrostatic interactions involved a salt bridge with the receptor key residue Lys 202 and a H bond with Ser 206 (Figure 5).

The analysis was further extended to identify any compounds from the initial XP hits (136 compounds), SP hits (1532 compounds), that showed considerable interactions with the key residues and possess a good XP/Glide score. This was done through the docked poses filter tool on maestro. The filtration comprised a set of criteria that must be all fulfilled including contact with Leu 142 within 5 Å, π-π contact with Tyr 145 within 5 Å, H bond with Lys 202 within 2.5 Å, maximum H bond distance of 2.5 Å, minimum H bond–donor angle of 90 and H bond–acceptor angle of 60, and contact with Phe 12, Val 16, and Leu 20. The XP/SP score cutoff was set to −5.0 to ensure that other interactions are not negatively impacting the overall binding.

The shortlisted compounds were docked again in the GLP-1 allosteric site using extra precision and their binding energies were calculated to further assess their binding affinity. All compounds from the XP docking that had a score of −3.3 (crystal ligand XP score) and below with a negative binding energy are considered as final hits. This resulted in identifying 58 more hits that exhibited the spotted interactions together with fulfilling the XP score threshold and were found to be in the range −7.4 to −3.7 (Table 3). The MM-GBSA ranged from −52.77 for hit 14 to −13.02 for hit 48 which indicates favorable interactions between those hits and the GLP-1 allosteric site in terms of binding energies (see Supplementary Materials).

Hit 12 exhibited similar interactions to hit 1 including salt bridge with Lys 202 and H bond with Ser 206. However, unlike hit 1, hit 12 established the sought π-π stacking with Tyr 145. The hit also maintained the vdW interactions with Leu 142, Phe 12, Val 16, and Leu 20 (Figure 6).

Interestingly, Ampicillin (hit 44) was identified as a hit molecule, being able to establish the key interactions. Ampicillin bound the key residue Tyr 145 with its β-lactam oxygen through the H bond. Although this type of interaction with Tyr 145 was not mentioned in the binding requirements of the crystal structure, more investigations need to be performed to establish its feasibility and impact. Again, the salt bridge with Ser 202 and the H bond with Ser 206 were maintained for Ampicillin. Additional H bond with Asp 198, as well as the previously discussed interactions with Leu 142, Phe 12, and Val 16 were also noted (Figure 7).

In conclusion, the current docking studies not only identified new hits for GLP-1 positive allosteric modulation, but also shed the light about the importance of Ser 206 as a crucial residue with Lys 202 and to a lesser extent Glu 138, and Asp 198 that would enhance the affinity and selectivity of potential ligands.

The top hit molecules were selected based on a detailed visual analysis of the ligand–protein interaction profile. Hit compounds which had a favorable binding score and binding interactions with key residues in the binding pocket were selected for molecular dynamic simulation studies which were essential to validate the stability and strength of the binding interaction over a period of 100 nanoseconds. Of the 10 ligands used in MD simulation studies, 3 exhibited a robust and stable binding pattern within the allosteric binding pocket of the GLP-1 protein. Root mean square deviation (RMSD) plot analysis was used to measure the average displacement of atoms for a particular frame with respect to a reference frame. The analysis revealed a stable binding profile, as indicated by an average RMSD fluctuation range within 4 Å of the ligand with reference to the protein backbone throughout the simulation period. Figure 8 shows the RMSD plots for the top selected three compounds expressing the favorable binding profile with key residues in the binding pocket and also suited for drug repurposing: compound 3 (CNP0086660.2), compound 5 (CNP0039190.2), and compound 2 (CNP0549010.1). While a certain degree of fluctuation is anticipated within the simulation period of an allosteric binding pocket comprised the GLP-1 protein and its interacting peptide, the hit compounds appeared to maintain key interactions highlighted in Table 3 throughout a significant period of the 100 ns simulation period as depicted in Figure 9, Figure 10 and Figure 11. Of significance, the protein backbone in Figure 8 appears to express an expected degree of fluctuation but this too becomes stabilized within a tight range overtime. For further verification, the energy calculations were repeated for the hit compounds after molecular dynamic simulation runs using the final frame of the simulation trajectories. For the hit compounds, compound 3 and compound 2 and compound 5, the MM-GBSA calculations were found to be −39.22, −38.18, and −48.15, respectively; these results were found to be comparable to initial energy calculations using the optimal docking orientations.

Hit 3 of the macrolide class exhibits one of the most favorable binding profiles, interacting with a number of residues of the GLP-1 allosteric binding pocket as shown in Figure 9. Notably, GLP-1 receptor binding residues Lys 197 and Ser 206 form a stable hydrogen bond with the ligand for over 70% and 50% of the simulation period, respectively, in addition to a stable hydrogen bond formed with the GLP-1 peptide residue Tyr 19 which was stable for over 60% of the simulation period. All significant interactions are represented in Figure 9A,B which highlight the residues with the strongest ligand interactions stable throughout the entire simulation period. Figure 9C depicts all significant interactions displayed by the ligand and interacting residues occurring for over 10% of the simulation period. Macrolides are a class of broad-spectrum antibiotics which are granted FDA approval for the management of a wide range of infectious disease conditions. The approval status and safety profile of this class of medications renders it particularly amenable to subsequent repurposing studies.

MD simulation results for hit 5 (CNP0039190.2) in Figure 10 show significant interactions with the key binding residues of the allosteric pocket. Key binding residues Tyr 145 forms hydrophobic interactions in the form of π-π stacking with the compound whereas a complex interaction is observed with Lys 202 and the terminal COO– group of hit 5. The interactions with the GLP-1 receptor are observed for over 30% of the simulation period.

Figure 11 shows the binding profile of hit 2 (CNP0549010.1) with residues in the allosteric GLP-1 receptor binding pocket. Hit 2 forms moderate interactions with the key binding residues Tyr 145 and Lys 202 of the receptor binding pocket: hydrophobic and water bridge interaction with Tyr 145 and a complex hydrogen and water bridge interaction with Lys 202. Additionally, a complex hydrogen bond and water bridge interaction is observed between hit 2 and the GLP peptide residue Gln 23 for over 20% of the simulation period. While the strength of the interaction is relatively weaker compared to other hit compounds, this hit remains significant owing to its classification as a beta lactam compound. Beta lactam skeletons form the essential component of the most potent class of antibiotics clinically available. They are widely used and have an excellent safety profile making this compound a desirable target for further drug repositioning studies.

The hit compounds were grouped according to their natural product classification and were surveyed in the literature for their biological function to identify potential correlations with postulated GLP-1 modulation activity. The following classes have been identified and shortlisted from the literature given their favorable pharmacodynamic profile:

Novel scaffolds for potential GLP-1 positive allosteric modulation were identified from the final XP hits as well as from the original SP and XP docking shortlists after docked poses filtration and XP redocking and MM-GBSA calculations (Figure 12). Those final hits were structurally analyzed to identify common structural scaffolds. These scaffolds would serve as lead compounds to design and develop new GLP-1 allosteric modulators. The list encompassed a variety of compounds from different natural products classes including β-lactams, alkaloids, steroids, sesquiterpenoids, macrolides, and polyketides that possess diverse chemical scaffolds. Peptidomimetic scaffolds were the most frequent scaffolds among the identified hits (36 compounds out of 68 compounds). In these compounds, the peptide moiety was attached to various chemical entities including indole derivatives (hits; 1, 13, 14, 19, 20, 22, 27, 29, 32, 33, 40, 47, 49, 52, 54, 55, 62), coumarin derivatives (hits; 11, 12, 35, 39, 42), azines and benzazines derivatives (hits; 8, 10, 14, 38, 53, 59, 65, 68), and azoles derivatives (hits; 23, 45, 56). The peptide moiety was also found to be incorporated in a cyclic structure including hits; 14, 21, 22, 27, 52, 53, 55, 59 as well as a chain structure including hits; 1, 2, 5, 10, 11, 12, 13, 19, 20, 25, 29, 32, 33, 35, 37, 38, 40, 43, 44, 45, 46, 47, 49, 62, 64, 65, 67, 68. β-lactams which could also be considered as cyclic peptidomimetic structure were among the most frequently identified hits (hits; 2, 25, 44, 46, 67). A significant number of coumarin derivatives were also identified including hits; 36, 39, 41, 42, 66. Sulfonamide derivatives were represented by five hits including hits; 31, 34, 37, 56, 59. In conclusion, the structure analysis of the identified hits revealed common structural scaffolds including peptidomimetics, β-lactams, coumarins and sulfonamides, among which peptidomimetics possessed the most significant frequency especially in indole and coumarin cores.

Some interesting features were shown for the top ten hits from the ADMET predicted calculations that influence pharmacokinetic and pharmacodynamic properties (Table 4). Water solubility was highly variable among the hits, with the lowest for hit compound 9 at −4.41 log mol/L, which may indicate potential formulation challenges. Intestinal absorption is notably higher for compounds 5 to 9, exceeding 90%, which suggests that these molecules have a favorable bioavailability profile. On the other hand, poor absorption was observed with compounds 2 and 4, probably due to their low solubility and permeability. In the Caco-2 permeability test, compound 9 was identified as the most permeable, with a log Papp of 1.376, further supporting its potential for effective gastrointestinal absorption.

Compounds 8 and 9 may cross the blood–brain barrier considering log BB > 0 given the data for BBB permeability distribution, thus standing out as candidates for CNS-targeting therapies. However, the value for CNS permeability logs PS for all compounds remained low, reflecting restricted access to the tissues of the central nervous system. compound 5 showed the most favorable CNS permeability characteristics, recorded at −2.132 log PS. These results emphasize the importance of refining the distribution parameters during the lead optimization phase, depending on the target of the intended therapy.

Metabolic interactions with the cytochrome P450 enzymes are variable amongst the compounds. Notably, compound 8 was a substrate of CYP2D6 and CYP3A4, possibly suggesting potential for metabolic interactions. The remaining compounds interacted considerably less with the enzymes, and in most instances, one would not predict clinically significant drug–drug interactions. Compound 9 also inhibited CYP2C19 and may affect its metabolism; in vivo investigation is required.

The excretion analysis, judged from total clearance parameters, identified six exhibiting a high clearance value of 1.101 log mL/min/kg; this outcome indicates a possibly short half-life. However, compound 9 revealed negative clearance values, suggesting a possible accumulation and a long systemic exposure. Renal excretion studies revealed that compounds 6 and 8 are substrates for renal OCT2, indicating that active transport plays a primary role in their elimination processes.

Toxicity analyses revealed that all the hit compounds expressed no significant toxicity since none were AMES toxic nor inhibited the hERG channels. However, compounds 1, 2, 3, 5, and 10 were hepatotoxic, which should be well considered in further stages of development. The maximum tolerated dose was predicted; the highest tolerance was for compound 2, calculated as 1.50 log mg/kg/day, while the lowest was for 8.

The selected molecules in this study were diverse in their physicochemical, drug-likeness profiles, and medicinal chemistry attributes (Table 5). The MW of the compounds ranged between 264.36 and 410.44 g/mol, hence within the acceptable threshold for small molecules that would enhance permeability and good bioavailability. Lipophilicity was represented by log P values ranging from −1.40 to 3.37. This indicated an appropriate range of hydrophilic–lipophilic properties, supporting the solubility and membrane permeability of most compounds. The number of hydrogen bond acceptors and donors corresponded to Lipinski’s Rule of Five, showing that the structures would be in compliance with the criteria for drug-likeness. At the same time, rotatable bonds and molar refractivity values support these molecules for conformational flexibility and molecular interaction.

From a drug-likeness perspective, all the compounds passed Lipinski’s rule without violation and would thus be orally bioavailable. Veber and Egan showed some alerts due to high topological polar surface area (TPSA), above the threshold in some cases. This minor deviation does not disturb the overall good bioavailability score in the dataset, as the majority of the molecules showed bioavailability scores above 0.55. Although the Ghose and Muegge filters also occasionally violated other parameters, such as XLOGP3 and weight, these were molecule-specific and hence too weak to overpower their overall potential to be good candidates.

Medicinal chemistry analysis showed that none of the compounds had any PAINS alerts; hence, none of these compounds was likely to show false-positive activity upon high-throughput screening. Brenk alerts were seen for some of the molecules; these had structural features like phthalimide which may need further optimization. It was noticed that more than half of these molecules congregate under lead-likeness criteria, while a few exceptions took values above the threshold in rotatable bonds and/or molecular weight. Synthetic accessibility scores spanned a medium range from 3.92 to 5.58, reflecting the fact that it would not result in considerable problems.

In silico drug experiments using molecular docking to target GLP-1 receptor (PDB ID: 6VCB) were obtained from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB). PDB DOI: “https://doi.org/10.2210/pdb6VCB/pdb↗ (accessed on 20 June 2024)”. Natural products database was downloaded from the Coconut website: “https://coconut.naturalproducts.net/download↗ (accessed on 20 June 2024)”. All data generated or analyzed during this study are included in this published article.