Computational Prediction and Biochemical Analyses of New Inverse Agonists for the CB1 Receptor.

Human cannabinoid type 1 (CB1) G-protein coupled receptor is a potential therapeutic target for obesity. The previously predicted and experimentally validated ensemble of ligand-free conformations of CB1 [Scott, C. E. et al. Protein Sci. 2013 , 22 , 101 - 113 ; Ahn, K. H. et al. Proteins 2013 , 81 , 1304 - 1317] are used here to predict the binding sites for known CB1-selective inverse agonists including rimonabant and its seven known derivatives. This binding pocket, which differs significantly from previously published models, is used to identify 16 novel compounds expected to be CB1 inverse agonists by exploiting potential new interactions. We show experimentally that two of these compounds exhibit inverse agonist properties including inhibition of basal and agonist-induced G-protein coupling activity, as well as an enhanced level of CB1 cell surface localization. This demonstrates the utility of using the predicted binding sites for an ensemble of CB1 receptor structures for designing new CB1 inverse agonists.

Introduction

G-protein coupled receptors (GPCRs) are located in the cellular membrane and act as mediators for cell signaling, making them promising therapeutic targets. However, drug design for GPCRs has been challenging due to the paucity of validated 3D structures. Although these integral membrane proteins have long been difficult targets for structure determination efforts, recent methodological developments now provide structures for ∼33 of the 819 human GPCRs. Six of these, bovine rhodopsin,[1−4] human β2 adrenergic receptor,[5,6] human adenosine A2A receptor,[7−10] human muscarinic M2 receptor,[11] human neurotensin NTS1 receptor,[12] and human P2Y12 receptor,[13] have been crystallized in an active or partially active conformation, but only one, human β2 adrenergic receptor, has been cocrystallized with the full heterotrimeric Gs protein.[5]

The human cannabinoid 1 (CB1) is a GPCR with high therapeutic potential as a drug target. This receptor is activated by Δ9-tetrahydrocannabinol (THC) in marijuana and increases appetite for AIDS and cancer patients. The CB1-selective inverse agonist/antagonist rimonabant (also known as SR141716A)[14] is an antiobesity drug that was available in Europe and was in phase III of US FDA clinical trials but had to be withdrawn due to severe depressive effects. There is some indication that these side effects are caused by rimonabant antagonizing the CB1 receptor activity in brain,[15] underscoring the need for biased inverse agonists for CB1 that are devoid of these unacceptable side effects. Nonetheless, the intracellular pathways mediating the beneficial and unwanted effects have not yet been elucidated.

CB1 has not yet been crystallized; however our previous studies, using the GEnSeMBLE (GPCR Ensemble of Structures in Membrane Bilayer Environment) complete conformational sampling method, predicted the ensemble of ten low energy CB1 conformations[16,17] expected to play a role in binding various ligands and controlling function. Briefly, GEnSeMBLE carries out systematic sampling of trillions of seven-helix bundles by rotating and tilting the transmembrane (TM) helices starting with templates from previous calculations or X-ray experiments on other GPCRs. We consider it important to examine this full conformational space for CB1, because GPCRs are dynamic and flexible, enabling different ligands to bind to and stabilize quite different GPCR conformations. This, in turn, facilitates coupling to multiple types of intracellular proteins.[18−20]

For each of these ten CB1 conformations we used the hierarchical binding site prediction methods, DarwinDock and GenDock,[21−25] to predict the optimum structures for binding of rimonabant and seven other ligands. Our predicted ligand–receptor complexes are in good agreement with published site-directed mutagenesis data[26−35] and structure–activity relationship (SAR) data.[36,37]

We then used these predicted binding sites to design and identify 16 new ligands that we expected might act as inverse agonists upon binding to the CB1 receptor. Subsequently we tested five of these ligands experimentally and found that two compounds possess inverse agonist properties including reduced G-protein coupling and changes in cellular localization of the receptor. This provides new design strategies for developing more selective and potent CB1 inverse agonists through computational optimization of chemical functionalities of the compounds coupled with further pharmacological characterization.

Methods

We used our GEnSeMBLE method[22,38−41] to predict an ensemble of ten low energy CB1 receptor structures, described elsewhere.[16,17,41] These predictions started with four distinct templates [bovine rhodopsin,[42] turkey ß1 adrenergic receptor,[43] human ß2 adrenergic receptor,[44] and human adenosine A2A receptor[45]], optimized the helix shapes within the turkey ß1 adrenergic receptor template, sampled all (12)7 = 35 million rotations of each of the seven helices independently, followed by sampling of 13 trillion combinations of helix rotations and tilts to select by energy the best ten conformations of the seven-helix bundles. We then docked various ligands, finding the best GPCR conformation for each ligand. Previous papers[16,17,41] have analyzed thoroughly this ensemble for the wild-type and several mutant receptors, which were used successfully to predict (in advance of experiment) mutations that would bias the ensemble toward the fully inactive, the constitutively active, or the essentially fully active states.

Docking of Rimonabant to CB1

