Finding non-peptide GLP-1 receptor activators using computer-based screening and simulation

The resultant databases underwent three steps of virtual screening in Schrödinger’s Glide: high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP) docking [40]. Ten percent of the top-scoring compounds were retained in each phase, while the scoring hierarchy was maintained due to the protocol’s flexible docking approach, which included post-docking energy minimization. The final output following the XP docking stage included 276 and 36 compounds from the COCONUT and Marine databases, respectively. The XP scores of the top 10 hits were correlated with their binding energies (MM-GBSA) (Figure 3). The docking methodology was first verified by calculating the root mean square deviation (RMSD) of the crystal and docked poses of the co-crystallized ligand, resulting in 2.03 Å, indicating that Glide can predict binding to the GLP-1R binding site with high accuracy. PF-06882961 was also redocked successfully into its native receptor (6X1A) with RMSD 2.03 Å and GlideScore of −12.040. Cross-docking to non-native GLP-1R structures (6X19 and 6ORV) yielded lower RMSD values (1.15 Å and 1.19 Å), indicating higher pose agreement and demonstrating receptor pocket similarity as well as method robustness. GlideScores were slightly less favorable in cross-docked poses but were still indicative of strong binding (Table 1). These results demonstrate the strengths of ensemble-based docking methods and the limitations of one-pose redocking evaluation.

The results of the top hits from the COCONUT and Marine databases, as shown in Table 2 and Table 3, provide valuable insights into how well each compound binds to the GLP-1 receptor (PDB ID: 6X1A) through the descriptors of XP Score, MM-BBSA dG_Bind, and Shape Similarity. In the COCONUT database, the top hits displayed XP scores ranging from −13.835 to −12.532, indicating strong binding affinities to the GLP-1 receptor, with 1 exhibiting the highest XP score of −13.835. The MM-BBSA dG Bind values, indicative of binding affinity and measured in kcal/mol for these hits, ranged from −84.69 to −66.07, showing a significant correlation with the XP scores and further implying strong interactions between the ligands and the receptor. The Shape Similarity value reflects the geometric similarity of the hits compared to a reference ligand, with hits falling in the intermediate similarity range of 0.305–0.346 compared to the control compound (the co-crystallized agonist PF-06882961) (Shape Similarity = 1.000). In the Marine database, the mean XP scores for the top 10 hit compounds were between −12.80 and −10.190 kcal/mol. The best among them is 11, which achieved the highest binding affinity (−12.82 kcal/mol). The binding energies—MM-BBSA dG_Bind values for those compounds are between −100.43 and −52.04 kcal/mol, from which compound 16 exhibited the most favorable chemical affinity of −100.43 kcal/mol. As for the comparisons of Shape Similarity of Marine hits from 0.302 to 0.375, different compounds fit moderate to high degrees of similarity to the control compound. It should be noted that only one compound from the Marine database, 15, obtained a Shape Similarity score of 0.375, which further underpins its potential as a promising hit. The results confirm the existence of compounds within the COCONUT and Marine databases that possess strong binding affinities and some promising shape characteristics to the reference ligand, making them suitable for further examination as non-peptide GLP-1RAs. Favorable binding-based 2D and 3D docking interactions with key amino acid residues involved LYS-197, TRP-203, and PHE-381, which proved essential for the high-affinity interactions observed with the GLP-1R binding pocket (Figure 3). These results imply that the screened compounds are promising hits for targeting the GLP-1 receptor, enabling a more potent and selective interaction.

Analysis of the binding interactions of the identified ligands was based on identifying their interactions with the key residues identified in the crystal structure as crucial for binding the GLP-1 peptide, namely, LYS-197, TRP-203, and PHE-381. Common additional interactions were also searched to identify potential binding poses and important residues that might drive future drug design and development.

Strong π–cation interactions were noted for hit 1 with the pyridine nitrogen and TRP-203 and PHE-381. The indole ring further anchored the compound through a π–π interaction with TRP-33. An additional H bond was established between the acetate side chain and GLN-221.

