G protein-coupled receptors (GPCRs) constitute one of the largest protein families in the human genome and represent the largest class of drug targets (Hauser et al., 2017; Nieto Gutierrez and McDonald, 2018; Sriram and Insel, 2018). Typically, a GPCR detects molecules outside the cell with its extracellular part, and undergoes global conformational changes to transduce cell signaling. The molecular mechanism of GPCR activation has been intensively studied. In class A receptors, though the ligand binding sites are highly diverse, the movements of the cytoplasmic part during activation are similar, featuring the outward movement of TM6 and the distortion of TM7 to open the space for G protein binding (Rasmussen et al., 2011b; Latorraca et al., 2017). For the connector region, which links the extracellular part and the cytoplasmic part, the P^5.50-I^3.40-F^6.44 [Ballesteros–Weinstein numbering (Ballesteros and Weinstein, 1995)] motif is a well-known microswitch (Rasmussen et al., 2011a; Wacker et al., 2013) that involves residues on three transmembrane (TM) helices. However, this P-I-F motif is not universal—among the 286 human nonolfactory class A GPCRs recorded in GPCRdb (Kooistra et al., 2021), only 94 possess this motif (Figure 1A). Therefore, GPCRs may have other mechanisms to relay the movements of the extracellular part to the cytoplasmic part.

Cannabinoid receptor 1 (CB1) is a GPCR that is highly expressed in the central nervous systems, and its endogenous ligands are a class of lipid mediators called endocannabinoids. It is also the major target of plant cannabis and synthetic cannabinoid drugs. CB1 has been extensively studied for ligand discovery, and its activation mechanism has also generated much attention. Interestingly, none of the three residues of the P-I-F motif is conserved in CB1. Instead, the residues in the three positions are L286^5.50, V204^3.40, and L352^6.44. The structures of CB1 in the inactive (Hua et al., 2016; Shao et al., 2016) and active states (Hua et al., 2017; Krishna Kumar et al., 2019; Hua et al., 2020) exhibited movements of the extracellular part that were not observed in the other class A GPCRs (Figure 1B); TM1 and TM2 move inward during activation, shrinking the space of the orthosteric site (binding site of endogenous ligands), and the N-terminus forms a V-shaped loop and inserts into the orthosteric site in the inactive state, but lies straight outside the helix bundle in the active state. Furthermore, a remarkable conformational change of F200^3.36 and W356^6.48 was observed (Figure 1C). F200^3.36 adopts different rotamers—in the inactive state, it points toward TM5, and in the active state, it points toward TM1. Upon activation, W356^6.48 takes the space left by F200^3.36, and correspondingly, the relative positions of TM3 and TM6 slide by 6.8 Å (measured by movement of W356^6.48 Cα atom when TM3 is aligned). While this previously described conformational change of the “twin toggle switch” residues F200^3.36/W356^6.48 has been proposed to be the key for CB1 activation (Hua et al., 2017), its mechanism of action has not been investigated.

Molecular dynamics (MD) simulations have been widely used in the study of activation mechanisms of GPCRs (Dror et al., 2011; Miao et al., 2013; Kohlhoff et al., 2014; Yuan et al., 2014; Lee et al., 2015; Weis and Kobilka, 2018; Lovera et al., 2019; Wifling et al., 2019; Mattedi et al., 2020). For example, simulations of β2 adrenergic receptor (β2AR) produced the deactivation process and confirmed the key role of the P-I-F motif as a connector (Dror et al., 2011). Here, we performed microsecond MD simulations of CB1 in different conditions (Table 1) to investigate its special activation mechanism. We found that wild-type CB1, unlike β2AR, could not deactivate in MD simulations (three independent runs each for 1 µs). But when we introduced double mutation F200A/W356A to CB1, MD simulations captured the deactivation process, validating the conformational change of the “twin toggle switch” residues F200^3.36/W356^6.48 as a barrier step in the activation mechanism. Simulations with different types of ligands (four full agonists, a partial agonist, a neutral antagonist, and two inverse agonists) showed that ligand efficacy is correlated with the movement of F200^3.36 and W356^6.48. In addition, we analyzed the structures of other class A GPCRs and discovered the 3.36/6.48 “twin toggle switch” is not unique to CB1.