To obtain a diverse set of ligand conformations for docking, we generated a total of 83 ligand conformations for rimonabant. The molecular structure of rimonabant was constructed with Maestro software,[46] and a conformational search was performed with MacroModel software,[47] which generated 83 conformations of the ligand. Sytematic and extended torsional sampling options were used where the selected rotatable bonds are rotated 360° in 30° increments. Ligand conformations that fall within an energy window of 10 kcal/mol and an RMSD diversity of 0.5 Å were saved for subsequent steps. The conformational search was conducted with the OPLS 2005 force field[48] and in a dielectric of 80.37 to match water. Subsequently, we preformed two rounds of clustering for each ligand; in the first round we clustered ligands with a 2.0 Å diversity followed by another round of clustering with a 1.0 Å diversity. The Mulliken populations of each atom were calculated with Jaguar software[49] using Density Functional Theory (DFT) with the B3LYP functional and the 6-31G** basis set. This led to the selection of nine ligand conformations. An additional ligand conformation was constructed from existing coordinates from the crystallized rimonabant in methanol solvate[50] deposited in the Cambridge Structural Database.[51] Each ligand was minimized using the Surface Generalize Born (SGB) solvation model[52] for 100 steps or to a convergence threshold of 0.2 kcal/mol/Å RMS force with the MPSim program.[53] For each of the ten conformations of rimonabant, we used the DarwinDock and GenDock[21−25] methods to predict the optimum binding site for each of the ten lowest energy protein structures predicted by GEnSeMBLE for the CB1 receptor.

DarwinDock

The DarwinDock method aims at generating a complete set of poses for the binding pocket while using RMSD clustering of the poses to dramatically reduce the computational cost. To provide flexibility and space for the ligand to identify favorable binding sites, we replaced the seven bulky hydrophobic residues (FILMYVW) with alanines. The mutated residues are called alanized residues, and the mutated protein is called alanized protein. For each of the best 100 ligand poses in the alanized protein, we then dealanized the mutated residues back to their original hydrophobic identity and optimized their positions along with those of other residues in the binding site using SCREAM.[54] This leads to a unique set of optimized residue side chains for each of the 100 ligand poses. In this process we did not replace the W5.43 residue with alanine because we consider this tryptophan to be critical for ligand interaction based on site-directed mutagenesis data,[26] leaving its side chain in the form predicted by SCREAM for each of the ten protein conformations.

In the pose generation step of DarwinDock, a ligand pose is acceptable if it clashed or bumped the receptor residues at six positions or less. First, we used Dock6[55] to generate 5,000 ligand poses (without evaluating an energy) and clustered them into families, where every family member is within a 2.00 Å RMSD of each other. Then, we added 5,000 more ligand poses from Dock6[55] and reclustered. This procedure of adding 5,000 poses and reclustering is repeated until the number of new families generated is less than 2% of the total number of families in the preceding iteration. Typically, 45,000 poses were generated leading to 6000–9000 2.0 Å families. At this point DarwinDock scores the energies of one representative from each family, the family head, and selects the 10% of family heads with the lowest Dreiding energies.[56] All members of these respective families are then scored energetically. From this list of approximately 5000 poses, we select the lowest 50 by each of three criteria: lowest hydrophobic energy, lowest polar energy, and lowest total energy, giving at most 150 poses.

GenDock

GenDock was used to refine the 150 docked ligand–receptor poses generated by DarwinDock. In the SCREAM step,[54] “alanized” residues were replaced with the original hydrophobic residues but using the optimum side chain rotamers to avoid any clashes with the ligand and other protein side chains. Then the entire complex was minimized for 10 steps for each case to remove any bad contacts. Next, the receptor was neutralized, so that the acidic residues (aspartic acid and glutamic acid) each gained a proton, and the basic residues (lysine and arginine) each lost a proton. The resulting receptor–ligand complexes were minimized for 60 steps using the Dreiding III FF.[56] Then for each of the ten receptor structures and each of the ten ligand conformations (100 complexes total), we selected the lowest binding energy structure, including strain and ligand solvation, which we expected to best represent the binding affinity. The binding energy with strain and ligand solvation is defined as the energy difference between the complex and the sum of the receptor and ligand energies with ligand strain and ligand solvation included. The selected complex was minimized with the Dreiding III FF[56] with the LJ vdW term (Dreiding III-LJ FF) in vacuum for 50 steps or to an RMS force threshold of 0.5 kcal/mol/Å using the MPSim program.[53]

Validation of Predicted Binding Sites

Building and Docking Rimonabant Derivatives for SAR Studies

