Distinct activation mechanisms regulate subtype selectivity of Cannabinoid receptors.

Abstract

Design of cannabinergic subtype selective ligands is challenging because of high sequence and structural similarities of cannabinoid receptors (CB₁ and CB₂). We hypothesize that the subtype selectivity of designed selective ligands can be explained by the ligand binding to the conformationally distinct states between cannabinoid receptors. Analysis of ~ 700 μs of unbiased simulations using Markov state models and VAMPnets identifies the similarities and distinctions between the activation mechanism of both receptors. Structural and dynamic comparisons of metastable intermediate states allow us to observe the distinction in the binding pocket volume change during CB₁ and CB₂ activation. Docking analysis reveals that only a few of the intermediate metastable states of CB₁ show high affinity towards CB₂ selective agonists. In contrast, all the CB₂ metastable states show a similar affinity for these agonists. These results mechanistically explain the subtype selectivity of these agonists by deciphering the activation mechanism of cannabinoid receptors.

Untitled section

Subject terms: Computational biophysics, Membrane biophysics, Molecular medicine, Computational chemistry

Web Summary

MD simulations and Markov state modeling provide mechanistic insight into the subtype selectivity of agonists of cannabinoid receptors CB1 and CB2.

Article notes

Untitled section

Received 2022 Oct 18; Accepted 2023 Apr 24; Collection date 2023.

Author notes

Introduction

Cannabinoid receptors (CB₁ and CB₂) were first discovered as a target of phytocannabinoids in the last decades of the twentieth-century^1,2. Important physiological and psychological roles of these receptors were soon perceived by the elucidation of the endocannabinoid signaling system^3–5. CB₁ is majorly expressed in the central nervous system and has a role in appetite control, neuroprotection, and neurotransmission^6,7. CB₂ is majorly expressed in the immune system; targeted for immunomodulation and inflammation^6,8. These receptors belong to the lipid subfamily of class A GPCRs^9–17.

Class A GPCRs are recognized by conserved motifs which undergo structural changes during activation of the protein and help to transduce intracellular signaling by β-arrestins, and G-proteins^18–20. GPCR are targets of more than 30% of FDA-approved drugs^21–23. Similarly, academic labs and pharmaceutical companies have led drug discovery efforts to develop molecules targeting cannabinoid receptors^24–27. Earlier discovery campaigns resulted in molecules structurally similar to the known phytocannabinoids and endocannabinoids, which was followed by generation of chemically diverse molecules²⁸. Most of these ligands bind to both the receptors and lead to non-specific side effects^29–31; only a few designed ligands are selective to CB₁ or CB₂³². Specifically, CB₂ selective ligands can be used to explore the biology of CB₂, which has been described as “a cannabinoid receptor with an identity crisis"^33,34. Furthermore, CB₂ selective agonists hold a promise to be an important drug target for pain and inflammation without the psychoactive side effects caused by the binding to CB₁. However, an explanation of mechanistic underpinning for these ligands’ selectivity is still missing. Understanding subtype selectivity between CB₁ and CB₂ would be important for selective drug design^35–37.

Lack of subtype selectivity is also a major issue in drug design for other subfamilies of class A GPCRs. For example, the main targets of ergotamine are 5-HT_1B and 5-HT_1D receptors; however, its usage is very limited due to off-target actions on the other 5-HT subtype receptors, particularly on the 5-HT_2B receptor that causes life-threatening cardiovascular disorders³⁸. Structure-based approaches like docking studies or data-driven methods like machine learning have been employed to tackle the subtype selective drug design for class A GPCRs^39–43. Large-scale docking studies have been performed where ligands were docked into the static crystal structure or homology models to explain specific interactions crucial for subtype selectivity. For instance, important protein residue difference in TM5 of 5-HT_1B and 5-HT_2B was found to be an essential factor in drug design⁴⁴. However, docking does not take into account the large structural changes in the receptor binding pocket. Therefore, docking might not reveal the correct binding poses for ligands, especially when the binding pocket undergoes a significant volume change. Previous crystal structure studies show that selectivity between A₁-AR and A_2A-AR appears due to structural rather than sequence changes⁴⁵. On the contrary, machine learning models (e.g., LiCABEDS) have provided interpretability about the ligand scaffold that might lead to subtype selectivity. However, these ML models fail to consider receptor-ligand interactions, which leads to high false positive rate in their predictions^46,47.

In the last few years, crystal structures of inactive and active states of CB₁ and CB₂ were determined, which revealed ligand binding poses, and important structural changes during the activation of CB₁ and CB₂ (Supplementary Table 1). Comparing the antagonist-bound pose of inactive CB₁ and CB₂ shows the distinction between binding poses of antagonist of CB₁ and CB₂ (Supplementary Fig. 1a). However, structural comparison of agonist bound pose of unselective ligands of CB₁ and CB₂ shows that binding pocket residues have high structural and sequence similarity (Supplementary Fig. 1b, e). Docking performed using the experimentally determined structures did not reveal differences in the binding poses nor could determine subtype selectivity for CB₁ and CB₂^9,11. Furthermore, during the activation process, the CB₁ binding pocket undergoes large conformational changes due to the upward movement of the N-terminus, whereas CB₂ pocket retains its overall shape (Supplementary Fig. 1c, d). Considering these two factors, we hypothesized that CB₂ selective agonists conformationally select for the CB₂ pocket shape. As the pocket shape remains similar through out the activation process, these ligands have the higher possibility of binding to the CB₂ rather than CB₁, where the pocket shape and volume changes significantly during the activation thereby decreasing the overall binding affinity (Fig. 1a, b).

Cartoon representation of the hypothesis describing the subtype selectivity of the CB2 selective agonists.a According to the proposed hypothesis, binding pocket volume of CB1 (Color: Orange) changes during the activation process. CB2 selective agonists (Color : Blue) may bind into specific intermediate states of CB1 decreasing overall binding affinity. b Contrastingly, these ligands prefer binding pocket shape of CB2 (Color: Green), which retains this generalized shape during activation process. — Fig. 1Cartoon representation of the hypothesis describing the subtype selectivity of the CB₂ selective agonists.a According to the proposed hypothesis, binding pocket volume of CB₁ (Color: Orange) changes during the activation process. CB₂ selective agonists (Color : Blue) may bind into specific intermediate states of CB₁ decreasing overall binding affinity. b Contrastingly, these ligands prefer binding pocket shape of CB₂ (Color: Green), which retains this generalized shape during activation process.

Here, we compared the activation mechanism for CB₁ and CB₂ using molecular dynamics simulations to obtain their entire conformational ensemble. Using 700 μs of aggregate unbiased simulation data, we deciphered the similarity and distinction of activation mechanism between CB₁ and CB₂. Analyzing the data using markov state model, we found important allosteric communications between different structural motifs, which facilitates the receptor activation. By implementing deep learning based VAMPnets architecture, we discretize the conformational space into different metastable states⁴⁸. Docking of the CB₂ selective ligands shows the docking affinities of these ligands are similar for all the metastable states in CB₂ whereas these ligands bind to only specific metastable states in CB₁. Binding of CB₂ selective ligands favorably to only specific metastable states might explain the subtype selectivity of the ligands. This mechanistic understanding of ligand selectivity between CB₁ and CB₂ will aid in the design of new subtype-selective drugs for cannabinoid receptors.

Results and discussions

Structural comparison of CB₁ and CB₂ crystal structures

In the last few years, a plethora of CB₁ and CB₂ structures are reported in different conditions as shown in Supplementary Table 1. A comparison of these crystal structures shows major structural changes in the receptors’ extracellular, transmembrane, and intracellular regions during the activation. As both these receptors belong to the lipid subfamily of class A GPCRs, some of these changes are conserved between CB₁, and CB₂^9,11,12,16. For example, both the receptors undergo large conformational changes in intracellular TM6 towards the outward direction, as shown in Fig. 2. Projection of intracellular TM6 movement (R^3.50-K^6.35 distance) feature on one dimensional space shows that this feature can distinguish active and inactive structures for both CB₁ and CB₂ (Fig. 3 and Supplementary Fig. 2 and Supplementary Note 1). Similarly, the sidechain of conserved Y^7.53 in the NPxxY region of intracellular TM7 moves towards TM5 during the full activation of the receptors (Fig. 2). Projection of intracellular TM7 distance (I^5.54-Y^7.53 distance) also shows that this feature changes during activation of CB₁ and CB₂ (Fig. 3). Therefore, these two conserved intracellular features are used as metrices to judge whether a structure is active or inactive. However, other conformational changes are not conserved between the CB₁ and CB₂. Structural overlap of a representative active and inactive CB₁ and CB₂ structures show the following non-conserved changes in CB₁ and CB₂.

Comparison of inactive and active structures of CB1 and CB2 highlighting the important motions during activation.Inactive crystal structure of CB1 (a) (PDB ID: 5TGZ9, color: Orange) and CB2 (b) (PDB ID: 5TZY12, color: Green) are shown as cartoon representations in the central panel. Protein residues in conserved motifs are shown as sticks and highlighted inside block boxes. Major conformational changes during the activation are shown as separated box, where the active (CB1 PDB ID: 5XRA11, color: Blue; CB2 PDB ID: 6KPF16, color: Pink) and inactive conformations are superposed. Directions of the conformational changes during the activation are shown as black arrow. — Fig. 2Comparison of inactive and active structures of CB₁ and CB₂ highlighting the important motions during activation.Inactive crystal structure of CB₁ (a) (PDB ID: 5TGZ⁹, color: Orange) and CB₂ (b) (PDB ID: 5TZY¹², color: Green) are shown as cartoon representations in the central panel. Protein residues in conserved motifs are shown as sticks and highlighted inside block boxes. Major conformational changes during the activation are shown as separated box, where the active (CB₁ PDB ID: 5XRA¹¹, color: Blue; CB₂ PDB ID: 6KPF¹⁶, color: Pink) and inactive conformations are superposed. Directions of the conformational changes during the activation are shown as black arrow.

Conversed and unconserved feature movements captured in active and inactive PDB structures for CB1 and CB2.Distance and angle feature values that undergo conformational changes during the activation for CB1 (a) and CB2 (b) are shown as scatter points on 1-D line. Markers are colored based on the activation state of the protein (CB1 active: blue, CB1 inactive: orange, CB2 active: magenta, CB2 inactive: green). The shape of the marker is based on the type of the ligand and downstream signaling partner (Antagonist bound: circle, Agonist bound: triangle, Agonist and NAM bound: star, Agonist and G protein bound: rectangle). PDB structures of a receptor with same binding partner and activation state are distinguished with different outline color of the marker. — Fig. 3Conversed and unconserved feature movements captured in active and inactive PDB structures for CB₁ and CB₂.Distance and angle feature values that undergo conformational changes during the activation for CB₁ (a) and CB₂ (b) are shown as scatter points on 1-D line. Markers are colored based on the activation state of the protein (CB₁ active: blue, CB₁ inactive: orange, CB₂ active: magenta, CB₂ inactive: green). The shape of the marker is based on the type of the ligand and downstream signaling partner (Antagonist bound: circle, Agonist bound: triangle, Agonist and NAM bound: star, Agonist and G protein bound: rectangle). PDB structures of a receptor with same binding partner and activation state are distinguished with different outline color of the marker.