MD simulations were performed with GROMACS 5.1.2 (Abraham et al., 2015) using AMBER FF14SB force field (Maier et al., 2015). The crystal structures of CB1 were used as starting points—active state, PDB entry 5XRA (Hua et al., 2017) (apo or agonists/neutral antagonist-bound, systems #1–8 in Table 1), and the inactive state, PDB entries 5TGZ (Hua et al., 2016) (AM6538-bound, system #9 in Table 1) or 5U09 (Shao et al., 2016) (taranabant-bound, system #10 in Table 1). Protein processing was performed with the Protein Preparation Wizard tool in Maestro software of the Schrodinger platform v2015-4 (Schrödinger, 2015). In this step, the missing side-chain atoms were added, and the N- and C-termini were capped as acetyl and N-methylamide, respectively. Mutating back the four mutations in the CB1 crystal construct, adding N-terminal residues F102 and M103 for the active state model, and modeling the truncated third intracellular loop were done with the Residue and Loop Mutation tool. For CB1 double-mutant F200A/W356A (system 2 in Table 1), the alanine residues were generated manually by removing the side-chain atoms. Models of active/inactive CB1 were embedded into a pre-equilibrated lipid bilayer containing 180 POPC (1-palmytoil-2-oleoyl-sn-glycerol-3-phosphatidylcholine) molecules and 12,894 water molecules [originally generated by CHARMM-GUI (Jo et al., 2008)] using the membed tool in the GROMACS program. Sodium ions were added at 0.15 M in water, and chloride ions were added to neutralize the system. The constructed systems contained 154–155 POPC molecules and ∼12,300 water molecules, making the total number of atoms ∼63,000. Three-dimensional structures of ligands (except for crystallized AM6538 and taranabant) were generated using the LigPrep tool in Maestro, and their binding poses in CB1 were predicted by molecular docking performed with Glide v6.9 (Friesner et al., 2004; Halgren et al., 2004) in extra precision (Friesner et al., 2006). The topology files of all the ligands and the lipid POPC were generated using AmberTools in UCSF Chimera (Pettersen et al., 2004) version 1.10.2 and were converted to GROMACS format with the ACPYPE tool (Sousa da Silva and Vranken, 2012).

For each system, energy minimization was performed with position restraints on non–hydrogen atoms of the protein and ligand—first using the steepest descent algorithm until the maximum force was smaller than 1,000 kJ mol⁻¹·nm⁻¹, then using the conjugate gradient algorithm until the maximum force was smaller than 100 kJ mol⁻¹·nm⁻¹ or converged to machine precision. The initial velocity of each atom was generated at 310 K (2–5 times with different seeds), and the system was relaxed in NVT (constant Number of atoms, Volume, and Temperature) ensemble (310 K, Berendsen coupling scheme) for 300 ps (time step 1 fs). The system was then equilibrated with position restraints to non–hydrogen atoms of protein and ligand in NPT (constant Number of atoms, Pressure, and Temperature) ensemble (310 K, 1 atm semi-isotropic, both Berendsen coupling scheme) for 15 ns in total—first with force constant 1,000 kJ mol⁻¹·nm⁻² for 5 ns (time step 1 fs), then gradually reducing the restraint force following 800, 600, 400, 200, and 0 kJ mol⁻¹·nm⁻², each for 2 ns (time step 2 fs). Finally, productive simulations were performed in the NPT ensemble (310 K, Nose–Hoover coupling scheme; 1 atm semi-isotropic, Parrinello–Rahman coupling scheme) for 1 μs. For all the simulations (relaxing, balancing, and productive), the leap-frog algorithm was used as the integrator, and constraints were set on the bonds with hydrogen atoms using the LINCS method; the cut-off distance for Van der Waals and short-range electrostatic interactions were set to 10 Å; and long-range electrostatic interactions were computed using smooth particle mesh Ewald method.

Snapshots were taken every 1 ns for analysis. The reference structures were—active state, cryo-EM structure of CB1/AM841/G protein (PDB entry 6KPG [Hua et al., 2020)] and inactive state, crystal structure of CB1/AM6538 [PDB entry 5TGZ (Hua et al., 2016)]. The overall root mean square deviations (RMSD) were calculated based on the Cα atoms of residues 106–306 and 335–411 (the third intracellular loop was neglected).

Structures of CB1 in inactive and active states showed uncommon conformational changes. Our MD simulations of deactivation process of CB1 captured conformational transitions at different regions. Based on these results, we propose a model of the special activation mechanism of CB1 (Figure 6).

There is an intermediate state of CB1 that is similar to the active structure at the extracellular part, but similar to the inactive structure at the cytoplasmic part. This is supported by our MD simulations starting from apo CB1 in the active state: for wild-type CB1, we observed movements of only the cytoplasmic part; for CB1 double-mutant F200A/W356A, transitions to the inactive state occurred earlier at the cytoplasmic part than the extracellular part. In the intermediate state, cytoplasmic ends of TM3 and TM6 were close, occupying the space for G protein binding. Eventually, such conformation of CB1 was already obtained in the experiment recently—in the crystal structure of CB1 bound with CP55,940 and ORG27569 (an allosteric modulator that enhances the binding of CP55,940 but reduces its signaling (Price et al., 2005; Shao et al., 2019).

This activation mechanism of CB1 is uncommon. For most class A GPCRs, the extracellular parts adopt similar conformations when bound with agonists or antagonists (Hua et al., 2017). Therefore, receptors with large-scale movement around the orthosteric pocket, including CB1 and P2Y purinoceptor 12 (Jin Zhang et al., 2014; Kaihua Zhang et al., 2014), may have special mechanisms for the coupling of the orthosteric pocket and the G protein–binding site. Moreover, a continuous water channel is an important hallmark for GPCR activation discovered by MD simulations (Yuan et al., 2013; Yuan et al., 2014; Yuan et al., 2015). However, we did not observe such a water channel in our simulations. This may be because the orthosteric pocket of CB1 is highly hydrophobic thus few water molecules could enter. We speculate this type of “dry activation” occurs in all the lipid receptors with highly hydrophobic orthosteric pockets similar to CB1.

The conformational change of the “twin toggle switch” residues F200^3.36/W356^6.48 occurs during transition from the intermediate state to active state. In our MD simulations of wild-type CB1, though the intermediate state was not achieved due to limited simulation time, we observed fluctuations of F200^3.36/W356^6.48 side-chains accompanied with movements of the cytoplasmic part, corresponding to the difference between the intermediate state and active state. Furthermore, in the structure of CB1 stabilized in the intermediate state by CP55,940 and ORG27569-bound (Shao et al., 2019), F200^3.36/W356^6.48 adopts the “inactive” conformation. For the process of this conformational change, we identified three possible transient configurations of F200^3.36/W356^6.48: “reversed-inactive,” “pushed,” and “sliding.” Among these configurations, “pushed” occurred most frequently and were observed in all the runs that captured fluctuations of F200^3.36/W356^6.48, thus might be a major configuration during the transition between the intermediate state and active state. For the other two configurations, run 2 of the apo and run 2 of the THC-bound CB1 mainly sampled “reversed-inactive” (Figure 2C and Supplementary Figure S4A), while run 1 of the THC-bound and run 3 of the O-2050-bound CB1 sampled more “sliding” than “reversed-inactive” (Figure 4C, Supplementary Figure S5B). So, we cannot reach any solid conclusions to discriminate partial agonist against neutral antagonist based on the limited MD simulations. The whole process of conformational changes of the “twin toggle switch” residues is beyond the time scale of this study and might be obtained with accelerated MD simulations or super long time conventional MD simulations. Also, free energy calculation in future investigations would help to obtain the energy barrier of the conformational change of “twin toggle switch” residues and decide the path of activation.

The role of conserved W6.48 in class A GPCR activation has long been debated. W6.48 is called the “toggle switch” because it was proposed to adopt different side-chain rotamers in the inactive and active states of β2AR by Monte Carlo simulations (Shi et al., 2002). But the crystal structure of the active state β2AR showed no change of W6.48 rotamer (Rasmussen et al., 2011a). In most structures of class A GPCRs, either inactive or active, W6.48 adopts the same rotamer as in β2AR, contradicting W6.48 as the “toggle switch” of common activation mechanism. However, there are several receptors showing W6.48 with different side-chain rotamers: 5-hydroxytryptamine receptor 2A (Kim et al., 2020) and μ-opioid receptor (Koehl et al., 2018) have unusual W6.48 in the active state; muscarinic acetylcholine receptors M1 (Thal et al., 2016), M2 (Suno et al., 2018), and CB2 (Li et al., 2019) have unusual W6.48 in the inactive state. These suggested there are diverse activation mechanisms of class A GPCRs at the connector region. “Twin toggle switch” of residues 3.36/6.48 in CB1 and MC4R provides a special mechanism and expands our understanding of the activation of GPCRs.