We took the most stable predicted rimonabant-CB1 receptor complex shown in Figure and extracted the rimonabant. Then, we modified it with Maestro software to look like the new derivative, calculated the Mulliken charges at the B3LYP DFT level (with Jaguar software), minimized this configuration with 100 steps of conjugate gradients (or to a threshold of 0.2 kcal/mol/Å) using the MPSim program, and predicted the binding site using DarwinDock as described above. However, instead of docking the new ligands to all ten CB1 conformations, we docked only to the alanized CB1 receptor conformation wild-type 6 (WT6) because it produced the lowest energy (most stable) complex with rimonabant shown in Figure . We also docked a single ligand conformation derived from the docked rimonabant pose. We repeated these steps for seven derivatives and for the original rimonabant ligand as a control. For rimonabant, we did not alter the ligand structure but rather docked the optimized ligand conformation. The final complexes were chosen according to the lowest binding energy including strain and ligand solvation. For comparison of the experimental binding affinities with the computational energies, we used the binding energy, that is, the difference in energy between the receptor–ligand complex and the receptor and the ligand structures calculated separately. The ligand strain and ligand solvation were ignored since we were comparing energetics across different ligands based on the same rimonabant conformation.

Figure 1

(A) The predicted structures for inverse agonist rimonabant to WT6, the sixth best energetic conformation of the CB1 receptor. Rimonabant is anchored to CB1 by hydrogen bonds to W5.43 and K3.28. Hydrogen bond heteroatom distances are indicated in black. (B) Our predicted pharmacophore shows that in our GEnSeMBLE-derived CB1 structure rimonabant forms two hydrogen bonds, including one with W5.43, and have strong aromatic interactions with the receptor.

Free Energy and Binding Affinity Calculations

The experimental change in free energy upon ligand binding to the receptor, ΔG, is obtained from the pKi[57] in eq :

The pKi is the experimental binding equilibrium constant of the inhibitory ligand. Since the pKi describes the strength of the interaction between the receptor and ligand, we use it to derive ΔG. The ΔΔGExp, or the difference in the change in experimental free energy upon ligand binding for a given inverse agonist (ΔG) with respect to the corresponding value of rimonabant, is determined by eq :

This expression was used to compare experimental changes in binding energy with our predicted changes in binding energy for the series of ligands.

Computational Discovery of Novel CB1 Ligands Using the Predicted Ligand Binding Site To Suggest New Inverse Agonists

Our predicted rimonabant-CB1 binding site was used to identify new potential CB1 inverse agonists. The predicted pharmacophore (shown in Figure B), was used in a preliminary search over the 2 million compounds in PubChem[58] for ligands that are similar to rimonabant according to the Tanimoto coefficient,[58,59] an indicator of how similar two 2D ligands structures are to one another. The Tanimoto coefficient ranges from 0 to 1, with 0 being no resemblance between the molecules and 1 being identical molecules. We wanted ligands that are similar to rimonabant, but which could exploit nearby underused polar and aromatic residues in the predicted binding site. For example, K7.32 is located near rimonabant but does not form a hydrogen bond or salt bridge with it; thus we searched for ligands similar to rimonabant, with an appropriately placed functional group to create the new polar interaction. We also wanted a ligand that is commercially available but has not been previously tested with CB1, so that our prediction can be tested easily. Once we had identified a ligand (MSC1, described later) from PubChem that met the above criteria, we performed a search in PubChem to identify ligands similar to MSC1 (scoring >0.90 on the Tanimoto coefficient) or that have a 2D structure that is 90% similar to the new PubChem ligand MSC1 and satisfy specific constraints coming from the CB1 binding site. The 16 ligands found by using PubChem in this protocol are denoted by the acronym “MSC” followed by a number (MSC1–MSC16).

CB1 Expression and Membrane Preparation

HEK293 cells were maintained in Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum and 3.5% mg/mL glucose at 37 °C in 5% CO2. One day prior to transfection, cells were seeded at approximately 1 million cells/100 mm dishes. The cells were transiently transfected by the calcium phosphate precipitation method.[60] At 24 h post-transfection, the cells were harvested in phosphate buffered saline (PBS) containing mammalian protease inhibitor cocktail ((4–2-aminoethyl)benzene-sulfonyl fluoride, pepstatin A, E-64, bestatin, leupeptin, and aprotinin) (Sigma-Aldrich, St. Louis, MO) and lysed by nitrogen cavitation at 750 psi for 5 min. The lysate was spun at 500 g for 10 min at 4 °C, and the supernatant was subsequently spun at 100,000 g for 45 min at 4 °C. The membrane-containing pellet was resuspended in TME buffer (25 mM Tris-HCl, 5 mM MgCl2, and 1 mM EDTA, pH 7.4) containing 7% w/v sucrose.

Radioligand Binding Assay

Competition binding assays were performed as described previously[61] to determine the binding affinity of the test compounds to the receptor. Briefly, 6 μg of membrane preparation was incubated for 60 min in TME buffer containing 0.1% fatty acid-free BSA with a fixed concentration of tracer [3H]SR141716A (43 Ci/mmol, PerkinElmer Life Sciences (Boston, MA)) typically at its Kd which was determined from a saturation binding isotherm. At least nine concentrations of the unlabeled test compound (ranging between 100 pM and 32 μM) were used for the binding assays. Nonspecific binding was determined in the presence of unlabeled SR141716A (1 μM). The reaction was terminated by filtration with a Brandell cell harvester through Whatman GF/C filter paper, and the radioactivity was measured.

GTPγS Binding Assay

