Structure-Based Virtual Screening and Molecular Dynamics Simulation Assessments of Depsidones as Possible Selective Cannabinoid Receptor Type 2 Agonists

The discovery of natural drug metabolites is a leading contributor to fulfilling the sustainable development goal of finding solutions to global health challenges. Depsidones are a class of polyketides that have been separated from lichens, fungi, sponges, and plants and possess various bioactivities, including cytotoxic, antimicrobial, antimalarial, antituberculosis, acetylcholinesterase and α-glucosidase inhibition, and anti-inflammatory effects. Endocannabinoid receptors (CB1 and CB2) are G-protein-coupled receptors (GPCRs), and their activation mediates many physiological processes. CB1 is the dominant subtype in the central nervous system, while CB2 is mainly expressed in the immune system. The two receptors exhibit high heterogeneity, making developing selective ligands a great challenge. Attempts to develop CB2 selective agonists for treating inflammatory diseases and neuropathic pain have not been successful due to the high homology of the binding sites of the CB receptors. In this work, 235 depsidones from various sources were investigated for the possibility of identifying CB2-selective agonists by performing multiple docking studies, including induced fit docking and Prime/molecular mechanics–generalized Born surface area (MM-GBSA) calculations to predict the binding mode and free energy. Simplicildone J (10), lobaric acid (110), mollicellin Q (101), garcinisidone E (215), mollicellin P (100), paucinervin Q (149), and boremexin C (161) had the highest binding scores (−12.134 kcal/mol, −11.944 kcal/mol, −11.479 kcal/mol, −11.394 kcal/mol, −11.322 kcal/mol, −11.305 kcal/mol, and −11.254 kcal/mol, respectively) when screened against the CB2 receptor (PDB ID: 6KPF). The molecular dynamic simulation was performed on the compounds with the highest binding scores. The computational outcomes show that garcinisidone E (215) and paucinervin Q (149) could be substantial candidates for CB2 receptor activation and warrant further in vivo and in vitro investigations.

In our continual interest to explore new bioactivities of the reported natural metabolites, 235 were depsidone derivatives, and 217 of them were reported to have naturally originated mainly from lichens (Figure S1-S9) inn addition, 18 nornidulin semisynthetic derivatives ( Figure S10). They were investigated for their possible CB2R agonistic potential utilizing molecular docking and molecular dynamics simulations (Tables S1 and S2), with the hope that the findings could highlight new possible CB2R ligands that encourage further in vitro and in vivo research to prove this effect. Remarkably, no previous investigations have been published on the possible CB2R-modulating activity of this class of metabolites and its related derivatives.  Lichens are composite symbiotic organisms that consist of a fungus with one or more photosynthetic partners (e.g., cyanobacterium, green alga, or both) [13]. Lichens are encountered in various ecosystems, including drought, icebergs, hot deserts, rocky coasts, toxic heaps, continuous light, and prolonged darkness [14]. They can produce varied structures and unique metabolites, such as polyketides, and in particular, depsones, depsidones, dibenzofurans, depsides, and chromones [14,15].
In our continual interest to explore new bioactivities of the reported natural metabolites, 235 were depsidone derivatives, and 217 of them were reported to have naturally originated mainly from lichens (Figure S1-S9) inn addition, 18 nornidulin semisynthetic derivatives ( Figure S10). They were investigated for their possible CB2R agonistic potential utilizing molecular docking and molecular dynamics simulations (Tables S1 and S2), with the hope that the findings could highlight new possible CB2R ligands that encourage further in vitro and in vivo research to prove this effect. Remarkably, no previous investigations have been published on the possible CB2R-modulating activity of this class of metabolites and its related derivatives. Depsidones have demonstrated diverse bioactivities, including anti-proliferative, antioxidant, antimalarial, antibacterial, cytotoxic, antihypertensive, anti-inflammatory, antifungal, and aromatase, protein tyrosine phosphatase, protein kinase, and HIV-1 integrase inhibitory properties [15][16][17].
In our continual interest to explore new bioactivities of the reported natural metabolites, 235 were depsidone derivatives, and 217 of them were reported to have naturally originated mainly from lichens (Figures S1-S9) inn addition, 18 nornidulin semisynthetic derivatives ( Figure S10). They were investigated for their possible CB2R agonistic potential utilizing molecular docking and molecular dynamics simulations (Tables S1 and S2), with the hope that the findings could highlight new possible CB2R ligands that encourage further in vitro and in vivo research to prove this effect. Remarkably, no previous investigations have been published on the possible CB2R-modulating activity of this class of metabolites and its related derivatives.
The depsidone derivatives were screened by generating a grid from the CB2 receptor (PDB ID: 6KPF) and docking the compounds into the grid, generating different docking scores, including the XP GScore and induced fit docking score. In addition, MM-GBSA was used to calculate the binding free energy between the protein and the ligand. Compounds with the highest docking scores were then subjected to a molecular dynamic (MD) simulation.

