The CBRs play a critical role in several human physiological and pathological conditions. However, the CNS side effects of CB1R ligands may limit the therapeutic use of such agents if they cross the Blood-Brain Barrier (BBB). That is the case of the CB1R inverse agonists Rimonabant and Taranabant (Moreira and Crippa, 2009; Martín-García et al., 2010). To overcome the central effects, peripheral CB1R antagonists were developed (Chorvat, 2013; El-Atawneh et al., 2019; Quarta and Cota, 2020). Another option is to develop pure antagonists (An et al., 2020; Stasiulewicz et al., 2020). Agonists of the CBRs may be used to treat anxiety (Stasiulewicz et al., 2020) or as analgesics, anti-inflammatory, neuroprotective and anti-emetic compounds (An et al., 2020). Peripheral CB1R antagonists combined with CB2R agonists may be used for treating liver diseases (Mallat et al., 2011) and diabetic complications (Gruden et al., 2016). This dual activity may be useful in treating obesity, abolishing diabetes-induced albuminuria, inflammation, tubular injury, and renal fibrosis (Barutta et al., 2017). Combining CB1R antagonism with CB2R agonism in the brain is shown to have a synergistic effect on reward processing (Gobira et al., 2019). Another option is to design selective CB2R agonists to benefit from their nociception and neuroinflammation role without psychoactive effects (Hollinshead et al., 2013; Verty et al., 2015; Poleszak et al., 2020). CB2R selective agonists are investigated to treat pain, inflammation, arthritis, addictions, cancer besides their neuroprotective role (An et al., 2020).

CB2R could be targeted with other receptors to attenuate inflammation for several autoimmune and inflammatory conditions. The peroxisome proliferator-activated receptor (PPAR)-γ is a nuclear receptor that plays a crucial role in regulating lipid metabolism and glucose homeostasis. It associated with metabolic disorders, such as atherosclerosis, obesity, metabolic syndrome, dyslipidemias, type 2 diabetes, and cancer (Decara et al., 2020). PPARγ agonists have been shown to prevent inflammation, dermal fibrosis, and lipoatrophy in preclinical models of systemic sclerosis (SSc) (Wei et al., 2010). SSc is an orphan autoimmune multi-organic disease that affects the connective tissue. Dual CB2/PPARγ agonists such as VCE-004.8 and JBT-101 (Ajulemic acid, Lenabasum) have alleviated skin fibrosis and inflammation in SSc models (Rio et al., 2018; García-Martín et al., 2019). JBT-101 is in clinical trials for SSc (NCT03398837), dermatomyositis (NCT03813160), and cystic fibrosis (NCT02465450). Additionally, PPARγ agonists can suppress the pro-inflammatory cytokines associated with chronic diseases such as Inflammatory Bowel Disease (IBD).

IBD, including ulcerative colitis (UC) and Crohn’s disease (CD), has been considered one of the most prevalent GI diseases with accelerating incidence in newly industrialized countries. Yet it lacks effective drug targets and medications (Seyedian et al., 2019). As a lifelong disease, therapy aims to induce remission in the short term and maintain remission in the long term. New drugs have diverse mechanisms of action, targeting mainly the inflammation pathways. The current anti-inflammatory small molecules used to treat IBD are associated with several side effects (5-amino salicylate and its prodrugs such as Olsalazine and Balsalazide), with more severe toxicity (Azathioprine, Mercaptopurine, Methotrexate) or with known long term negative impacts of steroid hormones (glucocorticoids). Biological drugs are expensive, require more intensive medical attention in a clinic or at home (self-injections), and, in the case of TNFalpha antibodies, elicit resistance by immune system response (Torres et al., 2020). Although the mechanism by which PPARγ acts on the pathogenesis of IBD has not been clarified (Decara et al., 2020), natural and chemical PPARγ ligands have ameliorated the fibrotic process in preliminary clinical trials and experimental models of intestinal fibrosis (Vetuschi et al., 2018). Moreover, many studies showed the anti-inflammatory role of PPARγ activation in intestinal tissues in UC and CD (Decara et al., 2020).

