From the CB2 indexed molecules, ligands showing EC₅₀, Emax, and activity data were selected (total of 6356) and retrieved from the ChEMBL webserver as a .csv file. The “activity” category includes compounds with not only agonist activity, but antagonist, inverse agonists, and allosteric modulators as well. DataWarrior, an open-source data visualization software, was used to further analyze the ligand output which included discarding ligands without an agonist profile (−568 ligands), removing duplicates (−2159 ligands), and eliminating ligands with low activity (−773 ligands). This resulted in a final CB2 library of 2856 unique molecules that included a variety of chemotypes.

Ligands that have been indexed for TRPV1 activity within the ChEMBL database were selected and filtered in search of agonists in accordance with the reported EC₅₀ and E_max values and activity. The resulting 7,436 compounds were exported from the ChEMBL webserver as a .csv file and uploaded to DataWarrior. The selection of ligands with TRPV1 activity from the ChEMBL database included antagonists, inverse agonists, possible allosteric modulators, ligands with low activity, and duplicates which were all removed using DataWarrior. The final curated TRPV1 library contained 3,830 unique molecules with a variety of chemotypes.

The CB2 agonist resolved with the active hCB2 structure (AM12033) was used as an internal standard for CB2 docking. Three internal standards were used for hTRPV1: capsaicin in the VBP, AEA in the tunnel as observed from MD simulations, and CBD at the putative CBD site.

JWH133, which acts as an agonist at both CB2 and TRPV1, was used as a molecular basis for this additional screen to explore more unique scaffold options that may not be present in the CB2 or TRPV1 curated libraries. A JWH133 similarity library was curated using PubChem Biosays, 2021 which included compounds that shared >0.85 Tanimoto similarity index with JWH133, while also following Lipinsky’s rules of drug likeness (apart from xLogP values, which were set to −1 to 6 due to the lipophilicity of cannabinoid ligands). The JWH133 similarity library consisted of 5081 that were screened at all sites (CB2 and the three TRPV1 sites), and the output was analyzed to identify dual potential chemotypes.

Each of the curated libraries were exported as .sdf files and their conformations were optimized using the LigPrep module of the Maestro suite (Schrödinger, LLC, New York, NY, 2019). The Epik software was employed to predict pKa values in the pH range of 7.0 ± 0.5 and to return all chemically sensible structures in accordance with the Hammett and Taft methodology. All compounds were minimized using the OPLS3e force field as implemented in Maestro.

Molecular docking was performed using the HTVS Glide-dock module integrated in the Schrödinger package. The HTVS was conducted under the default setting, ensuring that high-energy ionization and tautomer states were removed, and the planarity of conjugated pi systems were enhanced. Ligands were docked flexibly, allowing for exploration of an arbitrary number of torsional degrees of freedom, in addition to the six spatial degrees of freedom spanned by the translational and rotational parameters. Up to 10 poses per compound state were generated and ligand poses that were generated in this way were run through a series of hierarchical filters to evaluate ligand interactions with the receptor. Docking score, glide gscore, glide emodel, ionization penalty, and topological polar surface area (TPSA) were used to select the docking poses in the output. The output from the HTVS contained the top 10% of the best scoring compound states and were analyzed for use in the extra precise (XP) screen via their docking scores.

Top scoring compounds from the HTVS were then studied through high-precision docking calculations which was performed using the XP Glide module. As with the HTVS protocol, 10 poses of the short-listed ligands were docked flexibly in their respective receptor site within the generated grids. A post-docking minimization was performed and the top 20% of the best scoring ligands were retained. XP Glide uses two key features that impact the XP Glide scoring: the recognition of structural motifs that provide large contributions to binding affinity and the application of large desolvation penalties to ligand and protein polar and charged groups wherever appropriate. To accomplish this, the sampling algorithm and scoring functions have been simultaneously optimized in XP. Ligands making it through the XP screen were organized by their docking scores and analyzed for ligand/receptor interactions. Selected ligands for each receptor were investigated through manual docking based on the automatic docking score, binding mode, as well as reported activity.

Manual Docking Identification of Potential PAINS Off-Targets Evaluation. Selected compounds were subjected to manual docking at CB2 and TRPV1 for further investigation of key interactions. Docking at CB2 was performed following the protocols previously reported by us for cannabinoid and related GPCRs (Morales et al., 2017a). In the case of TRPV1, select ligands were positioned within the respective binding site with steric clashes being removed via ligand and/or receptor adjustment using a graphical interface. Minimization of the ligand and surrounding 6 Å of residues (due to complex size) was performed using Prime version 19.3 (Schrödinger Inc.) with the OPLS3e forcefield in an implicit membrane.