Hit 5 was found to bind the active site with a similar binding pose to hit 1. The pyridine nitrogen interacted through π–cation interactions with TRP-203 and PHE-381. The polycyclic aromatic system established several π–π interactions with TRP-33, and the dioxolane oxygen interacted with an H bond as acceptor with GLN-221. The same interactions were completely observed between hit 9 and the same previously mentioned residues. Hit 11, however, lacked the interaction with PHE-381 while keeping the same interaction profile with TRP-203 and GLN-221. The missed interaction was compensated for by additional H bonds with SER-31 and LEU-32. On the other hand, hit 13 lacked the interaction with PHE-381 without additional bonding (Figure 3).

In summary, the binding profile of the top hits exhibited consistent binding with the key residues TRP-203 and GLN-221 through π–cation or π–π interactions and H bond as acceptor, respectively. The binding analysis also identified TRP-33 as an important binding residue that would act as a gatekeeper at the entrance of the binding site and that is preferably bound through a π–π interaction with a planar aromatic system. Moreover, the high values of the docking scores (−13) and binding energy calculations (−85) validated the significance of the drawn binding pattern and residues.

In mutagenesis studies, the activity of the non-peptide control compound was found to be reduced by alanine mutations of many residues, including ARG-380, GLU-373, PHE-385, LYS-197, TRP-33, PHE-230, LEU-141, TRP-203, PHE-381, LEU-384, and ARG-190 [9]. With this knowledge, we have used MD simulations over a period of 100 and 500 nanoseconds to evaluate the stability of ligand-target complexes, the binding poses from the docking calculations, and the interaction pattern between each ligand with key GLP-1R amino acids. The studies were performed for top select candidates from virtual screening, namely compounds 1, 5, and 9, and for the control compound as a reference. To assess the stability of the protein, we evaluated the RMSD values of Cα atoms (Figure 4).

The proteins in the complexes of the control and compounds 1 and 5 stabilized at around 200 ns, while the protein in the 9-GLP-1R complex stabilized at around 250 ns. The proteins in all complexes remained stable for the rest of the simulation time with no significant structural deviations, indicating that they had reached equilibrium. Fluctuations from average Cα-atom RMSD values were within an acceptable range of 1–3.5 Å, which is expected since the GLP-1R protein is a large, multi-domain, G-protein coupled structure. The average protein Cα atoms RMSD values were in the range of 4.7–6.5 Å. RMSF analysis shows that the ligand binding site in the TM region remains stable throughout the simulation for all the complexes (Figure 5).

To assess ligand stability, RMSD values for ligand heavy atoms were measured with respect to the protein (Figure 6). None of the average ligand heavy-atom RMSD values is significantly larger than the average protein Cα atom RMSD value. This finding indicates that compounds 1, 5, 9, and the control were all bound efficiently within the binding site, which was also confirmed visually (Supplementary Materials Figures S1–S3, residues that interact for 45% or more of the simulation time are shown). Based on RMSD analysis and visual inspection, compound 5 seems to have changed its orientation after 200 ns has passed, but it remains stable. Similar to the control, compound 9 remains stable throughout the simulation time. Compound 1 shows some fluctuations within the binding side until the pyridoindole group is in an optimum position to form H-bonds with GLN-221, and the alcohol group is in an optimum position to form H-bonds with THR-298 and ARG-299. Compound 1 then binds stably to these residues for the rest of the simulation time. Except for compound 1, all ligands had no major structural deviations from their starting docking conformations. Average RMSD values for ligand heavy atoms ranged from 3.2–5 Å. Compounds 9 and 1 had the closest average RMSD values to the control, while compound 5 had the highest RMSD value.