Recent investigations suggest that serotonin (5-HT) can influence the development and severity of inflammation within the gut, particularly in the setting of IBD. 5-HT influences every major function inherent to the gut, including motility, secretion, blood flow, and sensation (Coates et al., 2017). Alterations in its receptor activity in disease conditions may result in many problematic symptoms, including abdominal pain, diarrhea, or constipation (Coates et al., 2017). The 5-HT4 receptor (5-HT4R) mediates enteric neuron survival and neurogenesis of adult mice (Liu et al., 2009). It promotes the reconstruction of an enteric neural circuit leading to the recovery of the defecation reflex in the distal gut (Matsuyoshi et al., 2010). 5-HT4R activation maintains motility in healthy colons of mice and guinea pigs and reduces inflammation in colons of mice with colitis (Spohn et al., 2016). PPARγ and 5-HT4R agonists may be combined with CB2R as a potential therapy for IBD (Turcotte et al., 2016). A peripheral CB2R agonist (Olorinab) reached phase II trials for abdominal pain in CD (NCT03155945) and irritable bowel syndrome (NCT04043455).

Computational methods allow us to examine options for designing or discovering multitargeting candidates in a reliable, fast, and low-cost manner (Sliwoski et al., 2014; Zhang et al., 2017). Screening candidates for binding against several targets to find single MTA differs from designing compounds based on conjugated pharmacophores by merging/fusing/linking molecules (Morphy and Rankovic, 2005; Zhou et al., 2019), which could take longer to synthesize and might increase the molecular weight and affect the drug-likeness properties.

Our research combines ligand and structure-based methods. Our algorithm for solving complex combinatorial problems, the 'Iterative stochastic elimination algorithm’ (ISE) (Stern and Goldblum, 2014; El-Atawneh and Goldblum, 2017), has been applied in recent years to molecular discovery (Zatsepin et al., 2016; Da’adoosh et al., 2019; El-Atawneh et al., 2019), including one example of multitargeting modeling: modeling the properties of molecules that may be remotely loaded to nanoliposomes and the properties that enable them to be stable inside the nanoliposomes, in a biological fluid (Cern et al., 2017). Molecules that had high scores in both loading and stability models were chosen. For any discovery of MTAs, virtual screening (VS) by separate ligand-based models is performed in sequential order.

After finding top candidate ligands, it is helpful to examine the structural aspects, since our classifications are based on physicochemical properties and not on structural elements. Molecules with similar properties might have different structures and sizes. Thus, we dock the top candidates to the target protein if such a structure has been reported. Structures of CB2R were deposited recently in the Protein Data Bank (PDB), one with a bound antagonist (PDB code 5ZTY) (Li et al., 2019) and the other with an agonist (PDB code 6KPC) (Hua et al., 2020), which makes structure-based design feasible (Tuccinardi et al., 2006; Cichero et al., 2011). CB2R shares 44% sequence identity and 68% similarity with CB1R in the transmembrane regions (Munro et al., 1993). The antagonist-binding pockets in both receptors are quite distinct, while the agonist-binding pockets in CB1R and CB2R, including side-chain rotamers, of the key residues involved in ligands interactions are almost identical (Li et al., 2019; Hua et al., 2020), which might be the source of cross-reactivity between their ligands and difficulty in attaining selectivity. There are also CB1R and PPARγ structures, with agonists and antagonists in both. Yet, there is no published atomic-level structure of 5-HT4R, but ligand-based modeling for 5-HT4R with ISE is possible due to its many known ligands.

All molecules were prepared by the “Molecular Database Wash” (v. 2011.10) (Molecular Operating Environment, 2021). This includes hydrogen adjustment, removing minor components, determining the protonation state, enumeration of ionization states, and tautomer forms. Mutagenic and reactive molecules (based on calculated descriptors by MOE) were removed from the learning sets.