Protein and Ligand Preparation
The low-energy ionization and tautomeric states of 3D structures were generated in 2D using LigPrep. The generated three-dimensional structures were docked into the CB2 receptor (PDB-ID: 6KPF) [18] prepared via the protein preparation wizard tool, where the H bonds were optimized, and the geometry minimized. Missing hydrogen atoms and the correct ionization states were added to confirm the correct formal charges and force field treatment assignment.
The binding interactions consisted of a pi-pi stacking with Phe183 and hydrogen bonds with Ser285 at a distance of 1.76 Å and with Lys278 at a distance of 1.98 Å. Simplicildone J (10) had different interactions with the amino acid residues. Figure 5 shows the 2D and 3D views of these interactions, which included an H bond through its hydroxyl group with Leu182 at a distance of 1.94 Å, as well as ionic binding interactions with Ile110.
The binding interactions of mollicellin Q (101) are shown in Figure 6, and included a pi-pi stacking at a distance of 4.92 Å with Phe183, a hydrogen bond with Thr114 at a distance of 2.15 Å, and an ionic bond with His95.
Lobaric acid (110) also interacted through its aromatic ring, forming a pi-pi stacking interaction with Phe87 at a distance of 5.47 Å, as well as a hydrogen bond with Lys 109 through its hydroxyl group and two ionic bonds with His95 and Phe117 ( Figure 7). Paucinervin Q (149) formed a hydrogen bond with Ser285; it also formed several pi-pi stacking interactions with Phe183, at a distance of 5.13 Å, and with Phe87, Phe94, and His95 ( Figure 8). Figure 9 illustrates the binding interactions of boremexin C (161) with the residues, including a pi-pi stacking with Phe87, at a distance of 5.38 Å, and with Phe183 and a hydrogen bond with Leu182 at a distance of 1.78 Å. A pi-pi stacking interaction at a distance of 5.18 Å and with Phe183 at a distance of 5.10 Å, and hydrogen bonding interactions with Tyr25 at a distance of 1.67 Å and Leu182 with mollicellin P (100) are shown in Figure 10. The binding interactions of these compounds with the protein were examined. To validate these findings, the native agonist E3R was redocked alongside the depsidone ligands, and the docking positions were examined. Figure 4 represents the binding interactions of the native agonist E3R with the CB2 receptor in 3D and 2D views. The binding interactions consisted of a pi-pi stacking with Phe183 and hydrogen bonds with Ser285 at a distance of 1.76 Å and with Lys278 at a distance of 1.98 Å. Sim-plicildone J (10) had different interactions with the amino acid residues. Figure 5 shows the 2D and 3D views of these interactions, which included an H bond through its hy-droxyl group with Leu182 at a distance of 1.94 Å, as well as ionic binding interactions The binding interactions of mollicellin Q (101) are shown in Figure 6, and included a pi-pi stacking at a distance of 4.92 Å with Phe183, a hydrogen bond with Thr114 at a distance of 2.15 Å, and an ionic bond with His95. Lobaric acid (110) also interacted through its aromatic ring, forming a pi-pi stacking interaction with Phe87 at a distance of 5.47 Å, as well as a hydrogen bond with Lys 109 through its hydroxyl group and two ionic bonds with His95 and Phe117 (Figure 7). Paucinervin Q (149) formed a hydrogen bond with Ser285; it also formed several pi-pi stacking interactions with Phe183, at a distance of 5.13 Å, and with Phe87, Phe94, and     The Prime/MM-GBSA equation was used to calculate the binding free energies (DGbind) of each ligand based on the docking complex. The more negative values indicate a stronger binding (Table 1). Table 1. Docking results of in silico screening against CB2 receptor (PDB-ID:6KPF) and binding  The Prime/MM-GBSA equation was used to calculate the binding free energies (DGbind) of each ligand based on the docking complex. The more negative values indicate a stronger binding (Table 1).