Control-GLP-1R complexThe control compound forms a hydrogen bond with GLN-221. More importantly, it forms strong and stable π–cation interactions with three key amino acids, TRP-33, TRP-203, and PHE-381, for 87%, 80%, and 49% of the simulation time, respectively (Figure 7A). These key interactions stabilize the complex and contribute to the known agonist activity of the control compound. The average binding free energy for the control-GLP-1R complex was −97.30 kcal/mol.Compound 1-GLP-1R complex

Compound 1 forms H-bonds with GLN-221, THR-298, and ARG-299 (Figure 7B), and the former two can stabilize the complex since they last for more than half the simulation time. It also forms a π–π stacking interaction with TRP-214 and a π–cation interaction with TRP-203. Although TRP-203 is important for the activity of non-peptide GLP-1R agonists, the interaction of this residue with compound 1 is not strong. The average binding free energy for the compound 1-GLP-1R complex was −97.87 kcal/mol, similar to the control ligand. This could indicate that compound 1 is a good binder.

Compound 5 forms a H-bond with GLU-138 (Figure 7C). It also forms π–π stacking interactions and π–cation interactions with two key amino acids, TRP-33 and TRP-203, with the latter lasting for more than half the simulation time. Such strong interactions with TRP-203 can stabilize the complex and contribute to compound activity, potentially as a GLP-1R agonist. The average binding free energy for the compound 5-GLP-1R complex was −92.50 kcal/mol, a value close to that of the control ligand. These findings indicate that compound 5 may be a potential GLP-1R agonist.

Compound 9 forms a H-bond with ASP-198 through a water bridge (Figure 7D). It also forms π–π stacking interactions with important amino acids TRP-203 and PHE-381 for more than half the simulation time. Stable π–cation interactions are also observed between compound 9 and the key amino acids TRP-33 and TRP-203, lasting for 89% and 93% of the simulation time, respectively. The average binding free energy for the compound 9-GLP-1R complex was −102.78 kcal/mol, the lowest value among all four complexes. Compound 9 seems the most promising as a GLP-1R agonist.

The average MM-GBSA ΔG_bind values obtained from the MD simulation frames of the control and the three compounds were low, comparable, and indicative of strong binding. The affinity ranking of compounds based on MM-GBSA ΔG_bind values often correlates with the ranking observed in experimental settings. Accordingly, the binding affinity of the ligands can be predictively ranked as 9 (−102.78 kcal/mol) > 1 (−97.87 kcal/mol) > control (−97.30 kcal/mol) > 5 (−92.50 kcal/mol). Compound 9 outperformed compounds 1, 5, and the control in terms of predicted binding affinity.

Results from the MD simulations and the binding free energy calculations indicate that compound 9 is the most promising since it demonstrates the highest predicted affinity. For most of the simulation time, this compound maintains stable interactions with three key amino acids that are important for non-peptide agonist activity. Its binding profile in terms of energy, interactions, and stability within the binding site is similar to, and sometimes better than, that of the control, which further supports its potential as a promising non-peptide GLP-1R agonist. On the other hand, compound 1 has favorable energetics and seems to strongly bind to GLP-1R through different amino acids than the other three compounds. This warrants further investigation into its activity to determine whether such binding behavior will also result in GLP-1R activation. Compounds 1 and 9 have comparable ligand heavy-atom average RMSD values to the control and remain stable in the binding site. Compound 5 maintains interactions with one important amino acid for more than half the simulation time and has a favorable binding free energy. All three hits interact with at least one residue for more than 50% of the simulation time, show protein stability, and exhibit minimum structural fluctuations for at least half of the simulation time.

Despite not having direct experimental or clinical evidence concerning the exact molecules, a complete review of specific compounds from the COCONUT and CMNPD (marine) databases reveals a diverse array of structural classes with potential anti-diabetic properties (Table 4 and Table 5). Among the compounds originating from coconut, CNP0593098.1 (3,4,5,6-tetradehydroyohimbine) stands out. No research has directly addressed this derivative, but anti-diabetic activity has been demonstrated for its parent compound, yohimbine, in animal models via α2-adrenoceptor antagonism, resulting in increased insulin secretion and improved glucose tolerance. Neither this derivative nor other derivatives of yohimbine have, however, been able to show GLP-1 agonist activity.