GTPγS binding assays were performed as described previously.[62] Briefly, 6 μg of membranes was incubated for 60 min at 30 °C in GTPγS binding assay buffer (50 mM Tris-HCl, pH 7.4, 3 mM MgCl2, 0.2 mM EGTA, and 100 mM NaCl) with the unlabeled 2-arachidonoylglycerol (2-AG) (at least nine different concentrations were used ranging between 100 pM and 1 μM), 0.1 nM [35S]GTPγS (1250 Ci/mmol; PerkinElmer Life Sciences, Boston, MA), 10 μM GDP (Sigma, St. Louis, MO), and 0.1% (w/v) BSA in the absence and presence of 10 μM test compounds. The effect of the compound on inhibiting the level of basal GTPγS binding was evaluated in the absence of agonist. Nonspecific binding was determined with 10 μM unlabeled GTPγS (Sigma, St. Louis, MO). After rapid filtration through Whatman GF/C filters the radioactivity trapped in the filters was determined by liquid scintillation counting.

Experimental Ligand and GTPγS Binding Data Analysis

Data are presented as the mean ± SE or the mean with the corresponding 95% confidence limits from at least three independent experiments. The Ki values of the test compounds were calculated by nonlinear regression using Prism 6.0 (Graphpad Software Inc., San Diego, CA).[61]

Confocal Microscopy

HEK293 cells expressing the CB1 receptor C-terminally fused to GFP were seeded onto 35 mm glass-bottomed dishes (MatTek, MA) precoated with poly-d-lysine. Cells were treated with 10 μM MSC compounds or 1 μM rimonabant for various lengths of time and then washed three times with PBS, followed by fixation with 4% paraformaldehyde for 10 min at room temperature. Cells were mounted in Vectashield mounting medium (Vector Laboratories, CA) and visualized using a Leica TCS SP2 confocal microscope (Leica Microsystems, Wetzler, Germany). Images were collected from at least three independently transfected cell dishes and processed for presentation in figures using Adobe Photoshop 6.0 (Adobe Systems, San Jose, CA).

Results and Discussion

Predicting the Binding Site of the Inverse Agonist Rimonabant to CB1

The predicted lowest energy pose for the CB1-rimonabant complex is shown in Figure A, and the corresponding pharmacophore is shown in Figure B. We predict that rimonabant is anchored by hydrogen bonds to W5.43 and K3.28 [Residues are numbered according to the Ballesteros–Weinstein scheme.[63]], and indeed previous site-directed mutagenesis data indicate that mutations of W5.43 and K3.28 to alanine have the largest effect in decreasing experimental binding affinity upon mutation to alanine with >1000-fold and 17.2-fold, respectively.[26,28] The components of the predicted binding energies for this complex are shown in Table , showing contributions from each residue in the binding site. The most important residue is W5.43 with a predicted binding contribution to rimonabant of −9.36 kcal/mol, consistent with the decrease by a factor of >1000 upon mutation to alanine.[26] Indeed our predicted pharmacophore for rimonabant (Figure B) has the significant polar and hydrophobic contacts with the CB1 receptor expected for a strongly bonding ligand.

Table 1

Cavity Analyses for the Predicted Complex of Rimonabant to WT CB1

residuea	no.b	VdWc	Coulombc	H-bondc	NonBondc,d
TRPe	5.43	–3.146	–2.085	–4.129	–9.360
MET	6.55	–6.043	–0.756	0	–6.799
LYS	3.28	–0.833	–1.410	–3.322	–5.566
PHE	3.36	–3.841	–0.054	0	–3.896
VAL	3.32	–2.904	–0.363	0	–3.267
THR	5.47	–1.912	–0.048	0	–1.960
LEU	6.51	–2.100	0.193	0	–1.907
TRP	6.48	–1.844	0.008	0	–1.837
THR	3.33	–1.889	0.093	0	–1.797
ILE	2.56	–1.449	0.000	0	–1.449
LEU	7.43	–1.382	0.119	0	–1.263
PHE	7.35	–0.959	–0.174	0	–1.134
SER	3.35	–1.385	0.302	0	–1.083

Only residues with energies with contributions stronger than 1.0 kcal/mol are shown.

Residues are numbered according to the Ballesteros–Weinstein scheme.[63]

All energies are in units of kcal/mol.

The nonbonding (NonBond) energy in the far right column is the sum of the van der Waals (VdW) energy, the electrostatic (Coulomb) energy, and hydrogen bond (H-bond) energy.

Residues forming hydrogen bonds are highlighted in bold.

According to the cavity analysis in Table , K3.28 has the third largest binding interaction energy (−5.57 kcal/mol) with rimonabant due to its hydrogen bond with the amide carbonyl of the ligand. This interaction is also supported by a mutagenesis study.[28] The same experiments with the rimonabant derivative VCHSR, an analogue of rimonabant with hydrocarbons replacing the amide and piperidine groups, showed that the binding affinity remained unaffected upon mutation of K3.28 to alanine. Our study shows that K3.28 interacts with the polar atoms of the amide group or piperidine ring, including the carbonyl, which agrees with the mutagenesis study on this ligand that lacks these polar atoms.[28]