(1) For CB₁, in the extracellular region, N-terminus shows large conformational differences in agonist bound structure as compared to the antagonist bound structures (Fig. 2a). Agonist and antagonist bound poses are also different for CB₁^11,13. Agonist binding lifts the N-terminus above the orthosteric pocket. The N-terminus distance (M^N-term-D^2.50 distance) from the ligand binding pocket clearly shows the difference in the agonist bound structures as compared to the antagonist bound structures (Fig. 3a). As agonist binding leads to the activation of CB₁, active and inactive structures have distinct N-terminus feature values. The exception to this rule is the agonist and Negative Allosteric Modulator (NAM) bound structure (PDB ID: 6KQI¹⁴), where the structure remains inactive because of NAM binding. On the other hand, in the active and inactive structure, N-terminus always remain above the orthosteric pocket; therefore, N-terminus feature values do not distinguish between active and inactive crystal structures for CB₂ (Figs. 2b, 3b).

(2) In CB₁, due to the extended structure of the antagonist, TM1 moves outwards towards the membrane in the inactive structures (Fig. 2a), whereas extracellular TM1 conformation of CB₂ stays relatively similar in the active and inactive state in the protein (Fig. 2b). The higher variation in extracellular TM1 movement between active and inactive structures in CB₁ as compared to CB₂ is shown in (Fig. 3a, b).

(3) Another major difference in activation of CB₁ and CB₂ is observed in toggle switch residue movement. For CB₁, twin toggle switch residues (W^6.48 and F^3.36) undergo translational changes (Fig. 2a). During the activation, W^6.48 moves towards the TM5, which leads to the displacement of F^3.36 towards TM2. The displacement is captured by the difference in the center of mass of the sidechain of these two residues in the z-direction (Supplementary Fig. 2). In the case of CB₂, F^3.36 does not undergo conformational change (Fig. 2b). Instead, agonist binding to CB₂ has shown to rotate the toggle switch residue W^6.48¹⁶. This rotational motion is captured by calculating the χ₂ angle of W^6.48. Projection of translational and rotational movement features of toggle switch show that former feature could capture the distinction between active and inactive structures for CB₁ and latter for CB₂ (Fig. 3a, b).

(4) For CB₂, inactive structures in intracellular TM5, the sidechain of Y^5.58 faces towards the membrane and creates more space for TM6 to move inside in the inactive structure. In the active structure, Y^5.58 side chain (I^2.43-Y^5.58 distance) moves inside the transmembrane domain and faces towards the TM2 (Figs. 2b, 3b). No such conformational change is observed for CB₁. Y^5.58 always remains inside the transmembrane domain (Figs. 2a, 3a).

Thermodynamic relevance of structural features in cannabinoid receptors’ activation

All extracellular (N-terminus and TM1 movements), transmembrane (Toggle Switch translation and rotation), and intracellular features (TM5, TM6 and TM7 movements) that distinguish active and inactive states of the receptors identified from different crystal structures may not be thermodynamically and kinetically relevant (Fig. 3). Experimentally determined side chain conformations occasionally do not have good resolutions. Hence, it may over or under emphasize some of the conformational changes. Furthermore, features that do not distinguish between active and inactive state for respective proteins (e.g., Toggle switch rotation for CB₁; TM1 movement for CB₂) may obtain distinct stabilized values in intermediate states during activation. Therefore, to check the thermodynamic relevance of each feature in CB₁ and CB₂ conformational ensemble, we compared the apo (without ligand bound) and holo (with ligand bound) simulation of both proteins.

Compairing ~ 20 μs of unbiased MD simulation for agonist bound active and antagonist bound inactive states of both CB₁ and CB₂ (Supplementary Table 1), it is observed that the conformational ensemble of every feature remains stable (Supplementary Fig. 3). It reinforces that the features observed in the crystal structures are thermodynamically relevant. The features that distinguish between CB₁ and CB₂ activation in experimentally determined structures are stabilized with different feature values.

The presence of intermediate states was judged by comparing the conformational ensemble of important features obtained from apo simulations with the ligand-bound active and inactive holo simulation (Supplementary Fig. 3). Simulation details are shown in Supplementary Tables 2 and 3. Three metrices are used to judge the feature importance : (1) sum of the absolute value of the mean difference between apo and holo simulations ∣E(Holo_inactive) − E(apo)∣ + ∣E(Holo_active) − E(apo)∣, where E(Feature) represents the mean of the feature from simulation. (2) Kullback-Leibler divergence (K-L divergence) of apo simulation with respect to agonist bound active simulation (3) K-L divergence of apo simulations with respect to antagonist bound inctive simulations⁴⁹. Here, features were normalized for better comparison using these metrices. As expected, all the features which show distinction in crystal structures for both proteins have larger values as shown in Tables 1 and 2. This indicates that these features undergo large changes during activation. On the other hand, features that were found to be only important for CB₁ or CB₂ show higher values for the metrices for only the respective protein. The only exception to this rule is the extracellular TM1 movement for CB₂, where ligand-bound active and inactive conformational ensembles show overlapping distribution, whereas apo activation conformational changes much wider distribution of feature values (Supplementary Fig. 3), which may reveal the newer conformations. Therefore, using our simulation, we observe changes in the conformational features of the proteins that are not deemed important from static crystal structures. Overall, for CB₁, features that undergoes conformational changes during activations are N-terminus distance, extracellular TM1 movement, toggle switch translational movement, intracellular TM6 and TM7 movements. Contrastingly for CB₂, conformational changes are observed in extracellular TM1 movement, toggle switch rotational movement, intracellular TM5, TM6 and TM7 movements.

Table 1Comparisons of conformational feature values between apo (without ligand) and holo (with agonist and antagonist) simulations of CB₁ using mean difference and K-L divergence.

Features	∣E(Holo_inactive) − E(apo)∣ + ∣E(Holo_active) − E(apo)∣	K-L divergence (Holo_active and apo)	K-L divergence (Holo_inactive and apo)
N-terminus Movement	0.4	2.1	1.3
Extracellular TM1 movement	0.2	1.9	1.1
Toggle Switch translational movement	0.3	1.3	1.8
Toggle Switch rotational movement	0.1	0.3	0.4
Intracellular TM5 movement	0.0	0.7	0.9
Intracellular TM6 movement	0.3	1.4	2.2
Intracellular TM7 movement	0.4	0.3	1.7

Table 2Comparisons of conformational feature values between apo (without ligand) and holo (with agonist and antagonist) simulations of CB₂ using mean difference and K-L divergence.

Features	∣E(Holo_inactive) − E(apo)∣ + ∣E(Holo_active) − E(apo)∣	K-L divergence (Holo_active and apo)	K-L divergence (Holo_inactive and apo)
N-terminus Movement	0.1	2.1	0.7
Extracellular TM1 movement	0.3	0.9	1.1
Toggle Switch translational movement	0.1	1.5	1.7
Toggle Switch rotational movement	0.2	0.8	2.4
Intracellular TM5 movement	0.5	0.4	1.5
Intracellular TM6 movement	0.2	1.0	1.7
Intracellular TM7 movement	0.5	0.1	2.9

Allosteric communications between CB₁ and CB₂ structural features

In the active state, conformational changes in the intracellular part of GPCR proteins help to transduce signals by secondary messenger (G protein and β-arrestin)^20,50. These intracellular changes are allosterically linked with extracellular and transmembrane regions of protein¹⁹. To observe how the conserved and non-conserved changes in CB₁ and CB₂ are allosterically connected, we predicted whether the change in one feature leads to the change in the other features. To estimate this, the absolute difference of P(Feature_active∣Condition_active) and P(Feature_active∣Condition_inactive) was determined for every combination of features. These conditional probabilities were measured by discretizing the features based on a threshold that divides the features into active and inactive. The threshold values are taken as an average of active and inactive crystal structures. For the conditional probability calculations, weight of each frame was updated using equilibrium probability distribution estimated using MSM with lag time 25 ns (More details in Method section).

Based on the calculations for CB₁, it shows that N-terminus movement is highly correlated with all the extracellular (TM1 movement), transmembrane (Toggle switch translational movement), and intracellular features (TM6 and TM7 movement), showing the importance of the N-terminus movement in activation (Fig. 4a and Supplementary Fig. 4a). Projection of MSM weighted free energy landscapes show that with the N-terminus inside the pocket, other features remain mostly in the stabilized inactive state (Supplementary Fig. 5a–d). On the other hand, extracellular TM1 movement is not strongly coupled with transmembrane and intracellular features (Fig. 4a and Supplementary Fig. 4b). Free energy projection of TM1 movement with respect to intracellular and transmembrane features shows those features are able to obtain active and inactive forms without depending on the TM1 value (Supplementary Fig. 5e–g). Similarly, toggle switch translational movement also shows relatively smaller values of coupling for intracellular TM6 movements and TM7 movements (Fig. 4a and Supplementary Fig. 4c). However, the MSM-weighted free energy landscape can explain this relatively smaller value. The free energy plots show an L-shaped landscape between the toggle switch and intracellular features depicting the change to be sequential (Supplementary Fig. 5h–i). Thus, changes in the toggle switch can lead to changes in the intracellular features. Lastly, both conditional probability difference and free energy landscape show the intracellular features are highly coupled with each other, which depicts that the TM6 movement will lead to the TM7 movement (Fig. 4a, Supplementary Figs. 4d, e and 5j). Based on the calculations for CB₁, it shows that N-terminus movement is highly correlated with all the extracellular (TM1 movement), transmembrane (Toggle switch translational movement), and intracellular features (TM6 and TM7 movement), showing the importance of the N-terminus movement in activation (Fig. 4a and Supplementary Fig. 4a). Projection of MSM weighted free energy landscapes show that with the N-terminus inside the pocket, other features remain mostly in the stabilized inactive state (Supplementary Fig. 5a–d). On the other hand, extracellular TM1 movement is not strongly coupled with transmembrane and intracellular features (Fig. 4a and Supplementary Fig. 4b). Free energy projection of TM1 movement with respect to intracellular and transmembrane features shows those features are able to obtain active and inactive forms without depending on the TM1 value (Supplementary Fig. 5e–g). Similarly, toggle switch translational movement also shows relatively smaller values of coupling for intracellular TM6 movements and TM7 movements (Fig. 4a and Supplementary Fig. 4c). However, the MSM-weighted free energy landscape can explain this relatively smaller value. The free energy plots show an L-shaped landscape between the toggle switch and intracellular features depicting the change to be sequential (Supplementary Fig. 5h, i). Thus, changes in the toggle switch can lead to changes in the intracellular features. Lastly, both conditional probability difference and free energy landscape show the intracellular features are highly coupled with each other, which depicts that the TM6 movement will lead to the TM7 movement (Fig. 4a, Supplementary Figs. 4d, e, and 5j).