In silico calculation of ADME properties. A set of 34 physico-chemical descriptors was computed using QikProp version 3.5 integrated in Maestro (Schrödinger, LLC, New York, United States). The QikProp descriptors are shown in Supplementary Tables S2, S3. The 3D conformations used in the calculation of QikProp descriptors were generated using LigPrep as previously detailed.

Identification of Potential PAINS. In the search of potential candidates, it is crucial to avoid the presence of potential promiscuous moieties or PAINS (pan-assay interference compounds) (Baell and Holloway, 2010; Capuzzi et al., 2017). Therefore, the selected molecules were subjected to a PAINS identification study using the swissADME webserver (Daina et al., 2017).

Off-Targets evaluation. XP Glide docks at potential off-target receptors including cannabinoid-related GPCRs such as CB1 (Shao et al., 2016; Hua et al., 2020), GPR55 (Kotsikorou et al., 2013; Lingerfelt et al., 2017), GPR18 (Sotudeh et al., 2019) and TRP channels such as TRPV2 (Pumroy et al., 2019), TRPV3 (Singh et al., 2018; Zubcevic et al., 2018), TRPA1 (Suo et al., 2020; Zhao et al., 2020), and TRPM8 (Diver et al., 2019; Yin et al., 2019). For this purpose, the cited available structures, whether crystal, cryoEMs, or models previously developed in our group, have been used. Results of these additional dockings can be found in Supplementary Tables S4, S5.

A chemical library of compounds with >0.85 Tanimoto similarity index with JWH133 was curated and screened at CB2 and TRPV1 using the methods described above (workflow depicted in Supplementary Figure S2A). Analysis of docking interactions of top-ranked XP results from the CB2 site and the TRPV1 sites revealed seventeen common ligands between the CB2 site and the VBP and CBD sites (Table 1). No common ligands were identified between the CB2 site and the TRPV1 tunnel. The selection of novel potential CB2/TRPV1 chemotypes includes key structural features and ligand-receptor site interactions at both targets, as well as the absence of previously reported activity at these receptors. This strategy allows for prioritization of molecularly diverse and novel compounds. To ensure the VBP and CBD sites were explored equally, two ligands were selected that targeted CB2 and the VBP (57756957 and 151332252), two ligands were selected that targeted CB2 and the CBD site (59824268 and 123533625), and one ligand was selected that targeted CB2, the VBP, and the CBD site (153641693), resulting in the selection of five ligands for further investigation via manual docking and pharmacokinetic profiling (Supplementary Table S2). The selected chemotypes have not been yet explored at CB2/TRPV1 and their reported activity is not significant, providing novel opportunities for the investigation of the endocannabinoid system.

Among the five selected candidates, 3-(4-methylbenzyl)-chromane 59824268 presented a better druggability profile (Supplementary Table S2) being therefore prioritized for future in vitro testing as CB2/TRPV1 dual modulator. As shown in Supplementary Table S2, candidates 57756957, 123533625, 153641693, and 151332252 exhibit HERG values that fall outside the range of approved drugs. Docking studies of 59824268 at CB2 and TRPV1 are shown in Figure 1. At CB2 this chromane sits a bit higher than JWH133 in the binding crevice being stabilized by hydrophobic and aromatic interactions with residues F2.64, F2.61, F2.57, F183 and F7.35. Regarding TRPV1, 59824268 orients itself in a way similar to JWH133 in the pocket maintaining overlap with the aromatic ring. While JWH133 appears to have primary interactions with Y584, 59824268 has interactions with Y584 in addition to Y632 and F639, further stabilizing the chromane in this pocket.

The second strategy for the identification of dual compounds is based on a HTVS of reported CB2 and TRPV1 agonists. CB2 agonists indexed in the ChEMBL database have been retrieved and studied in the three known TRPV1 binding sites as detailed in section 2.4. Likewise, TRPV1 ligands indexed in the ChEMBL database have been retrieved and studied in the CB2 binding site. Following the workflow depicted in Supplementary Figure S2B, five candidates were selected for further analysis at each receptor (Table 2). Reported activity at the known target, docking score (Table 2) and druggability profile (Supplementary Table S3) led us to select compounds 1288208, 1288239 (TRPV1 virtual screening) 1508577 and 1508215 (CB2 virtual screening).

1288208 passed the screening as a potential modulator of TRPV1 at two sites: the VBP and tunnel. While there is argument for the elimination of this ligand due to lack of site specificity, it was selected for exactly this reason. With the abundance of ligands reported to modulate TRPV1, and the variability in reported and putative binding locations, a ligand that shows the potential for interaction at multiple locations, both putative and confirmed, within the channel is worthy of further study to better understand why this is. In the VBP, the headgroup of 1288208 forms H-bonds with R557 and S512 via the backbone and hydroxy group, both reachable from within the tunnel. The α,β-unsaturated ketone oxygen H-bonds with Y511, and the addition of the naphthyl moiety at the tail end of the ligand provides pi-stacking capabilities farther up in the VBP with F543 and F591 (Figure 2A).