Other compounds, as found from COCONUT, include CNP0402650.0 and CNP0311770.0 (chromen derivatives). These belong to structural classes with an established reputation for antioxidant activities, inhibitory effects on carbohydrate-hydrolyzing enzymes (like α-glucosidase), and insulin-sensitizing effects. Despite this, no studies have researched these specific molecules for these functions. Other compounds from coconut, such as CNP0542406.3, CNP0510864.0, and CNP0294111.0, are also characterized by their structural cores of indole, pyridine, and benzofuran, which are recurrent in anti-diabetic studies, mainly for their oxidative stress reduction, insulin sensitizing, and pancreatic β-cell protecting effects; once again, no literature records exist for these particular compounds.

Continuing in the same vein, marine-derived compounds cataloged under CMNPD appear to harbor equally enticing avenues. Some of the quinazoline alkaloids such as CMNPD5314 (fumiquinazoline B) and CMNPD27343 (versiquinazoline B) resemble structurally known compounds that inhibit α-glucosidase and protein tyrosine phosphatase 1B (PTP1B), two well-studied anti-diabetic targets. Other marine substances, comprising CMNPD2041 (waixenicin A), CMNPD24592 and CMNPD24593 (speradines D and E), and CMNPD26002 (a gliotoxin derivative), are diterpenes and alkaloids that are known, at least in structurally related forms, for anti-inflammatory and metabolic regulating actions, including PTP1B inhibition and AMPK activation.

Marine-derived epipolythiodioxopiperazines, in addition, activate the Nrf2/HO-1 pathway and thus would confer cytoprotective and antioxidant effects that preserve β-cell function. Even though there are no direct studies regarding anti-diabetic or GLP-1 activities for specific compounds such as CMNPD27661, CMNPD4014, CMNPD9270, and CMNPD22795, the structural classes to which they belong, including polyketides, prostaglandin derivatives, and benzofurans, are increasingly associated with redox modulation, enzyme inhibition, and generalized metabolic benefits.

In summary, compounds from both coconut and marine natural product databases present strong potential as anti-diabetic agents from the standpoints of structural features, mechanistic plausibility, and class-based pharmacological evidence. Compounds from coconut are notably important for insulinotropic and β-cell-protective effects, while those marine ones would provide additional scaffolds for enzyme and redox pathways involved in glucose homeostasis. Although bioactivity data on these particular molecules are absent, the confluence of evidence available for their structurally related counterparts renders them promising agents in multi-target anti-diabetic drug discovery and clearly worthy of intensive evaluation.

The literature-supported findings affirm the pharmacophore hypotheses and in silico prediction proposed in the current investigation. They promote the contention that natural product scaffolds, selected via rational filters combined with docking and shape-similarity methods, constitute a rich source of non-peptidic GLP-1 receptor modulators.

In Table 6, the water solubility of all compounds ranges from −3.851 to −2.791 log mol/L, with 3 exhibiting the highest solubility and compound 6 showing the lowest solubility. Despite this difference, most compounds have moderate to low solubility, with 5 below average. With respect to the Caco-2 permeability scale, 1, 2, 5, 6, 8, and 9 can be stated to have good membrane permeability (>0.5 log Papp), allowing easy passage across the intestinal epithelial barriers; while 3, 4, 7, and 10, get below this threshold, showing less permeability. Intestinal absorption from the human GI tract is generally high, above 95% for compounds 1, 2, 5, 6, 8, and 9, implying good oral bioavailability. On the contrary, 3 is poorly absorbed (3.62% absorption), with subpar values for 4 and 10 (48%). Nearly all the candidates are P-glycoprotein substrates, except for 3 and 7, which probably increase retention inside the cell.

