Mechanistic origin of partial agonism of tetrahydrocannabinol for cannabinoid receptors

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.

Cannabinoid receptor 1 (CB 1 ) is a therapeutically relevant drug target for controlling pain, obesity, and other central nervous system disorders. However, full agonists and antagonists of CB 1 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 CB 1 . 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 CB 1 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.
Cannabinoid receptor 1(CB 1 ) belongs to the family of class A G-protein-coupled receptors (GPCRs) (1,2), which modulates diverse cellular signaling processes via intracellular G-proteins (3) and β-arrestins (4). CB 1 receptors were first discovered in the last decade of the 20th century as a target of plant cannabinoid molecules (5). Due to its ubiquitous presence in physiological processes, CB 1 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 CB 1 for the treatment of pain, obesity, and inflammation (6)(7)(8). CB 1 receptor agonists such as MDMB-fubinaca (3) and inverse agonists such as rimonabant (9,10) have been shown to modulate the receptor activity significantly. However, these designed agonists and antagonists of CB 1 also exhibit dangerous side effects. For instance, fubinaca, also known as "zombie drug," caused the hospitalization of thousands of patients in New York (11). Rimonabant had to be withdrawn from the market due to its psychological side effects such as depression and anxiety (12). 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 (13). Simulation revealed important differences in binding site and pathway between CB 1 and CB 2 , 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 (14). Dronabinol is a synthetic form of tetrahydrocannabinol (THC) (Fig. 1), the main psychoactive compound in marijuana, which binds to the CB 1 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 (15). 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 CB 1 receptor.
Structural studies of CB 1 have revealed that toggle switch residue (TRP356 6.48 and PHE200 3.36 ) movement by the agonist is crucial for the activation of CB 1 (Fig. 2B) (3,(16)(17)(18)(19)(20)(21). Kumar et al. (3) 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 (3). 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 CB 1 (22) due to the large size of the pocket as compared with the ligand volume (18). 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 CB 1 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-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.
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) (23)(24)(25). 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 (18). 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
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 CB 1 (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 TRP356 6.48 . The new position of the TRP356 6.48 is different compared with active state of CB 1 . This intermediate position TRP356 6.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 CB 1 Active and inactive structures of CB 1 reveal that orthostericbinding site volume undergoes a large change upon activation as compared with other class A GPCRs (18). Comparison of active (PDB ID: 5XRA (17)) and inactive (PDB ID: 5TGZ (18)) 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 agonistbound crystal structure of CB 1 or CB 2 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 agonistbinding region in a similar conformation as other agonists (17,18). 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 THCbinding mechanism in CB 1 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 TYR275 5.39 -Cα atom on TM5. THC is predicted to enter CB 1 -binding pocket through the space between the space of TM1, TM2, and N-loop ( Fig. S2) (16). From MD simulations, THC is further predicted to diffuse inside the pocket and stabilized in antagonistbinding pose where distance between THC(C1') and TYR275 5.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 antagonistlike 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) (17). 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 (13,26,27) (28). Relative free energy calculation using alchemical analysis shows that mutations of PHE174 2.61 and PHE177 2.64 decrease the stability (ΔΔG PHE174 2:61 ALA ¼ 1:62± 0:28 kcal/mol and ΔΔG PHE177 2:64 ALA ¼ 3:45±0:56 kcal/mol) of the antagonist pose 1 suggesting that these mutations might decrease the THC binding by destabilizing the antagonistbound pose (29). 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 CB 1 ; 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 SER383 7.39 same as other agonists. Furthermore, THC forms extensive hydrophobic interactions with amino acid residues of TM3 (VAL196 3.32 , LEU193 3.29 ), TM5 (TRP279 5.43 ), TM7 (PHE379 7.35 ), N-loop (MET103 N−loop ), and ECL2 (PHE268 ECL2 ) as shown in Figure 4D and Fig. S3C.
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 agonistbinding 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 CB 1 . 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 (30). Furthermore, adding a dimethyl group in the first carbon of the alkyl chain is hypothesized to increase the interaction with toggle switch TRP356 6.48 (18). 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 agonistlike 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 CB 1 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 0 to 4 0 (Fig. 5A). In this conformation (macrostate 3 0 ), 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 0 and 3 (2.4 ± 0.5 μs) is more accessible compared with the 3 0 to 4 0 (13.2 ± 3.5 μs) (Figs. 5B and S6). The binding of THC to the agonistlike pose shows the orientation of alkyl chain may flip between macrostate 4 and 4 0 as they are energetically favorable. However, THC has more tendency to stay in the orthogonal conformation (macrostate 4) compared with macrostate 4 0 . 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.
THC is predicted to rotate the toggle switch TRP356 6.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 CB 1 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 TRP356 6.48 movement of the receptor toward the TM5, which consequently leads to outward movement of TM6 (18). It was hypothesized that THC behaves as partial agonist due to the lack of interaction with toggle switch residues (3). However, we observed that binding of THC leads to the rotation of TRP356 6.48 (Figs. 6A and S7). The dihedral angle of TRP356 6.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 TRP356 6.48 leads to breakage of aromatic interactions with PHE200 3.36 , and it moves toward TM2 similar to active structure (Fig. 6B). Comparison of representative structures from partially active state 1 with inactive CB 2 shows that toggle switch TRP356 6.48 /TRP258 6.48 has similar rotamaric conformation as CB 2 inactive structure (Fig. S8). An inverse agonist of CB 2 , MRI2687, which was predicted to retain the toggle switch residue to similar conformation (31), acts as a partial agonist for CB 1 . This in turn supports our prediction that THC stabilizes the TRP 6.48 in intermediate state to partially activate the receptor.
Although, THC is predicted to rotate toggle switch TRP356 6.48 for subset of structures in agonist-binding position, we noticed another stable state around inactive pose of TRP356 6.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 TRP356 6.48 , we calculated the probability TRP356 6.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 TRP356 6.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 TRP356 6.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, TRP356 6.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 TRP356 6.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.
Comparison of partial agonist and agonist-bound simulation ensemble using data-driven approach We have performed 17 μs MD simulations of agonist (AM11542) bound CB 1 (Experimental procedures section) and compared the result with the ensemble with the partial agonist bound to the agonist-like pose (32). 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 (18). 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. CB 1 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, TRP356 6.48 -ASP163 2.50 , LYS343 6.35 -ARG214 3.50 , and TYR153 2.40 -TYR397 7.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 activelike 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 (33)(34)(35). Using this network, we estimated the shortest pathway between the toggle switch and the NPxxY region using Dijkstra's algorithm (36). 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.

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 CB 1 . 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 (30). 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 (37): 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 (13). 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 CB 1 .
Our results also show that the THC can rotate the toggle switch TRP356 6.48 conformation when binding to the agonistlike poses. The partially active conformation of CB 1 matches well with CB2 inactive structure, which helps to explain yinyang relation for some ligands between CB 1 and CB 2 . Recently published NMR study on another class A GPCR, AA2AR, has shown that partial agonist rotates TRP246 6.48 in a distinct conformation compared with full agonists or Figure 7. Movement of intracellular TM6 due to TRP356 6.48 rotation. A, MSM weighted free energy landscape to capture toggle switch TRP356 6.48 (TM6) χ 2 angle and intracellular TM6 movement of CB 1 . Intracellular TM6 movement is measured ARG214 3.50 -Cα (TM3) and LYS343 6.35 -Cα (TM6). B, box plot to show the probabilities of intracellular TM6 movement with TRP356 6.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 (Experimental procedures 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.
antagonists (38). 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 (39)(40)(41)(42).
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 CB 1 . The findings of the paper are crucial for guiding the chemical modification for cannabinoid ligands for future drug development.

System preparation
Inactive structure of CB 1 (PDB ID: 5TGZ (17)) 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 (43). 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 (44). Thermostabilized mutant residues are replaced with original residues using tleap (17). The modified protein structure is embedded in POPC bilayer using CHARMM-GUI (45). 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 (46,47). THC is added to MD system using packmol (48) to generate the starting structure. To perform agonist (AM11542) bound simulation, active-like agonist-bound crystal structure (PDB ID: 5XRA (18)) 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 (49). 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 (50) 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 (51). Simulations were performed with Amber FF14SB (52) and GAFF (46) 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 (27,(53)(54)(55) including the ligand-binding process (13,(56)(57)(58)(59). 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 (60)(61)(62) 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 T ij of transition probability matrix calculates the probability of transition using the equation T ij ¼ C ij = P n j¼1 C ij where C ij is the count of jump between the i and j, and C i 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 (63) 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) (64,65) to find the slowest components. Tic components were shown to correlate well with important features for THC binding (Fig. S18, A-C). The transformed data were clustered using kmeans clustering algorithm. For MSM building, lag time was chosen by finding out logarithmic convergence of process timescales computed from MSM eigenvalues (process timescale t ¼ − τ log e λ ; λ is the eigenvalue) (Fig. S19). To find out the optimal number of clusters and optimal tic variance, MSM was subjected to VAMP2 (64) 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 (G i = −k b T log e π 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 (66, 67) was implemented. TPT calculates the flux between two macrostates (which consists of one or multiply MSM states) using MSM transition matrix (68). The flux (F AB ) is given by the equation where q þ j is probability that state j will reach macrostate B before A (69). From the flux calculation, we computed MFPT using equation MFPT = τπ A /F AB . where π A is the probability the system was in macrostate A. TPT calculation was performed with Pyemma (63) 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')-TYR275 5.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 (70) and MDtraj (71). For the trajectory visualization and analysis, VMD (72) 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 (73). 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 (PHE174 2.61 ALA and PHE177 2.64 ALA). 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 (29).

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 s ij ¼ P j k¼1 T ik . If R lies between s i;j and s i;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 TYR275 5.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 (71). 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 (74). 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 (TRP356 6.48 ) and NPxxY region (TYR397 7.53 ) using Dijkstra's algorithm implementation of the Networkx python package (75,76).

Docking study
Docking study was performed using Auto dock vina (77) 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 agonistbound 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 (78). 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. Funding and additional information-D. S. acknowledges support from NSF Early CAREER Award (MCB-1845606).
Conflict of interest-The authors declare that they have no conflicts of interest with the contents of this article.