Pharmacophore-based Virtual Screening From Phytocannabinoids for Anti-obesity Activity

Search for new pharmacological alternatives for obesity is based on the design and development of compounds that can aid in weight loss so that they can be used safely and effectively over a long period while maintaining their function. The endocannabinoid system is related to obesity by increasing orexigenic signals and reducing satiety signals. Cannabis sativa is a medicinal plant of polypharmaceutical potential that has been widely studied for various medicinal purposes. The in silico evaluation of their natural cannabinoids (also called phytocannabinoids) for anti-obesity purpose stems from the existence of synthetic cannabinoid compounds that have already presented this result, but which did not guarantee patient safety. In order to nd new molecules from C. sativa phytocannabinoids, with the potential to interact with the pharmacological target cannabinoid receptor 1, a pharmacophore-based virtual screening was performed, including the evaluation of physicochemical, pharmacokinetic, toxicological predictions and molecular docking. The results obtained from the ZINC 12 database pointed to Zinc 69 (ZINC33053402) and Zinc 70 (ZINC19084698) molecules as promising anti-obesity agents. Molecular Dynamics (MD) studies discloses that both complexes were stable by analyzing the RMSD (Root Mean Square Deviation) values, and the binding free energy values demonstrate that the selected structures can interact and inhibit their catalytic activity. obtained a new hydrogen now with Ser383, also important in antagonist This docking study showed relevant amino acids from the active site within 10 Å: Phe102, Met103, Ile105, Ile119, Ser123, Ile169, Phe170, Phe174, His178, Leu193, Val196, Thr197, Phe268, Leu276, Val279, Val356, Leu359, Met363, Ala380, Ser383, Met384, C386 e Leu387.


Introduction
Obesity is considered a global epidemic disease, is a public health problem that has been growing annually, according to the World Health Organization (WHO). It may lead to comorbidities such as diabetes mellitus, cardiovascular disease, stroke, hypertension and certain types of cancer (endometrial, breast, ovarian, prostate, liver, gallbladder, kidney and colon) [1,2].
Currently, the search for new targets for the treatment of obesity motivates research worldwide. An example of this was rimonabant, a cannabinoid antagonist that was marketed in Europe, but it was discontinued because presented several side effects on the central nervous system (CNS) [3].
The endocannabinoid system involves the presence of the following components: cannabinoid receptors (rCB), endogenous cannabinoids (endocannabinoids), and enzymes involved in endocannabinoid synthesis and degradation [4].
Cannabinoid receptors belong to the G protein superfamily and are distributed throughout the human body [5] rCB1 primarily regulate central nervous system (CNS) neurotransmission, but are also present peripherally in gastrointestinal organs, adipose tissue, myocardium, vascular endothelium, peripheral nerves, liver, pancreas, and skeletal muscle [6]. It has been studied as a target for various diseases such as pain, anxiety, multiple sclerosis, Huntington's and Parkinson's diseases, and even for obesity.
The rCB1 are known for having endogenous neuromodulatory ligands that belong to the family of lipids, which present different chemical structure than the exogenous tetrahydrocannabinol (THC) ligand.

Design and optimization of molecules
The design of the chemical structure of the 104 phytocannabinoid molecules presented by Pertwee [15] was performed in the ChemSketch v2007 2.1 program and their structures were geometrically optimized in the HyperChem v8.0 program, using the semi-empirical method RM1 (Recife Model 1), in order to get structures with minimal conformational energy.