On the permeability of the blood–brain barrier (BBB), all compounds show a log BB value under the +0.3 threshold, thereby implying low or just moderate penetration into the CNS. Compounds 1, 2, and 8, albeit close to this threshold, are considered moderate. CNS Permeability (log PS) values are uniformly low, as all compounds apart from compound 8 (−1.894) are below the −2.0 mark, confirming little CNS exposure and so little risk of CNS-associated side effects.

Compounds 1, 2, 6, and 8 are substrates of CYP3A4. Compounds 3 and 4 are substrates of CYP2D6, which may impede their metabolic stability. Regarding CYP450 inhibitory activity, the scenario changes since 1, 2, 4, 5, 6, and 8 inhibit at least one isoform, interfering with the metabolism of some co-administered drugs. Compound 4 inhibits three CYP enzymes (i.e., CYP1A2, CYP2C9, and CYP3A4) and can, therefore, be considered at greater risk for metabolic interaction. Compounds 3, 7, and 10 show little or no CYP inhibition and therefore seem to possess a clean metabolic profile.

Total clearance is from 0.099 to 1.189 log mL/min/kg. A greater clearance above 1 is for 1, 2, 5, and 8, which means that they are quickly eliminated systemically. Compound 7 has a very low clearance value and therefore remains exposed longer. Only compounds 5 and 8 are substrates of renal OCT2, suggesting that renal excretion may occur, a factor that should be kept in mind when considering dosage for patients with renal impairment.

Toxicological profiling is suspected to be a risk for several candidates. Compounds 1, 2, and 8 are AMES-positive, which points toward mutagenicity. The maximum tolerated dose varies among compounds 3, 4, 6, 7, and 10, and is generally high. Hepatotoxicity seems to be the leading cause for concern, observed in all except compounds 5 and 7, thus supporting the presumption that liver safety is a top priority. It is interesting to note that no compounds are considered to inhibit hERG I channel, thus leaving cardiotoxicity out of consideration.

Most of the compounds (1, 2, 5, 6, and 8) exhibit appreciable absorption; however, compound 3 is poorly soluble and absorbed, and thus requires structural modification. None of the compounds are presumed to cross the BBB or to promote significant CNS penetration. Compounds 3, 7, and 10 exhibit minimal CYP inhibition, whereas compounds 4 and 6 demonstrate high metabolic liability. Most compounds were cleared efficiently, with compound 7 being an outlier. In terms of toxicity, compound 5 exhibits a more favorable profile, being free from AMES toxicity and hepatotoxicity. These profiles are used to lay the groundwork for prioritizing experimental validation.

The physicochemical and drug-likeness profiling of the tested compounds revealed overall favorable properties with respect to oral drug development guidelines (Table 7). Most compounds exhibited molecular weights below the 500 Da threshold, except for the control (543.59 Da), aligning with Lipinski’s rule of five. LogP values ranged from 0.68 to 3.22, indicating moderate lipophilicity, and most compounds had acceptable numbers of hydrogen bond donors and acceptors. Polar surface area (TPSA) was within reasonable values (mainly < 140 Ų), which are suitable for good membrane permeability, except for 3 (152.08 Ų). Most molecules met more than one drug-likeness filter (Ghose, Veber, Egan, and Muegge), but 3 always diverged and were not capable of crossing all filters and were more flexible and polar. PAINS and Brenk alerts were low across the dataset, suggesting a low risk of assay interference or structure risk. Synthetic accessibility scores ranged from 2.68 to 4.75, which translates to moderate to excellent synthetic viability. The data in total give support to the validity of these molecules as potential drug-like hits, with 3 and the control when it comes to polarity and complexity.

For our computational studies, we utilized the Schrödinger’s Maestro molecular modeling software (Versions 2024-4 and 2025-1) [32,67,68]. Our computational work was conducted on a desktop workstation equipped with an Intel^® Core™ i7-10700F Processor (Intel, Santa Clara, CA, USA), running the Linux Ubuntu 22.10 Operating System, and featuring an RTX 5000 graphics card (RDIA 12990-iau-2023-iau-R-3-1-HW: P.O. 6947 License key: 03cb87b8-723c-4fec-9b8c-8a58137d7a76).