In our docked pose of rimonabant with the WT receptor, the F3.36 residue has a significant interaction with the ligand, including a van der Waals component of −3.90 kcal/mol. In binding assays, the F3.36A mutation decreased the binding affinity of this inverse agonist by 20-fold,[27] or 1.78 kcal/mol, for the CB1 receptor. Interestingly, the F3.36L mutation decreased the binding affinity by only 2-fold, indicating the importance of the bulky hydrophobic residue for receptor–ligand interaction.[27]

We predicted that M6.55 has the second strongest interaction (−6.8 kcal/mol) with CB1, which is almost entirely hydrophobic (Table ). However, experiments show that mutating M6.55 to alanine decreases the binding affinity to CB1 by only 3-fold[35] or by 0.65 kcal/mol in energy. One possible explanation is that the M6.55A mutation alters the binding site and pose for rimonabant in a way that other residues become available to compensate for the lost interaction with methionine. We observed similar cases in our studies of binding of ligands to CCR5.[64]

Analysis of the Conformational Ensemble of CB1 Structures

A very important advantage of the GEnSeMBLE approach is that we have an ensemble of ten low energy structures to which rimonabant is allowed to bind. We showed for the adenosine A3 receptor that neither the selective agonists nor the selective antagonists prefer the lowest energy apoprotein structures.[24] Also, for CCR5 we showed that the antagonists all prefer to bind to structures that are not the lowest energy for the apoprotein.[64] Similarly for CB1, the most favorable binding is with WT6, the sixth best conformation for the apo-CB1 structure, not with the lowest energy apo-CB1 conformation (WT1) (Table ). Since the ligands do not preferentially bind to the lowest energy conformation, these docking results support the concept that GPCRs are highly dynamic structures that can sample many conformations, providing an opportunity for different ligands that can selectively bind to different conformations, and perhaps lead to different function. The ten CB1 conformations used in this study range in backbone RMSD from 0.4–2.1 Å. Of course the protein conformation may change when complexed with the ligand.

Table 2

Comparison of Receptor Conformations Selected by the Inverse Agonist Rimonabant

energy rank	WT conf no. bound to rimonabantb	binding energy with strain and rimonabant solvationa
1	WT6	–59.10
2	WT1	–57.24
3	WT9	–55.44
4	WT5	–54.39
5	WT2	–50.97
6	WT7	–50.12
7	WT4	–50.05
8	WT10	–49.59
9	WT8	–49.35
10	WT3	–48.63

The complexes are ranked according to best binding energy with strain and ligand solvation included.

Receptor conformation numbers that do not contain the conserved R3.50 and D6.30 salt bridge are in bold.

Surprisingly, the WT6 CB1 structure, to which the inverse agonist rimonabant preferentially binds, lacks the R3.50 and D6.30 ionic lock, which is generally believed to be important for preventing activation. Experimental data have shown that rimonabant inactivates the CB1 receptor,[14,37] while our calculations find that the inverse agonist has the best binding energies with activated CB1 conformations. Indeed of the top four rimonabant-CB1 complexes, we find that three conformations (WT6, WT1, and WT5) do not form an ionic lock to the conserved R3.50 and D6.30.

On the other hand, several studies have suggested that a R3.50 to D6.30 salt bridge is not necessary for maintaining the inactive conformation. For example, Audet and Bouvier[65] showed that the R3.50 and D6.30 salt bridge is broken in all but three of the 36 crystallized GPCR structures bound to antagonists and inverse agonists analyzed. It should be noted that many crystallized GPCR structures, such as the human adenosine A2A receptor, contain an inserted T4 lysozyme for stabilization, which may impact the conformation of the cytocellular ends of the helices and thus prevent the formation of the ionic lock.[45]

In our previous work with the human CB1 receptor, we showed that a different residue, R2.37, plays an important role in preventing receptor activation.[16,17,41] We found that this residue forms an “ionic lock” via a salt bridge interaction with D6.30 to stabilize the inactive conformation. Thus, the traditional R3.50 and D6.30 salt bridge was not essential. Indeed our calculations found that the rimonabant inverse agonist did not prefer to bind to the conformations that contained this signature contact as shown in Table . The rimonabant-bound CB1 conformation WT6 resembles the predicted structure of a constitutively active mutant, T3.46A/R2.37A, which also lacks the R3.50 and D6.30 ionic lock.[16,17] The two structures have a Cα-RMSD of 1.2 Å, which is smaller than most crystal structure resolutions[10] suggesting that the two conformations are very similar. For comparison, the difference between the WT receptor conformations WT2 (which contains the R3.50+D6.30 salt bridge) and WT6 (without that salt bridge) is 1.8 Å.

Analysis of Structure–Activity Relationship (SAR) with Rimonabant Derivatives

To validate our proposed binding site, we docked seven known derivatives of rimonabant (Figure A) to the CB1 conformation and plotted the calculated binding energies against the experimental binding affinity, pKi, to find the correlation between the two data sets (Figure ). For each ligand, we found the same binding site as for rimonabant (Figure B). Figure A shows a comparison between the pKi of the inverse agonists from Hurst et al.[36] and our calculated binding energies. We observe a very high correlation of 93.4%, indicating good agreement between the trends observed in experimental binding affinities and our own calculated binding energies. This strongly supports the selected binding site.