Induced Fit Docking Analysis
In reality, proteins undergo backbone or sidechain movements upon ligand binding, making it difficult to assume a rigid receptor and model the real binding process [19]. Induced fit binding is one of the most challenging factors in drug design. The IFD scores of the best binding ligands are listed in Table 2. According to the obtained IFD scores, all compounds show comparable results to the prepared ligand, which indicates good binding. Paucinervin Q (149), garcinisidone E (215), and lobaric acid (110) had better IFD scores than the prepared ligand, which is an indication of a better interaction with the CB2 receptor.

QM/MM (Quantum Mechanics/Molecular Mechanics) Analysis
After induced fit docking, the best position was chosen and introduced to the QM/MM computation. The QM method was used for the protein's active site, while the rest of the proteins were treated by the MM method. This offers the advantage of obtaining accurate results by avoiding the demands of the computational work for calculating the QM for a large number of atoms [20]. The quantum mechanics method offers the advantage of accurately calculating the internal energies and HOMO/LUMO values [19]. The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) are located at the outer boundaries of the molecule's electrons. The ionization potential is related to the HOMO, while the electron affinity is linked to the LUMO [21]. The energy gap between the HOMO and LUMO in the native agonist was the lowest, indicating higher chemical reactivity [21]. As shown in Table 2, the QM/MM energy of the native agonist was higher than the tested depsidone derivatives, which indicates that they had a higher potential than the native agonist.