We obtained a remarkable comprehensive natural products (NPs) database that contains 695,133 NPs from the Coconut website https://coconut.naturalproducts.net/↗ (accessed on 15 November 2024) [35]. In addition, we downloaded over 46,000 compounds from the Comprehensive Marine Natural Products database (CMNPD), which was retrieved from the online server (https://www.cmnpd.org/↗, accessed on 18 November 2024) [36].

The collected structures were filtered using Schrödinger’s Canvas based on their physicochemical descriptors, including MW, logP, number of rings, total heavy atoms count, overall charge, number of H-donors, number of H-acceptors, number of rotatable bonds, total charge, and QED score. The final database filtration ensured compliance with the rule of five [69] and maintained a natural product-like (NPL) score < 2 [70]. The NPL score, introduced by Ertl et al. is a numerical measure that quantifies how much a chemical structure resembles a natural product compared to a synthetic molecule. The NPL score ranges between −5 (if the compound is more similar to a synthetic compound) and 5 (if the compound is more similar to a natural product). A score of 0 or higher is often used to distinguish between natural and synthetic molecules. The higher the score, the more similar the molecule is to natural product structures. It is used to prioritize molecules for screening in drug discovery programs [71]. Given the big number of compounds screened in our research, we used a cut-off of 2 for NPL scores to exclude compounds that are more natural product like and would have higher probability of being orally inactive.

Our investigation was carried out using Schrödinger’s Maestro’s form screening tool. We started by adding a default RMSD of 0.30 Å, specifically for non-hydrogen atoms, to the filtered database in order to decrease its energy using the Optimized Potentials for Liquid Simulations 3 (OPLS3) force field. Following this optimization, shape screening was carried out with the crystal ligand structure serving as a guide. We used a number of different pharmacophore types in conjunction with a volume-scoring algorithm to precisely evaluate the compounds. We ranked the compounds according to their shape similarity score, setting a limit of 0.3 to further refine our virtual screening database. The evaluation was based on pharmacophore features [72].

The crystal structure of GLP1 in complex with the non-peptide agonist PF-06882961 (PDB ID: 6X1A) was obtained from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) (Figure 8) [73].

This protein structure was prepared for docking using Maestro’s Protein Preparation Workflow. To fine-tune the ionization states, the preparation and minimization process were carried out at a pH 7.4. We added polar hydrogens and eliminated extraneous water molecules from the structures. Using the OPLS3 force field to optimize the receptors was the last stage. Under the OPLS3 force field, the ligand-protein complexes were optimized and minimized using a preset RMSD value of 0.30 for non-hydrogen atoms [74]. Then, receptor grids with a radius of 1.00 A van der waals (vdW) were made using the center of the bound ligand. In this procedure, a vdW radius was applied and the threshold for partial charges was set at 0.25 [75]. The binding sites were contained within a grid box of 20 Å using default parameters and without any constraints

The LigPrep tool (Versions 2024-4), integrated into the virtual screening workflow of Maestro, was used to prepare the filtered ligands. The three-dimensional structures of the ligands were created by the addition of hydrogen atoms, followed by the generation of the most probable ionization states at pH 7 ± 2 using Epik. The geometry of the ligands was then optimized with the OPLS3 force field, which included tautomer generation, desalting, and producing a maximum of 32 isomers per ligand [67,76]. The resulting conformations served as the starting input structures for the virtual screening workflow.

The molecular docking protocols were evaluated in order to confirm Maestro Glide’s accuracy in predicting the ligand’s docking poses for the protein under study [77,78,79]. The library screening criteria were used to re-dock the cognate ligand into the GLP1 receptor.