Prediction of biological activity and pharmacophore derivation
Prediction of biological activity was performed on the PASS (Prediction of Activity Spectra for Substances) server [17] (http://www.way2drug.com/PASSOnline/) to identify those molecules that have CB1 receptor antagonist activity. From the 32 best results, then, the pharmacophoric pattern was derived in the PharmaGist server [18] (http://bioinfo3d.cs.tau.ac.il/PharmaGist/) to obtain the biological activity regions of the cluster [19].

Virtual screening
Performed on the ZINCPharmer web server (http://zincpharmer.csb.pitt.edu/) [20] which identi es molecules by the similarity of hydrophobic, aromatic, hydrogen donor and acceptor regions, positive and negative ions. Molecules selected from the ZINCpharmer web server have the same pharmacophoric regions previously predicted for phytocannabinoids. Compounds that combine with a well-de ned pharmacophore will serve as potential leaders for the possible discovery of rCB1 acting drugs through molecular modeling.

Physicochemical properties
Physicochemical properties were identi ed on the PubChem server (http://www.ncbi.nlm.nih.gov/pccompound) for phytocannabinoids and on the ZINC server (http://zinc.docking.org/) for the screened molecules on the ZINCpharmer server. These data allow the assessment, by Lipinski's Rule of Five [21], of the oral absorption of each molecule through its molecular weight, log P, and the number of hydrogen donors and acceptors.

Prediction of pharmacokinetic and toxicological properties of phytocannabinoids and Zinc molecules
To predict the pharmacokinetic properties of absorption, distribution, metabolism and excretion (ADME), the QikProp software [22] was used to calculate the parameters: stars, percent human oral absorption (% HOA), cellular permeability in Caco2 (pCaco-2) and MDCK (pMDCK), human serum albumin binding (plogKhsa), Van der Waals surface area (PSA), and blood-brain barrier permeability (pLogBB). Toxicological properties were determined by nephrotoxicity, hepatotoxicity, carcinogenicity and mutagenicity predictions in the DEREK NEXUS 2.1 software [23], using as a prediction lter for humans.

Molecular Docking Simulation
Docking analysis was performed using GOLD 5.7.1 (Genetic Optimization for Ligand Docking) software to predict the binding of molecules with better physicochemical, pharmacokinetic and toxicological pro les to deal with rCB1. The crystallographic structure of rCB1 complexed with an antagonist (taranabant) was obtained from the Protein Data Bank (PDB) server under the PDB ID: 5U09, with a resolution of 2.6 Å. Validation of binding at the binding site, that is, reproduction of the ligand conformation occurs by redocking calculating Root Mean-Square-Deviation (RMSD), which must have been less than 2 Å. The highest score anchors are those whose chemical structures have the most favorable binding free energy to bind to the binding site [24].
Protein preparation involved the addition of hydrogens, the extraction of water molecules, as well as the deletion of all ligands present in the structure, with ligand A09 being taranabant. Ligand A09 was selected in the binding site de nition step for the program to detect the binding cavity, restricting the selection of atoms to solvent and keeping the surface accessible. The phytocannabinoids and screening molecules evaluated (Zinc molecules) were pre-optimized in HyperChem under the RM1 semi-empirical method and inserted into the following binding site coordinates: x = 21.66, y = 3.21 and z = -9.06.

Molecular Dynamics Simulation
The MD simulation for systems was conducted using the AMBER 18 program [25], from the coordinates obtained in the molecular docking. The electrostatic potential of each of the compounds was calculated with the Gaussian 09 program [26], with the HF/6-31G* [27]. The PROPKA 3.0 [28] was used to compute the protonation states of the amino acid residues at neutral pH. The general Amber force elds GAFF [29] and MMFF99SB [30], were used for the parameterization of the ligands and the amino acid residues of the enzyme, respectively. The compounds chosen were solvated using a cubic water-box with the explicit solvation model TIP3P [31]. The distance of 12.0 Å between the cell wall and the solvated atoms of the system was used, and Clions were added to neutralize the overall charge of the complexes, whereas the MD minimization, heating and simulation stages were performed in the SANDER module [32] of AmberTools 18 was used in the MD calculations. Water molecules, hydrogens, ions, and the whole protein structure were minimized by four minimization steps in order to eliminate possible stereo spatial con icts between the atoms. Subsequently, the systems were gradually heated considering the NVT ensemble from zero (0) to 298 K for a period of 25 ps, with the protein atoms subjected to a restriction constant of 25 kcal/mol.Å 2 , the system was equilibrated by 250 ps considering the NPT ensemble (p=1.0 atm and T = 298 K). The Langevin algorithm [33], was used to regulate the system temperature, while the particle mesh Ewald (PME) method, was used to calculate long-range electrostatic interactions using a 9.0 Å cutting radius. In the course of the MD simulation, the bonds involving hydrogen atoms were restricted using the SHAKE algorithm [34], and the equations of motion were integrated every 2.0 fs using the Verlet algorithm [35]. After minimization and equilibrations of the system, 40 ns of MD was produced with an integration time of 2.0 fs and considering the NPT ensemble at 298 K and 1.0 atm of pressure for each complex. Finally, the stability of the systems was performed out by using the root mean square deviation (RMSD) in relation to the initial structure, calculated with the PTRAJ program, present in the AMBER 18 program, and the trajectories were visualized with the Chimera 1.14 program [36]. The calculation of the interaction energy between the ligand and the macromolecule was performed with the MMPBSA.py program, implemented in the AMBER 18 program.

Free binding energy calculations
To evaluate the a nity between the ligands and the binding free energy of the analyzed compounds was determined using the molecular mechanics-generalized Born surface area (MM-GBSA) and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methods. In both methods, the last 10 ns of MD simulations, which contains 1,200 frames were used to perform the binding free energy calculations. The snapshots were carefully extracted by using the CPPTRAJ module and submitted to the MMPBSA.py module [37] of AmberTools 18. The nal frames of trajectory were selected for binding free energy calculations due to the conformational stability of protein-ligand systems analyzed by the RMSD plots.
The binding free energy values were computed by using MM-PB(GB)/SA approaches [38] implemented on the MMPBSA.py module [39].

Chemical structure and biological activities prediction of phytocannabinoids
Cannabinoid molecules evaluated in this study were designed based on the work of Pertwee in the Handbook of Cannabis [15]. Of the 104 initial phytocannabinoids, only 41 showed a prediction of activity for rCB1 antagonism, and their chemical structures are presented in Figure 1. The PASS server has been used to predict the biological activity for the initial 104 molecules, in terms of P a (probability of being active) and P i (probability of being inactive).

Pharmacophoric pattern of phytocannabinoids with the prediction of antagonist activity on rCB1
Determination of regions on molecule structures that are most likely related to their biological activity, i.e. derivation of their pharmacophoric patterns, was performed on the PharmaGist server [19] with the 32 best-sorted molecules in the previous step. The alignment of such molecules presented the best score of 15.993 with the following points in the pharmacophoric pattern: 01 hydrogen-accepting region (HBA) (yellow sphere) and 04 hydrophobic regions (HPO) (green spheres) ( Figure 2).
Sharma and collaborators [40] de ned as a characteristic pharmacophore of peripheral-acting rCB1 antagonists a four-point model with the following regions: 1 HBA, 1 HPO and 2 aromatics (AR). In the same study, other pharmacophoric hypotheses were evaluated, including a ve-point model with: 1 HBA, 2 HPO and 2 AR. Wang and collaborators [41] also obtained two 5-point pharmacophoric patterns where the rst comprised 1 HBA, 2 AR and 2 HPO; while the second 1 HBA, 1 AR and 3 HPO. The information shows us that there is no de ned pharmacophoric pattern for cannabinoid antagonists.
The pharmacophoric pattern obtained in this study resembles Wang's second model, except for the aromatic region. However, it is noteworthy that the pharmacophoric patterns obtained in the literature, concerning cannabinoid antagonists, were built using individual molecules as a reference, and do not refer to the alignment of various antagonists. The proposition of a novel pharmacophore may prove useful in the search for safer cannabinoid antagonists for anti-obesity activity.

Virtual Screening
In order to perform a pharmacophore-based screening, we used our previously described pharmacophore (1 HBA and 4 HPO), by using the purchasable database from the ZINCPharmer web server and 1000 hits were obtained as results, with RMSD ranging from 0.40 to 0.88Å [20].
Subsequently, these hits were submitted to biological activity prediction on the PASS server. Results obtained showed that 78 molecules ( Figure 3) were favorable to cannabinoid antagonist activity, besides presenting indicatives of anti-obesity activity, like Rimonabant, as shown in Table 2.

Analysis of physicochemical properties
Physicochemical properties were evaluated according to Lipinski's Rule of Five (RO5), [21] which represents a fundamental guideline for getting to know the oral absorption of a drug candidate, regarding its permeability and solubility, as well as of what is necessary for pharmacological action. According to this, has a good prediction for oral absorption those molecules that present results with at most a single violation of the evaluated parameters/properties, as follows: molecular weight (g/mol) (MW) ≤ 500, log P ≤ 5; hydrogen bond donors (HD) ≤ 5 and hydrogen bond acceptors (HA) ≤ 10.
The data required in the analysis of 32 phytocannabinoids (the same selected for the pharmacophore derivation) was obtained from the Pubchem server. All of the 32 molecules evaluated at this stage, except molecule 14, showed results indicating good oral absorption (Table 3). Bold results indicate violations of RO5.
According to the results presented in Table 3, only phytocannabinoid 14 violated RO5, with molecular weight and log P above the evaluated parameters, which indicates that its physicochemical pro le does not favor its oral administration. The other molecules are in agreement with RO5 and also with the physicochemical characteristics of rimonabant [42] in which the only violation is related to log P (equal to 6.11).
Concerning the ability of a molecule to present good CNS permeability, the physicochemical property that is most related to this issue is the lipophilicity, which may be expressed by log P, i.e. the partition coe cient of a given molecule calculated in the system 1-octanol/water. With this in mind and taking into account, a generic analysis in the context of RO5, values of log P ≤ 4 indicates that a molecule may be inactive in the CNS, and log P > 4 indicates that may be active in the CNS. [27] Given this information, all phytocannabinoids presented in Table 3 have physicochemical potential to permeate the blood-brain barrier, an expected result for these metabolites that are known for their action in the CNS.
For the 78 Zinc molecules obtained from the virtual screening procedure, their physicochemical properties (Table 4) were collected from the ZINC server itself (http://zinc.docking.org/). All screened molecules were either within Lipinski's RO5. Lee and collaborators [24] presented the following physicochemical pro le of cannabinoid antagonists: MW ≤ 500, log P ≤ 7; HD ≤ 3 and HA ≤ 6. It is noted that log P increases due to the fat-soluble pro le of molecules that has CNS permeability characteristics, and the number of hydrogen donors and acceptors decreases. Evaluating phytocannabinoids through this prism it can be observed that molecule 14 would continue with 2 violations and that molecules 12, 13, 14 and 16 would violate the stipulated limit for log P. Evaluating Zinc molecules according to the physicochemical properties of cannabinoid antagonists, no molecule has more than a single violation [44]. It is noteworthy that the physicochemical pro le of cannabinoid antagonists is inherent in the molecules that permeate the blood-brain barrier, a fact related to the adverse effects of this class of molecules and unwanted in this study [42].

Pharmacokinetic properties analysis
The pharmacokinetic properties of phytocannabinoids are listed in Table 5 and were analyzed according to the parameters indicated in the QikProp manual. The evaluated parameters include: stars (similarity to known drugs: high: 0-2; medium: 3; and low: > 4), percentage of human oral absorption (% HOA) (high: > 80%; average: 25-80%; and low: < 25%), binding to human serum albumin (plogKhsa) (between -1.5 (low) and 1.5 (high)), pCaco-2 (intestinal cell) permeability (good: > 500 nm/sec; and poor: < 25 nm/sec) and pMDCK (renal cells) (good: > 500 nm/sec; and poor: < 25 nm/sec); logBB (low permeation in blood-brain barrier when < -1 and easy permeation in blood-brain barrier when > -1) and Van der Waals surface area (PSA > 70: does not exceed BBB; PSA < 70 exceeds BBB) [22]. Prediction of binding to human serum albumin indicates the ability of the molecule to be absorbed and bioavailable to interact with its target receptor. Evaluating molecules from Table 5 in terms of logKhsa, one can see that none showed low aggregation, 09 and 22 showed medium aggregation and others showed results that indicate high aggregation to serum albumin protein, a similar result to rimonabant.
Blood-brain barrier permeation was assessed by logBB (low permeation when < -1 and easy permeation when > -1) [45] and PSA (low brain permeation when > 70 Å) [46]. Rimonabant showed a prediction for high permeability in the blood-brain barrier, a factor associated with its side effects. According to logBB values obtained, only molecule 20 showed predictive results for low CNS permeation. On the other hand, according to PSA values obtained, molecules 03, 08 and 09 would have low CNS permeability. All other molecules tend to be well absorbed in oral administration and permeate the CNS due to their values of PSA < 70 [47]. The best results, therefore, were for molecules 03, 08 and 09, which had PSA> 70, despite medium cell permeability; and for molecules 01, 04, 06, 07, 15, 18, 21-25 and 27-40, with high cell permeability and prediction for action in the CNS.
Furthermore, we also predicted and evaluated the pharmacokinetic properties of the 78 molecules selected from the virtual screening using the ZINC database, as can be seen in Table 6.

Analysis of toxicological properties
Results obtained in the DEREK software were analyzed considering the following toxicological endpoints: mutagenicity, carcinogenicity, peroxisomal proliferation, hepatotoxicity, genotoxicity, neurotoxicity and reproductive effects. This software generates its output results using different con dence levels for each corresponding endpoint: CERTAIN (active), PROBABLE (probable), PLAUSIBLE / EQUIVOCAL (evidence is not strong / toxicity is not certain), DOUBTED (inactive but doubtful due to lack of experimental evidence) and IMPROBABLE / INACTIVE (unlikely) [49].
Rimonabant presented alerts for carcinogenicity (EQUIVOCAL) and HERG channel inhibition (PLAUSIBLE). When analyzing phytocannabinoids no alerts were red about their genotoxicity, neurotoxicity or reproductive effects endpoints. Molecule 24 presented human toxicity alerts regarding nephrotoxicity (equivocal) and hepatotoxicity (equivocal), referring to the styrene or derivative toxicophore. Molecule 32 showed alert for hepatotoxicity (plausible for the p-alkylphenol or derivative toxicophore). Such plausible and equivocal results do not indicate strong evidence of these toxicities and can only be con rmed by experimental evaluation [47,49].
Given the results obtained for the toxicological properties of phytocannabinoids, molecules 29 and 35 presented the best safety pro le (no toxicity alerts) but did not present satisfactory results for activity and ADME predictions. Thus, considering the predictions of biological activity and physicochemical properties, molecules that allow an additional in silico evaluation to elucidate their mode of interaction concerning our study objective, phytocannabinoids 01, 06, 07, 15 and 37 were considered, which present results for the prediction of toxicity even safer than those of Rimonabant.
The 78 molecules that were selected from the ZINC database were also evaluated in the DEREK software.
All alerts have been red as PLAUSIBLE, and their evidence of uncertainty can be con rmed or ruled out by experimental evaluation [47,49]. The results are summarized below: The best pro le Zinc molecules, considering the predictions of biological activity (anti-obesity activity potential >0.4) and physicochemical, pharmacokinetic and toxicological properties (safer than rimonabant), which were selected for in silico evaluation to better comprehend their interaction modes towards our target of interest are molecules 06, 17, 20, 38, 67, 69 and 70.

Molecular Docking Analysis
Molecular docking analysis was based on the results of Shao and collaborators [48] which conducted CB1 receptor molecular anchoring experiments in inactive conformation with PDB ID: 5U09. Validation by redocking [49] has been performed with a RMSD of 1.269 Å, as shown in Figure 4, when superimposing the originally complexed ligand taranabant from protein 5U09 (in green), and the best pose obtained from molecular docking (in blue). This result shows greater reliability when compared to the RMSD of Loo and collaborators [50], which was 2.74 Å for the best pose of antagonist AM115424 on the same protein.
The rCB1 structure has seven transmembrane domains (TM) connected through three extracellular loops (ECL1-3) and three intracellular loops (ICL1-3). The N-terminal domain has a V-shaped loop formed by residues 99-111 that acts as a lid and is, therefore, responsible for favoring or restricting access of a binder to the pouch from the extracellular side [51].
Filling the orthosteric pouch at the opening of TM1 and TM7 acts by blocking the entry of endocannabinoids, contributing to this blocking the amino acids Ile119 (TM1), Phe381 and Met384 (TM7). Interactions with amino acid residues of these two helices should favor the inactive conformation of the receptor. Leu387 belongs to the membrane access channel and favors the closing of the entrance channel, along with the amino acids mentioned above. TM1 is formed by residues Asn112 to Ser144, while TM7 is formed by residues Asn372 to Arg400. Arg214 and Asp338 favor canonical ion blockade. Lys192, Phe200 and Tyr275 indirectly contribute to conformational equilibrium, even without interacting with the antagonist. THC agonist binding is favored by Phe174, Leu193 and Ser383 [48].
The interactions that occur with the residues Phe170 and Phe174, from TM2, favor the inclination in relation to TM1 (Arg150-Val179), while the interaction with Lys192 (from TM3) decreases the a nity of agonists and stabilizes the conformation of ECL1, from the domain N-terminal and extracellular part of TM2 and TM3, [51] in addition to moving the receptor to its inactive state [41] .
Murineddu and collaborators [52], when docking antagonist AM6538 into the rCB1 structure (PDB ID: Finally, proceeding as described in materials and methods, we present our docking results, revealing potential interactions of Rimonabant (Fig 5a) and the ve selected phytocannabinoids that showed best physicochemical, pharmacokinetic, toxicological pro le, among those with Pa >0.4 in the PASS prediction, which were: 01 (Fig 5b), 06 (Fig 5c), 07 (Fig 5d), 15 (Fig 5e) and 37 (Fig 5f). Although they did not present a pharmacokinetic pro le favorable to peripheral action, knowing their way of binding to rCB1 allows obtaining additional parameters for analyzing the interaction of the virtual screening molecules on the ZINC server. The rCB1 orthosteric pocket shows hydrophobic characteristics since this receptor is activated by lipids [30], which explains the fact that the vast majority of occurring interactions are hydrophobic. All phytocannabinoids interacted with the N-terminal domain (mostly hydrophobic bonds), which is related to the restriction of membrane access by ligands from the extracellular side [51]. All had a radius of interaction with the active site below 10 Å, none showed interaction with Lys192, but molecule 37 (Fig. 5f) made an important hydrogen bond with Ser383 [52].
The molecule with the smallest radius of interaction with active site residues was 15 (Fig. 5e), which hinders for an agonist to enter by reducing the access channel, although it has interacted only with seven amino acid residues of three transmembrane domains (including Met384). Rimonabant (Fig. 5a) performed hydrogen bonding between Met103 and the HD region of the rimonabant. Its GoldScore was not surpassed by phytocannabinoids, the closest result being that of molecule 37 (Fig. 5f). The molecule that presented the greatest interaction with the transmembrane domains was molecule 06 (Fig. 5c), displaying fourteen interactions with eleven amino acid residues (including Phe170, Met384 and Leu387). However, it had the lowest GoldScore of the ranking (55.79).
The Zinc molecule 17 (Fig. 6b) interacted with ve transmembrane domains, including N-terminal, ECL2 and TM7, its GoldScore was close to 80, performing interactions with more amino acid residues from the active site (12) than rimonabant (10), performing interactions of hydrogen between Ser383 and hydrogen donor region, and between Pro269 hydrogen acceptor region. The Zinc molecule 20 (Fig. 6c)  GoldScore of 74.79 was slightly lower than that of phytocannabinoid 37 (Fig. 5f), and it presented two key interactions with Phe268 and a hydrogen bond with Ser383, which reduces endocannabinoid agonist binding and favors greater interaction strength and stability.

Page 29/45
The Zinc molecule 67 (Fig. 6e) had a higher GoldScore than phytocannabinoid 37 (Fig. 5f), but no higher than rimonabant (Fig. 5a). It performed 13 interactions with the active site in six transmembrane domains, with two hydrogen bonds between Ser383 and the hydrogen donor and acceptor region.
The Zinc molecule 69 (Fig. 6f) had the highest number of interactions with the active site (18) in ve transmembrane domains, with GoldScore greater than 80. It performed hydrogen bonds between Ser383 and the hydrogen acceptor region, and between Pro269 and region hydrogen donor. The Zinc molecule 70 ( Fig. 6g) was the one that had the best result in docking. Its GoldScore of 87.29, surpassing that of rimonabanto ( Fig. 5a)

Molecular Dynamics Analysis
The docking poses of the Zinc 69 and Zinc 70 compounds were selected as the initial conformations for the MD simulations. For the systems evaluated ( Fig. 7 the RMSD values of each system tend to equilibrium after 10 ns simulation time and the average values are less than 2.81 ± 0.55 Å, suggesting that the structures were stable throughout the MD simulation As show in Figure 8, the B-factor plot shows the loop region formed by the Val367-Lys373 residues (represented in red), that presented greater uctuation to the Zinc 69 compound. In addition, there is another exible region in the protein at the N-terminal domain.
The results of the binding free energy calculations of the two analyzed compounds with the MM-PBSA/GBSA methods are shown in Table 7. Table 7 shows that the formation of the two complexes formed by the enzyme and the ligands is favorable. Furthermore, comparing the binding free energy values (MM-GB(PB)SA) shows that the Zinc 69 ligand presented the lowest energy values when compared to the Zinc 70. The energy contribution of each residue for the complexes formed by the Zinc 69 and Zinc 70 obtained by the MM-GBSA method are shown in Figure 9. In this gure on the right side is show the residue contributions to binding energies using the plugin CHEWD [53] with Chimera software. From this, favorable interactions are at the blue of the color scale, that is, those that contribute to the stabilization of ligand into complex, while the red color represents interaction by residue with positive values and corresponds to unfavorable interactions values. Figure 9 (A and B) was used to explore the relative energy contribution to binding free energy decomposition of the complexes. Figure 9A revealed the residues that most signi cantly contribute to the total binding energy and, therefore, to the stabilization of the complex are Phe102, Met103, Phe170, Val196, Phe200, Phe268, Trp279, Trp356, Leu359 and Cys386. Figure 9B shows some of the residues that contribute to the stabilization of the enzyme-Zinc 70 complex.

Conclusion
Molecular modeling studies performed in this work allowed us to disclose that the selected phytocannabinoids evaluated here, have shown to be potential rCB1 antagonists, according to their predicted biological activities by the PASS server. Moreover, docking results showed a broad possibility of establishing key interactions between these molecules and the receptor, considering GoldScores values greater than the original ligand taranabant from the crystallographic structure. It is worth to note that the phytocannabinoid that would most likely present a strong tendency to act as a rCB1 antagonist, with interest in the anti-obesity activity, is molecule 37 since it presented a relevant hydrogen bond with Ser383 that is essential for antagonists' activity. However, all have a prediction of action in the CNS, pro le that is related to the triggering of adverse reactions inherent to the pharmacological target studied.
When comparing to another virtual screening campaigns previously found in the literature, our virtual screening performance showed promising results, because it was able to achieve potential molecules with positive predictions for rCB1 antagonism and also indicated their potential to possess anti-obesity activity, with Pa > 0.4. Their corresponding docking results were also signi cant since it was observed similar interactions to the reference and GoldScores values greater than 80.00, with the highest result for the Zinc 70 molecule surpassing rimonabant, and this may suggest binding a nities as intense as those observed for phytocannabinoids. To sum up, the best pro le screened Zinc molecules to act as rCB1 antagonist in peripheral antiobesity activity were 69 e 70 (for presenting key binding for the blockade of anandamide endocannabinoid agonist and the high GoldScores). MD simulations indicate that the selected ligands are stable, with little variation of the RMSD values along the trajectory, and the prediction of the binding energy with the MM-GB(PB)SA methods showed that the formation of the complexes are favorable. Besides, the decomposition of the binding free energy revealed that the complexes formed between the Zinc 69 and 70 form interactions with residues Phe102, Met103, Phe170, Val196, Phe268 and Trp356, which are important for catalytic activity.
In addition, these molecules have a possible factor favorable to the safety of the use of a peripheral level cannabinoid antagonist, due to its pharmacokinetic pro le, which may overcome the barriers of adverse effects inherent in this pharmacological class. Declarations