Allosteric dependances of the important structural feature movements during the activation.Absolute differences of conditional probabilities (P(FAi∣FAj) and P(FAi∣FIj)) of structural features for CB1 (a) and CB2 (b). FAi denotes the ith feature is in active state. FIi denotes the ith feature is in inactive state. Errors in the conditional probability calculations were calculated based on bootstrapping. 200 bootstrap samples were selected with 80% of the total number of trajectories. Bar plots and error bars represent mean and standard deviation of the 200 samples, respectively. Each sample is shown on the plot as black dot. — Fig. 4Allosteric dependances of the important structural feature movements during the activation.Absolute differences of conditional probabilities ( $P (F_{A}^{i} ∣ F_{A}^{j})$ and $P (F_{A}^{i} ∣ F_{I}^{j})$ ) of structural features for CB₁ (a) and CB₂ (b). $F_{A}^{i}$ denotes the ith feature is in active state. $F_{I}^{i}$ denotes the ith feature is in inactive state. Errors in the conditional probability calculations were calculated based on bootstrapping. 200 bootstrap samples were selected with 80% of the total number of trajectories. Bar plots and error bars represent mean and standard deviation of the 200 samples, respectively. Each sample is shown on the plot as black dot.

For the CB₂, all the intracellular features are highly coupled with the extracellular TM1 movement (Fig. 4b and Supplementary Fig. 6a). MSM weighted free energy landscape shows that when intracellular features remains in the inactive states, extracellular TM1 can have a wider range of movement. Whereas, in active state, extracellular TM1 movement is restricted (Supplementary Fig. 7a–d). From the experimentally obtained structures, it has been shown that toggle switch residue rotates when the agonist molecule binds to the receptor. Therefore, it has been hypothesized that toggle switch rotation caused by agonists leads to intracellular changes in the CB₂ receptor^15,16. However, our calculations show that toggle switch rotation has minimal effect on the intracellular movement (Fig. 4b and Supplementary Fig. 6b). From the free energy projection, it is also clear that intracellular features can adopt both active and inactive forms irrespective of toggle switch rotation feature value (Supplementary Fig. 7e–g). Agonist-bound CB₂ structure (PDB ID: 6KPC) remains in an inactive state as there is no intracellular movement observed, which supports our finding that there is a lack of correlation between toggle switch rotation and intracellular movement¹⁶. For CB₂, intracellular TM5 movement is shown to be important for canonical GPCR activation. The free energy landscape shows that L shaped landscape between intracellular TM5 movement with respect to the intracellular TM6 and TM7 movement, which depicts that TM5 movements may be the cause of the activation of the GPCR, not the toggle switch rotation (Supplementary Fig. 7h, i). Lastly, similar to CB₁, intracellular TM6 and TM7 movement are highly correlated with each other (Fig. 4b and Supplementary Fig. 7j). Therefore, these analyses show that most of these features (Except for TM1 movement for CB₁; toggle switch rotation in CB₂), which are physically far away in protein topology, are allosterically influence conformational changes of other features resulting in protein activation.

Activation mechanism for CB₁ and CB₂

Previous sections show the important conformational changes for CB₁ and CB₂ during the activation. However, individual feature movements or projections of a few features do not reveal the activation mechanism of the receptors and intermediate states involved in the process. To better understand the activation mechanism, we clustered the protein conformational ensemble using neural network architecture, VAMPnets⁴⁸. VAMPnets architecture discretizes the protein conformational space into metastable states by building a coarse-grain Markov state model (MSM). The number of metastable states is decided by minimizing the error in implied timescale with a condition that each metastable state will have minimum population of 4% of entire ensemble (Supplementary Figs. 8, 9, 10, and 11). Using these criteria, we built a six-state model for both CB₁ and CB₂. These six states are projected on the tic1 and tic2 dimensions, which are linear combinations of the slowest features (Figs. 5 and 6). Here, two of the metastable states are presented as active and inactive based on the projection of the active and inactive experimentally determined structures onto tics. Except for the active and inactive states, there are four intermediate metastable states obtained from the model, named I1, I2, I3, and I4 for both CB₁ and CB₂ (Figs. 5 and 6). For CB₁, I1 and I2 are structurally similar to the inactive state of the receptor (Supplementary Table 4 and Supplementary Fig. 12). Whereas, I3 and I4 are structurally similar to active state (Supplementary Table 4 and Supplementary Fig. 12). The structural differences between different metastable states are calculated using K-L divergence based on the closest heavy atom distances between all residue pairs with respect to the inactive state. In Figs. 5 and 6, the calculated K-L divergence is shown as a color bar on the representative protein structure, where the red represents structural divergence. Furthermore, normalized values of the important features (that have been discussed in the earlier section) are calculated for each metastable state as shown in Figs. 5 and 6. Transition kinetics between the metastable states are shown as mean free passage time (MFPT) if there is a direct transition. Using the above analyses, it is observed that from the CB₁ inactive state, both I1 and I2 transitions are possible. I1 transitions are kinetically more favorable compared to I2. Inactive to I1 transition involves the movement of intracellular TM1 region as shown in K-L divergence analysis and bar plot (Fig. 5). In the case of inactive to I2 and I1 to I2, the transition involves the movement of the N-terminus in the upward direction. Movement of the N-terminus is a slow timescale process as these movements are related to the first eigenvectors of MSM (Supplementary Fig. 13). Intracellular and transmembrane regions of I1 and I2 metastable states remain similar to inactive states. Conformations of the I2 metastable state are similar to the agonist, and NAM-bound crystal structure (PDB ID: 6KQI¹⁴), where the structure remains in the inactive state in spite of the binding of the agonist in the pocket (Supplementary Fig. 14). I1 to I3 transition involves rearrangement of the intracellular and transmembrane features as N-terminus starts to move out of the pocket. K-L divergence analysis and bar plot show that the major difference occurs in N-terminus, toggle switch region, intracellular TM6, and TM7 region (Fig. 5 and Supplementary Fig. 12). This step is the the kinetically slowest process during the activation. Once the intracellular and transmembrane rearrangement happens in I3, the rearrangement of structural rearrangement happens between I3 & I4 and I3 & active state of the receptor. In the case of the I3 to I4 transition, N-terminus moves slightly further upward the pocket, and further intracellular features movement happens (Fig. 5). During I3 to active state transition, the N-terminus moves out of the binding pocket, which TM1 to come inside the pocket; however, intracellular features remain relatively the same (Fig. 5). Between these two transitions, the later transition is kinetically faster as I3 to I4 transition involves N-terminus movement (Fig. 5).

Mean first passage time and important structural changes in CB1 activation ensemble.Representative structures from six CB1 metastable states obtained from VAMPnets are projected on tic 1 and tic 2 spaces. Cartoon representations of receptors are colored based on K-L divergence. Red color is representing higher K-L divergence. K-L divergence of each metastable state was calculated compared to inactive metastable state. Transition between metastable states are shown as arrows. Thickness of the arrows represents mean free passage time (MFPT) of the transition which was calculated using transition path theory (TPT). Red dashed line is representing the slowest transition. The unit of MFPT is in microseconds. Bar plots per metastable state shows the normalized values of important structural features, where [1], [2], [3], [4], [5] represent N-terminus movement, extracellular TM1 movement, toggle switch translational movement, intracellular TM6 movement, intracellular TM7 movement respectively. Errors in bar plots were calculated based on bootstrapping. 20 bootstrap samples were prepared by selecting 1000 frames from each metastable states. Bar plots and error bars represent mean and standard deviation of the 20 samples, respectively. Each sample is shown on the plot as black dot. — Fig. 5Mean first passage time and important structural changes in CB₁ activation ensemble.Representative structures from six CB₁ metastable states obtained from VAMPnets are projected on tic 1 and tic 2 spaces. Cartoon representations of receptors are colored based on K-L divergence. Red color is representing higher K-L divergence. K-L divergence of each metastable state was calculated compared to inactive metastable state. Transition between metastable states are shown as arrows. Thickness of the arrows represents mean free passage time (MFPT) of the transition which was calculated using transition path theory (TPT). Red dashed line is representing the slowest transition. The unit of MFPT is in microseconds. Bar plots per metastable state shows the normalized values of important structural features, where [1], [2], [3], [4], [5] represent N-terminus movement, extracellular TM1 movement, toggle switch translational movement, intracellular TM6 movement, intracellular TM7 movement respectively. Errors in bar plots were calculated based on bootstrapping. 20 bootstrap samples were prepared by selecting 1000 frames from each metastable states. Bar plots and error bars represent mean and standard deviation of the 20 samples, respectively. Each sample is shown on the plot as black dot.

Mean first passage time and important structural changes in CB2 activation ensemble.Representative structures from six CB2 metastable states obtained from VAMPnets are projected on tic 1 and tic 2 spaces. Cartoon representations of receptors are colored based on K-L divergence. Red color is representing higher K-L divergence. K-L divergence of each metastable state was calculated compared to inactive metastable state. Transition between metastable states are shown as arrows. Thickness of the arrows represents mean free passage time (MFPT) of the transition which was calculated using transition path theory (TPT). Red dashed line is representing the slowest transition. The unit of MFPT is in microseconds. Bar plots per metastable state shows the normalized values of important structural features, where [1], [2], [3], [4], [5] represent extracellular TM1 movement, toggle switch rotational movement, intracellular TM5 movement, intracellular TM6 movement, intracellular TM7 movement respectively. Errors in bar plots were calculated based on bootstrapping. 20 bootstrap samples were prepared by selecting 1000 frames from each metastable states. Bar plots and error bars represent mean and standard deviation of the 20 samples, respectively. Each sample is shown on the plot as black dot. — Fig. 6Mean first passage time and important structural changes in CB₂ activation ensemble.Representative structures from six CB₂ metastable states obtained from VAMPnets are projected on tic 1 and tic 2 spaces. Cartoon representations of receptors are colored based on K-L divergence. Red color is representing higher K-L divergence. K-L divergence of each metastable state was calculated compared to inactive metastable state. Transition between metastable states are shown as arrows. Thickness of the arrows represents mean free passage time (MFPT) of the transition which was calculated using transition path theory (TPT). Red dashed line is representing the slowest transition. The unit of MFPT is in microseconds. Bar plots per metastable state shows the normalized values of important structural features, where [1], [2], [3], [4], [5] represent extracellular TM1 movement, toggle switch rotational movement, intracellular TM5 movement, intracellular TM6 movement, intracellular TM7 movement respectively. Errors in bar plots were calculated based on bootstrapping. 20 bootstrap samples were prepared by selecting 1000 frames from each metastable states. Bar plots and error bars represent mean and standard deviation of the 20 samples, respectively. Each sample is shown on the plot as black dot.