Figure 2

Rimonabant derivatives used in SAR study. (A) Rimonabant and the seven derivatives used in the SAR study. The portions of the ligand that are different from rimonabant are colored in red. The derivatives in the first two rows were used in the Hurst et al., 2006 study,[36] and the derivatives in the third row were used in the D’Antona et al., 2006 studies.[37,66] (B) Rimonabant and seven derivatives docked to the CB1 complex. Rimonabant has carbon atoms colored in cyan. The other derivatives have carbon atoms colored in gray.

Figure 3

Comparison of experimental binding affinities and computationally calculated energies between receptor and inverse agonists. (A) Plot showing the correlation of the experimental binding affinity, pKi, for the six inverse agonists in Hurst et al.[36] versus our calculated binding energy (red circles). This leads to R2 = 0.93, indicating that the error bar for our calculations is ∼2 kcal/mol, while the range is 18 kcal/mol. (B) Plot showing the correlation of the ΔΔGexp for seven previously published inverse agonists with respect to that of rimonabant[36,37] versus the calculated binding energy. R2 = 0.89, indicating that the error bar for our calculations is ∼2 kcal/mol, while the range is 14 kcal/mol.

For the inverse agonists from D’Antona et al.,[37,66] we compared these calculated energies with the ΔG determined from eq . The ΔG values are within 1 kcal/mol of each other since the energies are −9.33 kcal/mol for AM-281, −10.0 kcal/mol for AM-251, and −9.70 kcal/mol for rimonabant. Similarly, both the calculated cavity and binding energies for the three ligands are approximately within 1 kcal/mol. The binding energies are −67.4 to −69.0 kcal/mol, which are very similar and agree with the experimentally observed negligible changes in free energy upon binding. This excellent correlation relies on the predicted energies that are based purely on enthalpy and do not include entropy. However, the entropy change upon binding of all the ligands should be similar since they all have the same number of rotatable bonds.

Figure B combines the results from Figure A with those regarding the inverse agonists from D’Antona et al.[37] for a single comparison of our calculated binding energies and the binding affinities from the two experiments. The respective ΔΔGExp values from eq are plotted against the calculated binding energies for all eight ligands, leading to an 89.4% correlation. The excellent agreement between our predicted energies and the experimental binding affinity shown in Figure B provides confidence that our predicted binding site for rimonabant is reasonable.

Discovery of Novel CB1-Targeting Inverse Agonists

With a reliable rimonabant binding pharmacophore for CB1, which is supported by site-directed mutagenesis and SAR data, we aimed to design a new, more selective CB1 inverse agonist. After the informed search in PubChem,[58] described in the Methods section, we selected Zinc08587042, which we refer to as MSC1 (Figure ), a ligand with a 68% similarity to rimonabant according to the 2D Tanimoto coefficient.[67] Rather than relying exclusively on PubChem’s search algorithm, which looks for molecules based on a 0.90 Tanimoto similarity score with respect to rimonabant, we used the predicted CB1 binding site and searched for a ligand that maintained key structural aspects of rimonabant and additionally took advantage of underused residues in the binding site. For example, we wanted to find a ligand having the same amide and pyrazole functional groups to maintain the hydrogen bonds with K3.28 and W5.43, respectively. Yet, we wanted to replace the rimonabant piperidine ring with a phenyl ring and a substituted polar group to improve the interactions with the aromatic residues that we found nearby and to gain an additional hydrogen bond with K7.32, which is within the rimonabant binding site. Furthermore, we wanted to replace the methyl group of the pyrazole group with a polar substitute to form a hydrogen bond with S7.39. MSC1 is attractive because it contains an acetylphenyl group, which we predicted would reach into the aromatic pocket of the extracellular end, and a triazole ring, which would replace a methyl group with a nitrogen atom. We also wanted a small molecule that was commercially available and that had not been previously tested with CB1 in bioactivity assays. MSC1 met all of the above pretesting criteria.

Figure 4

Proposed CB1 inverse agonist MSC1 (Zinc08587042). Structure comparison of rimonabant (top) with MSC1 (bottom). MSC1 has a 2D structure that is 68% similar to rimonabant. Portions of MSC1 highlighted in red are different from the corresponding groups in rimonabant.

PubChem also identified 15 other small molecule ligands that have a 2D structural similarity Tanimoto score of 0.90 with MSC1 or have 2D structures that are 90% similar to that of MSC1 (Figure ). We docked these 16 ligands to the ensemble of predicted CB1 structures, found their corresponding binding energies, and predicted their respective pKi values based on Figure . The predicted binding energies are given in Supporting Information (SI) Table S1. Of the 16 MSC small molecules, 14 compounds, MSC1– MSC6, MSC8–MSC11, and MSC13–MSC16 were commercially available.

Figure 5

Proposed CB1-selective inverse agonists based on PubChem[58] similarity search with MSC1. PubChem identified 15 ligands that have 2D structures similar to MSC1.