The standard descriptors we calculated for building the models are the 2-dimensional (2D) descriptors by QuaSAR- MOE (v.2011.10) with 186 descriptors. The complete descriptors list is given at (http://www.cadaster.eu/sites/cadaster.eu/files/challenge/descr.htm). Descriptors with low variance (Smialowski et al., 2010), or highly correlated descriptors (Pearson correlation coefficient > 0.9), were excluded, using the Knime platform (v. 4.0.1) (Berthold et al., 2008) to exclude out of two highly correlated descriptors the one which has greater similarity to other descriptors. We have also tested the performance of 3D descriptors for CB2R (see results and discussion).

Our generic ISE algorithm has been applied to many problems related to drug discovery and has been presented in reviews, with details of the mathematical and statistical criteria to distinguish between two activities based on physicochemical properties (descriptors) of known active vs. inactive compounds (Stern and Goldblum, 2014; El-Atawneh and Goldblum, 2017). For each model, five cross-validations were performed (James et al., 2013), with 4 out of the five-folds producing the model, and the fifth fold was used as a test set. We include some of the main details of model construction and screening in Supplementary Data section 1.1.

The “Atom-pair” fingerprints for the active molecules were generated using RDKit toolkit (RDKit, 2018) (in Knime platform v. 4.0.1) (Berthold et al., 2008). The “Tanimoto similarity coefficient” (Tc) for the fingerprints is based on the CDK toolkit.

The two structures of CB2R were downloaded from the PDB (5ZTY (Li et al., 2019) and 6KPC (Hua et al., 2020)), and prepared by the “Protein Preparation Wizard” (Schrödinger Suit 2019-3) (Madhavi Sastry et al., 2013). For 5ZTY, we allowed C-OH rotations of SER90, THR114, TYR190; for 6KPC, we allowed such rotations of TYR25, SER90, THR114, TYR190, and SER285 for the grid construction. Alanine (ALA) scan was performed to assign the critical residues in the binding site of the two structures for 23 residues detected by PDBsum (Laskowski, 2009). The screened molecules were prepared using “LigPrep” (Schrödinger Release, 2018), with default settings, except the chirality option that was set to “Generate all combinations” for the Enamine database (5,024,833 entries were generated). Molecular docking was performed with Glide HTVS and SP (Richard A. Friesner et al., 2006).

In the docking analysis, we examined the geometric character of binding by requiring the docked molecules to be in contact with residues that were found to be “hot spots” by performing a virtual ALA scan.

To find multitargeting candidates for the different indications, we performed hierarchical VS. First, focusing on the CBRs, we screened the Enamine database (DB) through the different CBR activity models, considering desired activity, i.e., of CB2R agonists, and the unwanted activity as anti targets. Molecules with a positive index pass the model, and those with a negative score are considered to fail. We found 241,260 CB2 selective agonists (about 11% of the dataset); those molecules passed the CB2R agonist model and did not pass the CB2R antagonist model. They also did not pass the CB1R agonist and antagonist models. Adding the CB1R agonists or antagonists to CB2 agonists, we found many less candidates (63,735 and 324, respectively), as shown in Figure 1. Raising the index cutoff above 0.0 reduces these numbers.

Looking for additional activities of the selective CB2R agonists, we screened those 241,260 candidates through the PPARγ and 5-HT4R agonist models (Figure 1). To avoid anti-targets we screened the same set by the antagonist models of PPARγ and 5-HT4R. This yielded 374 CB2R and PPARγ agonists, and 14,008 candidates for CB2R and 5-HT4R agonism with no antagonism at any of the three receptors. We found 28 candidate agonists for simultaneously hitting all the three targets of CB2R, PPARγ, and 5-HT4R. All the mentioned hit sets are internally diverse, as well as being diverse (by Tanimoto criteria) towards the actives used for model construction: comparisons yield a low average Tanimoto coefficient of Tc ≤ 0.4 (Supplementary Table S3).

The structures of CB2R (6KPC (Hua et al., 2020) with an agonist and 5ZTY (Li et al., 2019) with an antagonist) have similar binding pockets and binding residues (Li et al., 2019; Hua et al., 2020) (Supplementary Table S4). Similarity is also observed between the CB2R and CB1R binding pockets (Li et al., 2019). This creates an obstacle to distinguishing between agonist and antagonist activity for the CB2R if we consider docking alone. We examined the binding residues in both structures by applying a virtual ALA scan (Schrödinger Suit 2019-3) (Madhavi Sastry et al., 2013) for 23 residues in the binding site (Supplementary Table S4). AM12033 (6KPC- CB2R agonist) has 19 interactions, mainly with hydrophobic and aromatic residues and 3 H-bonds, with LEU 182 and SER285. AM10257 (5ZTY- CB2R antagonist) has 16 interactions with no H-bonds (as shown in PDBsum (Laskowski, 2009)).

The calculated stability for the 23 residues (by virtual ALA scan) does not differ dramatically between 6KPC and 5ZTY. The considered contacts in the 6KPC agonist structure in order to suggest more successful docked ligands are: hydrogen bonding with LEU182 and SER285, and Van der Waals (VDW) interactions with the following: TYR25, PHE87, PHE91, PHE94, ILE110, PHE183, TYR190, LEU191, TRP194, LEU262, MET265, PHE281.

We compared ligand (ISE) and structure-based (docking) methods by performing VS of the Enamine DB (2,159,632 compounds) for CB2R. ISE screening is extremely fast compared to docking (Figure 4). A positive index in screening by the CB2R agonist model was assigned to 241,260 molecules. We pick molecules with higher indexes and better EF values to improve the quality of our candidates, thus resulting in fewer molecules. For example, with a high index cutoff ≥ 0.7, 41,102 molecules pass, and the EF equals 54. That EF is only 17 at a lower index cutoff >0.0 (for 241,260 molecules). Docking was applied to the ISE candidates with a positive index: SP docking to the 6KPC structure found 238,718 molecules with docking scores of 6.6 to −12.8 kcal/mol. Filtration was based on docking scores ≤ −9 kcal/mol and hydrogen bonds with LEU182 and SER285, to a final set of 131 candidates.

Docking to CB2R was performed in two stages with the same 6KPC structure. First, HTVS docking was executed for the whole Enamine DB. The docked poses have a docking score range from 10.4 to −12.5 kcal/mol. Molecules with docking scores of less than −9 kcal/mol were further docked by the SP protocol (130,358 molecules). Most of these molecules (130,080) passed SP with a 5.7 to −12.9 kcal/mol docking score. By picking those with a score better than −9 kcal/mol and hydrogen bonds with LEU182 and SER285, only 73 molecules remain. Ten out of the 73 docking hits have positive ISE index scores. Only nine molecules are shared between the two SP screenings. Both sets are diverse from the known active CB2R agonists, and from each other (average Tc ∼0.3).

Multitargeting by ISE could be based on molecules with known activities on two or more targets. One publication mentions the construction of such a database, but it is not accessible (Chen et al., 2017). It is highly unlikely that enough molecules will be found to enable ISE modeling. Therefore, in the main spirit of ISE, each “variable” (in that case, a target, with many ligands as its “values”) requires separate model construction. Screening and scoring through any single model reduce the molecular library size by 10-fold or more. In HTS, it is common to discover 1 out of 1,000 molecules tested for activity. However, that is a real activity in vitro, while we only suggest candidates for in vitro testing, which may include false positives. Therefore their numbers are much larger.

As we add more targets and anti-targets, the number of candidates decreases: we found, among our ∼2.1 million screened molecules, only 374 candidates for combined (simultaneous) CB2R and PPARγ agonism, which may be tested for SSc (Wei et al., 2010), dermatomyositis, cystic fibrosis, and IBD (Decara et al., 2020). Adding 5-HT4R agonists reduces that number to 28, while CB2R and 5-HT4R agonists that could be valuable for IBD have 14,008 candidates. The much larger number of shared molecules that could hit CB2R and 5-HT4R (compared to sharing between CB2R and PPARγ) reflects the fact that both are aminergic GPCRs of the A family with 27% sequence similarity, as calculated by blastp (McGinnis and Madden, 2004), and may have a greater chance for ligand cross-reactivity (Yang et al., 2021). PPARγ belongs to a different family of cytoplasmic nuclear receptors. Moreover, only 60 molecules are shared between PPARγ and 5-HT4R agonists (without screening through CB2R models).

Screening by ISE models has already succeeded in achieving “scaffold hopping” (Zatsepin et al., 2016; Da’adoosh et al., 2019; El-Atawneh et al., 2019) due to the use of physicochemical properties rather than of structures. Even in those cases of greater similarity among the actives (agonists of PPARγ (0.52) and of 5-HT4R (0.5), Supplementary Table S1), the top screened candidates are varied among themselves, i.e., Tc = 0.4 for the 28 multitargeted agonists of CB2R/PPARγ/5-HT4R. That is also the case of screened molecules vs. actives in the learning sets (all results in Supplementary Table S3).

The main substructure elements presented in Figure 2 may help to understand how it is possible that a single molecule binds to different binding sites: the amine moieties—frequently two amines in a molecule—are singly charged, and the first protonation reduces the pKa of the other amine. Amine protonation prevails in four out of the five multitargeted sets, except for CB2R/CB1R agonists in which a negative charge on the oxygen of the amides may have a leading role. It is also clear from the difference between the coupling of CB2R agonists with either CB1R agonists or antagonists, that it is possible to separate between these multitarget pairs. It would still be impossible to suggest a synthesis of multitargeted compounds based on these major fragments, but it is easy to pick molecules that contain these fragments for each multitargeted alternative by requiring to include these substructures with their statistical weight as in Figure 2 or even better, as in Supplementary Figure S1. None of these moieties resemble the structures of known cannabinoid ligands (classical, non-classical, amino-alkylindoles, and those with the eicosanoid group).

Structures of CB2R have been recently deposited in the PDB (Li et al., 2019; Hua et al., 2020) and enable to perform structure-based studies—docking, pharmacophore, and molecular dynamics. The similarity between CB2R agonist/antagonist complexes and CB1R and CB2R structures make it challenging to design ligands with high selectivity (Hua et al., 2020). Docking is considered a time-consuming approach, as shown in Figure 4. Screening by docking has been shown to be much less reliable statistically than our ligand-based approach for CB2R agonism. Our ISE models screen molecules based on their properties and not on structural elements. That may result in top screened molecules having similar properties but different sizes and volumes, which may or may not be accommodated by the targets. Some of these molecules might not fit into binding sites and will be rejected. The results of our CB2R modeling confirm our preferable sequence of actions: ligand-based modeling should be followed by structure-based testing, which is better than structure-based docking alone.

Virtual ALA scan was used in this and other of our studies for picking “hot spots”—the main residues that contribute to the binding of smaller or larger ligands (i.e., including protein-protein interactions). Those “hot spots” determine the region of the grids for screening by docking and provide the initial geometric criteria that are applied prior to considering the docking scores. In ALA scan, we replace a larger side chain (of 18 amino acids, except for GLY and ALA) with a shorter one. We do not however apply any minimization or dynamics to that change, which positions a methyl group in the Cβ position, with tetrahedral angles vis-à-vis Cα, in place of a longer side chain, leaving some “void”. No other side chain position is modified around the virtually mutated one. This protocol is due to our wish to discover molecules that replace an existing ligand/protein with an exact similar conformation of side chains in the protein target, as in the PDB, in order to promote competition. That is clearly not the case with genetically mutated ALA scan. In that in vitro experiment, other side chains could change their conformations in the vicinity and more remote from the ALA mutated position. In vitro ALA scan may even change conformations of the main protein chain. Therefore, it is rewarding if mutagenesis studies support some of our results such as for PHE87, PHE91, PHE94, HIS95 (Li et al., 2019), and TRP194 (Zhang et al., 2011). TYR190 mutation to Ile resulted in a loss of ligand recognition and function (McAllister et al., 2002).

This is a theoretical study, which includes statistics (AUC, EF) that clarify what are the chances for discovering multitargeted actives. Naturally, the next step is to pick top candidates from each set for biochemical experiments. Our multitargeting results also suggest which multitargeting sets have a greater chance to be experimentally confirmed. Previously, we published our theoretical predictions and experimental validations of the binding of 8 molecules out of 15 predicted candidates (picked by ISE modeling from a library of 1.8 million) (El-Atawneh et al., 2019). Finally, only in vitro testing of candidates predicted by each method in silico will confirm or refute the VS results conducted by ISE and docking approaches.