For CB₂, intermediate states I1, I2, and I3 are structurally and kinetically similar to the inactive state of the receptor, whereas I4 is structurally similar to the active state (Supplementary Table 5; Fig. 6 and Supplementary Fig. 15). K-L divergence shows no major structural changes between the intermediate I1 and I2 & inactive metastable states. Therefore, these states can interconvert with a faster timescale between each other (Fig. 6). The inactive and I3 states are separated in tic one space. Intracellular TM6 and TM7 features are correlated in the tic 1 space (Supplementary Fig. 16). In this transition, TM6 moves out slightly, whereas Y^7.53 moves towards the TM5 (Supplementary Table 5 and Supplementary Fig. 15). From I3 to I4 and active state movement, this TM6 movement leads the Y^5.58 to move inside the transmembrane region. Both of these transitions are kinetically slower as it involves rearrangement of the intracellular regions of the CB₂.

MFPT calculations based on transition path theory (TPT) clearly show that the transition between inactive to active state in CB₂ is much faster than CB₁. A Kinetic monte Carlo (kMC) study on MSM based transition matrices also verifies that inactive to active transition is faster for CB₂ transition (Supplementary Fig. 17a, b). kMC analysis also indicates that transition between inactive, I1, and I2 is kinetically much faster for CB₁, compared to I1 to I3 transition. I1 to I3 transition leads to the protein moving from the inactive to the active region. Similarly, the kMC study verifies that inactive, I1, I2, and I3 are kinetically closer to each other for CB₂. Therefore, the transitions between the states are relatively faster. Once the transition happens from I3 to active or active like I4, it remains stabilized in active-like proteins. The above analyses reveal intermediate states during CB₁ and CB₂ activation, by which the sequential conformational change happens that leads the protein from the inactive to the active state, and the distinction between changes between CB₁ and CB₂.

We also performed long-timescale MD simulations from the representative structure of every metastable state of CB₁ and CB₂ (12 systems in total) to observe the stability of each state. Simulation length for each system was approximately 600 ns (Supplementary Note 2). RMSD and RMSF analysis (Supplementary Figs. 18 and 20) on simulation data shows that RMSD of each metastable state remains within 4 Å for both CB₁ and CB₂, where RMSD for CB₁ metastable states are lower than the CB₂. This can be explained by the timescale of the interstate transition between CB₁ and CB₂. Transition path theory (TPT) (Figs. 4, 5) and kinetic monte carlo (kMC) simulation (Supplementary Fig. 17) on the MSM weighted simulation show that CB₁ interstate transitions are much slower compared to CB₂. Therefore, CB₂ structures can fluctuate more rapidly compared to CB₁, resulting higher RMSD. We observe the same phenomena with time independent component (tIC) projection of long timescale MD simulations. The tIC projections were calculated from the distances between the residue pairs calculated based on residue-residue contact score (RRCS) analysis (Method section). The tIC projections for CB₂ show that the simulations started in different metastable states are overlapping with each other within MD timescale (Supplementary Fig. 19b), whereas simulations of different metastable states of CB₁ remain separated due to larger kinetic barriers (Supplementary Fig. 19a).

We further compared how the conformational changes in the distinct metastable states lead to different allosteric communication networks. Allosteric pipelines in different metastable states were obtained following the analyses described in Bhattacharya and Vaidehi⁵¹. From these analyses, major distinctions in the allosteric pipelines were observed for CB₁ in the inactive and inactive intermediate states (I1 and I2) compared to active and active intermediate states (I3 and I4) (Supplementary Fig. 21). In inactive, I1, and I2 states, extracellular TM1, and TM2 are allosterically connected to ECL2 via N-terminus. The downward movement of the N-terminus can be explained to be responsible for these allosteric pipelines as there is more interaction between N-terminus and ECL2. There are fewer pipelines connecting extracellular and intracellular domains. In contrast, major pipelines are observed connecting these two domains in active and active intermediate states. Especially all these three states have an allosteric pipeline through TM6, which indicates the outward movement of TM6 helps to create an allosteric path. On the other hand, for CB₂, no such clear distinction is observed in allosteric pipelines between active like (active, I4) and inactive like states (Inactive, I1, I2, and I3). Most of the metastable states contains extracellular to intracellular communications (Supplementary Fig. 22). This contrasting network pattern in CB₂ compared to CB₁ can be explained by the lack of the major extracellular conformational changes in CB₂. Hence, these analyses show how conformational changes and the consequent change in allosteric communications facilitate receptor activation.

Explaining selectivity of Cannabinoid receptors

Ligand binding to a receptor can be explained using two distinct mechanisms, i.e., conformational selection and induced fit⁵². In the induced fit mechanism, ligand binding leads to conformational changes in the receptor. Without the ligand, it is kinetically and thermodynamically hard to reach the bound conformation of the receptor. On the other hand, in conformational selection mechanism, the receptor can adopt many possible conformations. However, the ligand preferentially binds to a particular protein conformational state.

As active and inactive binding pocket volumes of experimentally obtained CB₂ structures are relatively similar, contrasting to the dissimilarity between CB₁ structures, we hypothesized that the CB₂ pocket shape remains similar throughout the activation process. If so, then the binding of CB₂ selective agonists would demonstrate a conformational selection mechanism in preference to this generalized binding pocket conformation seen during CB₂ activation. On the contrary, the CB₁ binding pocket shape changes with the movement of the N-terminus during the activation. Therefore, these selective ligands can only bind to certain metastable states for CB₁ which might have similar pocket characteristics as of CB₂, decreasing the overall affinity. To test this hypothesis, we first calculated the binding pocket volumes of the metastable states during CB₁ and CB₂ activation. Our analyses show that CB₁ metastable states have different binding pocket volumes based on the position of the N-terminus (Fig. 7a). As the N-terminus remains completely inside the pocket during inactive and I1 states, agonist binding volume is comparatively smaller to the other four metastable states. In the case of CB₂, the binding pocket volume does not change during the activation due to a lack of N-terminus motion (Fig. 7b). Here our results illustrate how CB₂ maintains a competent conformational state for selective ligand binding throughout all stages of activation.

Binding pocket volume and docking calculations for metastable states of cannabinoid receptors.Binding pocket volumes for each metastable state of CB1 (a) and CB2 (b) are plotted against N-terminus distance as a scatter plot. Docking affinity (kcal/mol) of the four CB2 selective agonists for each metastable state of CB1 (c) and CB2 (d) are plotted as a bar plot. Docking and volume calculations were performed on 100 structures from each metastable state selected based on the MSM probabilities. Colors for CB1 metastable states change from orange to blue gradually from inactive to active state. Colors for CB2 metastable states change from green to magenta gradually from inactive to active state. For each protein structure, docking affinity is obtained as an average of top three dock poses. Error shown in docking calculations is the standard error which the standard deviation of docking affinity distribution obtained from 100 structures divided by 1002. — Fig. 7Binding pocket volume and docking calculations for metastable states of cannabinoid receptors.Binding pocket volumes for each metastable state of CB₁ (a) and CB₂ (b) are plotted against N-terminus distance as a scatter plot. Docking affinity (kcal/mol) of the four CB₂ selective agonists for each metastable state of CB₁ (c) and CB₂ (d) are plotted as a bar plot. Docking and volume calculations were performed on 100 structures from each metastable state selected based on the MSM probabilities. Colors for CB₁ metastable states change from orange to blue gradually from inactive to active state. Colors for CB₂ metastable states change from green to magenta gradually from inactive to active state. For each protein structure, docking affinity is obtained as an average of top three dock poses. Error shown in docking calculations is the standard error which the standard deviation of docking affinity distribution obtained from 100 structures divided by $\sqrt[2]{100}$ .

To find whether CB₂ selective agonist binding affinity depends on the binding pocket shape, we docked CB₂ selective ligands into CB₁ and CB₂ metastable states (described in the method sections). Four CB₂ selective ligands (JWH-133, HU-308, JWH-015, AM-1241) are selected from the literature³², which shows at least ten-fold selectivity (Supplementary Fig. 23). Docking results of CB₂ selective ligands to CB₁ show that docking affinity differs for different metastable states (Fig. 7c). For inactive and I1 states where the binding pocket volume is small, the docking is hard due to the steric hindrance of the N-terminus. In the other four states, docking results show variational docking affinity, with I2 states having the highest affinity for all the CB₂ selective ligands compared to the active state. Therefore, in the CB₁ active state, these ligands have a less binding probability. This differential binding affinity is explained by the differential binding pose of the ligands in I2 and active states. In I2, there is a larger gap between TM2 and ECL2, which allow the ligands to bind between the space of ECL2, TM1, and TM2 (Supplementary Fig. 24a, b). In the active state, the ligands bind deep in the pocket due to the movement of the toggle switch. Toggle switch movement creates a space for ligand to bind (Fig. 8a, c, e, g).

Docking poses of CB2 selective ligands in metastable states of cannabinoid receptors.Differences in most stable docking poses between active (color: Blue) and I2 (color: Brown) metastable states for CB1 for four different CB2 selective ligands (JWH-133 (a), HU-308 (c), JWH-015 (e), AM1241(g)). Differences in most stable docking poses between active (color: Magenta) and inactive (color: Green) metastable states for CB2 for four different CB2 selective ligands (JWH-133 (b), HU-308 (d), JWH-015 (f), AM1241(h)). — Fig. 8Docking poses of CB₂ selective ligands in metastable states of cannabinoid receptors.Differences in most stable docking poses between active (color: Blue) and I2 (color: Brown) metastable states for CB₁ for four different CB₂ selective ligands (JWH-133 (a), HU-308 (c), JWH-015 (e), AM1241(g)). Differences in most stable docking poses between active (color: Magenta) and inactive (color: Green) metastable states for CB₂ for four different CB₂ selective ligands (JWH-133 (b), HU-308 (d), JWH-015 (f), AM1241(h)).