To validate the robustness of the molecular docking protocol, a cross-docking experiment was conducted on GLP-1R. Three high-resolution crystal structures of GLP-1R in complex with different small-molecule ligands were retrieved from the Protein Data Bank (PDB): 6X1A, 6X19, and 6ORV. These structures represent slightly different conformations of the receptor–ligand binding domain. PF-06882961 was extracted from its native complex (PDB: 6X1A) and docked into all three receptor grids using the Glide Extra Precision (XP) docking protocol in Maestro. The results were calculated using the Ligand RMSD tool and GlideScore (Versions 2024-4).

Three distinct methods were employed to structure the virtual screening docking workflow: High-Throughput Virtual Screening (HTVS), Standard Precision (SP), and Extra Precision (EP). Therefore, flexible docking was utilized in conjunction with post-docking minimization for each approach, where each compound generated three poses and retained up to 10% of the highest-scoring compounds; this filtration criterion is acceptable when screening large databases [32]. A more rigorous hit identification stage was then applied to the initial filter through XP docking, MM-GBSA calculation, and lastly 500 ns MD simulations for the final hits. During the docking procedure, no additional filters or constraints were applied. The actual parameters included a vdw radius scaling factor of 0.80 and a partial charge cut-off of 0.15. The ligands and hit molecules were ranked using the XP score. For binding free energy calculations in molecular mechanics-generalized Born surface area (MM-GBSA), the virtual screening pipeline estimated the binding affinity of top-scoring hits using the Prime MM-GBSA module (Figure 9).

Molecular dynamics (MD) simulations were conducted for compounds 1, 5, 9, and the control compound, using the Desmond Module in the Schrödinger suite [80,81,82]. Desmond’s system builder tool was used to prepare the protein–ligand complexes for MD simulations. Since GLP-1R is a transmembrane (TM) protein and our ligands mainly bind in the TM region, a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipid bilayer membrane was added at 300 K [83]. The simple point charge (SPC) water model was used for an orthorhombic box with dimensions of 10 Å × 10 Å × 10 Å [84]. A 20 Å distance from the ligand covered the region for the exclusion of ion and salt placement. The number of required counter ions was calculated and subsequently added to neutralize the system. Sodium and chloride ions were used to maintain a 0.15 M salt concentration, and the OPLS5 force field was selected for running the simulations [85]. The environmental conditions included a temperature of 300 K, a pressure of 1 bar, and the NP γT ensemble were set [86,87]. The membrane model system was relaxed before running 100–500 ns simulations. Coordinates were recorded every 500 ps, and 1000 frames were generated for the trajectory. The temperature was controlled using the Nose–Hoover chain coupling scheme with a 1 ps coupling constant. The pressure was controlled using the Martyna–Tobias–Klein coupling scheme with a 2 ps coupling constant. A time step of 2 fs was selected for RESPA (Reference System Propagator Algorithm), and a 9 Å cutoff radius was set for Coulombic short-range interactions [87]. Desmond’s Simulation Interaction Diagram tool was used to analyze data generated from the MD simulations [88]. The RMSD values were evaluated to determine how stable the ligand–protein complexes are, and root-mean square fluctuation (RMSF) was used to assess individual residue fluctuations [89].

The binding free energy was calculated for frames generated from MD simulations for compounds 1, 5, 9, and the control. Schrödinger’s MM-GBSA continuum solvation method and the thermal_mmgbsa script were utilized [90]. A step size of five was used, and average ΔGbind values were reported. These calculations provide more thorough and accurate insights into the binding energetics for these compounds and can help rank them based on predicted affinities [91,92].

To predict the ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties as well as drug-likeness characteristics for the selected putative agonists used as modeling input, we utilized the pkCSM web server (http://biosig.unimelb.edu.au/pkcsm/prediction↗ (accessed on 13 December 2024)) [92]. Eight molecular descriptors were generated based on the ADMET descriptors identified in the potential hits. Additionally, Swiss ADME (www.swissadme.ch/↗ (accessed on 13 December 2024)) was employed for the in silico calculations of physicochemical parameters, medicinal chemistry-like properties, and drug-likeness attributes [93,94].