Molecular Dynamic Simulation
Molecular dynamic simulations utilize Newtonian physics in order to study atomic motions and investigate the behavior of molecules in action [22]. A variety of biomolecular processes can be captured during a molecular dynamic simulation, including ligand binding, protein folding, and conformational changes [23]. The simulation is carried out under precisely controlled conditions. Molecular dynamic simulations were performed on the native agonist and the best binding agonists to validate the docking scores. The MD outcome of the native agonist was used as a control. The root mean square deviation (RMSD) may be utilized to study the conformational stability of a structure during the simulation by measuring the average change in the displacement of atoms concerning a reference [24]. The average change in the displacement of boremexin C (161) ( Figure 11B) and garcinisidone E (215) ( Figure 11C) concerning a reference agonist (PDB: E3R) ( Figure 11A) was measured using the RMSD to predict the stability of the complex. The RMSD progression of the CB2 protein (left Y-axis) aligned on the reference frame backbone and the RMSD of the ligand (right Y-axis) are both shown in the plot. The RMSD plot of the native agonist ( Figure 11A) reveals that there is high fluctuation, more than 1-3 Å, which indicates that the CB2 receptor and native agonist complex are unstable and there was large conformational change in the protein during the simulation. However, the RMSD plots of 161 ( Figure 11B) and 215 ( Figure 11C) show that both complexes were stable and the slight fluctuations that occurred at the time of the simulation were within the acceptable range (1-3 Å).  Figure 12A represents a scheme of the binding interactions of the native agonist (PDB ID: E3R) and the amino acid residues of the CB2 receptor (PDB ID: 6KPF). Throughout the simulation, docked positions were maintained. The native agonist interacted through a hydrogen bond with Ser285, which was maintained for over 90% of the simulation time (100.00 ns), as well as a pi-pi stacking interaction with Phe87 (90%). Other interac-  Figure 12A represents a scheme of the binding interactions of the native agonist (PDB ID: E3R) and the amino acid residues of the CB2 receptor (PDB ID: 6KPF). Throughout the simulation, docked positions were maintained. The native agonist interacted through a hydrogen bond with Ser285, which was maintained for over 90% of the simulation time (100.00 ns), as well as a pi-pi stacking interaction with Phe87 (90%). Other interactions included a pi-pi stacking interaction with Phe183 (77%) and a hydrogen bond with Leu182 (40%). The interactions are also presented as stacked bars and categorized into three types: hydrogen bonds, hydrophobic interactions, and water bridges ( Figure 12B). The major contacts between the ligand and the CB2 receptor included hydrophobic interactions with Phe87 and Phe183 for over 80% of the trajectory time. The protein-ligand interactions also included hydrogen bonding interactions with Ser285 and Leu182, which were maintained for over 90% of the simulation. The binding interactions of boremexin C (161) with the CB2 receptor are presented in Figure 13A and include a hydrogen bond between the carbonyl oxygen and His96, which lasted for 84% of the simulation time. In addition, the hydroxyl formed two water-bridged hydrogen bonds with Glu181 and Asp24, which lasted for 61% and 33% of the trajectory time, respectively. Figure 13B displays the residue contacts of 161 with the CB2 receptor in the form of stacked bars. As shown in Figure 13A, the hydrogen bond interaction with His95 was maintained for more than 80% of the trajectory time. Boremexin C (161) interacted through both a hydrogen bond and a water-bridged interaction with Glu18. Several hydrophobic interactions were formed with various residues, including Phe91, Phe94, and Phe183. Figure 14A is a detailed 2D schematic representation of garcinisidone E (215) interactions with CB2 receptor residues that occurred for more than 30.0% of the simulation time in the selected trajectory (0.00 through 100.00 nsec). Molecular interactions included two pi-pi stacking interactions with Phe183 and Phe87 for almost 45% and 31% of the simulation time, respectively. The scheme also shows an intramolecular H bond between the hydroxyl group and the carbonyl oxygen, which increased the rigidity of the compound and occurred for 69% of the simulation time. The contacts of boremexin C (161) and the CB2 receptor are illustrated as stacked bars, which are normalized over the course of the trajectory and were classified into three main types: hydrophobic, hydrogen bonds, and water-bridged interactions ( Figure 14B). These contacts consisted mainly of hydrophobic interactions and a few hydrogen bonds. The hydrophobic interactions showed different degrees of maintenance throughout the simulation: Phe281 (80%), Phe183 (70%), and Phe117 (50%). The binding interactions of boremexin C (161) with the CB2 receptor are presented in Figure 13A and include a hydrogen bond between the carbonyl oxygen and His96, which lasted for 84% of the simulation time. In addition, the hydroxyl formed two water-bridged hydrogen bonds with Glu181 and Asp24, which lasted for 61% and 33% of the trajectory time, respectively. Figure 13B displays the residue contacts of 161 with the CB2 receptor in the form of stacked bars. As shown in Figure 13A, the hydrogen bond interaction with His95 was maintained for more than 80% of the trajectory time. Boremexin C (161) interacted through both a hydrogen bond and a water-bridged interaction with Glu18. Several hydrophobic interactions were formed with various residues, including Phe91, Phe94, and Phe183. Figure 14A is a detailed 2D schematic representation of garcinisidone E (215) interactions with CB2 receptor residues that occurred for more than 30.0% of the simulation time in the selected trajectory (0.00 through 100.00 nsec). Molecular interactions included two pi-pi stacking interactions with Phe183 and Phe87 for almost 45% and 31% of the simulation time, respectively. The scheme also shows an intramolecular H bond between the hydroxyl group and the carbonyl oxygen, which increased the rigidity of the compound and occurred for 69% of the simulation time. The contacts of boremexin C (161) and the CB2 receptor are illustrated as stacked bars, which are normalized over the course of the trajectory and were classified into three main types: hydrophobic, hydrogen bonds, and water-bridged interactions ( Figure 14B). These contacts consisted mainly of hydrophobic interactions and a few hydrogen bonds. The hydrophobic interactions showed different degrees of maintenance throughout the simulation: Phe281 (80%), Phe183 (70%), and Phe117 (50%).