To contrast the binding of ligands in different states, we performed MD simulations on the ligand bound systems of representative bound structures from intermediate state 2 and active state of the CB₁. Each system (in total 8 systems) was run for approximately 200 ns (Supplementary Note 3). RMSD analyses show that ligand RMSD remains within 2.5 Å of the initial predicted docking pose, depicting the stability of the predicted pose (Supplementary Figs. 25 and 26; Supplementary Table 6). To show the difference in binding positions of the ligands, we computed the stable dynamic contacts and corresponding interaction energies between the ligand and receptor in both the metastable states of the receptors. Calculations of the dynamic contacts show that there are major differences in stable contacts between the two states of CB₁ (Supplementary Fig. 27). In active state of CB₁, the ligands are bound deeper inside the pockets. Hence, in the active state, the ligands are forming major interactions with the TM5 residues. The interaction energy contributions for the residues reveals the important residues in different metastable state (Supplementary Fig. 28). K192^3.92 contributes to CB₁ I2 state binding for every ligand, whereas W279^5.43 contributes distinctively to CB₁ active state binding for JWH-133, JWH-105 and AM1241.

On the other hand, the docking affinities for these ligands are similar to each metastable state for CB₂ as the receptor pocket shape remains relatively similar (Figs. 7d, 8b, d, f, h). Further MD simulations were performed for holo (ligand bound) active and inactive states of CB₂. RMSD calculations reveals the ligand stability of the docking predicted poses in both the metastable states (Supplementary Figs. 25 and 26; Supplementary Table 6). As the ligand binding poses are relatively similar in distinct metastable states of CB₂, the differences in number of stable contacts and corresponding energy contributions are also small compared to CB₁ (Supplementary Fig. 29 and Supplementary Table 7). Therefore, these analyses show that how the difference in pocket shape due to distinct activation mechanism of CB₁ and CB₂ may lead to ligand selectivity towards CB₂.

Conclusions

CB₂ selective agonists are emerging as potential drug targets for treating inflammation. A number of CB₂ selective agonists have entered clinical trial phase^36,37. However, none of these molecules has been approved as a drug by FDA. There are two of the major drawbacks in developing CB₂ selective drugs. First, high lipophilic nature of the CB₂ selective agonists are not suitable to be drug target. Second, some CB₁ activity may retain in those molecules which leads to off-target side-effect. In this work, we mechanistically explain the selectivity of CB₂ selective agonists by hypothesizing that reason for selectivity lies in the difference of activation mechanism. Here, we have studied the activation mechanism of both CB₁ and CB₂ by with combinations of using molecular dynamics, Markov state model and neural network.

Our results showed three major distinctions in conformational changes in both receptors. (1) In the extracellular region due to the movement of the N-terminus, we observe a large change in the binding pocket volume during the activation of the CB₁, whereas, for CB₂ there are hardly any changes in the binding pocket. (2) For CB₁, toggle switch W^6.48 undergoes translational changes during the activation mechanism, which lead to the changes in the intracellular matrices. While, for CB₂, toggle switch W^6.48 only rotates in the same position. Our analysis also shows that toggle switch rotation does not effect the movement intracellular movements for CB₂. (3) Intracellular TM5 movement of CB₂ affects the TM6 and TM7 movement critical for receptor activation. On the other hand, no major conformational change is observed in intracellular TM5 of CB₁.

VAMPnets analysis further reveals the metastable states and sequential transition of the structural features during the activation. During the CB₁ activation, initial extracellular movements of N-terminus, TM1 lead to the change in further change in the transmembrane and intracellular features. Conversely, due to lack of the extracellular changes, the activation of CB₂ is kinetically faster compared to the CB₁.

Further analyses show that due to the lack of the extracellular movement during CB₂ activation, binding pocket volume remains relatively same for CB₂. Docking of metastable states show that CB₂ selective agonists show the similar high binding affinities for each metastable states of CB₂, depicting the preferential binding of these ligands in the generalized CB₂ binding conformation. Whereas, these ligands can only be selectively bind to few metastable states in CB₁ due to the larger conformational change in the CB₁ binding pocket. This explains lack of overall binding affinity of these ligands towards CB₁. Overall, we provide mechanistic knowledge about cannabinoid receptor selectivity by studying the activation mechanism. This study will guide to design new CB₂ selective agonists with better druggability profile by acting as virtual screening criteria for these ligands.

Methods

System preparation

To perform apo (without ligand) simulation, we started molecular dynamics simulations from inactive (CB₁ inactive state PDBID: 5TGZ⁹; CB₂ inactive state PDBID: 5ZTY¹²) and active (CB₁ active state PDBID: 5XRA¹¹; CB₂ active state PDBID: 6KPF¹⁶) state structures of CB₁ and CB₂. Non-protein residues were deleted from the CB₁ and CB₂ crystal structure. The thermostabilized mutations were mutated back to original residue using tleap. Hydrogen atoms were added to protein amino acid residues using reduce command of AMBER package⁵³. To add consistency in the residue numbering, 4 residues (CB₁ residue number 99-103) were added in the truncated N-terminus for CB₁ active structure. Intracellular region between TM5 and TM6 was truncated for both CB₁ and CB₂. Terminal end of N-terminus and C-terminus & unconnected residues of TM5 and TM6 were neutralized with ACE and NME residue. Both CB₁ and CB₂ was embedded in POPC bilayer with 150 mM NaCl solution in both extracellular and intracellular direction. TIP3P water model was used for MD system as this water model gives better performance in ref. ⁵⁴. AMBER forcefield was selected to parameterize the MD systems as the this forcefield was used to perform simulations for other GPCRs and all the necessary parameters for molecules in MD system are available in this forcefield^55–57. AMBER ff14SB and lipid17 forcefield was used for protein and lipid parameterization^53,58. Agonist bound active (CB₁ agonist: AM11542, CB₂ agonist: AM841) and antagonist inactive (CB₁ antagonist: AM6538, CB₂ antagonist: AM10257) holo simulation were also performed. GAFF forcefield parameters for the ligands are obtained using antechamber^59,60. Similar modifications were made in the proteins using tleap and embedded in bilayer using CHARMM-GUI⁶¹.

Simulation details

Atomistic unbiased molecular dynamics simulation was performed with AMBER v18 MD engine⁶². Before running the production MD run, prepared MD systems were minimized and equilibrated. Minimization was performed for 15,000 steps; where the systems were minimized for 5000 steps with gradient descent algorithm and followed by conjugate gradient algorithm for the rest 10,000 steps. Minimized systems was heated and pressurized sequentially in NVT and NPT ensemble, respectively. Heating step was performed for 2 ns where system was heated from 0 K to 300 K in NVT ensemble. Anisotropic pressure was applied to the heated systems to reach the constant pressure of 1 bar with NPT ensemble. Pressurization step was also performed for 2 ns. Berendsen thermostat and barostat was used to fix the temperature and pressure⁶³. During the heating and pressurization, protein backbone for each system was restrainted using a spring force with a spring constant of 10 kcal/mol/Å². Subsequently, restraint was removed and each system was equilibrated for 46 ns for collecting data. Numerical integration was performed using verlet algorithm to update velocity and coordinate. Time step for each numerical integration is 2 fs. Shake algorithm was implemented by constraint bond movement of the hydrogens to bonded to heavy atoms⁶⁴. Periodic boundary condition was applied to each MD system. Particle mesh Ewald method (PME) was implemented for long-range electrostatic calculation⁶⁵.

Featurization of protein conformational ensemble

To perform adaptive sampling and markov state model (MSM) building (Discussed in the following sections), protein conformations in MD system need to be represented by the descriptors that can capture protein conformational changes. To find the descriptors, we have implemented residue-residue contact score (RRCS) analysis^66,67. Using RRCS method, a score was assigned to each residue pair distances for a particular protein conformation which are more than one helical tern away. This score was calculated based on the every atom pair distances between two residues following the same algorithm of Zhou et al.⁶⁶. To find the major structural changes during activation for CB₁ and CB₂, we compared the change in RRCS between the active and inactive structure. α-carbon distances of the residue pairs, where the RRCS is changing more than 3, were used as a feature to describe a protein conformation (Supplementary Tables 8, 9). As shown in Supplementary Fig. 30, these descriptors can capture the extracellular, transmembrane and intracellular changes for both CB₁ and CB₂.

Adaptive sampling

Previous experimental and computational studies depict that GPCR activation occurs in a microsecond to millisecond timescale^20,68–70. However, each numerical integration time step for MD simulation is of femtosecond timescale⁷¹. Therefore, capturing the entire activation process by running a single long trajectory is time consuming process. Furthermore, to capture the thermodynamic and kinetics of the conformational changes better sampling is needed in high energy transition states⁷². Thus, to capture the entire protein conformational ensemble, adaptive sampling protocol was implemented^72–75. Adaptive sampling is a well established sampling technique used for studying ligand binding^76–78, protein conformational change^{67,68,79–81} and ligand selectivity⁸². First, conformations obtained from MD simulations are expressed in terms of collective variables. Second, collective variable space are clustered into different states using k-means clustering. Third, based on the population of the each state, MD frames are selected from least populated clusters to start simulation for the next round. GPCRs consist of the multiple allosterically coupled conversed motifs distributed in extracellular, transmembrane and intracellular regions. These motifs structurally rearrange themselves when protein conformational change happens from inactive to active. Therefore, finding suitable adaptive sampling matrices to minimize the amount of simulation time is difficult. Here, we used our preexisting knowledge of the structural changes to select CVs and modify the CV based on the requirement. For initial round of sampling, we selected the extracellular TM1-TM6 distance and intracellular TM3-TM6 distance as a collective variable for adaptive sampling. For both CB₁ and CB₂, sampling was performed from both active and inactive structures until energy landscape was connected. Next, 24 and 20 distinct features (Extracellular distances, Toggle switch movements, Intracellular TM6 and TM7 movements) were selected for CB₁ and CB₂ to improve the speed of the sampling (Supplementary Tables 10,11). To check whether simulations able to observe the transition between the inactive to active states in the RRCS descriptor space (as discussed in previous section), we performed dimensionality reduction using time independent component analysis. On the projection of slowest two tic components, transition between active and inactive state was observed in CB₂ conformational ensemble. However, there was a disconnect between active and inactive minima in tic projection of CB₁. Therefore, additional rounds of sampling were performed from the clusters that are closest from each other in the two disconnected regions to observe the transition between active and inactive minima is observed for CB₁. In total, 490 and 278 μs simulations were performed to capture apo activation process for CB₁ and CB₂ (Supplementary Fig. 34a, b). Adaptive sampling for holo systems is performed on the 24, 20 structurally known features for CB₁ and CB₂ holo systems for approximately 20 μs.

Markov state models

