Advertisement

Mechanistic origin of partial agonism of tetrahydrocannabinol for cannabinoid receptors

  • Soumajit Dutta
    Affiliations
    Department of Chemical and Biomolecular Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
    Search for articles by this author
  • Balaji Selvam
    Affiliations
    Department of Chemical and Biomolecular Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
    Search for articles by this author
  • Aditi Das
    Affiliations
    Department of Comparative Biosciences, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Department of Bioengineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Cancer Center at Illinois, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
    Search for articles by this author
  • Diwakar Shukla
    Correspondence
    For correspondence: Diwakar Shukla
    Affiliations
    Department of Chemical and Biomolecular Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    Cancer Center at Illinois, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    National Center for Supercomputing Applications, University of Illinois, Urbana, Illinois, USA

    Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

    NIH Center for Macromolecular Modeling and Bioinformatics, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
    Search for articles by this author
Open AccessPublished:February 25, 2022DOI:https://doi.org/10.1016/j.jbc.2022.101764
      Cannabinoid receptor 1 (CB1) is a therapeutically relevant drug target for controlling pain, obesity, and other central nervous system disorders. However, full agonists and antagonists of CB1 have been reported to cause serious side effects in patients. Therefore, partial agonists have emerged as a viable alternative as they can mitigate overstimulation and side effects. One of the key bottlenecks in the design of partial agonists, however, is the lack of understanding of the molecular mechanism of partial agonism itself. In this study, we examine two mechanistic hypotheses for the origin of partial agonism in cannabinoid receptors and predict the mechanistic basis of partial agonism exhibited by Δ9-Tetrahydrocannabinol (THC) against CB1. In particular, we inspect whether partial agonism emerges from the ability of THC to bind in both agonist and antagonist-binding poses or from its ability to only partially activate the receptor. We used extensive molecular dynamics simulations and Markov state modeling to capture the THC binding in both antagonist and agonist-binding poses in the CB1 receptor. Furthermore, we predict that binding of THC in the agonist-binding pose leads to rotation of toggle switch residues and causes partial outward movement of intracellular transmembrane helix 6 (TM6). Our simulations also suggest that the alkyl side chain of THC plays a crucial role in determining partial agonism by stabilizing the ligand in the agonist and antagonist-like poses within the pocket. Taken together, this study provides important insights into the mechanistic origin of the partial agonism of THC.

      Keywords

      Abbreviations:

      CB1 (cannabinoid receptor 1), GPCR (G-protein-coupled receptor), K-L (Kullback–Leibler), kMC (kinetic Monte Carlo), MD (molecular dynamic), MSM (Markov state model), NAM (negative allosteric modulator), PCA (principal component analysis), TI (thermodynamics integration), THC (tetrahydrocannabinol), TM6 (transmembrane helix 6)
      Cannabinoid receptor 1(CB1) belongs to the family of class A G-protein-coupled receptors (GPCRs) (
      • Hilger D.
      • Masureel M.
      • Kobilka B.K.
      Structure and dynamics of GPCR signaling complexes.
      ,
      • Rosenbaum D.M.
      • Rasmussen S.G.
      • Kobilka B.K.
      The structure and function of G-protein-coupled receptors.
      ), which modulates diverse cellular signaling processes via intracellular G-proteins (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ) and β-arrestins (
      • Ibsen M.S.
      • Finlay D.B.
      • Patel M.
      • Javitch J.A.
      • Glass M.
      • Grimsey N.L.
      Cannabinoid CB1 and CB2 receptor-mediated arrestin translocation: Species, subtype, and agonist-dependence.
      ). CB1 receptors were first discovered in the last decade of the 20th century as a target of plant cannabinoid molecules (
      • Zou S.
      • Kumar U.
      Cannabinoid receptors and the endocannabinoid system: Signaling and function in the central nervous system.
      ). Due to its ubiquitous presence in physiological processes, CB1 is an important drug target for the potential treatment of a variety of diseases. In the last 30 years, several synthetic molecules have been designed to target the CB1 for the treatment of pain, obesity, and inflammation (
      • Pertwee R.G.
      Pharmacology of cannabinoid receptor ligands.
      ,
      • Pertwee R.G.
      Cannabinoid pharmacology: The first 66 years.
      ,
      • Pertwee R.G.
      • Ross R.
      Cannabinoid receptors and their ligands.
      ). CB1 receptor agonists such as MDMB-fubinaca (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ) and inverse agonists such as rimonabant (
      • MacLennan S.
      • Reynen P.
      • Kwan J.
      • Bonhaus D.
      Evidence for inverse agonism of SR141716A at human recombinant cannabinoid CB1 and CB2 receptors.
      ,
      • Pertwee R.G.
      Inverse agonism and neutral antagonism at cannabinoid CB1 receptors.
      ) have been shown to modulate the receptor activity significantly. However, these designed agonists and antagonists of CB1 also exhibit dangerous side effects. For instance, fubinaca, also known as ”zombie drug,” caused the hospitalization of thousands of patients in New York (
      • Adams A.J.
      • Banister S.D.
      • Irizarry L.
      • Trecki J.
      • Schwartz M.
      • Gerona R.
      “Zombie” outbreak caused by the synthetic cannabinoid AMB-FUBINACA in New York.
      ). Rimonabant had to be withdrawn from the market due to its psychological side effects such as depression and anxiety (
      • Sam A.H.
      • Salem V.
      • Ghatei M.A.
      Rimonabant: From RIO to ban.
      ). While these designed agonists and antagonists have failed to meet the drug safety guidelines, alternate approaches (e.g., allosteric modulator, partial agonist) can be explored for designing therapeutics. We recently studied the binding of a negative allosteric modulator (NAM), sodium ion (Na+), to cannabinoid receptors (CBs) using molecular dynamics simulation (
      • Dutta S.
      • Selvam B.
      • Shukla D.
      Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors.
      ). Simulation revealed important differences in binding site and pathway between CB1 and CB2, which can be exploited to design a selective NAM drug. Similarly, a partial agonist dronabinol has been used as an appetite stimulant drug for AIDS patients and an antiemetic drug for chemotherapy (
      • Badowski M.E.
      • Yanful P.K.
      Dronabinol oral solution in the management of anorexia and weight loss in AIDS and cancer.
      ). Dronabinol is a synthetic form of tetrahydrocannabinol (THC) (Fig. 1), the main psychoactive compound in marijuana, which binds to the CB1 as a partial agonist and affects the endocannabinoid signaling pathway. THC has been shown to demonstrate positive effects for treating Huntington’s disease, Parkinson’s disease, and Alzheimer’s disease (
      • Walther S.
      • Halpern M.
      Cannabinoids and dementia: A review of clinical and preclinical data.
      ). Although THC has the potential to become a valuable drug for several diseases, this drug is still banned by the FDA due to its side effects. Therefore, molecular-level understanding of THC and the mechanism by which it partially activates the cannabinoid receptors will inform the design of potential partial agonist drugs targeting CB1 receptor.
      Figure thumbnail gr1
      Figure 12D and 3D representation of THC structure. A, 2D representation. B, 3D representation. Numbering of carbon atoms and rings are mentioned in the figure. THC, tetrahydrocannabinol.
      Structural studies of CB1 have revealed that toggle switch residue (TRP3566.48 and PHE2003.36) movement by the agonist is crucial for the activation of CB1 (Fig. 2B) (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ,
      • Hua T.
      • Li X.
      • Wu L.
      • Iliopoulos-Tsoutsouvas C.
      • Wang Y.
      • Wu M.
      • Shen L.
      • Brust C.A.
      • Nikas S.P.
      • Song F.
      • Song X.
      • Yuan S.
      • Sun Q.
      • Wu Y.
      • Jiang S.
      • et al.
      Activation and signaling mechanism revealed by cannabinoid receptor-Gi complex structures.
      ,
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ,
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ,
      • Shao Z.
      • Yin J.
      • Chapman K.
      • Grzemska M.
      • Clark L.
      • Wang J.
      • Rosenbaum D.M.
      High-resolution crystal structure of the human CB1 cannabinoid receptor.
      ,
      • Shao Z.
      • Yan W.
      • Chapman K.
      • Ramesh K.
      • Ferrell A.J.
      • Yin J.
      • Wang X.
      • Xu Q.
      • Rosenbaum D.M.
      Structure of an allosteric modulator bound to the CB1 cannabinoid receptor.
      ,
      • Li X.
      • Shen L.
      • Hua T.
      • Liu Z.-J.
      Structural and functional insights into cannabinoid receptors.
      ). Kumar et al. (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ) proposed that due to the smaller size of THC as compared with the full agonists, there is less interaction between toggle switch residues and ligand. Furthermore, they proposed that the THC-binding position with downward facing alkyl chain can shift the toggle switch residues and activate the receptor (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ). However, the crystal structures and docking studies only provide static interactions of the ligand in the binding pose. These studies do not reveal the mechanism of ligand conformational switching inside the binding pocket. To date, this proposed hypothesis has not been rigorously examined from either the experimental or computational approaches. Therefore, it is difficult to obtain structural understanding of the origin of partial agonism exhibited by THC. A recent study employing metadynamics simulations predicted that a partial agonist GAT228 binds in multiple positions inside the ligand-binding pocket of CB1 (
      • Saleh N.
      • Hucke O.
      • Kramer G.
      • Schmidt E.
      • Montel F.
      • Lipinski R.
      • Ferger B.
      • Clark T.
      • Hildebrand P.W.
      • Tautermann C.S.
      Multiple binding sites contribute to the mechanism of mixed agonistic and positive allosteric modulators of the cannabinoid CB1 receptor.
      ) due to the large size of the pocket as compared with the ligand volume (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). This prediction is also consistent with the distinct poses of agonist and antagonist in the binding pocket (Fig. 2A). Therefore, it is likely that partial agonist THC might also be stabilized in both the agonist and antagonist-like pose inside the binding pocket (Equations 2 and 3), and only the subset of THC molecules bound to CB1 in the agonist-like pose activates the receptor (Equation 2). This phenomenon would decrease the maximum response of the secondary messenger. The first hypothesis is shown as Equations (1), (2), (3) where R represents the receptor in apo form, RA∗∗ is the agonist bound active state; RPA∗∗ and [RPA] represent partial agonist bound active and inactive state, respectively. In RPA∗∗ and [RPA], partial agonist binds in agonist and antagonist-bound poses, respectively. Agonist and partial agonist are represented as A and PA.
      Figure thumbnail gr2
      Figure 2Pictorial representation of both the hypotheses. A and B are representing Hypothesis 1. Superposition of active (PDB ID: 5XRA (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ), color: pink) and inactive (PDB ID: 5TGZ (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ), color: cyan) structures of CB1 is shown in all panels of (A). In the top and bottom panel of (A), antagonist (AM6538, color: yellow) and agonist (AM11542, color: silver) are shown as stick and mesh representation. N-loop is colored differently (Active: red, Inactive: Blue). TM6 and TM7 are not shown in the cartoon representation for better visualization of ligand poses. In (B) toggle switch residues (TRP3566.48 and PHE2003.36) are shown as sticks. Superposition of active (PDB ID: 3SN6 (
      • Rasmussen S.G.F.
      • DeVree B.T.
      • Zou Y.
      • Kruse A.C.
      • Chung K.Y.
      • Kobilka T.S.
      • Thian F.S.
      • Chae P.S.
      • Pardon E.
      • Calinski D.
      • Mathiesen J.M.
      • Shah S.T.A.
      • Lyons J.A.
      • Caffrey M.
      • Gellman S.H.
      • et al.
      Crystal structure of the β2 adrenergic receptor–Gs protein complex.
      ), color: pink), partially active (PDB ID: 6CSY (
      • Masureel M.
      • Zou Y.
      • Picard L.P.
      • Westhuizen E.V.D.
      • Mahoney J.P.
      • Rodrigues J.P.G.L.M.
      • Mildorf T.J.
      • Dror R.O.
      • Shaw D.E.
      • Bouvier M.
      • Pardon E.
      • Steyaert J.
      • Sunahara R.K.
      • Weis W.I.
      • Zhang C.
      • et al.
      Structural insights into binding specificity, efficacy and bias of a β2AR partial agonist.
      ), color: green), and inactive (PDB ID: 2RH1 (
      • Cherezov V.
      • Rosenbaum D.M.
      • Hanson M.A.
      • Rasmussen S.G.F.
      • Thian F.S.
      • Kobilka T.S.
      • Choi H.-J.
      • Kuhn P.
      • Weis W.I.
      • Kobilka B.K.
      • Stevens R.C.
      High-resolution crystal structure of an engineered human β 2 -adrenergic G protein–coupled receptor.
      ), color: cyan) structures of β2-AR is shown in (C). Intracellular TM6 movement is highlighted in a separated box. CB1, cannabinoid receptor 1; TM6, transmembrane helix 6.

      Hypothesis I

      R+ARA
      (1)


      R+PAAgonistPoseRPA
      (2)


      R+PAAntagonistPose[RPA]
      (3)


      Canonical class A GPCR activation is characterized by intracellular TM6 movement, which facilitates the G-protein binding. A partial agonist, salmeterol, binds in the orthosteric pocket in the β2-adrenergic receptor and causes partial movement of TM6 compared with the full movement by an agonist (Fig. 2C) (
      • Rasmussen S.G.F.
      • DeVree B.T.
      • Zou Y.
      • Kruse A.C.
      • Chung K.Y.
      • Kobilka T.S.
      • Thian F.S.
      • Chae P.S.
      • Pardon E.
      • Calinski D.
      • Mathiesen J.M.
      • Shah S.T.A.
      • Lyons J.A.
      • Caffrey M.
      • Gellman S.H.
      • et al.
      Crystal structure of the β2 adrenergic receptor–Gs protein complex.
      ,
      • Masureel M.
      • Zou Y.
      • Picard L.P.
      • Westhuizen E.V.D.
      • Mahoney J.P.
      • Rodrigues J.P.G.L.M.
      • Mildorf T.J.
      • Dror R.O.
      • Shaw D.E.
      • Bouvier M.
      • Pardon E.
      • Steyaert J.
      • Sunahara R.K.
      • Weis W.I.
      • Zhang C.
      • et al.
      Structural insights into binding specificity, efficacy and bias of a β2AR partial agonist.
      ,
      • Cherezov V.
      • Rosenbaum D.M.
      • Hanson M.A.
      • Rasmussen S.G.F.
      • Thian F.S.
      • Kobilka T.S.
      • Choi H.-J.
      • Kuhn P.
      • Weis W.I.
      • Kobilka B.K.
      • Stevens R.C.
      High-resolution crystal structure of an engineered human β 2 -adrenergic G protein–coupled receptor.
      ). Therefore, we hypothesized that partial activation may happen if the partial agonist stabilizes the receptor in different conformation than the active structure. THC has smaller side chain compared with agonist AM11542 (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). Furthermore, absence of the dimethyl group at the first carbon (C1') of the alkyl chain decreases the interaction with toggle switch residues. Thus, we proposed that THC binding may also cause the partial outward movement of the TM6.

      Hypothesis II

      R+ARA
      (4)


      R+PARPA
      (5)


      The second proposed hypothesis is explained in Equations 4 and 5 where RPA∗ represents partial agonist-bound partially active form of the receptor. Other notations are similar to the reactions 1, 2, 3.
      We assessed the validity of these hypotheses by running extensive simulations of THC binding to CB1 (see Experimental procedures section). Using the Markovian property of molecular dynamic (MD) simulations, we built Markov state model (MSM) using the simulation data. From the eigenvalues and vectors of MSM, we predicted the timescale and the thermodynamics of the ligand-binding process. MSM weighted data predict that THC is stabilized in both antagonist and agonist-bound poses in agreement with our first hypothesis. The free energy barrier of ∼1 kcal/mol is estimated for THC to transition between the antagonist and agonist-bound pose, which implies that thermal fluctuations at the human body temperature could easily allow the THC to transition between these poses. In the agonist-bound pose, THC rotates the important toggle switch residue TRP3566.48. The new position of the TRP3566.48 is different compared with active state of CB1. This intermediate position TRP3566.48 leads to partial outward movement of TM6 suggesting that THC can only partially active the receptor according to the second hypothesis. Furthermore, comparison of the ensemble where THC bound in agonist pose with the agonist (AM11542)-bound simulation showed that major differences are in the movement and signaling of the activation microswitches, supporting our results that the partial agonist cannot lead to the full activation of the receptor. This mechanistic study predicts the reason behind the partial agonism behavior of THC compared with other agonists and will aid future drug development targeting cannabinoid receptors.

      Results

      THC is predicted to be stabilized in both antagonist and agonist-binding poses in the orthosteric pocket of CB1

      Active and inactive structures of CB1 reveal that orthosteric-binding site volume undergoes a large change upon activation as compared with other class A GPCRs (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). Comparison of active (PDB ID: 5XRA (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      )) and inactive (PDB ID: 5TGZ (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      )) structures also reveals that the agonist and antagonist bind in the different regions within the pocket (Fig. 2A). Agonist molecule binds in a region close to the TM5, whereas antagonist binds in the extended pocket formed by TM1 and TM2. In the inactive structure, the downward movement of the N-loop toward the binding pocket separates the agonist and antagonist-binding regions (Fig. 2A). However, partial agonist-bound crystal structure of CB1 or CB2 is not reported in the literature. Therefore, the partial agonist-binding position is not well documented. Preliminary docking studies in both inactive and active structures reveal that THC binds to the agonist-binding region in a similar conformation as other agonists (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ,
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). However, docking alone cannot infer the exact binding pose of a ligand, which undergoes a dynamic conformational change within the binding pocket. Therefore, we performed ∼143 μs of MD simulations to predict the THC-binding mechanism in CB1 starting with inactive structure.
      To capture the THC-binding process and upward movement of the N-loop, we projected the MD simulation data along the two metrics that characterize these motions (Figs. 3A and S1, A and B). MSM weighted free energy landscape plot (Fig. 3A) shows the predicted movement of THC molecule toward TM5 as indicated by the distance between THC-C1' atom (Fig. 1) and TYR2755.39-Cα atom on TM5. THC is predicted to enter CB1-binding pocket through the space between the space of TM1, TM2, and N-loop (Fig. S2) (
      • Hua T.
      • Li X.
      • Wu L.
      • Iliopoulos-Tsoutsouvas C.
      • Wang Y.
      • Wu M.
      • Shen L.
      • Brust C.A.
      • Nikas S.P.
      • Song F.
      • Song X.
      • Yuan S.
      • Sun Q.
      • Wu Y.
      • Jiang S.
      • et al.
      Activation and signaling mechanism revealed by cannabinoid receptor-Gi complex structures.
      ). From MD simulations, THC is further predicted to diffuse inside the pocket and stabilized in antagonist-binding pose where distance between THC(C1') and TYR2755.39(Cα) is 15 to 21 Å (Fig. 3B). We observed two stable local energy wells in the antagonist-binding pose. The free energy well away from agonist-bound pose is named antagonist-like pose 1, and second well is named antagonist-like pose 2. These two minima are separated by the activation barrier of 0.55 ± 0.43 kcal/mol. Superposition of predicted MD structures of THC bound in antagonist-like pose 1 and 2 with inactive structure shows that B and C ring of the tricyclic dibenzopyran group of THC binds in same position as arm 3 of the antagonist, AM6538 (Fig. 3, C and D) (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ). The aromatic ring of THC matches with pyrazole ring of the antagonist. However, THC alkyl chain (side chain) orientation varies between the two binding poses. In antagonist-like pose 1, it extends toward the conserved sodium-binding site (
      • Dutta S.
      • Selvam B.
      • Shukla D.
      Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors.
      • Katritch V.
      • Fenalti G.
      • Abola E.E.
      • Roth B.L.
      • Cherezov V.
      • Stevens R.C.
      Allosteric sodium in class A GPCR signaling.
      ,
      • Selvam B.
      • Shamsi Z.
      • Shukla D.
      Universality of the sodium ion binding mechanism in class AG-protein-coupled receptors.
      ) similar to the arm 1 of antagonist, whereas, in antagonist-like pose 2, it orients itself in the direction of agonist-binding site similar to the arm 2 of the antagonist (Fig. 3, C and D). In both the antagonist-like poses, THC forms stable polar interaction with SER3837.39 and hydrophobic interaction with N-loop (PHE102N−loop, MET103N−loop), TM1 (SER1231.39, ILE1191.35), TM2 (PHE1702.57, PHE1742.61, PHE1772.64), and TM7 (ALA3807.36) (Figs. 4, B and C, S2B and S3A). Previous experiments have shown the importance of PHE1742.61 and PHE1772.64 in the THC binding (
      • Shim J.-Y.
      • Bertalovitz A.C.
      • Kendall D.A.
      Identification of essential cannabinoid-binding domains.
      ). Relative free energy calculation using alchemical analysis shows that mutations of PHE1742.61 and PHE1772.64 decrease the stability (ΔΔGPHE1742.61ALA=1.62±0.28 kcal/mol and ΔΔGPHE1772.64ALA=3.45±0.56 kcal/mol) of the antagonist pose 1 suggesting that these mutations might decrease the THC binding by destabilizing the antagonist-bound pose (
      • Song L.F.
      • Lee T.-S.
      • Zhu C.
      • York D.M.
      • Merz K.M.
      Using AMBER18 for relative free energy calculations.
      ). In these poses, N-loop remains inside the pocket and restricts the movement of THC toward the agonist-binding region. The upward movement of N-loop allows the movement of THC inside the pocket and stabilizes it in the agonist-like pose (Figs. 3E and S16). THC has smaller alkyl side chain than cannabinoid-like full agonists of CB1; thus, first carbon of THC (C1') binds closer to TM5 compared with the agonist, AM11542. In agonist-binding pose, THC(OH) forms polar interaction with SER3837.39 same as other agonists. Furthermore, THC forms extensive hydrophobic interactions with amino acid residues of TM3 (VAL1963.32, LEU1933.29), TM5 (TRP2795.43), TM7 (PHE3797.35), N-loop (MET103N−loop), and ECL2 (PHE268ECL2) as shown in Figure 4D and Fig. S3C.
      Figure thumbnail gr3
      Figure 3Distinct stabilized poses of THC inside binding pocket of CB1. A, MSM weighted free energy landscape to capture THC binding and N-loop upward motion. THC binding distance is measured between THC-C1' and TYR2755.39-Cα (TM5) and N-loop upward motion is measured between MET103N−loop-Cα (N-loop) and ASP1632.50-Cα (TM2). B, one-dimensional free energy diagram depicting stabilized binding position of THC and activation barrier between them. C and D, superposition of inactive (PDB ID: 5TGZ (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ), color: cyan) structure of CB1 and MD snapshots from antagonist-like pose 1 (C) and 2 (D). E, superposition of active (PDB ID: 5XRA (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ), color: pink) structure of CB1 and MD snapshot from agonist-like pose. MD snapshots are shown in green color. Agonist (AM11542 (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      )), partial agonist (THC) and antagonist (AM6538 (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      )) are represented as sticks with silver, violet, and yellow color respectively. TM6 and TM7 are not shown in the cartoon representation for better visualization of ligand poses. MD, molecular dynamic; MSM, Markov state model; THC, tetrahydrocannabinol.
      Figure thumbnail gr4
      Figure 4Important interactions between protein residues and THC at different stabilized positions during binding (side view). A, a representative structure when THC enters the receptor through the space between N-loop, TM1 and TM2. B and C, representative structures from antagonist-like pose 1 and 2, respectively. D, a representative structure from agonist-like pose. Stable interactions were measured using GetContacts package. Protein structures are shown as cartoon representation (color: green). THC (color: violet) and interactive residues (color: green) are shown as stick. THC, tetrahydrocannabinol.
      To capture the timescale of this entire binding process, we ran kinetic Monte Carlo (kMC) simulations on MD data. kMC utilizes MSM transition probability matrix to find the probable pathway for binding (method section). In total, 150 μs long KMC trajectory reveals that entire binding process from unbound to bound poses takes approximately 100 μs (Fig. S4A). From the solution, THC is first stabilized in antagonist-bound pose in ∼50 μs. THC occupies antagonist pose I and II for approximately 30 μs and subsequently moves to the agonist-binding pose (Fig. S4A).

      THC chain orientation plays an important role in partial agonism

      Alkyl chain of THC plays an important role in binding and activation of CB1. Modification of Alkyl chain leads to change in binding affinity. For example, increasing the chain length of Δ8-THC (structural homolog of Δ9-THC) from five carbons to eight carbons increases the binding affinity from ∼40 nM to ∼8 nM (
      • Bow E.W.
      • Rimoldi J.M.
      The structure–function relationships of classical cannabinoids: CB1/CB2 modulation.
      ). Furthermore, adding a dimethyl group in the first carbon of the alkyl chain is hypothesized to increase the interaction with toggle switch TRP3566.48 (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). To characterize the importance of the alkyl chain of THC, we observed chain dihedral angle (C2-C3-C1'-C2') (Fig. 1) movement during binding. Positive dihedral is crucial to orient the alkyl chain of THC in orthogonal direction of aromatic group as in agonist-like pose (Fig. 3E). The free energy landscape of THC alkyl chain dihedral with respect to binding (Figs. 5A and S5) reveals that THC binds to CB1 in two different chain orientations. If hydrogen atom on C1 faces downward (or positive dihedral angle) during the binding, it moves to agonist-binding pocket with maximum free energy barrier of ∼2 ± 0.4 kcal/mol. This alkyl chain orientation enables THC to pressurize N-loop in upward direction to take orthogonal conformation similar to full agonist. However, if THC enters the receptor with negative dihedral, it is stabilized in the antagonist-like poses. High free energy barrier of ∼4 ± 0.4 kcal/mol is required for THC to move from macrostate 3′ to 4′ (Fig. 5A). In this conformation (macrostate 3′), tricyclic dibenzopyran group of THC is unable shift N-loop upward as the active structure. Hence, alkyl chain of THC shifts the population to macrostate 3, which is more accessible as the free energy barrier is lower than earlier transition. Applying transition path theory (TPT) on MSM states, we calculated effective timescales of these macrostate transitions. TPT provides mean free passage time (MFPT) for the transition between two macrostates by taking into account all possible pathways through the intermediate states. Calculated MFPT shows that transition between state 3′ and 3 (2.4 ± 0.5 μs) is more accessible compared with the 3′ to 4′ (13.2 ± 3.5 μs) (Figs. 5B and S6). The binding of THC to the agonist-like pose shows the orientation of alkyl chain may flip between macrostate 4 and 4′ as they are energetically favorable. However, THC has more tendency to stay in the orthogonal conformation (macrostate 4) compared with macrostate 4′. The conditional probability of THC to be in orthogonal conformation in the agonist-like pose is 0.9 ± 0.0. Our results show that THC is stabilized in agonist and antagonist-like poses that supports our first hypothesis.
      Figure thumbnail gr5
      Figure 5THC alkyl chain movement during binding to CB1. A, MSM weighted free energy landscape to capture THC binding and THC side chain dihedral. THC binding distance is measured between THC-C1' and TYR275-Cα (TM5) and THC sidechain dihedral is measured between C2,C3,C1',C2'. Manually defined macrostate regions are numbered in the figure. B, mean free passage time (MFPT) of transitions between the eight macrostates is shown. Each macrostate is represented by MD snapshot from the region. Different range of MFPT is shown with distinguished arrow thickness. Protein structures are shown as cartoon and THC molecules are shown as stick. TM6 and TM7 are not shown in the cartoon representation for better visualization of ligand poses. THC, tetrahydrocannabinol; TM6, transmembrane helix 6.

      THC is predicted to rotate the toggle switch TRP3566.48 in the agonist-binding pocket

      Although, we established the prediction that THC can be stabilized in different positions of the orthosteric-binding pocket, it is not clear how THC activates the receptor. Crystal structures of the CB1 receptor in active and inactive position depict that toggle switch residues play an important role in the activation of the receptor. Agonist molecule triggers the movement of toggle switch residue TRP3566.48 movement of the receptor toward the TM5, which consequently leads to outward movement of TM6 (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). It was hypothesized that THC behaves as partial agonist due to the lack of interaction with toggle switch residues (
      • Kumar K.K.
      • Shalev-Benami M.
      • Robertson M.J.
      • Hu H.
      • Banister S.D.
      • Hollingsworth S.A
      • Latorraca N.R.
      • Kato H.E.
      • Hilger D.
      • Maeda S.
      • Weis W.I.
      • Farrens D.L.
      • Dror R.O.
      • Malhotra S.V.
      • Kobilka B.K.
      • et al.
      Structure of a signaling cannabinoid receptor 1-G protein complex.
      ). However, we observed that binding of THC leads to the rotation of TRP3566.48 (Figs. 6A and S7). The dihedral angle of TRP3566.48 shifts from inactive conformation (χ2 dihedral angle between 60° and 120°) to new intermediate active conformation (−30° to 30°) as well as relatively less stable state (−120° to −60°) (Fig. 6, A and B). We referred energetically favorably accessible states as partially active state 1 and state 2, respectively. The rotation of TRP3566.48 leads to breakage of aromatic interactions with PHE2003.36, and it moves toward TM2 similar to active structure (Fig. 6B). Comparison of representative structures from partially active state 1 with inactive CB2 shows that toggle switch TRP3566.48/TRP2586.48 has similar rotamaric conformation as CB2 inactive structure (Fig. S8). An inverse agonist of CB2, MRI2687, which was predicted to retain the toggle switch residue to similar conformation (
      • Li X.
      • Hua T.
      • Vemuri K.
      • Ho J.H.
      • Wu Y.
      • Wu L.
      • Popov P.
      • Benchama O.
      • Zvonok N.
      • Locke K.
      • Qu L.
      • Han G.W.
      • Iyer M.R.
      • Cinar R.
      • Coffey N.J.
      • et al.
      Crystal structure of the human cannabinoid receptor CB2.
      ), acts as a partial agonist for CB1. This in turn supports our prediction that THC stabilizes the TRP6.48 in intermediate state to partially activate the receptor.
      Figure thumbnail gr6
      Figure 6TRP3566.48 rotation due to THC binding in Agonist-like pose. A, MSM weighted free energy landscape to capture THC binding and toggle switch TRP3566.48 χ2 angle. THC binding distance is measured between THC-C1' and TYR2755.39-Cα (TM5). Inactive and partially active (state 1 and 2) macrostates are marked when the THC is bound to the agonist bound pocket. B, mean free passage time (MFPT) of transitions between the inactive to partially active conformations of toggle switch residues is shown. Direction of the conformational change in TRP3566.48 and PHE2003.36 are shown via arrow. Different range of MFPT is shown with distinguished arrow thickness. Protein structures (color: green) are shown as cartoon. THC molecules (color: violet) and toggle switch residues (color: green) are shown as stick. C, box plot to show the probability of TRP3566.48 rotation with positive or negative THC dihedral in agonist-like pose. Blue and red boxes show the conditional probabilities when TRP3566.48 is in inactive and partially active pose. Data distribution in box plot is generated with 200 rounds of bootstrap sampling with 80% of total number of trajectories ( section). MSM, Markov state model; THC, tetrahydrocannabinol.
      Although, THC is predicted to rotate toggle switch TRP3566.48 for subset of structures in agonist-binding position, we noticed another stable state around inactive pose of TRP3566.48 (Fig. 6A). In this state, THC is likely unable to shift the toggle switch movement. To explain the reason for the two different orientation of TRP3566.48, we calculated the probability TRP3566.48 rotation with respect to the THC alkyl chain dihedral in agonist-like pose (Fig. 6C). Our calculations show that with negative alkyl chain dihedral, THC can move TRP3566.48 for only 4.0 ± 0.8% structures. With negative dihedral, the first carbon of the alkyl chain (C1') binds far from the toggle switch and cannot induce conformational change (Fig. 5B). However, binding with a positive dihedral leads to rotation of TRP3566.48 for 95.5 ± 0.8% of the structures due to more interactions.
      The changes in the toggle switch residue lead to the outward movement of TM6 as shown in the Figure 7A and Fig. S9. When the toggle switch remains in the inactive position, intracellular TM6 can move easily 2 to 3 Å either side inactive structure (low free energy barrier). For this case, TM3-TM6 distance is similar or less than the inactive structure with a probability of 0.74 ± 0.0 (Fig. 7B). However, TRP3566.48 rotation creates a torsion in TM6 and stabilizes the intracellular part of TM6 of the receptor in between inactive and active-like conformation with probability of 0.6 ± 0.0 (Fig. 7B). Comparison of representative structure from partially active state minima with β2-AR inactive structure shows that intracellular TM5, TM6, and TM7 match well (Fig. S10). Therefore, THC in partially activated state forms favorable interaction with TRP3566.48 and stabilizes in intermediate conformation compared with the inactive structure. These results support our second hypothesis that THC can only activate the receptor partially.
      Figure thumbnail gr7
      Figure 7Movement of intracellular TM6 due to TRP3566.48 rotation. A, MSM weighted free energy landscape to capture toggle switch TRP3566.48 (TM6) χ2 angle and intracellular TM6 movement of CB1. Intracellular TM6 movement is measured ARG2143.50-Cα (TM3) and LYS3436.35-Cα (TM6). B, box plot to show the probabilities of intracellular TM6 movement with TRP3566.48 in inactive and partially active condition. Data distribution in box plot is generated with 200 rounds of bootstrap sampling with 80% of total number of trajectories ( section). C and D, representative structure from partial active minima (color: green) is superimposed with inactive (color: cyan) (C) and active (color: pink) (D) structure of CB1 to depict the toggle switch movement and intracellular movement. MSM, Markov state model; TM6, transmembrane helix 6.

      Comparison of partial agonist and agonist-bound simulation ensemble using data-driven approach

      We have performed 17 μs MD simulations of agonist (AM11542) bound CB1 (Experimental procedures section) and compared the result with the ensemble with the partial agonist bound to the agonist-like pose (
      • Fleetwood O.
      • Carlsson J.
      • Delemotte L.
      Identification of ligand-specific G protein-coupled receptor states and prediction of downstream efficacy via data-driven modeling.
      ). Principal component analysis (PCA) was implemented on all possible interresidue distances to compare the ligand-bound ensembles of the partial and full agonists. PCA projection onto the first two components shows that the receptor explores different conformational spaces for different ligands (Fig. 8A). As PC1 and PC2 emphasize large amplitude motions, we calculated distance pairs with highest weight contributions for PC1 and PC2. For both agonist and partial agonist-bound ensembles, these distances correspond to the distances between truncated N-terminus and transmembrane domain (Figs. S11 and S12). Binding of agonist and partial agonist leads to outward movement of N-terminus from the orthosteric binding pocket compared with the inactive state (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      ). As the N-terminus moves out of the binding pocket when CB1 is in the bound state, it fluctuates more because of the reduced interactions with binding pocket residues. Hence, PC1 and PC2 capture these large amplitude motions for both the ensembles. We also performed Kullback–Leibler (K-L) divergence analysis for each interresidue distance to determine the major differences between these two ensembles. K-L divergence analysis shows that the higher divergence occurs in intracellular TM5, TM6, and TM7 (Figs. 8B and S13), depicting the difference in activation microswitches for agonist and antagonist. CB1 has three well-known microswitches: toggle switch movement, intracellular TM6 outward movement (DRY motif salt bridge breaking), and TM7 (NPxxY motif) movement (Fig. S14). Comparison of the agonist bound active-like structure and inactive structure shows that these microswitches rearrange themselves during the protein activation. To capture these motions of the toggle switch, TM6, and TM7, TRP3566.48-ASP1632.50, LYS3436.35-ARG2143.50, and TYR1532.40-TYR3977.53 distances were calculated. Projection of activation microswitch distances onto a 2-D scatter plot shows the difference in activation microswitch movement for partial agonist and agonist-bound ensembles (Fig. 8, C and D). Although the partial agonist affects the movement of these microswitches, it cannot lead to the full activation of the receptor as described earlier. In contrast, an agonist stabilizes the receptor in an active-like pose. Difference in the intracellular microswitch movement might happen due to the difference in allosteric communication pathway between different microswitches. Specifically, we studied how the signaling is transduced from the toggle switch residue to the NPxxY region and compared the pathway for agonist and partial agonist-bound ensemble. An allosteric communication pathway was determined by estimating a mutual information network between residue dihedral angles (
      • Bhattacharya S.
      • Vaidehi N.
      Differences in allosteric communication pipelines in the inactive and active states of a GPCR.
      ,
      • Vaidehi N.
      • Bhattacharya S.
      Allosteric communication pipelines in G-protein-coupled receptors.
      ,
      • Nivedha A.K.
      • Tautermann C.S.
      • Bhattacharya S.
      • Lee S.
      • Casarosa P.
      • Kollak I.
      • Kiechle T.
      • Vaidehi N.
      Identifying functional hotspot residues for biased ligand design in G-protein-coupled receptors.
      ). Using this network, we estimated the shortest pathway between the toggle switch and the NPxxY region using Dijkstra’s algorithm (
      • Kapoor A.
      • Martinez-Rosell G.
      • Provasi D.
      • de Fabritiis G.
      • Filizola M.
      Dynamic and kinetic elements of μ-opioid receptor functional selectivity.
      ). Results show that the shortest paths for signal transduction between toggle switch and NPxxY region are different for partial agonist and agonist (Fig. S15, A and B). For the agonist, we observed an allosteric signaling pathway between the NPxxY region and the toggle switch transduce through the outward movement of intracellular TM6, whereas, due to the lack of intracellular TM6 movement, allosteric signal transduction happens from a different pathway for partial agonist. Therefore, these analyses show how the activation signaling differs between partial agonists and agonists as the ligands affect the toggle switch differently.
      Figure thumbnail gr8
      Figure 8Comparison of agonist (pink) and partial agonist (green) bound ensembles. A, PCA projection of agonist and partial agonist bound simulation onto first two components. B, important differences between agonist and partial agonist bound ensemble were derived by computing the Kullback-Leibler (KL) divergence along all inter-residue features. K-L divergence values were averaged over each residue to determine the contribution of each residue. C, 2-D scatter plot showing intracellular TM6 movement with respect to toggle switch movement (agonist: pink; partial agonist: green). D, 2-D scatter plot projecting NPxxY movement with respect to toggle switch movement (agonist: pink; partial agonist: green). Inactive and active-like agonist bound structure was shown as highlighted points on 2-D scatter plot. PCA, principal component analysis; TM6, transmembrane helix 6.

      Discussion

      In this study, we proposed two hypotheses for the partial agonism of THC molecule for Cannabinoid receptor 1. Our first hypothesis is based on different binding position of the THC inside the orthosteric binding pocket. Second hypothesis states that THC may only able to move TM6 partially. To test the hypotheses, we performed unbiased molecular dynamics simulation for THC binding to CB1. Our results support both our hypotheses. Simulations show that during binding THC is stabilized in agonist-like and antagonist-like poses. While binding in the antagonist-like poses, the aromatic group of THC orients itself as arm 2 of the antagonist, whereas alkyl chain of THC can take the conformation of arm 1 and arm 3. During the binding process, alkyl chain orientation is shown to be important factor of determining the binding of THC to agonist pose. When THC enters the receptor with positive side chain dihedral angle (C2-C3-C1'-C2'), it leads to the upward movement of N-loop and favors the THC binding to the agonist-like pose, whereas, with negative dihedral angle, THC increases the free energy barrier for transition from antagonist-like binding pose to agonist-like binding pose. Therefore, chemical modification, which can stabilize the side chain in perpendicular direction to the aromatic ring, may increase the agonistic property of the ligand.
      Over the years, various effects on THC alkyl chain modifications on binding and functionality have been proposed. Δ9-THC and Δ8-THC have similar binding affinity and agonistic property. It has been shown that Δ8-THC analog with cis double bond between C1 and C2 (which fixes the side chain in perpendicular direction) with same side chain length increases the agonistic property of the ligand (
      • Bow E.W.
      • Rimoldi J.M.
      The structure–function relationships of classical cannabinoids: CB1/CB2 modulation.
      ). Our simulations also support this experimental observation. With cis bond, THC analog can move to the agonist-like pose with less free energy barrier (Fig. 4A) and therefore has more agonist property, which bolsters our claim in the first hypothesis.
      To generalize our first hypothesis that THC can be stabilized at both agonist and antagonist pocket, we performed docking of other partial agonists on representative CB1 structures of THC bound agonist-binding pose and antagonist-binding pose (Fig. S17). Five available partial agonists are selected for GPCRdb database (
      • Pándy-Szekeres G.
      • Munk C.
      • Tsonkov T.M.
      • Mordalski S.
      • Harpsøe K.
      • Hauser A.S.
      • Bojarski A.J.
      • Gloriam D.E.
      GPCRdb in 2018: Adding GPCR structure models and ligands.
      ): magnolol, AM4089, NMP-4, (S)-Δ3-THC, (R)-Δ3-THC. Docking studies reveal stable docked poses for partial agonists in both the binding pockets. In the antagonist-binding pose, we observed that the partial agonists extend downward toward Na+-binding site similar to THC (
      • Dutta S.
      • Selvam B.
      • Shukla D.
      Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors.
      ). These partial agonists bind with similar affinity in agonist and antagonist-binding pocket. These docking results show that our hypothesis for partial agonism may be universally valid for other partial agonists for CB1.
      Our results also show that the THC can rotate the toggle switch TRP3566.48 conformation when binding to the agonist-like poses. The partially active conformation of CB1 matches well with CB2 inactive structure, which helps to explain yin–yang relation for some ligands between CB1 and CB2. Recently published NMR study on another class A GPCR, AA2AR, has shown that partial agonist rotates TRP2466.48 in a distinct conformation compared with full agonists or antagonists (
      • Eddy M.T.
      • Martin B.T.
      • Wüthrich K.
      A2A adenosine receptor partial agonism related to structural rearrangements in an activation microswitch.
      ). These experimental observations support our results obtained computationally. Furthermore, our simulations reveal that partially active toggle switch movement consequently affects the intracellular side of the receptor and leads to the partial outward movement of TM6. In future study, this partial movement of the TM6 can be experimentally validated by DEER spectroscopy (
      • Jeschke G.
      DEER distance measurements on proteins.
      ,
      • Manglik A.
      • Kim T.H.
      • Masureel M.
      • Altenbach C.
      • Yang Z.
      • Hilger D.
      • Lerch M.T.
      • Kobilka T.S.
      • Thian F.S.
      • Hubbell W.L.
      • Prosser R.S.
      • Kobilka B.K.
      Structural insights into the dynamic process of β 2 -adrenergic receptor signaling.
      ,
      • Mittal S.
      • Shukla D.
      Maximizing kinetic information gain of Markov state models for optimal design of spectroscopy experiments.
      ,
      • Mittal S.
      • Shukla D.
      Recruiting machine learning methods for molecular simulations of proteins.
      ).
      Overall, our present computational study with MD simulations and MSM is competent with previous experimental results and also provides new predictions for partial agonism of THC for CB1. The findings of the paper are crucial for guiding the chemical modification for cannabinoid ligands for future drug development.

      Experimental procedures

      System preparation

      Inactive structure of CB1 (PDB ID: 5TGZ (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      )) is used as a starting structure for MD simulation for the THC binding. This structure does not contain first 98 residues in the N-terminus. However, previous study has shown that most of the N-terminal residues (1–103) can be deleted without abolishing the ligand binding for synthetic cannabinoids (
      • Fay J.F.
      • Farrens D.L.
      The membrane proximal region of the cannabinoid receptor CB1N-terminus can allosterically modulate ligand affinity.
      ). Truncated N-terminus and C-terminus and unconnected residues of TM5 and TM6 are capped with neutral terminal residues (acetyl and methylamide groups). Hydrogen atoms are added to protein amino acid residues using reduce command of AMBER package (
      • Case D.A.
      • Aktulga H.M.
      • Belfon K.
      • Ben-Shalom I.Y.
      • Brozell S.R.
      • Cerutti D.S.
      • Cheatham T.E.
      • Cisneros G.A.
      • Cruzeiro V.W.D.
      • Darden T.A.
      • Duke R.E.
      • Giambasu G.
      • Gilson M.K.
      • Gohlke H.
      • Goetz A.W.
      • et al.
      AMBER 2018.
      ). Thermostabilized mutant residues are replaced with original residues using tleap (
      • Hua T.
      • Vemuri K.
      • Pu M.
      • Qu L.
      • Han G.W.
      • Wu Y.
      • Zhao S.
      • Shui W.
      • Li S.
      • Korde A.
      • Laprairie R.B.
      • Stahl E.L.
      • Ho J.H.
      • Zvonok N.
      • Zhou H.
      • et al.
      Crystal structure of the human cannabinoid receptor CB1.
      ). The modified protein structure is embedded in POPC bilayer using CHARMM-GUI (
      • Jo S.
      • Kim T.
      • Iyer V.G.
      • Im W.
      CHARMM-GUI: A web-based graphical user interface for CHARMM.
      ). Salt concentration of 150 mM (Na+ and Cl−) is to neutralize the system. MD system is solvated using TIP3P water model. 3-D structure of THC is obtained from PubChem in sdf format. Forcefield parameters of THC are obtained using antechamber (
      • Wang J.
      • Wolf R.M.
      • Caldwell J.W.
      • Kollman P.A.
      • Case D.A.
      Development and testing of a general amber force field.
      ,
      • Wang J.
      • Wang W.
      • Kollman P.A.
      • Case D.A.
      Automatic atom type and bond type perception in molecular mechanical calculations.
      ). THC is added to MD system using packmol (
      • Martínez L.
      • Andrade R.
      • Birgin E.G.
      • Martínez J.M.
      PACKMOL: A package for building initial configurations for molecular dynamics simulations.
      ) to generate the starting structure. To perform agonist (AM11542) bound simulation, active-like agonist-bound crystal structure (PDB ID: 5XRA (
      • Hua T.
      • Vemuri K.
      • Nikas S.P.
      • Laprairie R.B.
      • Wu Y.
      • Qu L.
      • Pu M.
      • Korde A.
      • Jiang S.
      • Ho J.H.
      • Han G.W.
      • Ding K.
      • Li X.
      • Liu H.
      • Hanson M.A.
      • et al.
      Crystal structures of agonist-bound human cannabinoid receptor CB1.
      )) was selected. Necessary modifications were performed as discussed as for THC-binding simulation. Agonist structure was parameterized using GAFF forcefield. CHARMM-GUI software was used to embed the agonist bound structure in POPC bilayer with NaCl solution (150 mM) in extracellular and intracellular regions.

      Simulation details

      MD simulations were performed using AMBER18. The MD systems were subjected to minimization with gradient descent and conjugate gradient algorithm for 5000 and 10,000 steps, respectively. Minimized systems were slowly heated from 0 K to 10 K and from 10 K to 300 K in NVT ensemble to increase the temperature at desired level. Each step was done for 1 ns period. To control the pressure of the heated system at 1 bar, NPT ensemble was employed. During modulation of temperature and pressure, protein backbone Cα atoms are restrained with a spring force. Furthermore, total 50 ns of equilibration was performed in NPT ensemble for each system to maintain the temperature and pressure at desired level of 300 K at 1 bar without any restraint force. Production runs were also performed in NPT ensemble. The temperature and pressure of the systems were maintained by Berendsen thermostat and barostat (
      • Berendsen H.J.
      • Postma J.V.
      • van Gunsteren W.F.
      • DiNola A.
      • Haak J.R.
      Molecular dynamics with coupling to an external bath.
      ). Simulation timestep of 2 fs was used for MD simulation. As hydrogen atom is the lightest atom in the system, vibrational frequency of the hydrogen can be lower than the specified timescale. It can create instability in the system due to the large fluctuation of hydrogen atom. Therefore, SHAKE algorithm (
      • Kräutler V.
      • Van Gunsteren W.F.
      • Hünenberger P.H.
      A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations.
      ) was used to put restraint in the movement of hydrogen atom by implementing Lagrangian multiplier. For nonbonded force calculation, 10 Å cutoff distance was used. Periodic boundary condition was applied for all simulations. To consider the nonbonded long-range interactions, particle mesh Ewald method was implemented (
      • Essmann U.
      • Perera L.
      • Berkowitz M.L.
      • Darden T.
      • Lee H.
      • Pedersen L.G.
      A smooth particle mesh Ewald method.
      ). Simulations were performed with Amber FF14SB (
      • Maier J.A.
      • Martinez C.
      • Kasavajhala K.
      • Wickstrom L.
      • Hauser K.E.
      • Simmerling C.
      ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB.
      ) and GAFF (
      • Wang J.
      • Wolf R.M.
      • Caldwell J.W.
      • Kollman P.A.
      • Case D.A.
      Development and testing of a general amber force field.
      ) forcefield parameters.

      Adaptive sampling

      Ligand binding to GPCR is rare event compared with MD simulations timescale. To capture the entire THC-binding process, adaptive sampling technique is utilized. This approach has been shown to sample the conformational ensemble of a variety of biological systems (
      • Selvam B.
      • Shamsi Z.
      • Shukla D.
      Universality of the sodium ion binding mechanism in class AG-protein-coupled receptors.
      ,
      • Selvam B.
      • Mittal S.
      • Shukla D.
      Free energy landscape of the complete transport cycle in a key bacterial transporter.
      ,
      • Selvam B.
      • Yu Y.-C.
      • Chen L.-Q.
      • Shukla D.
      Molecular basis of the glucose transport mechanism in plants.
      ,
      • Kohlhoff K.J.
      • Shukla D.
      • Lawrenz M.
      • Bowman G.R.
      • Konerding D.E.
      • Belov D.
      • Altman R.B.
      • Pande V.S.
      Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways.
      ) including the ligand-binding process (
      • Dutta S.
      • Selvam B.
      • Shukla D.
      Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors.
      ,
      • Shukla S.
      • Zhao C.
      • Shukla D.
      Dewetting controls plant hormone perception and initiation of drought resistance signaling.
      ,
      • Chen J.
      • White A.
      • Nelson D.C.
      • Shukla D.
      Role of substrate recognition in modulating strigolactone receptor selectivity in witchweed.
      ,
      • Aldukhi F.
      • Deb A.
      • Zhao C.
      • Moffett A.S.
      • Shukla D.
      Molecular mechanism of brassinosteroid perception by the plant growth receptor BRI1.
      ,
      • Lawrenz M.
      • Shukla D.
      • Pande V.S.
      Cloud computing approaches for prediction of ligand binding poses and pathways.
      ). In this technique, data from one round of simulation are represented by conformational space of protein by using biological relevant features. The features are used to cluster the protein space using k-means clustering algorithm. The structures for next round of simulation are selected randomly from the cluster centers with lowest count of data points. This process is repeated in each round. For this case, intracellular and extracellular helical distances and THC-binding distance from TM5 were used as adaptive sampling feature matrices. Although adaptive sampling helps to parallelize the simulation by sampling from lower probability space, it affects the ensemble distribution and therefore brings sampling error in free energy calculation. To overcome this caveat, MSM (discussed below) is implemented to remove the sampling the bias from the simulation data. A total 143 μs of simulation data were collected using adaptive sampling protocol (Table S2). Agonist-bound simulations were also performed using adaptive sampling protocol using similar matrices as partial agonist THC-binding simulation. A total 17 μs of simulation data were collected for agonist-bound simulation.

      Markov state model

      The theory of MSM (
      • Bowman G.R.
      • Pande V.S.
      • Noé F.
      ,
      • Prinz J.-H.
      • Wu H.
      • Sarich M.
      • Keller B.
      • Senne M.
      • Held M.
      • Chodera J.D.
      • Schütte C.
      • Noé F.
      Markov models of molecular kinetics: Generation and validation.
      ,
      • Shukla D.
      • Hernández C.X.
      • Weber J.K.
      • Pande V.S.
      Markov state models provide insights into dynamic modulation of protein function.
      ) depicts that in a sequence of event, the transition probability of moving from one state to another depends only on the present state and not on the path the system has taken to be there. The trajectory of MD simulation follows the same principle. According to the verlet algorithm, the evolution coordinate and momenta of atoms in system only depend on the present state. Therefore, MD trajectory can be assumed to be Markovian and according to its nature, the probability distribution of the protein conformational ensemble can be calculated from the transition probability matrix between the states. Each element Tij of transition probability matrix calculates the probability of transition using the equation Tij=Cij/j=1nCij where Cij is the count of jump between the i and j, and Ci is the count of frame in state i. To make the transition probability matrix statistically significant, states with conformation are clustered together into microstates assuming that there is no large energy barrier in the same microstate. However, clustering of conformational space increases the memory of each state, which can invoke “non-markovinity” in our system. This issue can be overcome by increasing the lag time (τ) such that it preserves the Markov property and hence satisfies the equation p(t+τ)=p(t)T(τ) where p(t+τ) and p(t) are vectors representing the probability of the microstates. Pyemma (
      • Wehmeyer C.
      • Scherer M.K.
      • Hempel T.
      • Husic B.E.
      • Olsson S.
      • Noé F.
      Introduction to Markov state modeling with the PyEMMA softwar v1. 0.
      ) package was used to construct MSM. We selected 23 biologically important features to capture the THC binding and protein conformational changes (Table S1 and Fig. S22). Features were transformed into time-lagged independent components (tics) (
      • Noé F.
      • Clementi C.
      Kinetic distance and kinetic maps from molecular dynamics simulation.
      ,
      • Pérez-Hernández G.
      • Paul F.
      • Giorgino T.
      • De Fabritiis G.
      • Noé F.
      Identification of slow molecular order parameters for Markov model construction.
      ) to find the slowest components. Tic components were shown to correlate well with important features for THC binding (Fig. S18, AC). The transformed data were clustered using k-means clustering algorithm. For MSM building, lag time was chosen by finding out logarithmic convergence of process timescales computed from MSM eigenvalues (process timescale t=τlogeλ; λ is the eigenvalue) (Fig. S19). To find out the optimal number of clusters and optimal tic variance, MSM was subjected to VAMP2 (
      • Noé F.
      • Clementi C.
      Kinetic distance and kinetic maps from molecular dynamics simulation.
      ) scoring to measure the kinetic variance (Fig. S20). Based on transition matrix, MSM predicts the population (stationary probability, πi) of each clustered state, which is needed to calculate thermodynamics (Gi = −kbT logeπi) of the process. Using similar steps, optimal MSM was built for agonist-bound simulation with 65% tic variance and 25 ns of MSM lag time.

      Transition path theory

      To obtain the kinetics for transition between MSM macrostates, TPT (
      • Metzner P.
      • Schütte C.
      • Vanden-Eijnden E.
      Transition path theory for Markov jump processes.
      ,
      • Weinan E.
      • Vanden-Eijnden E.
      Towards a theory of transition paths.
      ) was implemented. TPT calculates the flux between two macrostates (which consists of one or multiply MSM states) using MSM transition matrix (
      • Meng Y.
      • Shukla D.
      • Pande V.S.
      • Roux B.
      Transition path theory analysis of c-Src kinase activation.
      ). The flux (FAB) is given by the equation FAB=iAjAπiTijqj+ where qj+ is probability that state j will reach macrostate B before A (
      • Noé F.
      • Schütte C.
      • Vanden-Eijnden E.
      • Reich L.
      • Weikl T.R.
      Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations.
      ). From the flux calculation, we computed MFPT using equation MFPT = τπA/FAB. where πA is the probability the system was in macrostate A. TPT calculation was performed with Pyemma (
      • Wehmeyer C.
      • Scherer M.K.
      • Hempel T.
      • Husic B.E.
      • Olsson S.
      • Noé F.
      Introduction to Markov state modeling with the PyEMMA softwar v1. 0.
      ) package. Here, the macrostates were defined manually. Five MSM states with highest raw count inside area of interest are defined as macrostates. For example, macrostate 4 in Figure 5A consisted of five MSM states with highest raw count inside the area where THC(C1')-TYR2755.39(Cα) distance is between 9 Å and 13 Å & THC dihedral angle is between 60° and 120°.

      Trajectory analysis

      All the feature calculation and data processing from the MD simulation were done CPPtraj (
      • Roe D.R.
      • Cheatham III, T.E.
      PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data.
      ) and MDtraj (
      • McGibbon R.T.
      • Beauchamp K.A.
      • Harrigan M.P.
      • Klein C.
      • Swails J.M.
      • Hernández C.X.
      • Schwantes C.R.
      • Wang L.-P.
      • Lane T.J.
      • Pande V.S.
      MDTraj: A modern open library for the analysis of molecular dynamics trajectories.
      ). For the trajectory visualization and analysis, VMD (
      • Humphrey W.
      • Dalke A.
      • Schulten K.
      VMD: Visual molecular dynamics.
      ) and pymol v2.1 (Schrödinger, LLC) software package were used. To calculate ligand and protein residue contact, GetContacts package (https://getcontacts.github.io/) was implemented.

      Alchemical calculation

      Alchemical calculations were performed on a reference frame selected from the antagonist-bound pose with thermodynamics integration (TI) technique (
      • Mey A.S.
      • Allen B.K.
      • Macdonald H.E.B.
      • Chodera J.D.
      • Hahn D.F.
      • Kuhn M.
      • Michel J.
      • Mobley D.L.
      • Naden L.N.
      • Prasad S.
      • Rizzi A.
      • Scheen J.
      • Shirts M.R.
      • Tresadern G.
      • Xu H.
      Best practices for alchemical free energy calculations [article v1.0].
      ). To prepare the system for TI calculations, alanine mutations were added to the system using tiMerge function of ParmEd software package. Two different systems were prepared for two mutations (PHE1742.61ALA and PHE1772.64ALA). AMBER PMEMD engine was used to carry out TI calculation. Simulation conditions were maintained same as the original simulation condition. For TI, 11 λ windows were prepared between 0 and 1 with a gap of 0.1. Softcore potential was used for van der Waals (vdw) and charge contribution for the dummy atoms. In each window, 2 ns simulation data were generated. Free energy estimation was performed with pyamber package (
      • Song L.F.
      • Lee T.-S.
      • Zhu C.
      • York D.M.
      • Merz K.M.
      Using AMBER18 for relative free energy calculations.
      ).

      Kinetic Monte Carlo simulation

      kMC simulation is a stochastic method to visualize the evolution of a system based on probability of the transition between different states. In this case, transition probability of different states was obtained from MSM. There are few steps to implement kMC in MSM weighted MD data. First, we generated a random number (R) between 0 and 1. Second, we described the possible transition event from i as discrete cumulation probability distribution. Therefore, the cumulative transition probability between i and j can be written as sij=k=1jTik. If R lies between si,j and si,j+1, then system transition happens from ith to jth state. Time required for each transition is taken to be the lag time considered for MSM. Same procedure is repeated for desired number of steps. To build a trajectory from kMC state evolution, we picked a random frame at each step from the chosen state. As we were interested to observe the ligand-binding dynamics, kMC simulation is started from a state where THC was 25 Å away from TYR2755.39(Cα).

      Data-driven modeling

      Data-driven analysis was performed to compare agonist and partial agonist-bound ensembles. For THC, 20,000 frames were selected from clusters where THC is bound to the agonist-bound pose based on cluster probability calculated using MSM. Similarly, 20,000 frames were selected from agonist-bound simulation based on cluster probability calculated using MSM for the agonist-bound simulation. All closest heavy carton distances were calculated using MDtraj (
      • McGibbon R.T.
      • Beauchamp K.A.
      • Harrigan M.P.
      • Klein C.
      • Swails J.M.
      • Hernández C.X.
      • Schwantes C.R.
      • Wang L.-P.
      • Lane T.J.
      • Pande V.S.
      MDTraj: A modern open library for the analysis of molecular dynamics trajectories.
      ). PCA was performed using sklearn package of python. Here, symmetrized version of K-L divergence was used to calculate the difference in distance distribution of two ensembles to find out major distinctions. Lastly, to find the activation path, normalized mutual information was calculated for dihedral angles for every residue pair using DihedralMutualInformation class of mdentropy package (
      • Hernández C.X.
      • Pande V.S.
      MDEntropy: Information-theoretic analyses for molecular dynamics.
      ). With this information, we developed a weighted undirected graph where each residue denotes a node. The edges of the graph were created between two residue pairs with the criteria that the residues are within 6 Å from each other and there is nonzero mutual information between two residues. From this graph, we calculated the shortest path between toggle switch (TRP3566.48) and NPxxY region (TYR3977.53) using Dijkstra’s algorithm implementation of the Networkx python package (
      • Hagberg A.
      • Swart P.
      • Chult D.S.
      ,
      • Dijkstra E.W.
      A note on two problems in connexion with graphs.
      ).

      Docking study

      Docking study was performed using Auto dock vina (
      • Trott O.
      • Olson A.J.
      AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading.
      ) software. 3-D structures of partial agonists are selected from PubChem in sdf format. Antechamber is used to convert the ligands in mol2 format and add partial changes. Then, Auto dock is used to convert the ligands into required pdbqt format. For the docking of the ligands in both agonist and antagonist pose, the box size is calculated from THC-binding structure in respective poses. To dock the ligand in antagonist and agonist-bound poses, MD predicted structures are used.

      Error analysis

      To determine error in our thermodynamics (free energy, conditional probability) and kinetics (TPT) calculations, we performed bootstrap analysis on MD data (
      • Sultan M.M.
      • Kiss G.
      • Pande V.S.
      Towards simple kinetic models of functional dynamics for a kinase subfamily.
      ). In each bootstrap sample, we randomly picked N trajectories, where N is equal to 80% of total number of trajectories. We kept the state labeling same as the original MSM (Fig. S21). For every bootstrap sample, MSM is computed to determine thermodynamics and kinetics. In total, 200 bootstrap samples are generated for error calculations.

      Code and data availability

      MSM feature matrix, final MSM object, codes, and pdb files used to generate figures can be obtained from https://github.com/ShuklaGroup/THC_binding.git.

      Supporting information

      This article contains supporting information.

      Conflict of interest

      The authors declare that they have no conflicts of interest with the contents of this article.

      Acknowledgments

      This research is performed as a part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993 ) the State of Illinois, and as of December, 2019, the National Geospatial-Intelligence Agency. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.

      Author contributions

      D. S. conceptualization; S. D. and B. S. formal analysis; D. S. funding acquisition; S. D. investigation; D. S. project administration; A. D. and D. S. supervision; S. D. writing—original draft; S. D., B. S., A. D., and D. S. writing—review and editing.

      Funding and additional information

      D. S. acknowledges support from NSF Early CAREER Award ( MCB-1845606 ).

      Supporting information

      References

        • Hilger D.
        • Masureel M.
        • Kobilka B.K.
        Structure and dynamics of GPCR signaling complexes.
        Nat. Struct. Mol. Biol. 2018; 25: 4
        • Rosenbaum D.M.
        • Rasmussen S.G.
        • Kobilka B.K.
        The structure and function of G-protein-coupled receptors.
        Nature. 2009; 459: 356-363
        • Kumar K.K.
        • Shalev-Benami M.
        • Robertson M.J.
        • Hu H.
        • Banister S.D.
        • Hollingsworth S.A
        • Latorraca N.R.
        • Kato H.E.
        • Hilger D.
        • Maeda S.
        • Weis W.I.
        • Farrens D.L.
        • Dror R.O.
        • Malhotra S.V.
        • Kobilka B.K.
        • et al.
        Structure of a signaling cannabinoid receptor 1-G protein complex.
        Cell. 2019; 176: 448-458.e12
        • Ibsen M.S.
        • Finlay D.B.
        • Patel M.
        • Javitch J.A.
        • Glass M.
        • Grimsey N.L.
        Cannabinoid CB1 and CB2 receptor-mediated arrestin translocation: Species, subtype, and agonist-dependence.
        Front. Pharmacol. 2019; 10: 350
        • Zou S.
        • Kumar U.
        Cannabinoid receptors and the endocannabinoid system: Signaling and function in the central nervous system.
        Int. J. Mol. Sci. 2018; 19: 833
        • Pertwee R.G.
        Pharmacology of cannabinoid receptor ligands.
        Curr. Med. Chem. 1999; 6: 635-664
        • Pertwee R.G.
        Cannabinoid pharmacology: The first 66 years.
        Br. J. Pharmacol. 2006; 147: S163-S171
        • Pertwee R.G.
        • Ross R.
        Cannabinoid receptors and their ligands.
        Prostaglandins Leukot. Essent. Fatty Acids. 2002; 66: 101-121
        • MacLennan S.
        • Reynen P.
        • Kwan J.
        • Bonhaus D.
        Evidence for inverse agonism of SR141716A at human recombinant cannabinoid CB1 and CB2 receptors.
        Br. J. Pharmacol. 1998; 124: 619-622
        • Pertwee R.G.
        Inverse agonism and neutral antagonism at cannabinoid CB1 receptors.
        Life Sci. 2005; 76: 1307-1324
        • Adams A.J.
        • Banister S.D.
        • Irizarry L.
        • Trecki J.
        • Schwartz M.
        • Gerona R.
        “Zombie” outbreak caused by the synthetic cannabinoid AMB-FUBINACA in New York.
        N. Engl. J. Med. 2017; 376: 235-242
        • Sam A.H.
        • Salem V.
        • Ghatei M.A.
        Rimonabant: From RIO to ban.
        J. Obes. 2011; 2011: 432607
        • Dutta S.
        • Selvam B.
        • Shukla D.
        Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors.
        ACS Chem. Neurosci. 2022; 13: 379-389
        • Badowski M.E.
        • Yanful P.K.
        Dronabinol oral solution in the management of anorexia and weight loss in AIDS and cancer.
        Ther. Clin. Risk Manag. 2018; 14: 643
        • Walther S.
        • Halpern M.
        Cannabinoids and dementia: A review of clinical and preclinical data.
        Pharmaceuticals. 2010; 3: 2689-2708
        • Hua T.
        • Li X.
        • Wu L.
        • Iliopoulos-Tsoutsouvas C.
        • Wang Y.
        • Wu M.
        • Shen L.
        • Brust C.A.
        • Nikas S.P.
        • Song F.
        • Song X.
        • Yuan S.
        • Sun Q.
        • Wu Y.
        • Jiang S.
        • et al.
        Activation and signaling mechanism revealed by cannabinoid receptor-Gi complex structures.
        Cell. 2020; 180: 655-665.e18
        • Hua T.
        • Vemuri K.
        • Pu M.
        • Qu L.
        • Han G.W.
        • Wu Y.
        • Zhao S.
        • Shui W.
        • Li S.
        • Korde A.
        • Laprairie R.B.
        • Stahl E.L.
        • Ho J.H.
        • Zvonok N.
        • Zhou H.
        • et al.
        Crystal structure of the human cannabinoid receptor CB1.
        Cell. 2016; 167: 750-762.e14
        • Hua T.
        • Vemuri K.
        • Nikas S.P.
        • Laprairie R.B.
        • Wu Y.
        • Qu L.
        • Pu M.
        • Korde A.
        • Jiang S.
        • Ho J.H.
        • Han G.W.
        • Ding K.
        • Li X.
        • Liu H.
        • Hanson M.A.
        • et al.
        Crystal structures of agonist-bound human cannabinoid receptor CB1.
        Nature. 2017; 547: 468-471
        • Shao Z.
        • Yin J.
        • Chapman K.
        • Grzemska M.
        • Clark L.
        • Wang J.
        • Rosenbaum D.M.
        High-resolution crystal structure of the human CB1 cannabinoid receptor.
        Nature. 2016; 540: 602-606
        • Shao Z.
        • Yan W.
        • Chapman K.
        • Ramesh K.
        • Ferrell A.J.
        • Yin J.
        • Wang X.
        • Xu Q.
        • Rosenbaum D.M.
        Structure of an allosteric modulator bound to the CB1 cannabinoid receptor.
        Nat. Chem. Biol. 2019; 15: 1199-1205
        • Li X.
        • Shen L.
        • Hua T.
        • Liu Z.-J.
        Structural and functional insights into cannabinoid receptors.
        Trends Pharmacol. Sci. 2020; 41: 665-677
        • Saleh N.
        • Hucke O.
        • Kramer G.
        • Schmidt E.
        • Montel F.
        • Lipinski R.
        • Ferger B.
        • Clark T.
        • Hildebrand P.W.
        • Tautermann C.S.
        Multiple binding sites contribute to the mechanism of mixed agonistic and positive allosteric modulators of the cannabinoid CB1 receptor.
        Angew. Chem. Int. Ed. Engl. 2018; 130: 2610-2615
        • Rasmussen S.G.F.
        • DeVree B.T.
        • Zou Y.
        • Kruse A.C.
        • Chung K.Y.
        • Kobilka T.S.
        • Thian F.S.
        • Chae P.S.
        • Pardon E.
        • Calinski D.
        • Mathiesen J.M.
        • Shah S.T.A.
        • Lyons J.A.
        • Caffrey M.
        • Gellman S.H.
        • et al.
        Crystal structure of the β2 adrenergic receptor–Gs protein complex.
        Nature. 2011; 477: 549-555
        • Masureel M.
        • Zou Y.
        • Picard L.P.
        • Westhuizen E.V.D.
        • Mahoney J.P.
        • Rodrigues J.P.G.L.M.
        • Mildorf T.J.
        • Dror R.O.
        • Shaw D.E.
        • Bouvier M.
        • Pardon E.
        • Steyaert J.
        • Sunahara R.K.
        • Weis W.I.
        • Zhang C.
        • et al.
        Structural insights into binding specificity, efficacy and bias of a β2AR partial agonist.
        Nat. Chem. Biol. 2018; 14: 1059-1066
        • Cherezov V.
        • Rosenbaum D.M.
        • Hanson M.A.
        • Rasmussen S.G.F.
        • Thian F.S.
        • Kobilka T.S.
        • Choi H.-J.
        • Kuhn P.
        • Weis W.I.
        • Kobilka B.K.
        • Stevens R.C.
        High-resolution crystal structure of an engineered human β 2 -adrenergic G protein–coupled receptor.
        Science. 2007; 318: 1258-1265
        • Katritch V.
        • Fenalti G.
        • Abola E.E.
        • Roth B.L.
        • Cherezov V.
        • Stevens R.C.
        Allosteric sodium in class A GPCR signaling.
        Trends Biochem. Sci. 2014; 39: 233-244
        • Selvam B.
        • Shamsi Z.
        • Shukla D.
        Universality of the sodium ion binding mechanism in class AG-protein-coupled receptors.
        Angew. Chem. Int. Ed. Engl. 2018; 130: 3102-3107
        • Shim J.-Y.
        • Bertalovitz A.C.
        • Kendall D.A.
        Identification of essential cannabinoid-binding domains.
        J. Biol. Chem. 2011; 286: 33422-33435
        • Song L.F.
        • Lee T.-S.
        • Zhu C.
        • York D.M.
        • Merz K.M.
        Using AMBER18 for relative free energy calculations.
        J. Chem. Inf. Model. 2019; 59: 3128-3135
        • Bow E.W.
        • Rimoldi J.M.
        The structure–function relationships of classical cannabinoids: CB1/CB2 modulation.
        Perspect. Med. Chem. 2016; 8: 17-39
        • Li X.
        • Hua T.
        • Vemuri K.
        • Ho J.H.
        • Wu Y.
        • Wu L.
        • Popov P.
        • Benchama O.
        • Zvonok N.
        • Locke K.
        • Qu L.
        • Han G.W.
        • Iyer M.R.
        • Cinar R.
        • Coffey N.J.
        • et al.
        Crystal structure of the human cannabinoid receptor CB2.
        Cell. 2019; 176: 459-467.e13
        • Fleetwood O.
        • Carlsson J.
        • Delemotte L.
        Identification of ligand-specific G protein-coupled receptor states and prediction of downstream efficacy via data-driven modeling.
        Elife. 2020; https://doi.org/10.7554/elife.60715.sa2
        • Bhattacharya S.
        • Vaidehi N.
        Differences in allosteric communication pipelines in the inactive and active states of a GPCR.
        Biophys. J. 2014; 107: 422-434
        • Vaidehi N.
        • Bhattacharya S.
        Allosteric communication pipelines in G-protein-coupled receptors.
        Curr. Opin. Pharmacol. 2016; 30: 76-83
        • Nivedha A.K.
        • Tautermann C.S.
        • Bhattacharya S.
        • Lee S.
        • Casarosa P.
        • Kollak I.
        • Kiechle T.
        • Vaidehi N.
        Identifying functional hotspot residues for biased ligand design in G-protein-coupled receptors.
        Mol. Pharmacol. 2018; 93: 288-296
        • Kapoor A.
        • Martinez-Rosell G.
        • Provasi D.
        • de Fabritiis G.
        • Filizola M.
        Dynamic and kinetic elements of μ-opioid receptor functional selectivity.
        Sci. Rep. 2017; 7: 11255
        • Pándy-Szekeres G.
        • Munk C.
        • Tsonkov T.M.
        • Mordalski S.
        • Harpsøe K.
        • Hauser A.S.
        • Bojarski A.J.
        • Gloriam D.E.
        GPCRdb in 2018: Adding GPCR structure models and ligands.
        Nucleic Acids Res. 2018; 46: D440-D446
        • Eddy M.T.
        • Martin B.T.
        • Wüthrich K.
        A2A adenosine receptor partial agonism related to structural rearrangements in an activation microswitch.
        Structure. 2021; 29: 170-176
        • Jeschke G.
        DEER distance measurements on proteins.
        Annu. Rev. Phys. Chem. 2012; 63: 419-446
        • Manglik A.
        • Kim T.H.
        • Masureel M.
        • Altenbach C.
        • Yang Z.
        • Hilger D.
        • Lerch M.T.
        • Kobilka T.S.
        • Thian F.S.
        • Hubbell W.L.
        • Prosser R.S.
        • Kobilka B.K.
        Structural insights into the dynamic process of β 2 -adrenergic receptor signaling.
        Cell. 2015; 161: 1101-1111
        • Mittal S.
        • Shukla D.
        Maximizing kinetic information gain of Markov state models for optimal design of spectroscopy experiments.
        J. Phys. Chem. B. 2018; 122: 10793-10805
        • Mittal S.
        • Shukla D.
        Recruiting machine learning methods for molecular simulations of proteins.
        Mol. Simul. 2018; 44: 891-904
        • Fay J.F.
        • Farrens D.L.
        The membrane proximal region of the cannabinoid receptor CB1N-terminus can allosterically modulate ligand affinity.
        Biochemistry. 2013; 52: 8286-8294
        • Case D.A.
        • Aktulga H.M.
        • Belfon K.
        • Ben-Shalom I.Y.
        • Brozell S.R.
        • Cerutti D.S.
        • Cheatham T.E.
        • Cisneros G.A.
        • Cruzeiro V.W.D.
        • Darden T.A.
        • Duke R.E.
        • Giambasu G.
        • Gilson M.K.
        • Gohlke H.
        • Goetz A.W.
        • et al.
        AMBER 2018.
        University of California, San Francisco, CA2018
        • Jo S.
        • Kim T.
        • Iyer V.G.
        • Im W.
        CHARMM-GUI: A web-based graphical user interface for CHARMM.
        J. Comput. Chem. 2008; 29: 1859-1865
        • Wang J.
        • Wolf R.M.
        • Caldwell J.W.
        • Kollman P.A.
        • Case D.A.
        Development and testing of a general amber force field.
        J. Comput. Chem. 2004; 25: 1157-1174
        • Wang J.
        • Wang W.
        • Kollman P.A.
        • Case D.A.
        Automatic atom type and bond type perception in molecular mechanical calculations.
        J. Mol. Graph. Model. 2006; 25: 247-260
        • Martínez L.
        • Andrade R.
        • Birgin E.G.
        • Martínez J.M.
        PACKMOL: A package for building initial configurations for molecular dynamics simulations.
        J. Comput. Chem. 2009; 30: 2157-2164
        • Berendsen H.J.
        • Postma J.V.
        • van Gunsteren W.F.
        • DiNola A.
        • Haak J.R.
        Molecular dynamics with coupling to an external bath.
        J. Chem. Phys. 1984; 81: 3684-3690
        • Kräutler V.
        • Van Gunsteren W.F.
        • Hünenberger P.H.
        A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations.
        J. Comput. Chem. 2001; 22: 501-508
        • Essmann U.
        • Perera L.
        • Berkowitz M.L.
        • Darden T.
        • Lee H.
        • Pedersen L.G.
        A smooth particle mesh Ewald method.
        J. Chem. Phys. 1995; 103: 8577-8593
        • Maier J.A.
        • Martinez C.
        • Kasavajhala K.
        • Wickstrom L.
        • Hauser K.E.
        • Simmerling C.
        ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB.
        J. Chem. Theory Comput. 2015; 11: 3696-3713
        • Selvam B.
        • Mittal S.
        • Shukla D.
        Free energy landscape of the complete transport cycle in a key bacterial transporter.
        ACS Cent. Sci. 2018; 4: 1146-1154
        • Selvam B.
        • Yu Y.-C.
        • Chen L.-Q.
        • Shukla D.
        Molecular basis of the glucose transport mechanism in plants.
        ACS Cent. Sci. 2019; 5: 1085-1096
        • Kohlhoff K.J.
        • Shukla D.
        • Lawrenz M.
        • Bowman G.R.
        • Konerding D.E.
        • Belov D.
        • Altman R.B.
        • Pande V.S.
        Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways.
        Nat. Chem. 2014; 6: 15
        • Shukla S.
        • Zhao C.
        • Shukla D.
        Dewetting controls plant hormone perception and initiation of drought resistance signaling.
        Structure. 2019; 27: 692-702
        • Chen J.
        • White A.
        • Nelson D.C.
        • Shukla D.
        Role of substrate recognition in modulating strigolactone receptor selectivity in witchweed.
        biorXiv. 2020; ([preprint])https://doi.org/10.1101/2020.07.28.225722
        • Aldukhi F.
        • Deb A.
        • Zhao C.
        • Moffett A.S.
        • Shukla D.
        Molecular mechanism of brassinosteroid perception by the plant growth receptor BRI1.
        J. Phys. Chem. B. 2019; 124: 355-365
        • Lawrenz M.
        • Shukla D.
        • Pande V.S.
        Cloud computing approaches for prediction of ligand binding poses and pathways.
        Sci. Rep. 2015; 5: 7918
        • Bowman G.R.
        • Pande V.S.
        • Noé F.
        An introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation. Vol. 797. Springer Science & Business Media, Berlin2013
        • Prinz J.-H.
        • Wu H.
        • Sarich M.
        • Keller B.
        • Senne M.
        • Held M.
        • Chodera J.D.
        • Schütte C.
        • Noé F.
        Markov models of molecular kinetics: Generation and validation.
        J. Chem. Phys. 2011; 134: 174105
        • Shukla D.
        • Hernández C.X.
        • Weber J.K.
        • Pande V.S.
        Markov state models provide insights into dynamic modulation of protein function.
        Acc. Chem. Res. 2015; 48: 414-422
        • Wehmeyer C.
        • Scherer M.K.
        • Hempel T.
        • Husic B.E.
        • Olsson S.
        • Noé F.
        Introduction to Markov state modeling with the PyEMMA softwar v1. 0.
        LiveCoMS. 2018; 1: 1-12
        • Noé F.
        • Clementi C.
        Kinetic distance and kinetic maps from molecular dynamics simulation.
        J. Chem. Theory Comput. 2015; 11: 5002-5011
        • Pérez-Hernández G.
        • Paul F.
        • Giorgino T.
        • De Fabritiis G.
        • Noé F.
        Identification of slow molecular order parameters for Markov model construction.
        J. Chem. Phys. 2013; 139015102
        • Metzner P.
        • Schütte C.
        • Vanden-Eijnden E.
        Transition path theory for Markov jump processes.
        Multiscale Model. Simul. 2009; 7: 1192-1219
        • Weinan E.
        • Vanden-Eijnden E.
        Towards a theory of transition paths.
        J. Stat. Phys. 2006; 123: 503
        • Meng Y.
        • Shukla D.
        • Pande V.S.
        • Roux B.
        Transition path theory analysis of c-Src kinase activation.
        Proc. Natl. Acad. Sci. U. S. A. 2016; 113: 9193-9198
        • Noé F.
        • Schütte C.
        • Vanden-Eijnden E.
        • Reich L.
        • Weikl T.R.
        Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations.
        Proc. Natl. Acad. Sci. U. S. A. 2009; 106: 19011-19016
        • Roe D.R.
        • Cheatham III, T.E.
        PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data.
        J. Chem. Theory Comput. 2013; 9: 3084-3095
        • McGibbon R.T.
        • Beauchamp K.A.
        • Harrigan M.P.
        • Klein C.
        • Swails J.M.
        • Hernández C.X.
        • Schwantes C.R.
        • Wang L.-P.
        • Lane T.J.
        • Pande V.S.
        MDTraj: A modern open library for the analysis of molecular dynamics trajectories.
        Biophys. J. 2015; 109: 1528-1532
        • Humphrey W.
        • Dalke A.
        • Schulten K.
        VMD: Visual molecular dynamics.
        J. Mol. Graph. 1996; 14: 33-38
        • Mey A.S.
        • Allen B.K.
        • Macdonald H.E.B.
        • Chodera J.D.
        • Hahn D.F.
        • Kuhn M.
        • Michel J.
        • Mobley D.L.
        • Naden L.N.
        • Prasad S.
        • Rizzi A.
        • Scheen J.
        • Shirts M.R.
        • Tresadern G.
        • Xu H.
        Best practices for alchemical free energy calculations [article v1.0].
        Living J. Comput. Mol. Sci. 2020; 2: 18378
        • Hernández C.X.
        • Pande V.S.
        MDEntropy: Information-theoretic analyses for molecular dynamics.
        J. Open Source Softw. 2017; 2: 427
        • Hagberg A.
        • Swart P.
        • Chult D.S.
        Exploring Network Structure, Dynamics, and Function Using Networkx. Proceedings of the 7th Python in Science Conference (SciPy2008), Pasadena, CA2008
        • Dijkstra E.W.
        A note on two problems in connexion with graphs.
        Numer. Math. 1959; 1: 269-271
        • Trott O.
        • Olson A.J.
        AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading.
        J. Comput. Chem. 2010; 31: 455-461
        • Sultan M.M.
        • Kiss G.
        • Pande V.S.
        Towards simple kinetic models of functional dynamics for a kinase subfamily.
        Nat. Chem. 2018; 10: 903-909