ADMET Properties
The ADMET properties, which include the absorption, distribution, metabolism, and excretion of the depsidone derivatives were examined using QikProp of the Schrodinger module and are listed in Table 3. This enables the prediction of many physiochemical properties and biological functions, which could help predict the use-

ADMET Properties
The ADMET properties, which include the absorption, distribution, metabolism, and excretion of the depsidone derivatives were examined using QikProp of the Schrodinger module and are listed in Table 3. This enables the prediction of many physiochemical properties and biological functions, which could help predict the use-

ADMET Properties
The ADMET properties, which include the absorption, distribution, metabolism, and excretion of the depsidone derivatives were examined using QikProp of the Schrodinger module and are listed in Table 3. This enables the prediction of many physiochemical properties and biological functions, which could help predict the usefulness of the drug and reduce failures throughout the discovery and development process. In addition to predicting the physiochemical properties of the compounds, other factors were also examined, including the molecular weight, metabolism, dipole moment, hydrogen bond donors and acceptors, CNS activity, human oral absorption, and the partition coefficient. It was found that most of the descriptors were within range except for a few, such as poor absorption and high metabolism.

Preparation of Protein
The CB2 receptor X-ray crystal structure (PDB ID:6KPF) [25] was downloaded from the protein data bank. Using the protein preparation wizard in Maestro Schrödinger [18], the protein was prepared by adding the missing hydrogens to the residues, correcting the metal ionization states, and removing water molecules beyond 5 Å. In addition, the correct charges were assigned, and the protein underwent restrained minimization using an OPLS4 force field.

Ligand Preparation
Schrödinger's LigPrep tool [26] was used for docking 235 depsidone derivates, and the OPSL3 forefield was utilized for the energy-minimized 3-dimensional structures. Subsequently, Epik generated all the applicable ionization states and tautomeric forms at a pH of 7.0 ± 0.2 and added hydrogens. Additionally, PROPKA-optimized H bonds and any water molecules beyond 3 Å away from the HET groups were taken away.

Grid Generation and Molecular Docking
The ligand in the co-crystallized structure of 6KPF was selected to define the grid box, and the binding site was defined using Glide's Receptor-Grid Generation tool [27]. The grid box had a length of 10 Å in each of the X, Y, and Z dimensions. The assigned grid box was used for docking the prepared ligands using the Ligand Docking Tool in the Schrödinger suite. The chosen protocol was the extra-precision (XP) protocol. The partial charge cut-off was 0.25, and the VdW radii scaling factor was set to 1.0. The co-crystallized ligand was redocked alongside depsidone derivatives to validate the docking method. The other settings were kept at their default. The G score (ranks the compounds), Emodel (ranks the conformers), and XP GScore were calculated to assess the docking scores. In order to approximate the ligand binding free energy and rank the positions of different ligands, the GlideScore was calculated. Glide uses Emodel to select the best position of the ligand and then uses the GlideScore to rank the best positions against each other. The XP Glide was used to rank the ability of the ligands to bind to a specific conformation of the protein receptor and is presented by the following equation: where E is the energy (calculated for each of the following descriptors); E coul is the Coulomb energy, E VdW is the van der Waal force, E bind is the energy that favors binding, E penalty is the penalty that disfavors binding, E hyd_enclosure is the hydrophobic enclosure, E hb_nn_motif is a special neutral-neutral hydrogen bond motif, E hb_cc_motif is a special charged-charged hydrogen bond motif, E PI is a pi-cation interaction, E hb_pair is a hydrogen bond pair, E phobic_pair is a lipophilic pair, and E desolv is the desolvation energy [28].
Prime was used to calculate the molecular mechanics-generalized born surface area (MM-GBSA) [27] for re-scoring the docked positions. The binding free energy was calculated using the following equation: where (Complex) is the protein-ligand complex, (Eprotein) is the free protein, and (Eligand) is the free ligands.