From adaptive sampling, we generate plethora of short trajectories that only capture the small changes in protein conformation. Markov state model was developed to connect the information from individual trajectories to by building transition network between distinct states^83–85. In markovian processes, the future step of the process should only depend on present states, not on the past states⁸⁶. Molecular dynamic simulations follow the same principle as the position and velocity of the future step is calculated from the energy and force calculation based on the current state. From MSM, we obtain stationary probability distribution of the protein ensemble and timescales of the slowest process by calculating transition probability matrix (TM). TM calculates probability of transition between discretized protein space by calculating the jump between different states after lag time. Lag time is the minimum time at which markovian property of the discretized states is valid. Building markov state model consist of four different steps: featurization, dimensionality reduction, clustering and hyperparameter optimization. Featurization of the protein conformational space was done from the distances calculated from RRCS (discussed in previous subsection). Dimensionality reduction was performed on the featurized space with time independent component analysis (tICA)⁸⁷. MSM captures the slowest timescale processes. Using tICA, we obtain new representations of features which is a linear combination of slowest features. Therefore, in the next discretization step, clustering algorithm can discretize the protein space which are distinct in slowest timescale. K-Means clustering algorithm was implemented in this case. The hyperparameters that needs to be optimized are lagtime at which markovian property of the discretized space is valid, number of tic dimensions on which clustering is performed and number of cluster components or microstates that can distinguish the slowest components. With a minimum lag time of 25 ns for both the systems, slowest timescales obtained from MSM converges, atleast in the log scale (Supplementary Fig. 31a, b). Hence, the markovian property is valid. For further testing of markovian property of models at lag time of 25 ns, Chapman-Kolmogorov test (C-K test) was performed (Supplementary Fig. 32a, b)⁸⁸. To optimize other two hyperparameters, the VAMP-2 score was calculated^89,90. For reversible process, which follows the detailed balance, VAMP-2 score is the sum of the squares of the highest eigenvalues⁹¹. MSM with the highest VAMP-2 score better captures the slowest timescale processes. To avoid overfitting for calculation of the VAMP-2 scores, 10-fold cross validation process was used. The hyperparameters which maximizes the VAMP-2 scores are considered for the final MSM estimation (Supplementary Fig. 31c, d). Final optimized MSM reweighed the conformational populations obtained from adaptive sampling (Supplementary Fig. 33a, b). Pyemma v2.5.6 software package was used for entire MSM building process⁹².

Trajectory analysis

Features such as distances, angles are calculated using python library MDtraj v1.9.8⁹³. VMD v1.9.3 software package was used for protein figure making and trajectory visualization⁹⁴. Ambertools CPPTraj v18.01 was used for frame selection and trajectory modulations⁹⁵. All the analysis codes (e.g. conditional probability calculation) were written in python programming language and matplotlib library was used for plotting the figures as a graph. POVME v3.0 software package was used for volume calculation of the binding pocket⁹⁶.

Metastable state estimatation using VAMPnets

VAMPnets is a deep learning architecture used for building coarse grain MSM^48,91,97. This architecture consists of two parallel deep learning networks which take the input feature values at time t and t+τ, where τ is the lagtime. This deep learning model is optimized by maximizing the VAMP-2 score. The output of the VAMPnets is the probability of the each frame belonging to a certain metastable state. VAMPnets was implemented using python library deeptime v0.4.3⁹⁸. RRCS calculated features were as an input. Lagtime for the VAMPnets is selected as 25 ns. The number of the metastable states for each system was selected by the two creteria : (1) Each metastable state should have atleast 4 percent of the total population (Supplementary Fig. 10, 11). (2) error bar in the timescale calculations is minimum (Supplementary Figs. 8, 9). Based on these creteria, the number of metastable states chosen for CB₁ and CB₂ is six.

Data-driven analysis

To capture the major structural changes between different macro states, closest heavy carbon atom distances was calculated between every residue pair of each frame. These calculations were performed on 1000 frames from each metastable state. These frames are selected based on the the probability of microstates belonging to the metastable states. Distance distribution between certain residue pair is compared to the identical residue pair distance distribution between distinct metastable states. Comparison is performed with symmetric K-L divergence analysis. To estimate the contribution of each residue, average K-L divergence of the all the residue pair distances was calculated that are associated with that residue⁴⁹.

Allosteric pathway calculation

To predict the allosteric pathway in different metastable states, similar procedure as explained in Bhattacharya and Vaidehi⁵¹. First, normalized mutual information for each residue pair was calculated per metastable state based on dihedral angle movement. These calculations was performed on the 1000 frames as discussed in Data-driven analysis section. Based on the residue pairs, where closest heavy atom are within 5 Å, an undirected acyclic graph was created. Nodes in the graph denote a residue and edges are between the two residues which are within 5 Å of one another. Edge weight represents the MI_max - MI_ij, where MI_ij is the mutual information between the two nodes i and j. Shortest allosteric path was calculated for the residue pairs that are atleast 12 Å away from each other. Among all the calculated paths, 500 paths with highest mutual information were selected. These paths were clustered into optimal number of clusters based on the procedure explained in Bhattacharya and Vaidehi⁵¹. Mutual information calculations have been performed with MDentropy v0.3.0 software package⁹⁹.

Docking calculations

To find the docking pose and the docking score of the CB₂ selective ligands, in distinct metastable states of CB₁ and CB₂, four CB₂ selective ligands were docked into binding pocket of distinct metastable states of CB₁ and CB₂. Four ligands were selected based on the criterion that those ligands show at least 10 fold more affinity for CB₂ than CB₁³². Ligand structures were downloaded from pubchem using sdf format which is converted into the pdb format using Charmm-gui ligand builder^61,100. Amber tools Antechamber is used for converting the pdb format to mol2 format^59,60. For docking, mol2 files were converted to pdbqt format. To represent each metastable state, 100 protein structures were selected from each metastable state proportional to the probability of microstates and converted to pdbqt format. Autodock-vina v1.2.3 software was used for docking calculations^101,102. To specify the docking center, center of mass of alpha-carbon of binding pocket residues for CB₁ and CB₂ was selected (Supplementary Fig. 35).

Transition path theory

