pubs.acs.org/jcim Article A Single Point Mutation Blocks the Entrance of Ligands to the Cannabinoid CB2 Receptor via the Lipid Bilayer Nil Casajuana-Martin, Gemma Navarro, Angel Gonzalez, Claudia Llinas del Torrent, Marc Gómez-Autet, Aleix Quintana García, Rafael Franco, and Leonardo Pardo* Cite This: J. Chem. Inf. Model. 2022, 62, 5771−5779 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: Molecular dynamic (MD) simulations have become a common tool to study the pathway of ligand entry to the orthosteric binding site of G protein-coupled receptors. Here, we have combined MD simulations and site-directed mutagenesis to study the binding process of the potent JWH-133 agonist to the cannabinoid CB2 receptor (CB2R). In CB2R, the N-terminus and extracellular loop 2 fold over the ligand binding pocket, blocking access to the binding cavity from the extracellular environment. We, thus, hypothesized that the binding pathway is a multistage process consisting of the hydrophobic ligand diffusing in the lipid bilayer to contact a lipid-facing vestibule, from which the ligand enters an allosteric site inside the transmembrane bundle through a tunnel formed between TMs 1 and 7 and finally moving from the allosteric to the orthosteric binding cavity. This pathway was experimentally validated by the Ala2827.36Phe mutation that blocks the entrance of the ligand, as JWH-133 was not able to decrease the forskolin-induced cAMP levels in cells expressing the mutant receptor. This proposed ligand entry pathway defines transient binding sites that are potential cavities for the design of synthetic modulators. 1. INTRODUCTION However, the depth of ligand penetration into the TM bundle The approximately 800 (∼450 for sensory and ∼350 for is considerable,14 so ligands must also engage stable nonsensory functions) members of the G protein-coupled interactions at the entrance of the binding site. Ligand binding receptor (GPCR) family are mainly formed by a conserved pathway simulations on β-adrenergic receptors have identified architecture of seven transmembrane domains (TMs).1 The these additional cavities that are transiently occupied before intracellular domains and/or the intracellular part of the TM arrival to the orthosteric binding site.15,16 They have been segments bind a small repertoire of signaling proteins (16 G named extracellular vestibules15 or entrances,17 or secon- proteins, four arrestins, or seven G protein-coupled receptor dary16,18 or metastable19 binding sites, or exosites.20 Cavities kinases),2−4 whereas the extracellular domains (ECDs) and/or like those described for β-adrenergic receptors have also been found for the M3 muscarinic receptor,21 the adenosine A2A the extracellular part of the TM segments recognize exogenous receptor,22 the histamine H4 receptor,23 the dopamine D2 and signals (e.g., 1 trillion olfactory5 or ∼1 000 bitter molecules6), D3 receptors,24 and olfactory receptors,25 among others. or endogenous signals that range in size from ions (e.g., Ca2+) Remarkably, while the orthosteric binding site is conserved to small molecules (e.g., glutamate, neurotransmitters) to among family members, these sites at the ECDs are highly peptides (e.g., opioids, endothelin, glucagon) or large proteins divergent and have become potential binding sites for synthetic (e.g., chemokines, glycoprotein hormones), or ∼200 000 modulators.26 Thus, they have been used in the design of synthetic ligands,7 or ∼475 drugs.8 Thus, the mechanism of allosteric modulators27 or bitopic ligands,20,28 as well as binding of this set of highly diverse extracellular ligands varies proposed to bind short bivalent ligands.29,30 considerably depending on the ligand and receptor,9 as Class A GPCRs activated by hormone-like signaling subsequently confirmed by the analysis of known inactive molecules derived from lipid species with long hydrophobic and active structures of GPCRs.10 The high druggability of membrane embedded GPCRs11 is due to the fact that their malfunction, with 435 disease- Received: July 15, 2022 associated mutations,12 commonly translates into pathological Published: October 27, 2022 outcomes13 and because the orthosteric binding site, a conserved pocket within the 7TM domain that optimally accommodates the electrostatic and steric properties of the ligand, is normally accessible from the extracellular space. © 2022 The Authors. Published by American Chemical Society https://doi.org/10.1021/acs.jcim.2c00865 5771 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 1. Ligand diffusion to CB2R. (a) Cluster of tunnels, as calculated with CAVER,60 from the orthosteric cavity toward the extracellular domain (blue) and toward the lipid bilayer via either TMs 2 and 3 (orange) or TMs 1 and 7 (green). The probabilities of ligand transportation (throughput) of each tunnel, calculated from the ensemble of structures collected during the aggregated 100 μs of MD simulation, are shown. (b) In three of the 100 MD simulations, the ligand diffused from the lipid bilayer to the receptor surface (10 snapshots per simulation are shown) and remained bound to either TMs 5 and 6 (one simulation in green) or TMs 1 and 7 (two simulations in orange and yellow). (c−f) Comparison of the trajectories (Figure S3) in which the ligand spontaneously bound the lipid-facing part of TMs 1 and 7 (orange in panel b, surface-bound) and a trajectory in which the ligand did not bind the receptor surface (surface-unbound). Detailed view of key amino acids located at the entrance of the TMs 1 and 7 tunnel (c). Evolution of the Phe2837.37 side chain and the χ1 dihedral angle along the surface-unbound (left) and surface-bound (right) trajectories (d). Representative structures obtained during the MD simulations (Figure S3) in which the cavity is closed (in white) or open (in blue). The conformations of Phe2837.37 in panel d and TM 1 in panel f that resemble these closed or open structures are shown in white or blue, respectively. Evolution of TM 1 and the distance between the top of TMs 1 (Cα atom of Thr341.33) and 7 (Cα atom of Val2777.31) along the surface-unbound (left) and surface-bound (right) trajectories (f). moieties31 possess a distinctive structural signature at the ECD, 2. MATERIALS AND METHODS in comparison with other class A GPCRs that are activated by 2.1. Initial CB2R Models. The inactive structure of CB2R polar ligands.32 In the crystal structures of S1P1,33 LPA1,34 was used (PDB id 5ZTY).39 Protonation states were assigned FFAR1,35 CRTH2,36 CB1,37,38 and CB239 receptors, the with PROPKA.48 The system was oriented by the Orientations extracellular N-terminus and extracellular loop 2 fold over of Proteins in Membranes (OPM) database49 and embedded the ligand binding pocket. Thus, the entrance of ligands to the in four different lipid bilayer boxes (9 × 9 × 9.5 nm), orthosteric binding site through the membrane bilayer has constructed using PACKMOL-memgen,50 containing 1- been proposed for S1P1,33,40 CRTH2,36 and CB137,41 receptors palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), cho- lesterol (CHL; 10:6 POPC/CHL ratio), water molecules via TMs 1 and 7, for the FFAR1 via TMs 3 and 4,35 for the (TIP3P), and monatomic Na+ and Cl− ions (0.15 M). JWH- MT1 melatonin receptor via TMs 4 and 5,42 for opsin via TMs 133 (JWH) was inserted into the membrane by replacing one 5 and 6,43,44 and for CB245 and PAR146 receptors via TMs 6 in six CHL molecules from both leaflets, leading to a 10:5:1 and 7. In contrast, some authors propose that hydrophobic POPC/CHL/JWH ratio. This substitution of CHL by JWH ligands of LPA134 and CB1 and CB247 receptors may reach the was randomly performed five times for each of the four orthosteric site from the extracellular environment. constructed membrane systems (see Figure S1). In this manuscript, we combine molecular dynamic (MD) 2.2. Molecular Dynamics Simulations. These 20 simulations and site-directed mutagenesis to study the binding combinations of CB2R models were simulated with GRO- process of a potent and selective CB2R agonist, JWH-133. The MACS2016.4.51 The amber14sb-ildn force field was used for proposed binding pathway is a multistage process consisting of the protein, solvent. and ions;52 a GROMACS adaptation of lipid14 for lipids;53,54 and the general Amber force field the ligand diffusing in the lipid bilayer to contact a lipid-facing (GAFF2) with HF/6-31G*-derived RESP atomic charges for vestibule, from which the ligand enters an allosteric site inside JWH-133.55 Molecular systems were subjected to 10 000 steps the TM bundle through a tunnel formed between TMs 1 and of energy minimization, followed by 20 ns of gradual relaxation 7, and finally moving from the allosteric to the orthosteric of positional restraints (corresponding to 100, 50, 25, and 10 binding cavity. kJ mol−1 nm−2) at protein backbone coordinates before the 5772 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article production phase to hydrate the cavities and allow lipids to by sequencing to select those expressing the correct DNA pack around the protein. After equilibration, five replicas of 1 sequence. μs unrestrained MD trajectory, of each of the 20 combinations, 2.5. cAMP Determination Assays. Determination of was generated at a constant temperature of 300 K using cAMP levels in HEK-293T cells transiently expressing CB2R (1 separate v-rescale thermostats for the receptor, ligand, lipids, μg of cDNA) or the mutant version of the receptor was and solvent molecules. A time step of 2.0 fs was used for the performed using the Lance-Ultra cAMP kit (PerkinElmer). integration of equations of motions. Bonds involving hydrogen Two hours before initiating the experiment, the medium was atoms were kept frozen using the LINCS algorithm. Lennard- substituted by a serum-free medium. Then, transfected cells Jones interactions were computed using a cutoff of 1.1 nm, and were dispensed in white 384-well microplates at a density of the electrostatic interactions were treated using PME with the 3000 cells per well and incubated for 15 min at room same real-space cutoff under periodic boundary conditions temperature with compounds, followed by 15 min of (see Figure S1). The analysis of the trajectories was performed incubation with forskolin and 1 h more with homogeneous using MDAnalysis;56 the stability of the membrane was time-resolved fluorescence (HTRF) assay reagents. Fluores- monitored (see Figure S2) using FATSLiM.57 Visualization cence at 665 nm was analyzed on a PHERAstar Flagship and image rendering were performed with PyMOL58 and microplate reader equipped with an HTRF optical module VMD,59 and the tunnels of CB2R, from the orthosteric binding (BMG Labtech). Data analysis was made based on the cavity, were calculated with the CAVER program.60 fluorescence ratio emitted by the labeled cAMP probe (665 2.3. Metadynamics Simulations. Metadynamics simu- nm) over the light emitted by the europium cryptate-labeled lations (see Figure S1) were performed using GRO- anti-cAMP antibody (620 nm). A standard curve was used to MACS2016.4 patched with PLUMED version 2.561.62 The calculate cAMP concentration. Forskolin-stimulated cAMP starting point was taken from one of the 100 unbiased MD levels were normalized to 100%. trajectories in which JWH-133 spontaneously bound the lipid- facing part of TMs 1 and 7 (orange in Figure 1b). The end 3. RESULTS point is the docking pose of JWH-133 to the orthosteric 3.1. Computational Model to Study the Pathway of binding site of CB2R previously reported,63 using the Ligand Entry to CB2R. The selected ligand is JWH-133, a structurally similar AM12033 ligand (PDB id 6KPC)47 as a potent selective CB2R agonist.64 We have chosen JWH-133 template. The distance between the center of mass at the over other agonists for its rigidity and for having only one starting and end points of JWH-133 was used as a collective heteroatom in its structure. JWH-133 (6.4) and cholesterol variable (CV). This CV was selected for simplicity and because (8.7) have similar computed logP values, which makes JWH- it can efficiently represent the dimensional space of the free- 133 a potential through-the-membrane binder. The model to energy of binding. In metadynamics simulations, a history- study the entry of JWH-133 into the orthosteric binding site dependent potential is applied along the CV, built as a sum of consisted of inactive CB2R embedded in a lipid bilayer box Gaussian kernels, with a specific width of 0.1, height of 0.48 (see Materials and Methods). Four different membranes were kJ/mol, and pace of 5000 steps (i.e., 0.1 ns). A light RMSD constructed, and for each of them, JWH-133 was positioned at constraint (with an energy constant KAPPA of 200 IU) was five different positions. These 20 combinations were each applied to ensure a correct end conformation of JWH-133 at subjected to five rounds of unbiased 1 μs MD simulation with the orthosteric site. We used well-tempered metadynamics, an aggregate sampling of 100 μs (see Methods and Figure S1). where the height of the applied Gaussian potentials is rescaled Data was collected every 10 ns. so that the system does not explore free-energy regions beyond 3.2. Tunnels in CB2R. To delineate the pathways from the the expected values. Thus, a bias factor of 20, accounting for orthosteric binding cavity of CB2R, buried in the TM bundle, the difference between the system temperature and the to the surrounding solvent, either the aqueous extracellular temperature of the CV, was set up. Convergence of the well- environment or the lipid bilayer, we used the ensemble of tempered metadynamics simulations was determined when structures collected during the MD simulations of CB2R. JWH-133 had explored all of the CV space and the free-energy Figure 1a displays three different clusters of tunnels calculated profile was constant at 10 ns intervals. Free-energy profiles with the skeleton search algorithms implemented in the were calculated using PLUMED sum_hills function. CAVER program.60 These are toward the extracellular domain 2.4. CB2R Mutants. Mutations were produced using the and toward the lipid bilayer via either TMs 2 and 3 or TMs 1 QuikChange Site-Directed Mutagenesis Kit. The cDNA for and 7. The tunnel route preference was evaluated using hCB2R, cloned into pcDNA3.1, was amplified using sense and throughput, which is the probability of ligand transportation antisense primers harboring the triplets for the desired point throughout the pathway (the higher the value, the greater the mutation (Pfu turbo polymerase was used). The forward and importance of the pathway). These values of throughput reverse primers are ccccagaagacagctTTTgctgtgttgtgcactc and (Figure 1a), which consider the conformational flexibility of gagtgcacaacacagcAAAagctgtcttctgggg for Val361.35Phe, cccca- CB2R, suggest that the transit of molecules through the lipidic gaagacagctATGgctgtgttgtgcactc and gagtgcacaacacagcCA- phase via TMs 1 and 7 is favored compared to TMs 2 and 3 or Tagctgtcttctgggg for Val361.35 Met, gtcaa- to the extracellular route. gaaggcctttTTCttctgctccatgctg and cagcatggagcagaa- 3.3. Ligand Diffusion to CB2R and Tunnel Opening. GAAaaaggccttcttgac for Ala2827.36Phe, and None of the JWH-133 agonists, embedded at different gtcaagaaggcctttATGtctgctccatgctg and cagcatggagcagaaCA- positions in the membrane, spontaneously bound the CB2R Taaaggccttcttgac for Ala2827.36Phe, respectively, in which orthosteric binding site during the aggregated 100 μs of MD mutated nucleotides are written in upper case characters. simulation, but in three simulations the ligand remained bound The nonmutated DNA template was digested for 1 h with to the receptor surface, mainly to amino acids in TMs 5 and 6 DpnI. PCR products were used to transform XL1-blue (one simulation) or TMs 1 and 7 (two simulations). Figure 1b supercompetent cells. Finally, positive colonies were tested shows these stable surface-bound ligand conformations with 5773 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article root-mean-square deviation (rmsd) values < 20 Å relative to the reference docked binding mode of JWH-133 (see Materials and Methods), while they have displaced >20 Å from their initial position. These results, together with the throughput values of tunnel route preference (see above), reinforce TMs 1 and 7 as the most feasible pathway of ligand entry to the orthosteric site of CB2R. Figure 1c shows the key Val361.35, Leu391.38, Phe2837.37, and Met2867.40 residues (Ballesteros and Weinstein numbering scheme65 as superscript) located at the entrance of TMs 1 and 7. Of note is the conformation of Phe2837.37, which in the gauche+ conformation, like that observed in the crystal structures of CB2R,39,47,66 closes the tunnel, whereas in the trans conformation, it opens it. Moreover, it has been shown that the N-terminus of the S1P1 receptor packs against ECL-2 in the active conformation, leading to an opening of the ligand access vestibule between TMs 1 and 7.67 Thus, we have quantified in Figure S3 the number of snapshots, during the 100 μs of MD simulation, in which Phe2837.37 adopts the trans conformation or/and the distance between the top of TMs 1 (Cα atom of Thr341.33) and 7 (Cα atom of Val2777.31) increases from the initial value of 14.5 Å to values larger than 15.5 Å. GPCRs are dynamic proteins that permit rapid small-scale structural fluctuations,68 and accordingly, both events can be freely observed, simultaneously or not, during the calculated trajectories (Figure S3). Importantly, the trajectory in which the ligand (orange in Figure 1b) spontaneously binds the lipid-facing part Figure 2. Well-tempered metadynamics simulations of the pathway of of TMs 1 and 7 (surface-bound, Figure S3) triggers or ligand entry to the orthosteric binding site of CB2R. (a) The ligand stabilizes the trans conformation of Phe2837.37 (Figure 1d), pathway from the initial position of JWH-133 (in orange, see Figure and the distance between TMs 1 and 7 increases to 16.5 Å 1b) to the final reference position bound to the orthosteric site of (Figure 1f) most of the simulation time. In contrast, in other CB2R (in black), obtained during the well-tempered metadynamics, is trajectories in which the ligand did not bind the receptor represented by dots (the dot color gradient, from red to blue, corresponds to simulation time, from the beginning to the end, surface and these events did not spontaneously occur (surface- respectively). Five independent replicas are shown. (b) The free- unbound, Figure S3), Phe2837.37 remains in the gauche+ energy profile of the five replicas (the last 10 snapshots before conformation (Figure 1d) and the distance between TMs 1 convergence are plotted). The collective variable (CV) is the distance and 7 decreases to values lower than the initial 14.5 Å (Figure between the center of mass of JWH-133 in the initial and final 1f). Figure 1e shows two selected snapshots of these surface- positions. Three energetic minima are observed when the ligand is unbound and surface-bound trajectories (see Figure S3) in bound to the lipid-facing part of the tunnel (green rectangle), to the which the tunnel between TMs 1 and 7 is closed and open, bundle-facing part of the tunnel (blue), or to the orthosteric site respectively. (pink). These three positions in the ligand pathway are depicted in 3.4. Ligand Access to the Orthosteric Binding Site. panel a by an ellipse with a color matching to the colors in panel b. The initial state of the binding process is on the right part of the graph The above results, and the experimental evidence (see and the final state on the left part. (c−e) Representative structures of Introduction), led us to use the ligand bound to TMs 1 and JWH-133 and the interacting side chains, along the ligand pathway, at 7 of CB2R to progressively sample the binding event. Because the lipid-facing part (c), the bundle-facing part (d), and the entry to the orthosteric binding site can have the largest orthosteric site (e). The color of JWH-133 corresponds to the energetic barrier of the process, 15,40 in some cases color of the minima, whereas the final reference position is shown in unsurmountable by unbiased MD simulations, we used well- black. tempered metadynamics.69 In this technique, a biasing potential is applied to permit the system to explore CB2R tunnel mainly formed by amino acids (Val361.35, energetically inaccessible regions for unbiased MD simulations Cys401.39, Phe872.57, Phe912.61, His952.65, Ala2827.36, and on reasonable time scales. We used as a collective variable of Met2867.40) in TMs 1, 2, and 7 (Figure 2d). In the final the distance between the center of mass of the ligand in the third minimum of the binding process, JWH-133 reaches the initial conformation and the reference docking pose at the orthosteric binding site. This computationally derived binding orthosteric binding site (see Materials and Methods). In the pose, using metadynamics, reproduced the main contacts free-energy profile of five independent simulations (Figure 2a between JWH-133 and the receptor (Phe872.57, Phe912.61, and b), three minima can be noticed in most trajectories. Phe942.64, His952.65, Thr1143.33, Phe1173.36, Phe183ECL2, These steps in the entry pathway were designated as lipid- Trp1945.43, Phe2817.35, and Ser2857.39) previously proposed63 facing, bundle-facing, and orthosteric. In the first minimum, (Figure 2e). the ligand is bound to the lipid-facing part of the CB2R tunnel 3.5. Experimental Validation of the Pathway of mainly formed by amino acids (Val361.35, Cys401.39, Leu431.42, Ligand Entry. To experimentally validate the proposed Ala2827.36, Phe2837.37, Met2867.40, Leu2897.43, Ile2907.44, and mechanism of ligand entry by site-directed mutagenesis, we Met2937.47) in TMs 1 and 7 (Figure 2c). In the second explored the amino acids involved in the transition from the minimum, the ligand is bound to the bundle-facing part of the lipid- to the bundle-facing part of the CB2R tunnel (Figure 5774 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article 2c,d). Val361.35 and Ala2827.36 are located midway between Table 1. Functional Properties of JWH-133 at Nonmutated these two minima, delimiting the tunnel between TMs 1 and 7. and Mutated CB2R Notably, three aromatic, Phe912.61, Phe942.64, and His952.65, pEC50a Emaxb side chains in TM 2 are in their environment (Figure 3a). nonmutated 7.5 ± 0.1 47.0 ± 3.5 Val361.35Phe 6.8 ± 0.2 60.2 ± 6.7 Ala2827.36Phe 6.5 ± 0.6 84.3 ± 8.4 Val361.35 Met 7.9 ± 0.1 56.7 ± 1.7 Ala2827.36 Met 7.8 ± 0.2 45.4 ± 4.2 a pEC50 (nM). bEmax (%), the maximum inhibition of forskolin- stimulated cAMP levels (normalized to 100%). These values were calculated using nonlinear regression analysis. Data are expressed as the mean ± SE of at least nine independent experiments performed in triplicate. the effect of JWH-133, as the decrease of forskolin-induced cAMP is of less magnitude in the mutant receptors (the effect is more noticeable in Ala2827.36Phe than in Val361.35Phe). The importance of these stable aromatic interactions was tested by comparing the results with those obtained using Val361.35 Met and Ala2827.36 Met mutants (Figure 3b and Table 1). The highly polarizable sulfur atom of the Met side chain can form with other aromatic side chains stronger interactions than Phe.70,71 Interestingly, these new mutations do not affect the decrease in forskolin-induced cAMP in a significant manner, thus suggesting they do not block the entrance of the JWH-133 ligand probably due to the higher flexibility of Met relative to Phe. Our experimental results point out that the highly stable aromatic cluster between Phe912.61, Phe942.64, His952.65, and Phe2827.36 (in the Ala2827.36Phe mutation) blocks the entrance of JWH-133 to the orthosteric binding site. 3.6. Ligand Access to the Binding Site in the Val361.35Phe and Ala2827.36Phe Mutant Receptors. To study the influence of this aromatic cluster, attained in the Val361.35Phe and Ala2827.36Phe mutants, during the process of Figure 3. A single point mutation blocks the entrance of the ligand. ligand entry, we performed well-tempered metadynamics using (a) Extracellular view of CB2R (ECLs are omitted for clarity) in which the same protocol as for nonmutated CB2R. As expected, the the position of Val361.35 (green surface) and Ala2827.36 (red), located between the lipid-facing (Figure 2c) and bundle-facing (Figure 2d) Val361.35Phe mutation hampers the ability of the ligand to minima, and Phe912.61, Phe942.64, and His952.65 side chains in TM 2 access the orthosteric site in one of the replicates and the (blue) are shown. (b) Decrease of forskolin-induced cAMP possibility to explore the minimum corresponding to the (normalized to 100%), in HEK-293T cells, upon stimulation of bundle-facing part of the CB2R tunnel in three other replicates nonmutated (black line), Val361.35Phe (green) and Ala2827.36Phe (Figure 3c). On the other hand, the Ala2827.36Phe mutation (red) mutations (mutation to Phe is shown as a solid line), and blocks the entrance of the ligand in four of the five replicates Val361.35 Met (green) and Ala2827.36 Met (red) mutations (mutation performed (Figure 3d) in agreement with our hypothesis and to Met is shown as a dotted line) with the JWH-133 agonist. (c,d) experimental results. Ligand pathways and free-energy profiles (the last 10 snapshots before convergence are plotted) of five replicas, obtained in well-tempered 4. DISCUSSION AND CONCLUSIONS metadynamics of JWH-133 entry to the Val361.35Phe (c) or Ala2827.36Phe (d) mutant receptors. The collective variable (CV) is Drug-target residence time is a key element as a biological the distance between the center of mass of JWH-133 in the initial effect requires binding to the target receptor.72,73 Residence (orange) and final (black) conformations. See legend of Figure 2 for time is related to the ligand−receptor binding kinetics and, additional details. thus, to the pathway used in ligand association and dissociation.74 MD simulations have become a convenient tool to estimate residence times75,76 and to analyze such Thus, we obtained the Val361.35Phe and Ala2827.36Phe mutant pathways.15,16,21−25,40,77 These simulations show that ligand versions of the CB2R with the main idea that the aromatic side binding is a multistage process in which the ligand forms key chains would form a stable aromatic cluster with Phe912.61. interactions with amino acids at the entrance of the binding Phe942.64, and/or His952.65, blocking the access of the ligand to site. Most of these simulations are for class A GPCRs that bind the orthosteric binding site. Dose−response experiments in polar ligands directly from the extracellular environment. HEK-293T cells expressing the CB2R stimulated with forskolin Here, we have studied the binding pathway of the potent and treated with JWH-133 showed reduced cAMP production JWH-133 agonist to CB2R by MD simulations. In CB2R, the (Figure 3b and Table 1), as expected for Gi-coupled receptors. N-terminus and ECL-2 that lacks the disulfide bridge to TM 3, Clearly, both Val361.35Phe and Ala2827.36Phe mutations impair characteristic of the GPCR family, fold over the ligand binding 5775 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article pocket39 blocking access to the orthosteric binding cavity from common in the GPCR field to perform mutagenesis assays to the extracellular environment. We, thus, propose that the study the influence in affinity or efficacy of amino acids located hydrophobic JWH-133 (calculated logP of 6.4) diffuses at the orthosteric binding pocket. Here, we aimed at decreasing through the bilayer leaflet to contact a lipid-facing cavity of the residence time of the ligand by blocking its entrance and, CB2R mainly formed by amino acids in TMs 1 and 7 (Figure thus, blocking the signaling response (decrease of forskolin- 4). In this point of the pathway, the ligand forms hydrophobic/ induced cAMP). The Ala2827.36Phe mutation, but not Ala2827.36 Met, allows a blockade of the effect of JWH-133 because the decrease of forskolin-induced cAMP is of less magnitude in the mutant than in nonmutated receptor (Figure 3b). We propose that a highly stable aromatic cluster between Phe912.61, Phe942.64, His952.65, and Phe2827.36 (in the Ala2827.36Phe mutation), which cannot be accomplished with Met2827.36, blocks the entrance of JWH-133 to the orthosteric binding site. Figure S4 shows a structure-based sequence alignment of the amino acids forming this pathway among the members of class A GPCRs that bind hormone-like signaling molecules derived from lipid species. Finally, JWH-133 moves to the most energetically stable position in the pathway, at the orthosteric binding site (Figure 4). The pathway of JWH-133 entry to CB2R defines two transient binding sites, in addition to the final location at the orthosteric cavity (Figure 4). These transient sites have been proposed to be potential binding sites for synthetic modulators (allosteric modulators or bitopic ligands), becoming important sites for drug discovery due to the high conservation of the traditional orthosteric binding site.18,19,79 To analyze this suggestion, we have superimposed our previously reported binding modes of a negative allosteric modulator63 and a bitopic ligand78 of CB2R. Clearly, the bitopic ligand binds at the orthosteric site and the identified transient lipid-facing site. Figure 4. The binding pathway of JWH-133 to the orthosteric site of Thus, the lipid-facing site in CB2R is analogous to the vestibule CB2R. JWH-133 (in orange) diffuses through the membrane to reach or secondary or metastable binding site defined for other a minimum in the free-energy profile in which JWH-133 (in green) GPCRs that bind polar ligands (see Introduction), which is key contacts a membrane vestibule or lipid-facing cavity. Subsequently, for ligand desolvation15,16,40 and for selectivity,20,28,77 but it JWH-133 (in blue) moves to a second allosteric minimum or bundle- might not allosterically modulate receptor activity. In contrast, facing cavity in the trajectory. Finally, JWH-133 (in pink) moves to the negative allosteric modulator cannabidiol binds at the the most stable minimum in the pathway at the orthosteric binding identified bundle-facing site, from which it elicits the allosteric site. Ligand positions during the pathway are shown by dots. RMSD modulation.63 Thus, the bundle-facing site is an allosteric site values (light gray from unbiased MD simulations, dark gray from well- in the CB2R, located near the orthosteric site, as has been tempered metadynamics) relative to the reference docked binding mode of JWH-133 at the orthosteric site and a representative free- suggested for muscarinic receptors.79 In conclusion, our energy profile obtained in well-tempered metadynamics are shown. findings have shown that ligand pathway simulations might Color bars indicate the different locations of the ligand. Comparison be a useful tool to identify and delineate cavities for the design of our previously reported binding mode of the negative allosteric of synthetic modulators of GPCRs. modulator cannabidiol (CBD, gray)63 and JWH-133 at the allosteric minimum or bundle-facing cavity (in blue), and the binding mode of 5. DATA AND SOFTWARE AVAILABILITY a bitopic ligand (compound 22, gray)78 and JWH-133 at the vestibule or lipid-facing cavity (in green) and the orthosteric binding site (in PACKMOL-Memgen, distributed with AmberTools, is free of pink). charge. MDAnalysis, GROMACS, and PLUMED are open source. FATSLiM and CAVER are licensed under the GNU General Public License. VMD is available to noncommercial aromatic contacts with Leu39 1.38, Leu421.41, Leu431.44, users under a distribution-specific license, and PyMOL is Phe2837.37, Met2867.40, and Ile2907.44 (Figure 2c). Importantly, commercial software with a paid license. the binding of the ligand to this lipid-facing cavity triggers and/ or stabilizes the trans conformation of Phe2837.37 (Figure 1d) and increases the distance between TMs 1 and 7 (Figure 1f), opening the tunnel between TMs 1 and 7 (Figure 1e). ■ * ASSOCIATED CONTENT sı Supporting Information Subsequently, JWH-133 moves to a second energetic The Supporting Information is available free of charge at minimum in the pathway, located in the bundle-facing part https://pubs.acs.org/doi/10.1021/acs.jcim.2c00865. of the CB2R tunnel (Figure 4), to mainly interact with Val361.35, Cys401.39, Phe872.57, Phe912.61, His952.65, Ala2827.36, Workflow of the computational process, heatmap and Met2867.40 (Figure 2d). We experimentally validated this representation of the cavity opening events, sequence entrance tunnel by mutating Ala2827.36 that is not located at alignment of class A GPCRs that bind hydrophobic-like the orthosteric binding site but midway between the lipid- and signaling molecules, and analysis of membrane stability the bundle-facing part of the CB2R tunnel (Figure 3b). It is (PDF) 5776 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article ■ AUTHOR INFORMATION Corresponding Author coupled receptor-effector macromolecular membrane assemblies (GEMMAs). Pharmacol. Ther. 2022, 231, 107977. (2) Rasmussen, S. G.; DeVree, B. T.; Zou, Y.; Kruse, A. C.; Chung, Leonardo Pardo − Laboratory of Computational Medicine, K. Y.; Kobilka, T. S.; Thian, F. S.; Chae, P. S.; Pardon, E.; Calinski, D.; Biostatistics Unit, Faculty of Medicine, Universitat Autònoma Mathiesen, J. M.; Shah, S. T.; Lyons, J. A.; Caffrey, M.; Gellman, S. Barcelona, 08193 Barcelona, Spain; orcid.org/0000- H.; Steyaert, J.; Skiniotis, G.; Weis, W. I.; Sunahara, R. K.; Kobilka, B. 0003-1778-7420; Email: leonardo.pardo@uab.cat K. Crystal structure of the beta2 adrenergic receptor-Gs protein complex. Nature 2011, 477 (7366), 549−55. Authors (3) Kang, Y.; Zhou, X. E.; Gao, X.; He, Y.; Liu, W.; Ishchenko, A.; Nil Casajuana-Martin − Laboratory of Computational Barty, A.; White, T. A.; Yefanov, O.; Han, G. W.; Xu, Q.; de Waal, P. Medicine, Biostatistics Unit, Faculty of Medicine, Universitat W.; Ke, J.; Tan, M. H.; Zhang, C.; Moeller, A.; West, G. M.; Pascal, B. Autònoma Barcelona, 08193 Barcelona, Spain D.; Van Eps, N.; Caro, L. N.; Vishnivetskiy, S. A.; Lee, R. J.; Suino- Gemma Navarro − Department of Biochemistry and Powell, K. M.; Gu, X.; Pal, K.; Ma, J.; Zhi, X.; Boutet, S.; Williams, G. Physiology, Faculty of Pharmacy and Food Sciences, J.; Messerschmidt, M.; Gati, C.; Zatsepin, N. A.; Wang, D.; James, D.; Universitat de Barcelona, 08028 Barcelona, Spain; Centro de Basu, S.; Roy-Chowdhury, S.; Conrad, C. E.; Coe, J.; Liu, H.; Lisova, S.; Kupitz, C.; Grotjohann, I.; Fromme, R.; Jiang, Y.; Tan, M.; Yang, Investigación en Red, Enfermedades Neurodegenerativas H.; Li, J.; Wang, M.; Zheng, Z.; Li, D.; Howe, N.; Zhao, Y.; Standfuss, (CIBERNED), Instituto de Salud Carlos III, 28031 Madrid, J.; Diederichs, K.; Dong, Y.; Potter, C. S.; Carragher, B.; Caffrey, M.; Spain Jiang, H.; Chapman, H. N.; Spence, J. C.; Fromme, P.; Weierstall, U.; Angel Gonzalez − Laboratory of Computational Medicine, Ernst, O. P.; Katritch, V.; Gurevich, V. V.; Griffin, P. R.; Hubbell, W. Biostatistics Unit, Faculty of Medicine, Universitat Autònoma L.; Stevens, R. C.; Cherezov, V.; Melcher, K.; Xu, H. E. Crystal Barcelona, 08193 Barcelona, Spain structure of rhodopsin bound to arrestin by femtosecond X-ray laser. Claudia Llinas del Torrent − Laboratory of Computational Nature 2015, 523 (7562), 561−7. Medicine, Biostatistics Unit, Faculty of Medicine, Universitat (4) Chen, Q.; Plasencia, M.; Li, Z.; Mukherjee, S.; Patra, D.; Chen, Autònoma Barcelona, 08193 Barcelona, Spain C. L.; Klose, T.; Yao, X. Q.; Kossiakoff, A. A.; Chang, L.; Andrews, P. Marc Gómez-Autet − Laboratory of Computational Medicine, C.; Tesmer, J. J. G. Structures of rhodopsin in complex with G- Biostatistics Unit, Faculty of Medicine, Universitat Autònoma protein-coupled receptor kinase 1. Nature 2021, 595 (7868), 600− 605. Barcelona, 08193 Barcelona, Spain (5) Bushdid, C.; Magnasco, M. O.; Vosshall, L. B.; Keller, A. Aleix Quintana García − Laboratory of Computational Humans can discriminate more than 1 trillion olfactory stimuli. Medicine, Biostatistics Unit, Faculty of Medicine, Universitat Science 2014, 343 (6177), 1370−2. Autònoma Barcelona, 08193 Barcelona, Spain; orcid.org/ (6) Dagan-Wiener, A.; Di Pizio, A.; Nissim, I.; Bahia, M. S.; 0000-0003-4472-5129 Dubovski, N.; Margulis, E.; Niv, M. Y. BitterDB: taste ligands and Rafael Franco − Department of Biochemistry and Molecular receptors database in 2019. Nucleic Acids Res. 2019, 47 (D1), D1179− Biomedicine, Faculty of Biology, Universitat de Barcelona, D1185. 08028 Barcelona, Spain; Centro de Investigación en Red, (7) Kooistra, A. J.; Mordalski, S.; Pandy-Szekeres, G.; Esguerra, M.; Enfermedades Neurodegenerativas (CIBERNED), Instituto Mamyrbekov, A.; Munk, C.; Keseru, G. M.; Gloriam, D. E. GPCRdb de Salud Carlos III, 28031 Madrid, Spain; orcid.org/ in 2021: integrating GPCR sequence, structure and function. Nucleic 0000-0003-2549-4919 Acids Res. 2021, 49 (D1), D335−D343. (8) Hauser, A. S.; Attwood, M. M.; Rask-Andersen, M.; Schioth, H. Complete contact information is available at: B.; Gloriam, D. E. Trends in GPCR drug discovery: new agents, https://pubs.acs.org/10.1021/acs.jcim.2c00865 targets and indications. Nat. Rev. Drug Discov. 2017, 16 (12), 829− 842. Author Contributions (9) Deupi, X.; Dolker, N.; Luz Lopez-Rodriguez, M.; Campillo, M.; N.C.-M., A.G., and L.P. devised the project concept and Ballesteros, J.; Pardo, L. Structural Models of Class A G protein- designed experiments. N.C.-M., C.L.T., and M.G.-A. con- coupled receptors as a Tool for Drug Design: Insights on tributed with computational simulations. G.N. and R.F. Transmembrane Bundle Plasticity. Curr. Top. Med. Chem. 2007, 7 (10), 999−1006. performed biochemical assays. A.Q.G. contributed to the (10) Venkatakrishnan, A. J.; Deupi, X.; Lebon, G.; Heydenreich, F. data analysis. N.C.-M., R.F., and L.P. wrote the paper with M.; Flock, T.; Miljus, T.; Balaji, S.; Bouvier, M.; Veprintsev, D. B.; contributions from all other authors. All authors contributed to Tate, C. G.; Schertler, G. F.; Babu, M. M. Diverse activation pathways the data analysis and have given approval to the final version of in class A GPCRs converge near the G-protein-coupling region. the manuscript. Nature 2016, 536 (7617), 484−7. Notes (11) Santos, R.; Ursu, O.; Gaulton, A.; Bento, A. P.; Donadi, R. S.; The authors declare no competing financial interest. Bologa, C. G.; Karlsson, A.; Al-Lazikani, B.; Hersey, A.; Oprea, T. I.; Overington, J. P. A comprehensive map of molecular drug targets. ■ ACKNOWLEDGMENTS We acknowledge the financial support from the Spanish Nat. Rev. Drug Discov. 2017, 16 (1), 19−34. (12) Zhou, Q.; Yang, D.; Wu, M.; Guo, Y.; Guo, W.; Zhong, L.; Cai, X.; Dai, A.; Jang, W.; Shakhnovich, E. I.; Liu, Z. J.; Stevens, R. C.; Ministry of Economy and Innovation with FEDER funds Lambert, N. A.; Babu, M. M.; Wang, M. W.; Zhao, S. Common (RTI2018-094204-B-I00, PID2019-109240RB-I00, PID2020- activation mechanism of class A GPCRs. eLife 2019, 8, e50279. 113430RB-I00). C.L.T. is recipient of a FPI fellowship (BES- (13) Schoneberg, T.; Liebscher, I. Mutations in G Protein-Coupled 2017-081872), and M.G.-A. and A.Q.G. acknowledge the Receptors: Mechanisms, Pathophysiology and Potential Therapeutic Approaches. Pharmacol. Rev. 2021, 73 (1), 89−119. Universitat Autònoma de Barcelona for their predoctoral grant. (14) Venkatakrishnan, A. J.; Deupi, X.; Lebon, G.; Tate, C. G.; ■ REFERENCES (1) Ferre, S.; Ciruela, F.; Dessauer, C. W.; Gonzalez-Maeso, J.; Schertler, G. F.; Babu, M. M. Molecular signatures of G-protein- coupled receptors. Nature 2013, 494 (7436), 185−94. (15) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Hebert, T. E.; Jockers, R.; Logothetis, D. E.; Pardo, L. G protein- Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Pathway and mechanism 5777 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. Pharmacology of Class A Lipid GPCRs: Opportunities and U.S.A. 2011, 108 (32), 13118−23. Challenges for Drug Discovery. Pharmaceuticals (Basel) 2022, 15 (16) Gonzalez, A.; Perez-Acle, T.; Pardo, L.; Deupi, X. Molecular (1), 12. Basis of Ligand Dissociation in beta-Adrenergic Receptors. PloS one (32) Gonzalez, A.; Cordomi, A.; Caltabiano, G.; Pardo, L. Impact of 2011, 6 (9), No. e23815. helix irregularities on sequence alignment and homology modelling of (17) Wang, C.; Jiang, Y.; Ma, J.; Wu, H.; Wacker, D.; Katritch, V.; G protein-coupled receptors. Chembiochem 2012, 13 (10), 1393− Han, G. W.; Liu, W.; Huang, X. P.; Vardy, E.; McCorvy, J. D.; Gao, 1399. X.; Zhou, E. X.; Melcher, K.; Zhang, C.; Bai, F.; Yang, H.; Yang, L.; (33) Hanson, M. A.; Roth, C. B.; Jo, E.; Griffith, M. T.; Scott, F. L.; Jiang, H.; Roth, B. L.; Cherezov, V.; Stevens, R. C.; Xu, H. E. Reinhart, G.; Desale, H.; Clemons, B.; Cahalan, S. M.; Schuerer, S. C.; Structural Basis for Molecular Recognition at Serotonin Receptors. Sanna, M. G.; Han, G. W.; Kuhn, P.; Rosen, H.; Stevens, R. C. Crystal Science 2013, 340, 610−614. structure of a lipid G protein-coupled receptor. Science 2012, 335 (18) Newman, A. H.; Battiti, F. O.; Bonifazi, A. 2016 Philip S. (6070), 851−5. Portoghese Medicinal Chemistry Lectureship: Designing Bivalent or (34) Chrencik, J. E.; Roth, C. B.; Terakado, M.; Kurata, H.; Omi, R.; Bitopic Molecules for G-Protein Coupled Receptors. The Whole Is Kihara, Y.; Warshaviak, D.; Nakade, S.; Asmar-Rovira, G.; Mileni, M.; Greater Than the Sum of Its Parts. J. Med. Chem. 2020, 63 (5), 1779− Mizuno, H.; Griffith, M. T.; Rodgers, C.; Han, G. W.; Velasquez, J.; 1797. Chun, J.; Stevens, R. C.; Hanson, M. A. Crystal Structure of (19) Fronik, P.; Gaiser, B. I.; Sejer Pedersen, D. Bitopic Ligands and Antagonist Bound Human Lysophosphatidic Acid Receptor 1. Cell Metastable Binding Sites: Opportunities for G Protein-Coupled 2015, 161 (7), 1633−43. Receptor (GPCR) Medicinal Chemistry. J. Med. Chem. 2017, 60 (10), (35) Srivastava, A.; Yano, J.; Hirozane, Y.; Kefala, G.; Gruswitz, F.; 4126−4134. Snell, G.; Lane, W.; Ivetac, A.; Aertgeerts, K.; Nguyen, J.; Jennings, A.; (20) Masureel, M.; Zou, Y.; Picard, L. P.; van der Westhuizen, E.; Okada, K. High-resolution structure of the human GPR40 receptor Mahoney, J. P.; Rodrigues, J.; Mildorf, T. J.; Dror, R. O.; Shaw, D. E.; bound to allosteric agonist TAK-875. Nature 2014, 513 (7516), 124− Bouvier, M.; Pardon, E.; Steyaert, J.; Sunahara, R. K.; Weis, W. I.; 7. Zhang, C.; Kobilka, B. K. Structural insights into binding specificity, (36) Wang, L.; Yao, D.; Deepak, R.; Liu, H.; Xiao, Q.; Fan, H.; efficacy and bias of a beta2AR partial agonist. Nat. Chem. Biol. 2018, Gong, W.; Wei, Z.; Zhang, C. Structures of the Human PGD2 14 (11), 1059−1066. Receptor CRTH2 Reveal Novel Mechanisms for Ligand Recognition. (21) Kruse, A. C.; Hu, J.; Pan, A. C.; Arlow, D. H.; Rosenbaum, D. Mol. Cell 2018, 72 (1), 48−59. M.; Rosemond, E.; Green, H. F.; Liu, T.; Chae, P. S.; Dror, R. O.; (37) Shao, Z.; Yin, J.; Chapman, K.; Grzemska, M.; Clark, L.; Wang, Shaw, D. E.; Weis, W. I.; Wess, J.; Kobilka, B. K. Structure and J.; Rosenbaum, D. M. High-resolution crystal structure of the human dynamics of the M3 muscarinic acetylcholine receptor. Nature 2012, CB1 cannabinoid receptor. Nature 2016, 540 (7634), 602−606. 482 (7386), 552−6. (38) Hua, T.; Vemuri, K.; Pu, M.; Qu, L.; Han, G. W.; Wu, Y.; Zhao, (22) Sabbadin, D.; Moro, S. Supervised molecular dynamics S.; Shui, W.; Li, S.; Korde, A.; Laprairie, R. B.; Stahl, E. L.; Ho, J. H.; (SuMD) as a helpful tool to depict GPCR-ligand recognition pathway Zvonok, N.; Zhou, H.; Kufareva, I.; Wu, B.; Zhao, Q.; Hanson, M. A.; in a nanosecond time scale. J. Chem. Inf. Model. 2014, 54 (2), 372−6. Bohn, L. M.; Makriyannis, A.; Stevens, R. C.; Liu, Z. J. Crystal (23) Wittmann, H. J.; Strasser, A. Binding pathway of histamine to Structure of the Human Cannabinoid Receptor CB1. Cell 2016, 167 the hH4R, observed by unconstrained molecular dynamics. Bioorg. (3), 750−762. Med. Chem. Lett. 2015, 25 (6), 1259−68. (39) Li, X.; Hua, T.; Vemuri, K.; Ho, J. H.; Wu, Y.; Wu, L.; Popov, (24) Thomas, T.; Fang, Y.; Yuriev, E.; Chalmers, D. K. Ligand P.; Benchama, O.; Zvonok, N.; Locke, K.; Qu, L.; Han, G. W.; Iyer, Binding Pathways of Clozapine and Haloperidol in the Dopamine D2 M. R.; Cinar, R.; Coffey, N. J.; Wang, J.; Wu, M.; Katritch, V.; Zhao, and D3 Receptors. J. Chem. Inf. Model. 2016, 56 (2), 308−21. S.; Kunos, G.; Bohn, L. M.; Makriyannis, A.; Stevens, R. C.; Liu, Z. J. (25) Xu, Z.; Guo, L.; Qian, X.; Yu, C.; Li, S.; Zhu, C.; Ma, X.; Li, H.; Crystal Structure of the Human Cannabinoid Receptor CB2. Cell Zhu, G.; Zhou, H.; Dai, W.; Li, Q.; Gao, X. Two entry tunnels in 2019, 176 (3), 459−467. mouse TAAR9 suggest the possibility of multi-entry tunnels in (40) Stanley, N.; Pardo, L.; Fabritiis, G. D. The pathway of ligand olfactory receptors. Sci. Rep. 2022, 12 (1), 2691. entry from the membrane bilayer to a lipid G protein-coupled (26) Hedderich, J. B.; Persechino, M.; Becker, K.; Heydenreich, F. receptor. Sci. Rep. 2016, 6, 22639. M.; Gutermuth, T.; Bouvier, M.; Bunemann, M.; Kolb, P. The (41) Krishna Kumar, K.; Shalev-Benami, M.; Robertson, M. J.; Hu, pocketome of G-protein-coupled receptors reveals previously H.; Banister, S. D.; Hollingsworth, S. A.; Latorraca, N. R.; Kato, H. E.; untargeted allosteric sites. Nat. Commun. 2022, 13 (1), 2567. Hilger, D.; Maeda, S.; Weis, W. I.; Farrens, D. L.; Dror, R. O.; (27) Kruse, A. C.; Ring, A. M.; Manglik, A.; Hu, J.; Hu, K.; Eitel, K.; Malhotra, S. V.; Kobilka, B. K.; Skiniotis, G. Structure of a Signaling Hubner, H.; Pardon, E.; Valant, C.; Sexton, P. M.; Christopoulos, A.; Cannabinoid Receptor 1-G Protein Complex. Cell 2019, 176 (3), Felder, C. C.; Gmeiner, P.; Steyaert, J.; Weis, W. I.; Garcia, K. C.; 448−458. Wess, J.; Kobilka, B. K. Activation and allosteric modulation of a (42) Stauch, B.; Johansson, L. C.; McCorvy, J. D.; Patel, N.; Han, G. muscarinic acetylcholine receptor. Nature 2013, 504 (7478), 101−6. W.; Huang, X. P.; Gati, C.; Batyuk, A.; Slocum, S. T.; Ishchenko, A.; (28) Medina, R. A.; Vazquez-Villa, H.; Gomez-Tamayo, J. C.; Brehm, W.; White, T. A.; Michaelian, N.; Madsen, C.; Zhu, L.; Grant, Benhamu, B.; Martin-Fontecha, M.; de la Fuente, T.; Caltabiano, G.; T. D.; Grandner, J. M.; Shiriaeva, A.; Olsen, R. H. J.; Tribo, A. R.; Hedlund, P. B.; Pardo, L.; Lopez-Rodriguez, M. L. The extracellular Yous, S.; Stevens, R. C.; Weierstall, U.; Katritch, V.; Roth, B. L.; Liu, entrance provides selectivity to serotonin 5-HT7 receptor antagonists W.; Cherezov, V. Structural basis of ligand recognition at the human with antidepressant-like behavior in vivo. J. Med. Chem. 2014, 57 (15), MT1 melatonin receptor. Nature 2019, 569 (7755), 284−288. 6879−84. (43) Hildebrand, P. W.; Scheerer, P.; Park, J. H.; Choe, H. W.; (29) Perez-Benito, L.; Henry, A.; Matsoukas, M. T.; Lopez, L.; Piechnick, R.; Ernst, O. P.; Hofmann, K. P.; Heck, M. A ligand Pulido, D.; Royo, M.; Cordomi, A.; Tresadern, G.; Pardo, L. The size channel through the G protein coupled receptor opsin. PLoS One matters? A computational tool to design bivalent ligands. 2009, 4 (2), No. e4382. Bioinformatics 2018, 34 (22), 3857−3863. (44) Tian, H.; Gunnison, K. M.; Kazmi, M. A.; Sakmar, T. P.; (30) Pulido, D.; Casado-Anguera, V.; Gomez-Autet, M.; Llopart, N.; Huber, T. FRET sensors reveal the retinal entry pathway in the G Moreno, E.; Casajuana-Martin, N.; Ferre, S.; Pardo, L.; Casado, V.; protein-coupled receptor rhodopsin. iScience 2022, 25 (4), 104060. Royo, M. Heterobivalent Ligand for the Adenosine A2A-Dopamine (45) Hurst, D. P.; Grossfield, A.; Lynch, D. L.; Feller, S.; Romo, T. D2 Receptor Heteromer. J. Med. Chem. 2022, 65 (1), 616−632. D.; Gawrisch, K.; Pitman, M. C.; Reggio, P. H. A lipid pathway for (31) Krishna Deepak, R. N. V.; Verma, R. K.; Hartono, Y. D.; Yew, ligand binding is necessary for a cannabinoid G protein-coupled W. S.; Fan, H. Recent Advances in Structure, Function, and receptor. J. Biol. Chem. 2010, 285 (23), 17954−64. 5778 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article (46) Bokoch, M. P.; Jo, H.; Valcourt, J. R.; Srinivasan, Y.; Pan, A. C.; (64) Huffman, J. W.; Liddle, J.; Yu, S.; Aung, M. M.; Abood, M. E.; Capponi, S.; Grabe, M.; Dror, R. O.; Shaw, D. E.; DeGrado, W. F.; Wiley, J. L.; Martin, B. R. 3-(1’,1’-Dimethylbutyl)-1-deoxy-delta8- Coughlin, S. R. Entry from the Lipid Bilayer: A Possible Pathway for THC and related compounds: synthesis of selective ligands for the Inhibition of a Peptide G Protein-Coupled Receptor by a Lipophilic CB2 receptor. Bioorg. Med. Chem. Lett. 1999, 7 (12), 2905−14. Small Molecule. Biochemistry 2018, 57 (39), 5748−5758. (65) Ballesteros, J. A.; Weinstein, H. Integrated methods for the (47) Hua, T.; Li, X.; Wu, L.; Iliopoulos-Tsoutsouvas, C.; Wang, Y.; construction of three dimensional models and computational probing Wu, M.; Shen, L.; Brust, C. A.; Nikas, S. P.; Song, F.; Song, X.; Yuan, of structure-function relations in G-protein coupled receptors. S.; Sun, Q.; Wu, Y.; Jiang, S.; Grim, T. W.; Benchama, O.; Stahl, E. L.; Methods in Neurosciences 1995, 25, 366−428. Zvonok, N.; Zhao, S.; Bohn, L. M.; Makriyannis, A.; Liu, Z. J. (66) Xing, C.; Zhuang, Y.; Xu, T. H.; Feng, Z.; Zhou, X. E.; Chen, Activation and Signaling Mechanism Revealed by Cannabinoid M.; Wang, L.; Meng, X.; Xue, Y.; Wang, J.; Liu, H.; McGuire, T. F.; Receptor-Gi Complex Structures. Cell 2020, 180 (4), 655−665. Zhao, G.; Melcher, K.; Zhang, C.; Xu, H. E.; Xie, X. Q. Cryo-EM (48) Sondergaard, C. R.; Olsson, M. H.; Rostkowski, M.; Jensen, J. Structure of the Human Cannabinoid Receptor CB2-Gi Signaling H. Improved Treatment of Ligands and Coupling Effects in Empirical Complex. Cell 2020, 180 (4), 645−654. Calculation and Rationalization of pKa Values. J. Chem. Theory (67) Liu, S.; Paknejad, N.; Zhu, L.; Kihara, Y.; Ray, M.; Chun, J.; Comput. 2011, 7 (7), 2284−95. Liu, W.; Hite, R. K.; Huang, X. Y. Differential activation mechanisms (49) Lomize, M. A.; Pogozheva, I. D.; Joo, H.; Mosberg, H. I.; of lipid GPCRs by lysophosphatidic acid and sphingosine 1- Lomize, A. L. OPM database and PPM web server: resources for phosphate. Nat. Commun. 2022, 13 (1), 731. positioning of proteins in membranes. Nucleic Acids Res. 2012, 40, (68) Nygaard, R.; Zou, Y.; Dror, R. O.; Mildorf, T. J.; Arlow, D. H.; Manglik, A.; Pan, A. C.; Liu, C. W.; Fung, J. J.; Bokoch, M. P.; Thian, D370−D376. F. S.; Kobilka, T. S.; Shaw, D. E.; Mueller, L.; Prosser, R. S.; Kobilka, (50) Schott-Verdugo, S.; Gohlke, H. PACKMOL-Memgen: A B. K. The Dynamic Process of beta(2)-Adrenergic Receptor Simple-To-Use, Generalized Workflow for Membrane-Protein-Lipid- Activation. Cell 2013, 152 (3), 532−42. Bilayer System Building. J. Chem. Inf. Model. 2019, 59 (6), 2522− (69) Wu, X.; Xu, L. Y.; Li, E. M.; Dong, G. Application of molecular 2528. dynamics simulation in biomedicine. Chem. Biol. Drug Des. 2022, 99 (51) Abraham, M. J.; Murtola, T.; Schulz, R.; Pall, S.; Smith, J. C.; (5), 789−800. Hess, B.; Lindahl, E. GROMACS: High performance molecular (70) Gomez-Tamayo, J. C.; Cordomi, A.; Olivella, M.; Mayol, E.; simulations through multi-level parallelism from laptops to super- Fourmy, D.; Pardo, L. Analysis of the interactions of sulfur-containing computers. SoftwareX 2015, 1−2, 19−25. amino acids in membrane proteins. Protein Sci. 2016, 25 (8), 1517− (52) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; 24. Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of (71) Gibbs, C. A.; Weber, D. S.; Warren, J. J. Clustering of Aromatic Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Amino Acid Residues around Methionine in Proteins. Biomolecules Theory Comput. 2015, 11 (8), 3696−713. 2022, 12 (1), 6. (53) Dickson, C. J.; Madej, B. D.; Skjevik, A. A.; Betz, R. M.; Teigen, (72) Guo, D.; Hillger, J. M.; IJzerman, A. P.; Heitman, L. H. Drug- K.; Gould, I. R.; Walker, R. C. Lipid14: The Amber Lipid Force Field. target residence time−a case for G protein-coupled receptors. Med. J. Chem. Theory Comput. 2014, 10 (2), 865−879. Res. Rev. 2014, 34 (4), 856−92. (54) Sandoval-Perez, A.; Pluhackova, K.; Bockmann, R. A. Critical (73) Sykes, D. A.; Stoddart, L. A.; Kilpatrick, L. E.; Hill, S. J. Binding Comparison of Biomembrane Force Fields: Protein-Lipid Interactions kinetics of ligands acting at GPCRs. Mol. Cell. Endocrinol. 2019, 485, at the Membrane Interface. J. Chem. Theory Comput. 2017, 13 (5), 9−19. 2310−2321. (74) Strasser, A.; Wittmann, H. J.; Seifert, R. Binding Kinetics and (55) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. Pathways of Ligands to GPCRs. Trends Pharmacol. Sci. 2017, 38 (8), A. Development and testing of a general amber force field. J. Comput. 717−732. Chem. 2004, 25 (9), 1157−74. (75) Mollica, L.; Theret, I.; Antoine, M.; Perron-Sierra, F.; Charton, (56) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, Y.; Fourquez, J. M.; Wierzbicki, M.; Boutin, J. A.; Ferry, G.; O. MDAnalysis: a toolkit for the analysis of molecular dynamics Decherchi, S.; Bottegoni, G.; Ducrot, P.; Cavalli, A. Molecular simulations. J. Comput. Chem. 2011, 32 (10), 2319−27. Dynamics Simulations and Kinetic Measurements to Estimate and (57) Buchoux, S. FATSLiM: a fast and robust software to analyze Predict Protein-Ligand Residence Times. J. Med. Chem. 2016, 59 MD simulations of membranes. Bioinformatics 2017, 33 (1), 133− (15), 7167−76. 134. (76) Potterton, A.; Husseini, F. S.; Southey, M. W. Y.; Bodkin, M. J.; (58) The PyMOL Molecular Graphics System, Version 1.3r1; Heifetz, A.; Coveney, P. V.; Townsend-Nicholson, A. Ensemble-Based Schrodinger, LLC, 2010. Steered Molecular Dynamics Predicts Relative Residence Time of (59) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular A2A Receptor Binders. J. Chem. Theory Comput. 2019, 15 (5), 3316− dynamics. J. Mol. Graph. 1996, 14 (1), 33−38. 3330. (60) Chovancova, E.; Pavelka, A.; Benes, P.; Strnad, O.; Brezovsky, (77) Xu, X.; Kaindl, J.; Clark, M. J.; Hubner, H.; Hirata, K.; J.; Kozlikova, B.; Gora, A.; Sustr, V.; Klvana, M.; Medek, P.; Sunahara, R. K.; Gmeiner, P.; Kobilka, B. K.; Liu, X. Binding pathway Biedermannova, L.; Sochor, J.; Damborsky, J. CAVER 3.0: a tool for determines norepinephrine selectivity for the human beta1AR over the analysis of transport pathways in dynamic protein structures. PLoS beta2AR. Cell Res. 2021, 31 (5), 569−579. Comput. Biol. 2012, 8 (10), No. e1002708. (78) Morales, P.; Navarro, G.; Gomez-Autet, M.; Redondo, L.; (61) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Fernandez-Ruiz, J.; Perez-Benito, L.; Cordomi, A.; Pardo, L.; Franco, Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. R.; Jagerovic, N. Discovery of Homobivalent Bitopic Ligands of the Commun. 2014, 185, 604−613. Cannabinoid CB2 Receptor. Chemistry 2020, 26 (68), 15839−15842. (62) consortium, P. Promoting transparency and reproducibility in (79) Kruse, A. C.; Kobilka, B. K.; Gautam, D.; Sexton, P. M.; enhanced molecular simulations. Nat. Methods 2019, 16 (8), 670− Christopoulos, A.; Wess, J. Muscarinic acetylcholine receptors: novel 673. opportunities for drug development. Nat. Rev. Drug Discov. 2014, 13 (63) Navarro, G.; Gonzalez, A.; Sanchez-Morales, A.; Casajuana- (7), 549−60. Martin, N.; Gomez-Ventura, M.; Cordomi, A.; Busque, F.; Alibes, R.; Pardo, L.; Franco, R. Design of Negative and Positive Allosteric Modulators of the Cannabinoid CB2 Receptor Derived from the Natural Product Cannabidiol. J. Med. Chem. 2021, 64 (13), 9354− 9364. 5779 https://doi.org/10.1021/acs.jcim.2c00865 J. Chem. Inf. Model. 2022, 62, 5771−5779