Induced Fit Docking
Induced fit docking (IFD) protocol was performed in order to efficiently model the flexibility of the receptor and the ligand using the induced fit docking function in Maestro v9.1. (Schrodinger, LLC) [29,30]. The complex was used to generate the centroid of the residues by picking up the ligand from the protein. The ligands were docked into the protein using Glide. The sidechains were trimmed automatically based on their B-factor, with receptor and ligand van der Waals forces of 0.50 each. Twenty positions were generated and used to test the protein's plasticity. Residues were refined within 5.0 Å of the ligand positions and subjected to prime sidechain prediction and energy minimization, and then the sidechains were optimized. Ligand structures and conformations were induced to fit each position of the receptor. Then, redocking using Glide SP was performed on the best structures. Lastly, the ligand was docked into the induced fit receptor, and an IFD score for each position was generated.

Molecular Dynamic Simulation (MD)
Desmond software in the Schrödinger package [31] was used to run the molecular dynamic simulation of the compounds with the highest docking scores. A solvated system was prepared by immersing the protein-ligand complex into an orthorhombic water box, which was set to extend 10 Å beyond the atoms in the complex. The system was neutralized by adding Na and Cl counter ions. The simulation time was set at 100 ns, and the simulation was performed at a constant temperature of 300 K and constant pressure of 1.01325 bar.

Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations
The geometries used for performing the QM/MM calculations were obtained from the induced fit docking protocol detailed above. The QM/MM calculations were run on the protein (PDB-ID: 6KPF), and the native agonist (PDB-ID: E3R) and the depsidone derivatives had the highest docking scores. The Q site program [32] was used to calculate the QM/MM by including the ligands and the residues in the interaction. The QM region was created by picking up the free ligands and sidechains. The DFT-B3LP method, which is a measure of electron density, was chosen to perform the QM calculations, and the QM system had a charge of −1. The OPLS-2005 force field was selected to treat the MM region, which consisted of the rest of the system. After the QM/MM run was completed, the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital LUMO, and energy gaps were calculated.

Prediction of ADMET Properties
Predicting the ADMET properties of the compounds, which are absorption, distribution, metabolism, and elimination, provided some vital information related to the drugs/molecules. All compounds were assessed using the Maestro Schrodinger QikProp module to identify promising molecules according to their bioavailability properties and ADMET profiling [33].

Conclusions
Natural metabolites originating from various sources, e.g., animals, plants, and microorganisms display high structural complexity and diversity. Taking into consideration that a high-throughput in vivo/in vitro investigation for assessing the bioactivities of natural metabolites is costly and time-consuming, quite effective in silico tools could assist as rapid, auspicious, and cost-efficient strategies to predict the natural metabolites-target association before experimental validity. CB2R is a possible therapeutic target for treating various illnesses, such as neurodegenerative diseases (e.g., Huntington's, Alzheimer's, and Parkinson's diseases and multiple and amyotrophic lateral sclerosis), as well as inflammation, pain, myocardial infarction, irritable bowel syndrome, and various cancer types. For the identification of potential CB2R agonists, 235 depsidone derivatives were virtually screened for their CB2R-binding potential using multiple virtual studies: extra precision docking, induced fit docking, and Prime/MM-GBSA calculation, as well as MD simulations. These methods proved the CB2R agonistic activity of simplicildone J (10), lobaric acid (110), mollicellin Q (101), garcinisidone E (215), mollicellin P (100), paucinervin Q (149), and boremexin C (161) and stability during the MD simulation studies. The variation in the binding affinity of the compounds based on the methods of calculation-for example, induced fit docking and QM/MM-was attributed to the difference in parameters and conditions used for each calculation; however, there were no significant differences amongthe results. The findings of this work reveal that 149 and 215 demonstrated the highest binding affinity to the CB2R. Therefore, they could be lead candidates for CB2 receptor activation, which warrants further in vivo and in vitro investigations.

Data Availability Statement:
The data presented in this study are available in the article.