Pharmacophore based virtual screening were performed using the Pharmit web interface ( http://pharmit.csb.pitt.edu/), which has a number of built-in databases that provide access to comprehensive data; Molprot (4,742,020 molecules), ChEMBL34 (2,264,112 molecules), ZINC (13,127,550 molecules), ChemDiv (1,456,120 molecules), ChemSpace (50,181,678 molecules), Enamine (4,117,328 molecules), MCULE (39,843,637 molecule), MCULE-ULTIMATE (126,471,502 molecules), NCI Open Chemical Repository (52,237 molecules), LabNetwork (1,794,286 molecules), and PubChem (103,302,052 molecules). The pharmacophore model was constructed using selected PDB code 5XRA from the RCSB Protein Data Bank ( http://www.rcsb.org/structure/5xra), with agonist AM11542. The model utilized a pharmacophore framework built on five features, adhering to the default parameters of the Pharmit server. The Pharmit filters were applied based on the Lipinski rule of five and Veber’s rule to refine the screening process and identify the most selective CB1 agonist. Figure 2 illustrates the multi-step virtual screening process used in this work.

The ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) filter is one of the phases of the drug discovery and development process. ^{13– 15} This allows researchers to determine potentially druggable drug-like properties and toxicology early on to assess which compounds will be effective drugs with safety potential. Thus, relative to virtual screening, such a filter can detect adverse properties that would fail at much later stages, such as insufficient absorption, excessive toxicity, and/or suboptimal bioavailability. Two of the trusted websites relied upon to evaluate such properties based upon lipophilicity, logP, solubility, etc., are SwissADME ( http://www.swissadme.ch) and pKCSM ( https://biosig.lab.uq.edu.au/pkcsm/prediction). Therefore, when these are filtered out early on, time and costs will be effectively saved down the line.

To estimate the binding affinities of the 433 highest-ranked compounds to CB1, molecular docking was carried out using AutoDock Vina integrated within the PyRx platform. The crystal structure of CB1 (PDB ID: 5XRA), with a resolution of 2.8 Å, ^{8, 16, 17} was retrieved from the Protein Data Bank. Before docking, water molecules and the native ligand were removed from the protein, followed by the addition of polar hydrogen atoms and Kollman charges. The prepared protein structure was then energy-minimized using UCSF Chimera and saved in PDBQT format.

Prior to docking, all 433 compounds underwent pre-optimization utilizing the Universal Force Field (UFF) in combination with the conjugate gradient algorithm. ¹⁸ Optimization parameters were set to 2000 total steps, with updates occurring every step, and the process was programmed to terminate once the energy difference fell below 0.01 kcal/mol. ¹⁹ After optimization, the compounds were converted to PDBQT format and docked at specific binding site coordinates (x = -42.052, y = -164.338, z = 306.631). Molecules demonstrating the lowest binding energies and root-mean-square deviation (RMSD) values below 2 Å were selected for subsequent analyses.

Molecular dynamics (MD) simulations were carried out on the top-ranked docking conformations using the GROMACS 2024.4 software package, employing the CHARMM27 all-atom force field. ²⁰ The topology files for the CB1 receptor were generated using the pdb2gmx tool, ²¹ while the ligand topologies were prepared through the CHARMM General Force Field (CGenFF) using the Param-Chem server. ^{22– 24} Each receptor-ligand complex was embedded in a triclinic simulation box and solvated with TIP3P water molecules, maintaining a minimum distance of 1.0 nm from the box boundaries. To neutralize the systems, Na ⁺ and Cl ^- ions were added accordingly.

Prior to the production runs, energy minimization was performed using the steepest descent algorithm for 50,000 steps to eliminate steric clashes and ensure system stability. Subsequently, two equilibration phases were conducted: first under the NVT ensemble to stabilize temperature, followed by the NPT ensemble to equilibrate pressure and density. Input files for both equilibration and production stages were customized to adjust parameters such as trajectory-saving intervals, energy monitoring, and other essential simulation settings. Finally, each ligand-receptor system underwent 100 ns of production MD simulation at a constant temperature of 300 K, pressure of 1 bar, and a time step of 2 femtoseconds.

The binding free energy was calculated via MM/GBSA (Molecular Mechanics Generalized Born Surface Area). ²⁵ This approach utilizes force fields of molecular mechanics and an implicit solvent approach to conclude stability and binding tendencies of molecular complexes. It is particularly useful for classifying ligands and understanding overall energetic contributions of molecular recognition. The binding free energy ΔGbind of a complex can be calculated as follows: ΔG bind = ΔE MM + ΔG solv – TΔS (Eq.1)

Where: ΔE MM = E MM complex − ( E MM receptor + E MM ligand ) (Eq.2) ΔG solv = G solv complex − ( G solv receptor + G solv ligand ) (Eq.3) G solv = G GB + G SA

Where TΔS, ΔE _MM and ΔG _sol are the conformational entropy upon binding, the changes of the gas phase molecular mechanism (MM) energy and the solvation free energy, respectively.

The pharmacophore model of the CB1 receptor agonist was generated based on the AM11542 agonist. The results are shown in Figure 3. A five-point model with 1 aromatic, 2 hydrogen bond acceptors and 2 hydrophobic regions was generated using the Pharmit web server. ²⁶

Subsequently, virtual screening was performed across 11 databases, identifying 433 compounds containing molecular groups that matched the pharmacophore mode.

Pharmacophore-based virtual screening was performed using information from the previous pharmacophore model. Approximately 300 million compounds from the Pharmit database were filtered. The criteria for filtering the library were as follows: molecular weight limits were set at 250-500 Da, the maximum number of hydrogen bond donors was less than 4, the maximum number of hydrogen bond acceptors was less than 9, the maximum number of rotatable bonds was less than 9, Log P was between 2 and 5, and polar surface area was less than 140 Å. The results are computed and classified according to different criteria such as energy minimization. The 433 top-ranked compounds are presented in Table 1.

Next, molecular docking analysis was performed using AutoDock Vina, implemented in the PyRx software. Compounds with a binding affinity of less than - 9.00 kcal/mol were selected for further analysis, while the others were excluded. Additionally, the predicted binding modes were required to have a root-mean-square deviation (RMSD) of less than 2Å when superimposed onto the native agonist (AM11542). The top 61 selected compounds are listed in Table 2, while the remaining candidates are provided in Table S1 (see supplementary data).

The ranked compounds from docking analysis were evaluated for potential toxicity, including an AMES toxicity test, an acute oral toxicity test in rats (LD ₅₀), a skin sensitization test and a maximum tolerated dose analysis. Highly toxic compounds are not considered in further studies. The top selected compounds are presented in Table 3. The other are presented in the supplementary data ( Table S2).

As described in Table 3, all selected compounds exhibited no Skin Sensitisation and no AMES toxicity. Additionally, the LD ₅₀ values for oral rat toxicity, ranging from 1.439 to 3.354 mol/kg, suggest moderate acute toxicity levels, consistent with safety margins suitable for therapeutic use. Moreover, the maximum tolerated dose in humans’ ranges from -1.12 to 1.059 Log mg/kg/day.

Physicochemical properties were evaluated using Lipinski’s rule of five, ²⁷ Ghose’s rule, ²⁸ Veber’s rule, ²⁹ Egan’s rule ³⁰ and Muegge’s rule, ³¹ with results detailed in Table 4. Based on these guidelines, it is suggested that for a compound to be effectively absorbed and administered orally, it must meet specific physicochemical parameters. These criteria serve as essential benchmarks for assessing the compound’s potential for bioavailability and oral absorption. Compounds with two or more violations are not considered in the further analysis (see Table S3 in supplementary data).

To predict the pharmacokinetic properties of absorption, distribution, metabolism and excretion (ADME), pkCSM web server was used to calculate the following parameters: water solubility (log mol/L), Caco-2 cell permeability, human intestinal absorption (HIA), blood-brain barrier (BBB) permeability, central nervous system (CNS) permeability and total clearance.

Water solubility (LogS) indicates the solubility of a compound in water at 25°C. Generally, water-soluble drugs are more readily absorbed than lipid-soluble ones. in vitro Caco-2 cell permeability is a crucial measure of drug absorption, with a compound considered to have high Caco-2 permeability if its value surpasses 8 x 10 ⁻⁶ cm/s. Within the pkCSM model, high Caco-2 permeability aligns with predicted values exceeding 0.90. The intestine typically serves as the primary site for drug absorption from orally administered solutions. A compound with absorbance below 30% is deemed poorly absorbed. An established gauge of blood-brain barrier (BBB) penetration is the log BB ratio, which reflects drug molecule concentrations in the brain and blood. Compounds with log BB > 0.3 exhibit high BBB permeability, while those with log BB < 1.0 show limited BBB distribution. Furthermore, central nervous system (CNS) permeability is a vital parameter for assessing the blood-brain permeability of a drug candidate, expressed as LogPS. Compounds with LogPS <-3 are considered incapable of penetrating the CNS.

Based on these ADME properties, 18 compounds were selected for further analysis. The results for these compounds are detailed in Table 5, while the other parameters are shown in Table S4 in the supplementary data.

Docking validation was performed using AutoDock 4.2. The results are presented in Table 6.

To identify the top CB1 agonist several criteria were taken into account including energy values (e.g. binding energy or Ki values), interactions with key amino acids in the CB1 binding site, number of hydrogen bonds and distances between hydrogen bonds. These interactions are essential for designing effective and selective CB1 agonists. ¹⁰ For instance, π-π interactions with aromatic residues (Phe296, Phe170, Phe174, Trp 279) and hydrogen bonds with Ser383 or Thr197 stabilize the agonist-receptor complex, while residues such as Phe379 and Asp366 play a key role in receptor activation through electrostatic interactions. ¹⁰ Based on these criteria three compounds were identified as CB1 receptor agonists: ZINC64438506, ZINC64438485 and PUBChem157251136.

ZINC64438506: The docking analysis of CB1 receptor and ZINC64438506 selected agonist is shown in Table 7 and Figure 4. The ZINC64438506 agonist was fixed in the CB1 binding pocket (cavity size = 2963 Å) through various type of interactions, including hydrogen bonds with key amino-acid residues SER A: 173 and SER A: 383, hydrophobic interactions with residues PHE A: 108, PHE A: 177, PHE A: 189, LEU A: 193, THR A: 197, PHE A: 268, TYR A: 275, LEU A: 276, TRP A: 279 and PHE A: 379 and π-π interaction with PHE A: 268. These interactions likely contribute to the compound’s low binding affinity (-13.11 Kcal/mol) and low inhibition constant (Ki = 0.244 nM) ( Table 6). The strong binding affinity may be attributed to the presence of short hydrogen bonds with SER A: 173 (3.34 Å) and SER A: 383 (51.79 Å).

ZINC64438485: The docking analysis of CB1 receptor and ZINC64438485 selected agonist is shown in Table 7 and Figure 5. The ZINC64438485 agonist was fixed in the CB1 binding pocket (cavity size = 2963 Å) through various type of interactions, including hydrogen bonds with key amino-acid residues THR A: 197 and SER A: 383, hydrophobic interactions with residues PHE A: 174, PHE A: 177, PHE A: 189, LEU A: 193, VAL A: 196 and ILE A: 271 and π-π interaction with PHE A: 174, HIS A: 178, PHE A: 268 and TRP A: 279. These interactions likely contribute to the compound’s low binding affinity (-13.07 Kcal/mol) and low inhibition constant (Ki = 0.262 nM) ( Table 6). The strong binding affinity may be attributed to the presence of short hydrogen bonds with THR A: 197 (1.94 Å) and SER A: 383 (3.08 Å).

PUBChem-157251136: The docking analysis of CB1 receptor and PUBChem-157251136 selected agonist is shown in Table 7 and Figure 6. The PUBChem157251136 agonist was fixed in the CB1 binding pocket (cavity size = 2963 Å) through various type of interactions, including hydrogen bonds with key amino-acid residues SER A: 173, HIS A: 178 and SER A: 383, hydrophobic interactions with residues PHE A: 174, PHE A: 177, LEU A: 193, THR A: 197, PHE A: 268, LEU A: 276, TRP A: 279 and PHE A: 379 and π-π interaction with PHE A: 268 and TRP A: 279. These molecular interactions likely contribute to the compound’s low binding affinity (-11.84 Kcal/mol) and low inhibition constant (Ki = 2.09 nM) ( Table 6). The strong binding affinity may be attributed to the presence of three hydrogen bonds with SER A: 173 (3.23 Å), HIS A: 178 (1.94 Å) and SER A: 383 (2.04Å).

While molecular docking analyses are enough to understand possible bindings based on ligand positioning and receptor-ligand interactions, it should be noted that such methodologies evaluate the flexibility of the ligand only while keeping the protein in a rigid form. Thus, in order to evaluate the best-docked candidates for binding pose stability and dynamics of protein conformation, molecular dynamics (MD) simulations were performed over 100 ns. The results of the MD simulations, including root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and Ligand-Receptor Interaction Plot analyses, can be found in Figures 7– 10. These studies provide insights into the dynamic behavior and stability of the protein–ligand complexes over time.

7.1 Root Mean Square Deviation (RMSD)

RMSD analysis was carried out for the protein backbone to have an idea of each protein-ligand complex’s structural stability during the simulation. ^{32, 33} The results are shown in Figure 7. As shown in Figure 7, mean RMSD values for the CB1_ZINC64438506, CB1_ZINC64438485 and CB1_PUBChem157251136 complex are 0.29 nm, 0.45 nm and 0.64 nm, respectively. The global RMSD values for CB1_ZINC64438506, CB1_ZINC64438485 are small which indicates that these ligands remain stable during the simulation of 100 ns and remained in the binding pocket of the CB1 receptor. The stability of the ligands can be attributed to several strong hydrogen bonds being mediated between these agonists and some of the key amino acids positioned within the binding pocket of the CB1 receptor.

7.2 Root mean square fluctuation (RMSF)

Root mean square fluctuation (RMSF) was used to determine the rigid and flexible regions of the CB1 receptor over the 100 ns of MD simulations. ¹⁴ RMSF has been term definition, instead of just RMSD values, so that the maximum range of motion of a bound ligand can be seen. RMSF is defined as a standard measure of deviation of a molecule from its initial position. ³⁴ Molecules and residues should not present a high value, which indicates a flexibility, and those that appears low value has a greater rigidity. The RMSF plot for all complexes ( Figure 8) indicates that most residues located in the CB1 receptor has a low RMSF value, indicating they were rigid and retained stability over the entire 100 ns MD simulation.

7.3 Radius of gyration (Rg)

The radius of gyration (Rg) is another important metric measured during MD simulations to determine spatial characteristics of a protein-ligand complex. The radius of gyration is the root mean square distance of all atoms making up a given structure from a relative center common point (the center of mass) and signifies extension versus folding of the given structure. Therefore, a lower Rg would imply stability and a more tightly folded structure, while a higher Rg would imply extension and flexibility of possible structures. In this work, the Rg profiles remained stable over the course of 100 ns of MD simulation for the complexes CB1_ZINC64438506 (mean = 4.70 nm; Figure 9) and CB1_ZINC64438485 (mean = 4.72 nm), while the complex CB1_PUBChem157251136 decreased in Rg (mean = 2.81 nm).

7.4 Hydrogen bonds

One of the main factors influencing the affinity of a molecule for the protein binding pocket is its ability to form and maintain hydrogen bonds with the binding site residues. The stability of the selected agonist was assessed by analyzing the hydrogen bonds between the ligand and the protein. The results are presented in Figure 10.

Compounds ZINC64438506, ZINC64438485, and PUBChem157251136 form two, two and three hydrogen bonds with the CB1 binding site, respectively, indicating strong and specific interactions with the protein. These results are consistent with the docking results.

The binding free energy of all complexes was calculated to revalidate the binding affinity obtained from molecular docking analysis. The results are presented in Table 8. A ΔG _bind value below -7 Kcal/mol indicates strong binding, a value between -5 and -7 kcal/mol suggests moderate binding, and a value between -2 and -5 kcal/mol corresponds to weak binding. ³⁵ The results show that all the agonists exhibit exceptionally low binding free energies, indicating a remarkable affinity for the target protein.

The results also highlight the crucial role of Van der Waals and electrostatic interactions in mediating protein-ligand binding throughout the molecular dynamics simulations ( Figure 11).

Quinazoline-2,4(1H,3H)-dione derivatives represent a highly promising class of heterocyclic compounds with broad therapeutic potential, particularly as anticancer, ^{36– 38} antibacterial, ^{39, 40} antihypertensive, ⁴¹ phosphodiesterase (PDE) 4 inhibition, ⁴² 5-HT3A receptor antagonist, ⁴³ anti-inflammatory, ⁴⁴ and anan up-and-comingtiviral agents. ⁴⁵ Their fused benzopyrimidinedione scaffold provides an excellent pharmacophore for targeting key biological pathways ( Figure 12). According to the article’s findings, quinazoline-2,4(1H,3H)-dione derivatives may be the first of a new class of CB1 agonists.

To synthesise the substituted quinazoline-2,4-dione and the highest-ranking agonists identified by virtual screening, namely ZINC64438506, PUBChem157251136 and ZINC64438485, several plausible synthesis strategies were developed based on established methodologies. ⁴³ These routes use various precursors such as 2-aminobenzoic acid, ⁴⁴ methyl 2-nitrobenzoate, 2-iodobenzamides ⁴⁶ or 2H-benzo [d][1,3]oxazine-2,4(1H)-dione ⁴⁷ ( Figure 13). The choice of starting material is guided by both its availability and synthetic feasibility. Depending on the substrate, the desired compounds can be obtained by catalytic transformations, particularly transition metal-catalysed couplings, or by intramolecular cyclisation reactions that form the characteristic quinazoline skeleton. The synthetic flexibility offered by these precursors allows the reaction conditions to be adjusted to optimise yield and purity, making them suitable candidates for further pharmacological evaluation and development.

Translated with DeepL.com (free version)

To synthesise substituted quinazoline-2,4-dione and the top-ranked agonists identified by virtual screening - namely ZINC64438506, PUBChem157251136 and ZINC64438485 - several plausible synthetic strategies have been devised, based on established methodologies. ⁴⁶ These routes use various precursors such as 2-aminobenzoic acid, ⁴⁷ methyl 2-nitrobenzoate ⁴⁸ which undergoes catalytic hydrogenation under mild and green conditions (1 atm of H ₂ in ethanol at room temperature), ^{49, 50} 2-iodobenzamides ⁵¹ or 2H-benzo [d][1,3]oxazine-2,4(1H)-dione ⁵² ( Figure 13). The choice of starting material is guided by both availability and synthetic feasibility. Depending on the substrate, the desired compounds can be accessed by catalytic transformations, including transition metal-catalysed couplings, or by intramolecular cyclisation reactions that form the characteristic quinazoline scaffold. The synthetic flexibility offered by these precursors allows reaction conditions to be fine-tuned to optimize yield and purity, making them suitable candidates for further pharmacological evaluation and development.