Transition path theory was implemented to calculate the effective timescale between metastable state transition^103,104. Effective timescale between two metastable states (A and B) is obtained from mean free passage time (1/k_AB) as shown in equation 1/k_AB=( $τ \sum_{i = 1}^{m} π_{i} (1 - q_{i}^{+}) / F$ , where F is the total flux between metastable state A to B and $q_{i}^{+}$ is the committer probability which is defined as the probability of microstate i reaching the metastable state B before A. A microstate is assigned to a metastable state based on the highest probability of belonging to each metastable state. Thus, six metastable states are defined as distinct microstates. TPT calculations were performed using Pyemma v2.5.6 software⁹².

Error calculation

To calculate errors for thermodynamic and kinetic quantities obtained from simulation, bootstrapping approach was followed^81,88. For a certain calculation, N numbers of random bootstrap samples were generated with 80% of total trajectories. Based on the standard deviation of a calculated quantity from bootstrapped samples, error bars were generated. For example, for error calculation in free energy plots, MSM was created for 200 times using 80% of the original data to find the standard deviation in the probability density.

Kinetic Monte Carlo simulation

Based on the transition matrix of MSM, kinetic monte carlo (kMC) simulation was implemented to observe the evolution of features with time^68,105. kMC simulation is a stochastic process where in a given time step, next transition is selected based on the probability of the all possible transitions. The algorithm for KMC simulation is as follows: (1) Simulation starts at a particular microstate (obtained from clustering; as discussed in earliar section). (2) From the MSM transition matrix, the probabilities of all possible transitions (P_i,j) from that microstate are obtained. As the sum of all transition probabilities is 1, a cumulative distribution (S) of possible transitions is created where $S_{i, j} = \sum_{k = 0}^{k = j} P_{i, k}$ . A random number (R) is generated between 0 to 1 to find the next possible transition. If R is between S_i,j and S_i,j+1, ith state to j + 1th state transition is selected and simulation moves forward by lag time (τ). This algorithm is followed iteratively. In this case, kMC simulation was initialized from the inactive state of both the proteins.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information

Untitled section

Acknowledgements

This research is 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. D.S. acknowledges support from NIGMS MIRA award R35GM142745. S.D. thanks Balaji Selvam for the help during the initial phase of the project and Austin T. Weigle for suggestions during manuscript writing.

Author contributions

D.S. and S.D. designed the research. S.D. performed simulations, analyzed data, and wrote the article. D.S. participated in discussing results. D.S. supervised the research. All authors reviewed the manuscript.

Peer review

Peer review information

Communications Biology thanks Shaoyong Lu and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Yun Lyna Luo and Gene Chong.

Data availability

All the necessary files (MSM object, bootstrap files, feature files) and trajectories (apo and holo simulation) that have been used for the analysis to generate figures can be obtained from https://uofi.box.com/s/jzooa0o27z1w9ha0h6va3i51ir7l38j4.

Code availability

All the necessary python scripts have been used to generate figures can be obtained from https://github.com/ShuklaGroup/Cannabinoid_activation.git.

Competing interests

The authors declare no competing interests.

Footnotes

Footnote Group

Supplementary information

The online version contains supplementary material available at 10.1038/s42003-023-04868-1.

References

Untitled section

References

1.Herkenham M, et al. Characterization and localization of cannabinoid receptors in rat brain: a quantitative in vitro autoradiographic study. J. Neurosci. 1991;11:563–583. doi: 10.1523/JNEUROSCI.11-02-00563.1991.
2.Munro S, Thomas KL, Abu-Shaar M. Molecular characterization of a peripheral receptor for cannabinoids. Nature. 1993;365:61–65. doi: 10.1038/365061a0.
3.Battista N, Tommaso MD, Bari M, Maccarrone M. The endocannabinoid system: An overview. Front. Behav. Neurosci. 2012;6:9. doi: 10.3389/fnbeh.2012.00009.
4.Mechoulam R, Parker LA. The endocannabinoid system and the brain. Ann. Rev. Psychol. 2013;64:21–47. doi: 10.1146/annurev-psych-113011-143739.
5.Lu H-C, Mackie K. An introduction to the endogenous cannabinoid system. Biol. Psychiatry. 2016;79:516–525. doi: 10.1016/j.biopsych.2015.07.028.
6.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. doi: 10.3390/ijms19030833.
7.Walsh KB, Andersen HK. Molecular pharmacology of synthetic cannabinoids: Delineating CB1 receptor-mediated cell signaling. Int. J. Mol. Sci. 2020;21:6115. doi: 10.3390/ijms21176115.
8.Malfitano AM, Basu S, Maresz K, Bifulco M, Dittel BN. What we know and do not know about the cannabinoid receptor 2 (CB2) Semin. Immunol. 2014;26:369–379. doi: 10.1016/j.smim.2014.04.002.
9.Hua T, et al. Crystal structure of the human cannabinoid receptor CB1. Cell. 2016;167:750–762.e14. doi: 10.1016/j.cell.2016.10.004.
10.Shao Z, et al. High-resolution crystal structure of the human CB1 cannabinoid receptor. Nature. 2016;540:602–606. doi: 10.1038/nature20613.
11.Hua T, et al. Crystal structures of agonist-bound human cannabinoid receptor CB1. Nature. 2017;547:468–471. doi: 10.1038/nature23272.
12.Li X, et al. Crystal structure of the human cannabinoid receptor CB2. Cell. 2019;176:459–467.e13. doi: 10.1016/j.cell.2018.12.011.
13.Kumar KK, et al. Structure of a signaling cannabinoid receptor 1-g protein complex. Cell. 2019;176:448–458.e12. doi: 10.1016/j.cell.2018.11.040.
14.Shao Z, et al. Structure of an allosteric modulator bound to the CB1 cannabinoid receptor. Nat. Chem. Biol. 2019;15:1199–1205. doi: 10.1038/s41589-019-0387-2.
15.Xing C, et al. Cryo-EM structure of the human cannabinoid receptor CB2-gi signaling complex. Cell. 2020;180:645–654.e13. doi: 10.1016/j.cell.2020.01.007.
16.Hua T, et al. Activation and signaling mechanism revealed by cannabinoid receptor-gi complex structures. Cell. 2020;180:655–665.e18. doi: 10.1016/j.cell.2020.01.008.
17.Wang X, et al. A genetically encoded f-19 NMR probe reveals the allosteric modulation mechanism of cannabinoid receptor 1. J. American Chem. Soc. 2021;143:16320–16325. doi: 10.1021/jacs.1c06847.
18.Rosenbaum DM, Rasmussen SGF, Kobilka BK. The structure and function of g-protein-coupled receptors. Nature. 2009;459:356–363. doi: 10.1038/nature08144.
19.Latorraca NR, Venkatakrishnan AJ, Dror RO. GPCR dynamics: Structures in motion. Chem. Rev. 2016;117:139–155. doi: 10.1021/acs.chemrev.6b00177.
20.Weis WI, Kobilka BK. The molecular basis of g protein–coupled receptor activation. Ann. Rev. Biochem. 2018;87:897–919. doi: 10.1146/annurev-biochem-060614-033910.
21.Hauser AS, Attwood MM, Rask-Andersen M, Schiöth HB, Gloriam DE. Trends in GPCR drug discovery: New agents, targets and indications. Nat. Rev. Drug Discov. 2017;16:829–842. doi: 10.1038/nrd.2017.178.
22.Odoemelam CS, et al. G-protein coupled receptors: Structure and function in drug discovery. RSC Adv. 2020;10:36337–36348. doi: 10.1039/D0RA08003A.
23.Yang D, et al. G protein-coupled receptors: Structure- and function-based drug discovery. Signal Transduc. Targeted Ther. 2021;6:7. doi: 10.1038/s41392-020-00435-w.
24.Pertwee RG. Pharmacology of cannabinoid receptor ligands. Curr. Med. Chem. 1999;6:635–664. doi: 10.2174/0929867306666220401124036.
25.Howlett AC. International union of pharmacology. XXVII. classification of cannabinoid receptors. Pharmacol. Rev. 2002;54:161–202. doi: 10.1124/pr.54.2.161.
26.Pertwee RG, Ross R. Cannabinoid receptors and their ligands. Prostaglandins, Leukotrienes Essential Fatty Acids (PLEFA) 2002;66:101–121. doi: 10.1054/plef.2001.0341.
27.Pertwee RG. Cannabinoid pharmacology: The first 66 years. Brit. J. Pharmacol. 2006;147:S163–S171. doi: 10.1038/sj.bjp.0706406.
28.Shafi A, Berry AJ, Sumnall H, Wood DM, Tracy DK. New psychoactive substances: A review and updates. Ther. Adv. Psychopharmacol. 2020;10:204512532096719. doi: 10.1177/2045125320967197.
29.Wiley JL, Marusich JA, Huffman JW. Moving around the molecule: Relationship between chemical structure and in vivo activity of synthetic cannabinoids. Life Sci. 2014;97:55–63. doi: 10.1016/j.lfs.2013.09.011.
30.Davis MP. Cannabinoids in pain management: CB1, CB2 and non-classic receptor ligands. Expert Opin. Investigational Drugs. 2014;23:1123–1140. doi: 10.1517/13543784.2014.918603.
31.Li X, Shen L, Hua T, Liu Z-J. Structural and functional insights into cannabinoid receptors. Trends Pharmacol. Sci. 2020;41:665–677. doi: 10.1016/j.tips.2020.06.010.
32.Pertwee R. Receptors and channels targeted by synthetic cannabinoid receptor agonists and antagonists. Curr. Med. Chem. 2010;17:1360–1381. doi: 10.2174/092986710790980050.
33.Rogers N. Cannabinoid receptor with an ‘identity crisis’ gets a second look. Nat. Med. 2015;21:966–967. doi: 10.1038/nm0915-966.
34.Soethoudt, M. et al. Cannabinoid CB2 receptor ligand profiling reveals biased signalling and off-target activity. Nat. Commun.8, 13958 (2017).
35.Li X, Shen L, Hua T, Liu Z-J. Structural and functional insights into cannabinoid receptors. Trends Pharmacol. Sci. 2020;41:665–677. doi: 10.1016/j.tips.2020.06.010.
36.Magham SV, krishnamurthy PT, Shaji N, Mani L, Balasubramanian S. Cannabinoid receptor 2 selective agonists and alzheimers disease: An insight into the therapeutic potentials. J. Neurosci. Res. 2021;99:2888–2905. doi: 10.1002/jnr.24933.
37.Whiting ZM, Yin J, de la Harpe SM, Vernall AJ, Grimsey NL. Developing the cannabinoid receptor 2 (CB2) pharmacopoeia: Past, present, and future. Trends Pharmacol. Sci. 2022;43:754–771. doi: 10.1016/j.tips.2022.06.010.
38.Peng Y, et al. 5-HT2c receptor structures reveal the structural basis of GPCR polypharmacology. Cell. 2018;172:719–730.e14. doi: 10.1016/j.cell.2018.01.001.
39.Katritch V, Kufareva I, Abagyan R. Structure based prediction of subtype-selectivity for adenosine receptor antagonists. Neuropharmacology. 2011;60:108–115. doi: 10.1016/j.neuropharm.2010.07.009.
40.Ostopovici-Halip L, Curpăn R, Mracec M, Bologa CG. Structural determinants of the alpha2 adrenoceptor subtype selectivity. J. Mol. Graphics Modelling. 2011;29:1030–1038. doi: 10.1016/j.jmgm.2011.04.011.
41.Sallander J, Wallinder C, Hallberg A, Åqvist J, de Terán HG. Structural determinants of subtype selectivity and functional activity of angiotensin II receptors. Bioorganic Med. Chem. Lett. 2016;26:1355–1359. doi: 10.1016/j.bmcl.2015.10.084.
42.Rataj K, et al. Fingerprint-based machine learning approach to identify potent and selective 5-HT2br ligands. Molecules. 2018;23:1137. doi: 10.3390/molecules23051137.
43.Ma C, Wang L, Xie X-Q. Ligand classifier of adaptively boosting ensemble decision stumps (LiCABEDS) and its application on modeling ligand functionality for 5ht-subtype GPCR families. J. Chem. Inform. Modeling. 2011;51:521–531. doi: 10.1021/ci100399j.
44.Rodríguez D, Brea J, Loza MI, Carlsson J. Structure-based discovery of selective serotonin 5-HT 1b receptor ligands. Structure. 2014;22:1140–1151. doi: 10.1016/j.str.2014.05.017.
45.Glukhova A, et al. Structure of the adenosine a1 receptor reveals the basis for subtype selectivity. Cell. 2017;168:867–877.e13. doi: 10.1016/j.cell.2017.01.042.
46.Ma C, Wang L, Yang P, Myint KZ, Xie X-Q. LiCABEDS II. modeling of ligand selectivity for g-protein-coupled cannabinoid receptors. J. Chem. Inform. Modeling. 2013;53:11–26. doi: 10.1021/ci3003914.
47.Hu B, et al. Three-dimensional biologically relevant spectrum (BRS-3d): Shape similarity profile based on PDB ligands as molecular descriptors. Molecules. 2016;21:1554. doi: 10.3390/molecules21111554.
48.Mardt A, Pasquali L, Wu H, Noé F. VAMPnets for deep learning of molecular kinetics. Nat. Commun. 2018;9:5. doi: 10.1038/s41467-017-02388-1.
49.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. 2021;10:e60715. doi: 10.7554/eLife.60715.
50.Maharana J, Banerjee R, Yadav MK, Sarma P, Shukla AK. Emerging structural insights into GPCR–β-arrestin interaction and functional outcomes. Curr. Opin. Struc. Biol. 2022;75:102406. doi: 10.1016/j.sbi.2022.102406.
51.Bhattacharya S, Vaidehi N. Differences in allosteric communication pipelines in the inactive and active states of a gpcr. Biophys. J. 2014;107:422–434. doi: 10.1016/j.bpj.2014.06.015.
52.Paul F, Weikl TR. How to distinguish conformational selection and induced fit based on chemical relaxation rates. PLOS Comput. Biol. 2016;12:e1005067. doi: 10.1371/journal.pcbi.1005067.
53.Case, D. et al. AMBER 2018. University of California, San Francisco.https://ambermd.org/doc12/Amber18.pdf (2020).
54.Mark P, Nilsson L. Structure and dynamics of the TIP3p, SPC, and SPC/e water models at 298 k. J. Phys. Chem. A. 2001;105:9954–9960. doi: 10.1021/jp003020w.
55.Li Y, Sun J, Li D, Lin J. The full activation mechanism of the adenosine a1 receptor revealed by gamd and su-gamd simulations. Proc. Natl. Acad. Sci. 2022;119:e2203702119. doi: 10.1073/pnas.2203702119.
56.Ding T, Karlov DS, Pino-Angeles A, Tikhonova IG. Intermolecular interactions in g protein-coupled receptor allosteric sites at the membrane interface from molecular dynamics simulations and quantum chemical calculations. J. Chem. Inform. Modeling. 2022;62:4736–4747. doi: 10.1021/acs.jcim.2c00788.
57.Li, Y., Yin, C., Liu, P., Li, D. & Lin, J. Identification of a different agonist-binding site and activation mechanism of the human p2y1 receptor. Scientific Rep.7, 13764 (2017).
58.Maier JA, et al. ff14sb: improving the accuracy of protein side chain and backbone parameters from ff99sb. J. Chem. Theory Comput. 2015;11:3696–3713. doi: 10.1021/acs.jctc.5b00255.
59.Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing of a general amber force field. J. Computat. Chem. 2004;25:1157–1174. doi: 10.1002/jcc.20035.
60.Wang J, Wang W, Kollman PA, Case DA. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modelling. 2006;25:247–260. doi: 10.1016/j.jmgm.2005.12.005.
61.Jo S, Kim T, Iyer VG, Im W. Charmm-gui: a web-based graphical user interface for charmm. J. Comput. Chem. 2008;29:1859–1865. doi: 10.1002/jcc.20945.
62.Salomon-Ferrer R, Götz AW, Poole D, Grand SL, Walker RC. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. explicit solvent particle mesh ewald. J. Chem. Theory Comput. 2013;9:3878–3888. doi: 10.1021/ct400314y.
63.Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984;81:3684–3690. doi: 10.1063/1.448118.
64.Krautler V, van Gunsteren WF, Hunenberger PH. A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 2001;22:501–508. doi: 10.1002/1096-987X(20010415)22:5<501::AID-JCC1021>3.0.CO;2-V.
65.Essmann U, et al. A smooth particle mesh ewald method. J. Chem. Phys. 1995;103:8577–8593. doi: 10.1063/1.470117.
66.Zhou Q, et al. Common activation mechanism of class a GPCRs. eLife. 2019;8:50279. doi: 10.7554/eLife.50279.
67.Bansal PD, Dutta S, Shukla D. Activation mechanism of the human smoothened receptor. Biophys. J. 2023;122:1400–1413. doi: 10.1016/j.bpj.2023.03.007.
68.Kohlhoff KJ, et al. Cloud-based simulations on google exacycle reveal ligand modulation of GPCR activation pathways. Nat. Chem. 2013;6:15–21. doi: 10.1038/nchem.1821.
69.Dror RO, et al. Activation mechanism of the β 2-adrenergic receptor. Proc. Natl. Acad. Sci. 2011;108:18684–18689. doi: 10.1073/pnas.1110499108.
70.Vilardaga J-P. Theme and variations on kinetics of GPCR activation/deactivation. J. Recep. Signal Transduc. 2010;30:304–312. doi: 10.3109/10799893.2010.509728.
71.Braun E, et al. Best practices for foundations in molecular simulations [article v1.0] Living J. Comput. Molecular Sci. 2019;1:5957. doi: 10.33011/livecoms.1.1.5957.
72.Zimmerman MI, Porter JR, Sun X, Silva RR, Bowman GR. Choice of adaptive sampling strategy impacts state discovery, transition probabilities, and the apparent mechanism of conformational changes. J. Chem. Theory Comput. 2018;14:5459–5475. doi: 10.1021/acs.jctc.8b00500.
73.Bowman GR, Ensign DL, Pande VS. Enhanced modeling via network theory: Adaptive sampling of markov state models. J. Chem. Theory Comput. 2010;6:787–794. doi: 10.1021/ct900620b.
74.Kleiman DE, Shukla D. Multiagent reinforcement learning-based adaptive sampling for conformational dynamics of proteins. J. Chem. Theory Comput. 2022;18:5422–5434. doi: 10.1021/acs.jctc.2c00683.
75.Shamsi Z, Cheng KJ, Shukla D. Reinforcement learning based adaptive sampling: REAPing rewards by exploring protein conformational landscapes. J. Phys. Chem. B. 2018;122:8386–8395. doi: 10.1021/acs.jpcb.8b06521.
76.Dutta S, Selvam B, Das A, Shukla D. Mechanistic origin of partial agonism of tetrahydrocannabinol for cannabinoid receptors. J. Biol. Chem. 2022;298:101764. doi: 10.1016/j.jbc.2022.101764.
77.Shukla S, Zhao C, Shukla D. Dewetting controls plant hormone perception and initiation of drought resistance signaling. Structure. 2019;27:692–702.e3. doi: 10.1016/j.str.2018.12.005.
78.Buch I, Giorgino T, Fabritiis GD. Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. 2011;108:10184–10189. doi: 10.1073/pnas.1103547108.
79.Betz RM, Dror RO. How effectively can adaptive sampling methods capture spontaneous ligand binding? J. Chem. Theory Comput. 2019;15:2053–2063. doi: 10.1021/acs.jctc.8b00913.
80.Doerr S, De Fabritiis G. On-the-fly learning and sampling of ligand binding by high-throughput molecular simulations. J. Chem. Theory Comput. 2014;10:2064–2069. doi: 10.1021/ct400919u.
81.Dutta S, Selvam B, Shukla D. Distinct binding mechanisms for allosteric sodium ion in cannabinoid receptors. ACS Chem. Neurosci. 2022;13:379–389. doi: 10.1021/acschemneuro.1c00760.
82.Pereira GV, et al. Degradation of complex arabinoxylans by human colonic bacteroidetes. Nat. Commun. 2021;12:459. doi: 10.1038/s41467-020-20737-5.
83.Husic BE, Pande VS. Markov state models: From an art to a science. J. American Chem. Soc. 2018;140:2386–2396. doi: 10.1021/jacs.7b12191.
84.Konovalov KA, Unarta IC, Cao S, Goonetilleke EC, Huang X. Markov state models to study the functional dynamics of proteins in the wake of machine learning. JACS Au. 2021;1:1330–1341. doi: 10.1021/jacsau.1c00254.
85.Shukla D, Hernández CX, Weber JK, Pande VS. Markov state models provide insights into dynamic modulation of protein function. Acc. Chem. Res. 2015;48:414–422. doi: 10.1021/ar5002999.
86.Prinz J-H, et al. Markov models of molecular kinetics: Generation and validation. J. Chem. Phys. 2011;134:174105. doi: 10.1063/1.3565032.
87.Schwantes CR, Pande VS. Improvements in markov state model construction reveal many non-native interactions in the folding of ntl9. J. Chem. Theory Comput. 2013;9:2000–2009. doi: 10.1021/ct300878a.
88.Sultan MM, Kiss G, Pande VS. Towards simple kinetic models of functional dynamics for a kinase subfamily. Nat. Chem. 2018;10:903–909. doi: 10.1038/s41557-018-0077-9.
89.Noé F, Nuske F. A variational approach to modeling slow processes in stochastic dynamical systems. Multiscale Modeling Simulation. 2013;11:635–655. doi: 10.1137/110858616.
90.McGibbon RT, Pande VS. Variational cross-validation of slow dynamical modes in molecular kinetics. J. Chemical Phys. 2015;142:124105. doi: 10.1063/1.4916292.
91.Noé F, Tkatchenko A, Müller K-R, Clementi C. Machine learning for molecular simulation. Ann. Rev. Phys. Chem. 2020;71:361–390. doi: 10.1146/annurev-physchem-042018-052331.
92.Scherer MK, et al. Pyemma 2: A software package for estimation, validation, and analysis of markov models. J. Chem. Theory Comput. 2015;11:5525–5542. doi: 10.1021/acs.jctc.5b00743.
93.McGibbon RT, et al. Mdtraj: A modern open library for the analysis of molecular dynamics trajectories. Biophys. J. 2015;109:1528 – 1532. doi: 10.1016/j.bpj.2015.08.015.
94.Humphrey W, Dalke A, Schulten K. Vmd: visual molecular dynamics. J. Mol. Graphics. 1996;14:33–38. doi: 10.1016/0263-7855(96)00018-5.
95.Roe DR, Cheatham III TE. Ptraj and cpptraj: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 2013;9:3084–3095. doi: 10.1021/ct400341p.
96.Wagner JR, et al. Povme 3.0: software for mapping binding pocket flexibility. J. Chem. Theory. 2017;13:4584–4592. doi: 10.1021/acs.jctc.7b00500.
97.Glielmo A, et al. Unsupervised learning methods for molecular simulation data. Chem. Rev. 2021;121:9722–9758. doi: 10.1021/acs.chemrev.0c01195.
98.Hoffmann M, et al. Deeptime: a python library for machine learning dynamical models from time series data. Machine Learning: Sci. Technol. 2021;3:015009.
99.Hernández CX, Pande VS. MDEntropy: Information-theoretic analyses for molecular dynamics. J. Open Source Softw. 2017;2:427. doi: 10.21105/joss.00427.
100.Kim S, et al. PubChem in 2021: New data content and improved web interfaces. Nucleic Acids Res. 2020;49:D1388–D1395. doi: 10.1093/nar/gkaa971.
101.Trott O, Olson AJ. AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2009;31:455–461. doi: 10.1002/jcc.21334.
102.Eberhardt J, Santos-Martins D, Tillack AF, Forli S. AutoDock vina 1.2.0: New docking methods, expanded force field, and python bindings. J. Chem. Inform. Modeling. 2021;61:3891–3898. doi: 10.1021/acs.jcim.1c00203.
103.Noé F, Schütte C, Vanden-Eijnden E, Reich L, Weikl TR. Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations. Proc. Natl. Acad. Sci. 2009;106:19011–19016. doi: 10.1073/pnas.0905466106.
104.E W, Vanden-Eijnden E. Transition-path theory and path-finding algorithms for the study of rare events. Ann. Rev. Phys. Chem. 2010;61:391–420. doi: 10.1146/annurev.physchem.040808.090412.
105.Shukla D, Meng Y, Roux B, Pande VS. Activation pathway of src kinase reveals intermediate states as targets for drug design. Nat. Commun. 2014;5:3397. doi: 10.1038/ncomms4397.

Associated Data

Supplementary Materials

Data Availability Statement

All the necessary python scripts have been used to generate figures can be obtained from https://github.com/ShuklaGroup/Cannabinoid_activation.git.

Studies

Conditions

Cannabinoids and active compounds

Related research

Abstract

Untitled section

Web Summary

Article notes

Untitled section

Author notes

Introduction

Results and discussions

Structural comparison of CB1 and CB2 crystal structures

Thermodynamic relevance of structural features in cannabinoid receptors’ activation

Allosteric communications between CB1 and CB2 structural features

Activation mechanism for CB1 and CB2

Explaining selectivity of Cannabinoid receptors

Conclusions

Methods

System preparation

Simulation details

Featurization of protein conformational ensemble

Adaptive sampling

Markov state models

Trajectory analysis

Metastable state estimatation using VAMPnets

Data-driven analysis

Allosteric pathway calculation

Docking calculations

Transition path theory

Error calculation

Kinetic Monte Carlo simulation

Reporting summary

Supplementary information

Untitled section

Acknowledgements

Author contributions

Peer review

Peer review information

Data availability

Code availability

Competing interests

Footnotes

Footnote Group

Supplementary information

References

Untitled section

References

Associated Data

Supplementary Materials

Data Availability Statement

Additional materials

Structural comparison of CB₁ and CB₂ crystal structures

Allosteric communications between CB₁ and CB₂ structural features

Activation mechanism for CB₁ and CB₂