Experimental Binding Affinity Data Suggest MSC1 and MSC3 Are CB1 Ligand Candidates

From the 16 compounds identified from the computational analyses, we chose five compounds based on the 2D structural similarity, Tanimoto score, and commercial availability. These are MSC1, MSC3, MSC5, MSC8, and MSC9. To evaluate experimentally the binding affinities of the compounds, we performed competition binding experiments using [3H] rimonabant as a tracer. MSC1, MSC3, and MSC9 bound the receptor with Ki values of 502 nM, 495 nM, and 4619 nM, respectively (Table ). These data indicate that while removal of the chloro group from the chlorophenyl ring of MSC1 does not impact receptor binding (MSC1 versus MSC3), the position of the methyl group in one of the biphenyl rings and the acetyl group in the amide-linked phenyl ring is critical for receptor binding (MSC3 versus MSC9). MSC5 and MSC8 failed to bind CB1 up to 32 μM (Table ). Unexpectedly, our experimental Ki values suggested a weaker affinity than rimonabant compared to the predictions. This may be due to the limitations of force fields to reliably predict relative binding affinities.

Table 3

Experiment Binding Results for Rimonabant and the MSC Compound Bound to the WT CB1 Receptor

compound	Kia (nM)
rimonabant	3.4 (2.7–4.2)
MSC1	502 (306–825)
MSC3	496 (211–1164)
MSC5	no detectable binding
MSC8	no detectable binding
MSC9	4619 (528–40380)

Ki values were determined from competition binding assays using [3H]rimonabant as tracer. Data are the median and corresponding 95% confidence limits of three independent experiments performed in duplicate.

MSC1 and MSC3 Inhibit Both Basal and 2-AG-Induced G Protein Coupling and Enhanced Cell Surface Expression of CB1

Since the experimental binding affinities to the CB1 receptor were higher for MSC1 and MSC3 than for MSC5, MSC8, and MSC9, we further tested MSC1 and MSC3 using [35S]GTPγS binding assays to evaluate changes in G-protein coupling. We assessed the effect of the compounds on the basal and the 2-AG-induced GTPγS binding. We chose 2-AG, the endogenous eicosanoid CB1 agonist for the assay since it showed binding affinity comparable to MSC1 and MSC3. Thus, MSC1 and MSC3 can potentially compete with 2-AG.

The addition of 10 μM of MSC1 and MSC3 resulted in substantial reduction in Emax values with 90 fmol/mg and 81 fmol/mg, respectively, compared with that in the absence of these MSC compounds (Emax = 126 fmol/mg) (Figure A). These compounds also exhibited small but not statistically significant shifts in EC50 values (785 nM for MSC1, 770 nM for MSC3, compared to 521 nM for the absence of any MSC compound). Thus, the extent of G protein coupling (Emax values) changes, though the amount of agonist needed to achieve this level of GTPγS binding (EC50 values) does not. Furthermore, Figure B showed that both MSC1 and MSC3 substantially decreased the basal [35S]GTPγS binding (67 fmol/mg) with specific binding of 50 fmol/mg and 47 fmol/mg, respectively.

Figure 6

Effect of MSC1 and MSC3 on the [35S]GTPγS binding to HEK293 cell membranes expressing the CB1 WT receptor. (A) Dose–response curves for 2-AG-induced [35S]GTPγS binding in the absence (●) or presence of 10 μM (■) MSC1 and MSC3. (B) Inhibition of basal [35S]GTPγS binding by MSC1 and MSC3. The levels of untransfected cells (no CB1) and rimonabant treated samples are shown for comparison. Data are presented as specific binding of [35S]GTPγS to the membrane. Nonspecific binding was determined in the presence of 10 μM unlabeled GTPγS. All data are the mean ± SE of at least three independent experiments performed in duplicate. (C) Cellular distribution of the WT receptor upon treatment with MSC compounds. HEK293 cells expressing the CB1 WT-GFP receptor were incubated with vehicle alone (0.03% DMSO), MSC1 (10 μM) or MSC3 (10 μM) for 6 h. Rimonabant-treated (1 μM) cell is shown for comparison. Scale bar, 10 μm.

Whereas prolonged agonist treatment can remove receptors from the cell surface by endocytosis, inverse agonists have been shown to enhance the cell surface localization, consistent with receptor inactivation.[68,69] CB1 exhibits some constitutive activity, and the majority of CB1 receptors are localized on intracellular vesicles even in the absence of agonist.[61,70,71] We tested if the MSC compounds can affect cellular localization of the receptor using confocal microscopy of cells expressing GFP-tagged CB1. Similarly to rimonabant treatment, upon treatment with 10 μM MSC1 or MSC3 the level of cell surface localization of the receptor was enhanced though the extent for MSC3 is less than MSC1 (Figure C). Collectively, although MSC1 and MSC3 exhibited somewhat lower potency and efficacy compared to rimonabant, these data suggest that they possess inverse agonist properties as evidenced by reduction in G-protein coupling and enhancement of CB1 cell surface expression. We took the MSC1-bound CB1 structure and relaxed it in an explicit lipid bilayer environment using 50 ns of molecular dynamics (MD) simulations. These MD simulations are performed under isothermal–isobaric conditions (NPT) at 310 K and 1 atm, using periodic boundary conditions with Particle-Mesh Ewald summation method for calculating long-range coulomb interactions. The temperature is controlled using Langevin dynamics, and the pressure is maintained using a Langevin-Hoover barostat. These simulations showed that binding of MSC1 to the CB1 receptor constrains TM helix 6 more through additional intracellular couplings preventing the G-protein from engaging the receptor, which is consistent with experimental findings of MSC1 ligand’s inverse agonism. These results offer opportunities for developing new CB1 inverse agonists through the optimization of chemical functionalities of the compounds and further pharmacological characterization.