AEA docking in the tunnel shows headgroup interactions with several residues including Y554, Y555, Y487, D708, and N438 (Supplementary Figure S4). Like AEA, the hydroxyl portion of the headgroup of 1288208 maintains interactions with Y554 and D708, while additional H-bonding between R491 and the amide oxygen is present. The inclusion of an α,β-unsaturated ketone mid-tail allows for more H-bonding via Y487 and N438 near the entrance of the tunnel. The naphthyl moiety at the end of the tail displays pi-pi interactions with both Y487 and Y445 (Figure 2B). The additional interactions of 1288208 could aid in the stability of the ligand in the tunnel from an external standpoint, allowing the headgroup more time in the tunnel, potentially triggering channel activation as previously hypothesized from MD simulations (Muller et al., 2020, 2021).

The other selected TRPV1 ligand, 1288239 (Figure 2C), shows headgroup interactions with R491 and Y554, like 1288208, with an additional interaction with E513. The ketone found midway down the tail of ligand H-bonds with N438 and Y487, again similar to 1288208. One feature that differentiates 1288239 from 1288208 is a biphenyl moiety in place of a naphthyl moiety. The lower ring of the biphenyl moiety has aromatic interactions with Y487 and F488, and both rings interact with Y445.

Because of their high potency at TRPV1, their interaction pattern at the CB2 orthosteric pocket and their optimal drug-like properties, compounds 1508577 and 1508215 were selected as potential candidates in the in silico search of dual ligands. Compounds like ACPA-OH and JWH133 also ranked at the top, however, since we are looking for unexplored dual chemotypes, they were not selected in this in silico study. Consistent with the hydrophobic nature of the CB2 orthosteric pocket, compound 1508577 is mainly stabilized by aromatic and van der Waals interactions. As displayed in Figure 2D, π-π stackings are stablished between the methoxybenzene group with W5.43 and F183 and the phenylacetamide group with F183, F7.35, F2.64 and H2.65. Moreover, the acetamide hydrogen engages with the backbone carbonyl oxygen of V182 in a H-bond while the central amide H-bonds with S7.39. Compound 1508215 orients similarly in the binding crevice stablishing aromatic π-π interactions between the central indole core and F183, F2.61 and F2.64, and the fluorobenzyl group with F2.57. In addition, the imidazolidinedione group H-bonds with the backbone carbonyl oxygen of V182.

In summation, from this approach, compounds 1288208, 1288239, 1508577, and 1508215 have been selected for future in vitro appraisal as dual CB2/TRPV1 agonists. Other compounds such as 1644371 could also be remarkable candidates for testing at TRPV1 due to its nanomolar agonist potency at CB2.

The selected hits (59824268, 1288208, 1288239, 1508577 and 1508215) have also been docked in related receptors in order to identify potential off-target effects. These molecules have been screened at CB1 and the cannabinoid-related GPCRs GPR55 and GPR18 in their active and inactive states. In addition, cannabinoid-related channels including TRPV2, TRPV3, TRPA1, and TRPM8 have also been assessed. As shown in Supplementary Table S5, by comparing docking scores to their reference orthosteric ligands we can conclude that at cannabinoid GPCRs compounds 1288208 and 1288239 might be more promiscuous showing high interaction energies at the GPR55 active and GPR18 inactive models. Moreover, compounds 1508577 and 1508215 may moderately act at CB1 whereas 1288208 and 1288239 were reported to lack binding affinity (Osman et al., 2010). 59824268 may be less selective with higher energies for the apo TRPV3 structure as well as both TRPA1 structures. 1288208 and 1288239 both show energies that are either comparable to or better than the reference ligand for each respective receptor, perhaps suggesting that the ethanolamide head group may be too promiscuous of a moiety to include when aiming to develop ligands for selective dual targeting. 1508577 shows variable activity across the TRP channels with comparable energies to the reference compounds of TRPV3 and TRPA1 in both states, with 1508215 displaying the potential for promiscuity at TRPA1.In light of these results, compounds 59824268, 1508577 and 1508215 could be prioritized as TRPV1/CB2 dual modulators. However, compounds with moderate activity at other cannabinoid targets could also be beneficial when targeting specific pathologies in which the ECS is involved or avoided when searching for more selective cannabinoid modulators.

Nonetheless, off-targets cannot be completely ruled out and not only cannabinoid-related but also other receptor families should be tested experimentally at further stages of this project.