Comparison of Our Predicted Binding Site to Previous Predictions

Our predicted rimonabant binding site in CB1 is different from previous computational studies that used the bovine rhodopsin template to create a CB1 homology structure for docking rimonabant.[26,28,72,73]Supporting Information Figure S1 compares the rimonabant pharmacophore published by Lange and Kruse[74] (Figure S1A) with the pharmacophore from our current study (Figure S1B). Both agree that rimonabant (Figure S1A) spans the width of the binding site with the chlorophenyl rings near TM5 and the piperidine ring near TM3. However, we predict that rimonabant (Figure S1B) lies parallel to the z-axis passing from the extracellular to intracellular side of the membrane, which is perpendicular to the previous models (Figure S1A). Although both agree that a hydrogen bond forms between K3.28 and the amide carbonyl of rimonabant, the earlier models suggest that K3.28 also participates in a salt bridge with D6.58,[26,28,72] which we do not observe. In addition, the previous models have the W5.43 residue sandwiched by the two chlorophenyl groups of rimonabant, leading to strong aromatic interactions. Furthermore, multiple other aromatic residues in the binding pocket (F3.25, F3.36, W4.64, Y5.39, W6.48) were found to participate in stacking with the ligand,[26,72,73] whereas in our predicted binding site, it is W5.43 that has the biggest impact on binding affinity, from the hydrogen bond with the pyrazole ring of rimonabant. In our predicted structure, we find that the F3.36 residue is sandwiched by the chlorophenyl groups rather than W5.43. In addition, a paper from Shim et al. based its CB1 model on the human ß2 adrenergic receptor[35] to identify a binding site similar to the one from Lange and Krause[74] but with the ligand lying perpendicular to our predicted structure, with the piperidine ring pointing to the TM1–2–7 binding pocket and the chlorophenyl rings pointing toward TM5, where they form aromatic stacks with W5.43 and W6.48.

One would expect our predicted rimonabant binding site to be different from the previously published ones since our predicted receptor structures are quite different. Specifically, we sample all rotations of all helices about their axes and helix tilt and sweep angles, a total of nearly 13 trillion receptor conformations. In contrast, the previous studies used homology models based on the bovine rhodopsin or human ß2 adrenergic receptor structure templates, so that no rotations or tilting of the helices were conducted. These three GPCRs have very similar global conformations, but we performed considerable sampling to obtain CB1 receptor conformations that were significantly different and more stable than the starting template.[16,17,41] The method we used to identify our final receptor-inverse agonist complex is also different. The study from Shim et al.[35] used the selection criterion that one of the two chlorophenyl rings should be near the TM3–5–6 hydrophobic site, so their final pose was chosen based on the interactions of rimonabant with residues for which there was site-directed mutagenesis data. However, the structures in our study were based on the binding energy of the complexes, with no extra conditions imposed. The previously suggested rimonabant binding pose most similar to ours came from another study in which the ligand was docked to a bovine rhodopsin structure-based homology model.[75] This predicted binding site is similar to the current model, with rimonabant sitting vertically in the binding pocket—the piperidine ring points toward the extracellular end and the chlorophenyl rings point toward the intracellular end. However, it does not lead to any hydrogen bonds with W5.43 and finds a hydrogen bond with Y5.39, that we did not observe.

Conclusions

Since many antiobesity drugs including rimonabant have been suspended from the market, there remains an enormous unmet need for compounds that reduce food intake. In an effort to develop novel inverse agonists for CB1, we predicted the binding sites and energies for rimonabant and structurally related molecules to our previously predicted and validated CB1 receptor conformations. Our calculated binding energies for the compounds are in good agreement with the previously reported mutagenesis[26,28] and SAR data.[36,37] Based on the predicted binding site of these inverse agonists, we identified 16 new potential CB1 inverse agonists. We tested five of these experimentally and found that two of them (MSC1 and MSC3) antagonized G-protein coupling and enhanced surface localization of the receptor suggesting they are inverse agonists. The MD simulations of MSC1:CB1 complex in explicit lipid bilayer environment show that the relaxed receptor exhibits additional intracellular couplings that will prevent G-protein coupling, which is consistent with the characterization of MSC1 as an inverse agonist. Although MSC1 and MSC3 are less potent than rimonabant, they provide new starting points for developing new candidates for antiobesity drugs that target CB1.