Introduction

Cannabinoid receptor 1 (CB₁), which is majorly expressed in the central nervous system (CNS) belongs to the class A G-protein coupled receptor (GPCR) family proteins.^1–4 GPCRs are expressed in the cellular membrane and help transduce chemical signals from the extra-cellular to the intracellular direction with the help of the downstream signaling proteins (G-proteins and β-arrestin).^5–7 In addition, GPCRs are the largest family of drug targets due to their substantial involvement in human pathophysiology and druggability.^8,9 Significant research efforts have been invested in the discovery of drugs targeting CB₁, which helps to maintain homeostasis in neuron signaling and physiological process.^10,11

Initial drug discovery efforts, especially the design of synthetic agonists, were based on modifying the scaffolds of phytocannabinoids (e.g., Δ⁹-Tetrahydrocannabinol, cannabinol) and endocannabinoids (e.g., Anandamide, 2-arachidonoylglycerol) (Figure 1).^12–14 The synthetic molecules, which maintain the aromatic, pyran, and cyclohexenyl ring of the most common psychoactive phytocannabinoid Δ⁹-THC, are known as classical cannabinoids (Figure S1).^15–17 However, the pharmacological potential of these molecules was diminished due to their psychological and physiological side effects (“tetrad” side effect).^18–20 One such example of a synthetic cannabinoid is 1,1-Dimethylheptyl-11-hydroxy-tetrahydrocannabinol (commonly known as HU-210), which is a Schedule I controlled substance in the United States.^21,22

Figure 1.

Apart from the canonical structures of synthetic cannabinoids, molecules with diverse scaffolds were also synthesized through structure-activity studies^23–25. However, these molecules also lacked any pharmacological importance due to psychological side effects.^26,27 Due to the diverse structures and psychological effects, these molecules became unregulated substitutes for traditional illicit substances. ²⁸ These synthetic cannabinoids belong to a class of molecules known as new psychoactive substances(NPS) as these molecules are not scheduled under the Single Convention on Narcotic Drugs (1961) or the Convention on Psychotropic Substances (1971). ^28,29 Synthetic cannabinoids make up the largest category of NPS molecules.^30,31 NPS creates a significant challenge for drug enforcement agencies, as they appeal to drug users seeking “legal highs” to avoid the legal consequences of using traditional drugs and to be undetectable in drug screenings.²⁷

The molecular structures of NPS synthetic cannabinoids consist of four pharmacophore components: linked, linker, core, and tail groups.^27,32 The core usually consists of aromatic scaffolds (e.g., indole, indazole, Carbazole, Benzimidazole) (Figure S2).²⁴ The tail and linker groups are connected to the core. In the tail group, long alkyl chain-like scaffolds are ubiquitous in most NPSs; however, molecules with bulkier cyclic chains (e.g., AB-CHMINACA) are also present.³² Frequently encountered scaffolds in linker groups are methanone, ethanone, carboxamide, and carboxylate ester groups.³³ The linker acts as a bridge between the core and the linked group. In the initial NPS synthetic cannabinoids, the linked group included polyaromatic rings; however, non-cyclic linked groups have also been identified in NPS recently.^24,32 Structural diversity in every component, while maintaining high binding affinity and potency for CB₁ make these molecules easier for drug manufacturers and harder to ban by drug enforcement agencies.^34,34–37

The use of NPS synthetic cannabinoids has been found to cause more physiological side effects than traditional cannabinergic “tetrad” side effects. ³⁸ These side effects include tachy-cardia, drowsiness, dizziness, hypertension, seizures, convulsions, nausea, high blood pressure, and chest pain.^38,39 For instance, Gatch and Forster have shown that the high concentrations of AMB-FUBINACA, the molecule which caused “zoombie outbreak” in New York, induced tremors.^40,41 A recent biochemical study has linked these discriminatory effects with the differential signaling of β-arrestin.³⁹ According to Finlay et al., NPS shows higher β-arrestin signaling compared to the classical cannabinoids, which has also been confirmed by other β-arrestin signaling studies.^39,42 However, mechanistic understanding of these differential downstream signaling effects between NPS and classical cannabinoids is still missing.

Mutagenesis studies have shown that the conserved NPxxY motif of CB₁ have a larger role in downstream β-arrestin signaling than G-protein signaling.^43,44 Recently published MDMB-FUBINACA bound CB₁-β-arrestin-1 complex structure also points out the importance of the unique triad interaction (Y397^7.53-Y294^5.58-T210^3.46) involving NPxxY motif in β-arrestin-1 signaling.⁴⁴ However, structural comparison of the classical cannabinoid (AM841) and NPS (MDMB-FUBINACA) bound active CB₁-G_i complex shows a conformationally similar NPxxY motif (Figure 2). ^45,46 In light of these experimental observations, it can be inferred that higher β-arrestin signaling stems from higher dynamic propensity of triad interaction formation for NPS-bound CB₁. We hypothesized that distinct orthosteric pocket interactions for NPS and Classical Cannabinoids causes differential allosteric modulation of intracellular dynamics that facilitate triad interaction.

Figure 2.

To study these distinct dynamic effects, we compared the (un)binding of the classical cannabinoid (HU-210) and NPS (MDMB-FUBINACA) from the receptor binding site. These molecules have nanomolar affinities. Obtaining the initial pathway of ligand unbinding from unbiased sampling will be computationally expensive. Therefore, a well-tempered metadynamics approach was used to sample the unbinding event, where a time-dependent biased potential is deposited for the faster sampling of the metastable minima along the pathway. ⁴⁷ However, a detailed characterization of the unbinding processes is only possible through the thermodynamics and kinetics estimation of intermediate states. Thus, a transition operator-based approach is needed, which helps to estimate the transition timescale between the states and the stationary density of each state. Estimation from these approaches usually depends on the equilibrium between the local states, which can only be maintained by reversible sampling. For high-affinity ligands like MDMB-FUBINACA and HU-210, reversible sampling is expensive as ligands move from high energy unbound states to lower energy bound states irreversibly. Hence, we implemented a transition operator approach named the transition-based reweighting analysis (TRAM) method, which can tackle this lack of local equilibrium between states by combining unbiased and biased approaches.⁴⁸ TRAM has been used different simulation studies for estimating thermodynamics and kinetics of high energy processes. Characterization of the dimerization kinetics and sodium ion translocation of µ-opioid receptors have been studied using TRAM.^49,50 Furthermore, TRAM have been utilized for thermodynamics and kinetics estimations of ligand (un)binding for kinase and Serine protease.^48,51 To implement TRAM for our study, extensive sampling of the (un)binding process of both ligands was performed using a combination of umbrella sampling and unbiased simulations from the pathway obtained using metadynamics (see Methods section). ⁵² We showed that TRAM can produce consistent kinetic estimation with less unbiased simulation data compared to traditional methods like the Markov state model.⁵³

Based on estimates of thermodynamics and kinetics, it was observed that both NPS and classical cannabinoids have similar unbinding pathways. However, their unbinding mechanisms differ due to the aromatic tail of the MDMB-FUBINACA compared to the alkyl side chain of HU-210. Furthermore, dynamic interaction calculations reveal a major difference with TM7 between NPS and classical cannabinoid. Specifically, the hydroxyl group in the benzopyran moiety of HU-210 forms much stronger polar interactions with S383^7.39 compared to the carbonyl oxygen of the linker group in MDMB-FUBINACA. MD simulations of other classical cannabinoids and NPS molecules bound CB₁ also support these significant interaction differences. The ligand binding effect in intracellular signaling was estimated by measuring the probability of triad formation in the intracellular region. NPS bound CB₁ shows higher probability of forming triad interaction compared to the classical cannbinoids, which supports the experimental observations of high β-arrestin signaling of NPS bound receptors. To validate that the triad formation is indeed caused by the binding pocket interaction differences between the two ligands, allosteric strength binding pocket residues and NPxxY motif was estimated with the deep learning technique, Neural relational inference (NRI).⁵⁴ NRI network revealed that binding pocket residues of NPS bound ensemble have higher allosteric weights for the NPxxY motif compared to classical cannabinoids. These analyses validate our hypothesis that the differential dynamic allosteric control of the NPxxY motif might lead to the β-arrestin signaling for different ligands. This study provides a foundation for additional computational and experimental research to enhance our understanding of the connection between ligand scaffolds and downstream signaling. This knowledge will assist drug enforcement agencies in proactively banning these molecules and inform policies that can protect individuals from the effects of abuse.

Results and Discussion

Metadynamics simulations capture the unbinding paths of NPS and classical cannabinoids

The representative classical cannabinoid and NPS selected for this study are HU-210 and MDMB-FUBINACA.²² Compared to Δ⁹-THC, HU-210 has an extra hydroxyl group in the C-11 position and a 1’,1’-Dimethylheptyl group instead of a pentyl side chain (Figure 1C and S1). MDMB-FUBINACA is a derivative of AB-FUBINACA, which was originally developed by Pfizer (Figure 1C).⁴⁵ These ligands binds to CB₁ receptor with nanomolar affinities (MDMB-FUBINACA K_i : 1.14 nM;⁵⁵ HU-210 K_i : 0.61 nM^14,56).

Metadynamics simulation is a biased sampling method and has been widely used in protein-ligand binding and unbinding studies, as preexisting knowledge of the pathway is not necessary for performing these simulations.^57–60 In metadynamics, a time-dependent biased potential is deposited into the sampling process for the ligand to get out of stable minima at a faster pace.⁶¹ In this work, two replicates of well-tempered metadynamics were performed to capture the unbinding pathway of HU-210 and MDMB-Fubinaca. The commonly used collective variables were selected for metadynamics simulations: (1) z component distance between the center of mass of ligand and residue in the ligand binding pocket (W356^6.48), and (2) Contact number with the ligand heavy atom and α carbon of all binding pocket residues (Equation 4).

The z-component distance was plotted against the RMSD of the ligands from the bound pose, which indicates that ligands follow a similar pathway for each replica (Figures S3A and S3C). It is also observed that the dissociation happens via the opening formed by TM2, TM3, ECL2, and N-terminus for both ligands (Figures 3A, 3B, 3C, and 3D). These observations indicate that the pathway may be the minimum free energy pathway for the ligand unbinding in CB₁. Previous metadynamics binding simulation of another cannabinoid ligand also points to a similar pathway.⁶⁰ Reweighted probability density obtained from metadynamics calculation shows one highly dense region in the pocket, depicting the stability of the bound pose of the ligands (Figures S3C and S3D). However, time dependent external force applied during the metadynamics makes the sampling in the orthogonal direction of the CVs less extensive. Thus, the biased simulation might not sample some protein-ligand interactions that helps to characterize intermediate states. To properly sample intermediate transition states during the underlying unbinding process, extensive unbiased simulations were performed from the initial pathway obtained from the limited sampling of the well-tempered metadynamics (discussed below).

Figure 3.

Comparison of thermodynamics and kinetics estimates from Markov state Model and Transition-based reweighting analysis method

MSM and TRAM are both postprocessing techniques for estimating the kinetics and thermodynamics of underlying physical processes observed in MD simulation. MSM is applied to reversible equilibrium simulations, whereas TRAM estimations can be obtained from multi-ensemble simulations (combination of biased and unbiased simulations). The MSM depends on the local equilibrium between the Markovian states, which is also known as detailed balance (Equation 5). However, reversible local sampling becomes challenging with short parallel trajectories when the free energy difference between two local Markovian states is high. In those cases, reversibility is still assumed by forcing the detailed balance when estimating the transition probability matrix.⁵³ This leads to the incorrect estimation of the unbinding kinetics due to limited sampling from the stable bound state to the high energy unbound states.⁴⁸ Refining the state discretization (i.e., increasing the number of states) may resolve the issue. However, refined state discretization sometimes decreases the statistically significant transition count between all states, decreasing the model certainty. TRAM was shown to solve this problem by combining biased and unbiased simulations (see Methods section). Biased simulations (e.g., replica exchange, umbrella sampling) help to enhance the local sampling, either by increasing the temperature for faster sampling or by fixing collective variables with biased potential to have better sampling in orthogonal directions. It has been shown that compared to MSM, kinetics predicted using TRAM from the combination of biased and unbiased simulations are more aligned with the experiment results.⁴⁸

As unbinding of ligands with high binding affinity (nanomolar) are being studied here, asymmetric transitions might be observed along the pathway. Therefore, we compared the use of MSM and TRAM in estimating the kinetics and thermodynamics of the (un)binding process. For TRAM, unbiased simulation and umbrella sampling were run from the clusters in conformational ensemble obtained from metadynamics (refer to the Methods section for more details). For MSM estimation, only unbiased simulations starting from the metadynamics pathway were considered.

For thermodynamics comparison, standard free energy was estimated for the ligands considering volume correction.⁶² TRAM and MSM predictions of standard binding free energy are within 0.6 kcal/mol of each other for each ligand (Figure 4A). Although the absolute binding free energy differs from the experimentally predicted value by approximately 3 kcal/mol, the relative estimated free energy (ΔΔG) values are also within 0.6 kcal/mol of the experimentally determined values. Therefore, it indicates that with sufficient sampling, both MSM and TRAM converge to the same predictions of relative free energy.

Figure 4.

We also compared the kinetics obtained from the MSM and TRAM. Kinetic measurements were performed with transition path theory (TPT), which uses transition probability matrix from MSM or Tram to estimate mean free passage time between different states (see Methods section). Estimated binding times using TRAM and MSM match perfectly for both ligands. The estimated dissociation times are within one order of magnitude with each other. These observations agree with the previously reported computational research, where experimentally comparable estimation of k_off rates were shown to be more challenging compared to k_on.⁶³

Further analyses were performed to compare these methods in the low unbiased data regime. The difference between the dissociation time of ligands was measured with different amounts of unbiased data. It is observed that even with only 25% of original, unbiased data, TRAM can predict the kinetics within an order magnitude of the kinetics estimated with full dataset (Figure 4D). On the other hand, error in MSM predicted kinetics more rapidly compared to TRAM with lesser amount of unbiased data. A similar trend can be observed for ΔΔG prediction (Figure S4). Therefore, TRAM provides better predictions of thermodynamics and kinetics when less amount of unbiased data.

Unbinding mechanism for new psychoactive substance

Although the binding position of the ligand and the overall binding pathway are similar for both the ligands, extensive biased and unbiased simulation analyzed by TRAM shows a significant difference in the unbinding mechanism of the ligands. To capture the unbinding pathway for MDMB-FUBINACA, we projected the TRAM weighted free energy landscape of the distance between the linked part of the ligand (Leucinate group) and TM5 with respect to the distance between the ligand tail group and TM7 (Figure 5A). The free energy landscape was divided into the non-overlapping intermediate macrostates to obtain better description of the unbinding process. In each macrostate, the contact frequency of ligand with binding pocket residues were calculated along with corresponding contact energies. A metastable minima is observed for macrostate representing the bound pose of the ligand depicting the stability of the ligand (Figure 5A). In the bound pose, the major interactions form between the aromatic (F170^2.57, W279^5.43, F268^ECL2) and hydrophobic (L193^3.29) residues of the binding pocket (Figures 6A, 6B, and S5B).

Figure 5.

Figure 6.

The free energy landscape shows two probable mechanisms for MDMB-FUBINACA unbinding from the bound pose. The two pathways are differentiated by whether the linked or tail part of MDMB-FUBINACA dissociates first. One of the pathways, aromatic tail part of MDMB-FUBINACA moves away from TM5 and form interaction with aromatic residues in TM2 (F170^2.57 and F174^2.61) (Figures 5A, 6A, 6B, and S5C). This leads to the formation of intermediate metastable states, which we characterize as macrostate Intermediate state 1 (I1). This metastable minimum observed from I1 macrostate might be unique to the FUBINACA family of NPS synthetic cannabinoids as this family has the aromatic ring in tail group, unlike the long alkyl chain in other common synthetic cannabinoids. Along with the aromatic residues of TM2, major interaction with F268^ECL2 is maintained in macrostate I1 (Figures 6A and 6B). K-L divergence analysis between bound and I1 macrostates show that only minor changes in the binding pocket residues, especially in TM2 are needed to accommodate MDMB-FUBINACA in this conformational state (Figure 6C). The interconversion timescale (MFPT) between the macrostates were obtained from the transition path theory. MFPT calculations show that both the timescales are similar with slightly higher timescales for the bound pose compared to I1 transition (20.6 ± 2.3µs) (Figure 6B). In this pathway, ligand moves from I1 metastable state to space between N-terminus, TM2, TM3, and ECL2 before dissociating from the receptor (Figure 6B). This region between the unbinding ensemble have been characterized as macrostate I3 (Figure 5A). Contact analysis show significant drop in ligand residue contacts with only aromatic residues in TM2, TM3 and ECL2 forming dominant interactions (Figures 6A and S5D). K-L divergence shows that ligand positioning in this particular regions causes relatively higher divergence on TM2 compared to I1 (Figure 6C). Kinetically, the transition from I1 to I3 (33.7 ± 3.1µs) is much slower compared to reverse transition (0.8 ± 0.0µs), validating the higher stability of the I1 compared to I3 macrostate (Figure 6B). According to the TPT analysis, breaking the aromatic interactions for complete dissociation of MDMB-FUBINACA requires ∼ 371.9 ± 40.2µs, making it the slowest step in this pathway (Figure 6B).

In the other possible unbinding pathway, orientation of MDMB-FUBINACA in the pocket doesnot change compared to the bound pose. The linked part of the ligand moves to space between N-terminus, TM2, TM3, and ECL2 (Figure 6B). We label this macrostate as I2. In this state, we observe stable polar interaction with K376^7.32 and hydrophobic interactions with aromatic and other hydrophobic residues (F177^2.64, F268^ECL2, P269^ECL2) (Figures 6A and S5A). However, free energy of this macrostate is higher than the bound pose, depicting higher entropic cost associated with this state. This can be shown by the higher intrastate RMSD of I3 compared to the bound pose (Figure S6). Transition timescale from the bound pose to the I2 (136.7 ± 13.9µs) is one order of magnitude higher compared to the reverse transition (9.8 ± 13.9µs) (Figure 6B). K-L divergence analysis also shows higher divergence in the extracellular region of TM2 and N-terminus compared to bound pose (Figure 6C). Dissociation of MDMB-FUBINACA from I2 to the bulk is faster compared to dissociation from I3 (246.6 ± 26.0µs) (Figure 6B). However, the overall kinetic barrier for dissociation from the binding pose for the both unbinding mechanisms are relatively similar.

Unbinding mechanism for classical cannabinoid

For capturing the classical cannabinoid (un)binding mechanism, distances from the two terminal scaffolds (cyclohexenyl and alkyl chain) to TM5 and TM7 were measured similar to the NPS (Figure 5B). The free energy landscape of the unbinding of the HU-210 shows the differences in the mechanism from MDMB-FUBINACA. Similar to MDMB-FUBINACA, the HU-210 unbinding land-scapes were also divided into non-overlapping macrostates. Macrostate representing HU-210 bound pose shows a metastable energy minimum. Comparing the bound macrostate interactions of MDMB-FUBINACA, classical cannabinoid HU-210 shows higher interactions with TM7 residues (S383^7.39, F379^7.35) (Figures 6A, 7A, S5B and S7B). Previous experimentally determined structures of classical cannabinoid bound CB₁ have pointed out these conserved polar interactions of the hydroxyl group at C-1 position with S383^7.39.^46,64 Although MDMB-FUBINACA also maintain this polar interaction with carboxylic oxygen, the interaction energy for the HU-210 is much higher, depicting the importance of this conserved residue in stabilizing classical cannabinoids (Figure 6B). Mutagenesis studies also support this difference in interaction with S383^7.39 between the classical cannabinoids with hydroxyl group (HU-210) and CB₁ ligands, which have carboxylic oxygen in the equivalent position (WIN-55,212-2).^65,66 Alanine mutation of S383^7.39 have shown to decrease the lig- and affinity and downstream efficacy of classical cannabinoids by orders of magnitude, while having minimal effect on WIN-55,212-2, which have carboxylic oxygen in the linked part as MDMB-FUBINACA.^65,66 Other major interactions (F170^2.57 and F268^ECL2) in the bound pose are common between the two ligands (Figures 6B and 7B).

A relatively weaker metastable state is observed when the ligand moves relatively deeper (closer to TM5) inside the binding pocket. The flexible alkyl chain of HU-210 allows the ligand to have this deeper position (Figures 7A and 7B). Protein-ligand interaction analysis in the macrostate representing this region (I1) show that hydroxyl group at C-11 forms a major polar interaction with H178^2.65 (Figure S7A). The bound and I1 macrostates are kinetically close, as indicated by the rapid interconversion between these states (Figure 7B). K-L divergence between the two states show the highest divergence in extracellular TM2 and TM7, where major interaction switch have happened (Figure 7C).

Figure 7.

Contrasting to MDMB-FUBINACA, only one pathway was discovered with classical cannabinoid cyclic scaffold departing from the receptor first. Major interactions that break when the ligand moves out of the binding pose to macrostate I2 is the polar interaction with S383^7.39 and hydrophobic interaction with aromatic F170^2.57 (Figures 7A, 7B, and S7C). Breaking of these bonds leads to larger kinetic barrier of approximately 39.5 ± 1.4µs (Figure 7B). In this macrostate, the HU-210 forms major interactions with aromatic residues F268^ECL2 and F177^2.64 and polar interactions with S173^2.60 and D176^2.63 (Figure 7B). From this pose, HU-210 either dissociates from the receptor or obtain another relatively weak stabilized state (I3) in the receptor. In I3, the alkyl chain of the ligand is flipped in the pocket and stabilized by aromatic residues in TM2, TM3 and ECL2 (Figure 7B). This transition from I2 to I3 (528.1 ± 32.8µs) kinetically much slower compared to the reverse transition (6.9 ± 0.0µs) (Figure 7B). From both I2 and I3 macrostates, the ligand can dissociate from the pocket and mean free passage time for these transitions appear to be in the milisecond timescale, which is one order magnitude higher compared to the MDMB-FUBINACA unbinding (Figure 6B). This phenomena supports the relatively high affinity of the classic cannabinoid HU-210 compared to the NPS MDMB-FUBINACA.

Allosterically controlled distinct downstream signaling between new psychoactive substances and classical cannabinoids

As discussed in the previous section, major interaction differences between are observed in TM7. The universality of this observation is measured by performing unbiased MD simulation (1 µs each) of other NPS (AMB-FUBINACA, 5F-AMP, CUMYL-FUBINACA) and classical cannabinoids (AMG-41, JWH-133, O-1317) bound CB₁ (Figure S8A). Average distance of carbonyl oxygen of NPS molecules’ linker group from S383^7.39 is compared to equivalent distance of hydroxyl group of classical cannbinoids’ benzopyran ring. Larger mean distance in case of all NPS-bound CB₁ supports the universality of the weaker interaction between TM7 and NPS molecules (Figure S8B). This variation in binding pocket pocket interactions might lead to differential allosteric control of the intracellular dynamics that facilitate triad interaction (Y397^7.53-Y294^5.58-T210^3.46) important for β-arrestin binding.

We adopted a data-driven deep learning network known as Neural relational inference (NRI) to validate our hypothesis of allosteric control. NRI network has an architecture of variational autoencoder. The encoder part of the network predicts the interactions between the residues from the trajectory dynamics, and the decoder predicts the trajectories from the interaction. With this network, we try to produce alpha carbon coordinates at t + τ from the coordinates at time t. In the process of regenerating the future coordinates, the latent space of the network learns the dynamic interactions between different residues in the protein. These interactions are calculated from the estimated posterior probability q(z_ij|x). In this work, we trained the network with the NPS, and classical cannabinoid-bound unbiased trajectories (Method Section). Here, we compared the allosteric interaction weights between the binding pocket and the NPxxY motif which involves in triad interaction formation. Results show that each binding pocket residue in NPS bound ensemble shows higher allosteric weights with the NPxxY motif, indicating larger dynamic interactions between the NPxxY motif and binding pocket residues(Figure S9).

The probability of triad formation was estimated to observe the effect of the difference in allosteric control. TRAM weighted probability calculation showed that NPS-bound CB₁ have the higher probability of triad formation (Figure 8A). Comparison of the pairwise interaction of the triad residues shows that interaction between Y397^7.53-T210^3.46 is relatively more stable in case of NPS-bound CB₁, while other two interactions have similar behavior for both systems (Figures S10A, S10B, and S10C). Therefore, higher interaction between Y397^7.53 and T210^3.46 in NPS-bound receptor causes the triad interaction to be more probable.

Figure 8.

Furthermore, we also compared TM6 movement for both ligand bound ensemble which is another activation metric involved in both G-protein and β-arrestin binding. Comparison of TM6 distance from the DRY motif of TM3 shows similar distribution for HU-210 and MDMB-FUBINACA (Figure 8B). Therefore, difference in NPxxY dynamics due to the differential allosteric control validates the dynamic effect of the NPS ligand on the NPxxY motif.

Conclusions

Synthetic cannabinoids were designed as a potential therapeutics to target cannabinoid receptors. However, major side effects of these ligands diminish their therapeutic potential. Although both classical cannabinoids and NPS synthetic cannabinoids have been abused as recreational drugs, later poses larger threats for the society due the chemical diversity of the NPS structures makes them harder to control from being abused. Furthermore, physiological studies have shown NPS targeting cannabinoid receptors lead to the dangerous physiological effects compared to “tetrad” side effects associated with classical cannabinoids. Previous studies have relate these side effects with the higher downstream β-arrestin signaling of NPS. Mutagenesis studies have shown NPxxY motif have larger role to play in β-arrestin signaling. In this work, we proposed that NPS and classical cannabinoid have distinct allosteric control on NPxxY motif when bound to orthosteric pocket of CB₁. In this hypothesis driven study, we compared ligand-protein interactions of NPS MDMB-FUBINACA and classical cannabinoid HU-210 for CB₁ by studying their unbinding mechanism and downstream signaling.

As both ligands are stable binders with nanomolar affinity, well-tempered metadynamics simulations were performed to obtain the initial unbinding pathway. These simulations were able to find similar pathway via the opening formed by TM2, TM3, N-terminus, and ECL2 which matches with previous metadynamics binding simulations for cannabinoid receptors. For the proper characterization of intermediate states, the unbinding processes were further extensively sampled using unbiased simulation and umbrella sampling.

Effectiveness of the post processing techniques TRAM and MSM were compared in predicting the ligands with binding affinity and kinetics. MSM predicts the kinetics and thermo-dynamics from the eigendecomposition of the transition probability matrix. MSM assumes that the local equilibrium is maintained between the Markovian states. However, with limited sampling, this creterion may not valid between local high and low energy states. TRAM tries to solve this issue by combining biased and unbiased simulation, where biased simulations enhances the local sampling to maintain the equilibrium. We observed that with sufficient data both methods performed in similar way in estimating the standard binding free energy. The relative free energy estimated by both methods match with the experimental result within 0.6 kcal/mol. With lesser amount of unbiased data, TRAM predictions of kinetics and thermodynamics remain more consistent than the MSM as the biased simulations help to maintain local equilibrium.

TRAM estimated thermodynamics helped to decipher the differences between the un-binding of NPS MDMB-FUBINACA and classical cannabinoid HU-210. First, for MDMB-FUBINACA, a larger conformational change is observed within the pocket. A metastable intermediate state is observed when the aromatic tail of FUBINACA flip inside the pocket and reorient itself close to the aromatic residues of TM2. It was observed that both linked part and tail part of the ligands can lead the dissociation of the ligand from the receptor. Second, for HU-210, conserved cyclic group leads the dissociation from the receptor. It supports previous simulation where the alkyl side chain of the ligand binds to the receptor first.¹⁷ Third, aromatic residues in the pocket (F268^ECL2, F170^2.57) form major interactions with both classical cannabinoid and NPS molecules. Major differences in protein-ligand interactions were observed in TM7. Stronger interactions were observed for the classical cannabinoid with TM7, especially polar interaction with S383^7.39 and hydrophobic interaction with F379^7.35 compared to NPS molecules.

Finally, we showed that the variation in binding pocket interaction leads to the distinct dynamic allosteric communications in the intracellular region. Allosteric communication strength was measured by the variational autoencoder (NRI). NRI network learns the dynamic interactions between residues in the latent space by learning to reconstruct the dynamics. Dynamic allostery measured by the posterior probability of VAE shows higher allosteric weights from the binding pocket residues to the NPxxY motif region for NPS molecules. Higher allosteric control for NPS-bound CB₁ increases the probability of triad interaction formation, resulting higher β-arrestin signaling for NPS. Therefore, these data driven computational study helps us to distinguish between the receptor-protein interaction, unbinding mechanism and downstream signaling NPS compared to other classical cannabinoids.

Methods

System Preparation

For NPS unbinding simulation, G-protein bound active structure (PDB: 6N4B⁴⁵) was selected as the initial structure. G-protein subunits and Non-Protein residues other than orthosteric ligand MDMB-Fubinaca were removed from the PDB structure file. Missing residues in ICL3 (21 residues, 314-334) and ECL2 (6 residues, 258-263) were modeled sequentially using Remodel protocol of Rosetta loop modeling.^67,68 In each step, the remodeled structure with least energy was further refined using kinematic closure protocol. ⁶⁹

Starting 108 residues from CB₁ N-terminus were also missing from the cryo-EM structure. However, it is not feasible to model the entire N-terminus because of two following reasons.⁶⁷ First, a proper template is not available for modeling N-terminus regions as most of the class A GPCRs do not contain large N-terminus.⁷⁰ Second, it is challenging to model these large numbers of residues accurately with template-free ab initio modeling because of the combinatorial expansions of conformational space. Therefore, the closest 20 residues (89-108) were modeled as membrane proximal regions of the N-terminus were shown to be important for CB₁ signaling by allosterically modulating ligand affinity. ⁷¹ Furthermore, Δ89CB₁ (CB₁ with first 88 residues truncated in N-terminus) have similar ligand binding affinity compared to CB₁ with full sequence.⁷² Modeling of this membrane proximal region was also performed Remodel protocol of Rosetta loop modeling. A distance constraint is added during this modeling step between C98^N−term and C107^N−term to create the disulfide bond between the residues.^71,73

The “Ligand Reader & Modeler” module of CHARMM-GUI was used for ligand (e.g., MDMB-Fubinaca) parameterization using CHARMM General Force Field (CGenFF).⁷⁴ The ligand bound receptor was embedded in the bilayer membrane and salt solution (extracellular and intracellular region) using CHARMM-GUI.⁷⁵ As CB₁ is majorly expressed in central nervous system, an average brain membrane composition of asymmetric complex membrane was selected. The membrane composition was obtained from Ingólfsson et al. and proportionally downsized according to our system (Table S1).⁷⁶ 150 mM NaCl salt solution with TIP3P water model was used in the extracellular and intracellular regions.⁷⁷ CHARMM36m forcefield was used to parameterize the protein, lipid, water, and ions.⁷⁸

For building the classical cannabinoid system, the modeled PDB structure was used. In this case, the orthosteric NPS ligand was removed, and a classical cannabinoid HU-210 was docked into the orthosteric pocket using Autodock Vina. ⁷⁹ The docked bound pose was selected based on best overalled structure of HU-210 to the experimentally determined crystal structure of another bound classical cannabinoid (Ligand: AM841, PDB:6KPG⁴⁶) (Figure S11). The classical cannabinoid-bound system was built with identical complex membrane composition, salt concentration, and forcefield with NPS bound system.

Other classical cannabinoids (AMG-41, JWH-133, and O-1317) and NPS (AMB-FUBINACA, CUMYL-FUBINACA, 5F-AMP) bound systems were also set up. These ligands are docked into the orthosteric pocket. Best docking poses were selected based on optimizing the distance between the hydroxyl group of classical cannabinoid (linker oxygen for NPS) to S383^7.39 and the furtherest tail atom distance to W279^5.43. These systems also have identical complex membrane composition, salt concentration, and forcefield with previously described systems.

System minimization and Equilibration

All ligand-bound systems are minimized and equilibrated before the production run. Ten thousand minimization steps were performed with the conjugate gradient method. Six sequential equilibration steps were carried out to stabilize the systems at 300 K temperature and 1 atm pressure for the production simulations. The systems were heated to 300 K in the NVT ensemble in the initial two stages. Each of these steps was performed for 250 ps. Langevin dynamics was used to control the temperature with additional damping and random force. The damping coefficient for the damping force was set as 1/picosecond. The Langevin dynamics was turned off for the hydrogen atoms. All the bonded hydrogen atoms were constrained with the SHAKE algorithm with default parameters of NAMD.⁸⁰ The integration time step for these two NVT ensemble equilibration was one femtosecond (fs). Harmonic constraints were used for fixing the temperature coupling of the protein residues with constraint scaling term set to 10 in the first NVT ensemble equilibration, followed by 5. Temperature coupling of lipid molecules was also restrained with harmonic force with a force constant equal to 5. CHARMM-GUI Selected Dihedral and Improper bonds were also restrained with an extra bonded term with a force constant of 500 in the first NVT ensemble equilibration, followed by 200. The non-bonded cutoff distance for the van der Waals interactions was set to be 12 Å with a switch distance of 10 Å, at which a switching function is turned on to truncate van der Waals interactions at the cutoff distance. Non-bonded interactions for three consecutively bonded atoms were excluded. The particle mesh Ewald method was implemented for electrostatics calculation with grid size 1 Å.⁸¹

The next four equilibrations were performed in the NPT ensemble, where pressure was fixed to 1 atm with Langevin piston pressure control. The barostat oscillation period was set to 100 fs with a damping time 50 fs. These four NPT ensemble equilibrations were performed for 250, 500, 500, and 500 ps, respectively. The integration timestep for these equilibration steps was increased to two fs. The constraint for temperature coupling on the protein residues was decreased gradually with the constraint scaling term for the four NPT ensemble simulations set to 2.5, 1.0, 0.5, and 0.1, respectively. Similarly, the restraints on temperature coupling on the lipid molecules were also decreased gradually with force constant for the four NPT ensemble simulations set to 2, 1, 0.2, and 0.0, respectively. Furthermore, the restraints on the CHARMM-GUI Selected dihedral and improper bonds were decreased gradually with force constant for the four steps set to 100, 100, 50, and 0.0, respectively.

Well-tempered metadynamics

Well-tempered metadynamics was implemented for finding unbinding pathways.^47,82 Simulations were performed with the Collective variables module (Colvars) of NAMDv2.14.^83,84 In metadynamics, a history-dependent biasing potential (V_meta(S, t)) is added to the Hamiltonian of the MD simulation, which discourages the system from revisiting configurations that have already been sampled.^61,85 The V_meta(S, t) is a sum of Gaussians deposited along the system trajectory in the CVs space (S = (S₁(r), S₂(r)..S_d(r)) as shown in Equation 1, where W, σ, τ are Gaussian height, width, and deposition time step, respectively. With a sufficiently long simulation, the bias potential estimates the underlying free energy along the CVs.^61,82 The Well-tempered metadynamics was introduced to increase the convergence of bias potential by decreasing the Gaussian height with time (Equation 2).⁶¹ In the Equation 2, ω is the bias deposition rate and ωτ is equivalent to constant gaussian height for well-tempered metadynamics. The free energy is estimated from the bias potential using Equation 3, where the ΔT is an user defined parameter.

In this work, we selected two collective variables (CVs) for NPS (MDMB-FUBINACA) and classical cannabinoid (HU-210) unbinding according to Mahinthichaichan et al..⁵⁹ The first collective variable is the z-direction distance from the converged toggle switch to the center of mass of all heavy carbon atoms of ligands. The second collective variable is the coordination number (CN) as defined in Equation 4 where d_ij represents the distance from ith atom of the ligand and the alpha carbon of jth residue in the binding pocket. Selected binding pocket residues for CN calculation are shown in the Supplementary table. The ΔT, gaussian height (ωτ), deposition time step (τ) are selected as 4200K, 0.4 kcal/mol and 100 ps, respectively. The gaussian width for the two CVs are set to be 0.5 and 0.1, respectively.

Umbrella Sampling and Unbiased Sampling

Umbrella sampling was performed along ligand distance from TM5 to capture the unbinding process (Figures S12A and S12B). ⁵² The unbinding pathway obtained from the metadynamics was clustered into 300 bins by dividing the selected distances from 5 to 35 Å. The center of each bin was used as the center of each window for umbrella sampling. Five independent structures were selected from each cluster to simulate five independent umbrella runs in each umbrella window. If a cluster does not contain any structure, starting structures for that window were selected from the closest clusters. A constant harmonic biased potential of 10 kcal/mol is used for each window. openMM v7.8 MD engine was used to run the umbrella sampling runs.⁸⁶ The temperature and pressure of the systems are controlled at 300K and 1 atm by the Langevin thermostats and Monte Carlo barostats. The integration timestep was chosen to be two fs. Movements of the containing Hydrogen atom were constainted using HBonds commands with SHAKE (or SETTLE for water) algorithm. The cutoff distance for Non-bonded interaction other than electrostatic interaction was set to 12 Å, with a switching potential at 10 Å to make the potential to zero smoothly at the cutoff. The particle weld method was used to calculate the long-range electrostatic. Each simulation was run for 20 ns.

Identical starting structures and simulation conditions (Thermostat, barostat, cutoff, electrostatic calculation method, integration timestep, and constraints on Hydrogen bond) were selected for unbiased simulations. openMMv7.8 simulation software was used to run simulations. Each trajectory was run for 100 ns. All the simulations were performed on the distributive computing facility folding@home. ⁸⁷

Markov State Model

Markov state model (MSM) is used to estimate the thermodynamics and kinetics from the unbiased simulation.^53,88 The prevalence of MSM as a post-processing technique for MD simulations was due to its reliance on only local equilibration of MD trajectories to predict the global equilibrium properties.^89,90 Hence, MSM can combine information from distinct short trajectories, which can only attain the local equilibrium.^91–93 This attribute of MSM helps to utilize the power of parallel supercomputers to run parallel short simulations to sample a process with a high energy barrier.

In MSM, the conformational ensemble is discretized into states, assuming that the interstate transition will be Markovian. According to the Markovian property, interstate transitions will only depend on the current state, not on the past state. In order to uphold the principle of Markovian transition, the structures of the discretized states must be conformationally and kinetically similar, allowing for rapid interconversion compared to interstate transitions.⁸⁸ These faster intrastate conversions ensure that Markovian states lose their memory of the previous interstate transition after a lag time (τ).

Another necessary criterion for MSM is the reversibility in the interstate transition, which is also known as the detailed balance.⁹⁴ The detailed balance assures the local equilibrium between the state population, as shown in Equation 5, where µ_i is the stationary density of ith Markovian states and P_ij represents the probability of transition between ith to jth at lag time. These transition probabilities form a transition matrix with N × N dimension where N is the total number of states. From the eigendecomposition of this transition matrix, biomolecular processes are discretized into orthogonal processes (orthogonal eigenvectors). The timescale of these processes is calculated based on eigenvalues. The highest eigenvalue is one; the corresponding left eigenvector of the transition matrix determines the stationary density. Eigenvectors with eigenvalues close to one denote the slow processes.

We can estimate the thermodynamics and kinetics of processes from the transition matrix, whose elements are transition probabilities. These probabilities are usually calculated from the transition counts (C_ij) using the maximum likelihood approach, where the likelihood function (Equation 6) is maximized with two constraints: detailed balance and normalization constraints (Σ_j p_ij = 1). The constrained optimization problem is usually solved iteratively, as shown in Prinz et al..⁵³

The following steps are taken for the practical implementation of the MSM from the MD data.^4,17,95–97

1. Each frame obtained from the MD simulation was featurized using features important for capturing the conformational ensemble. In this case, the unbinding process for each ligand was featurized using distances that characterize the ligand distances to the binding pocket and binding pocket conformational change. Specifically, all heavy atom distances from each the Cα carbon atom of all binding pocket residues were calculated (Figure S13). Additionally, all possible combinations of Cα carbon atom distances between all the binding pocket residues were included to capture the binding pocket motion. Feature calculations were performed with the python library MDTraj v1.9.8.⁹⁸ The total number features selected for MSM building of MDMB-FUBINACA and HU-210 are 297 and 288, respectively.

2. Dimensionality reduction was performed using time-independent component analysis (TICA).^99,100 We found the orthogonal projections (time-independent components) with TICA, which are linear combinations of the slowest features. In tIC space, two spatially close points are kinetically close. The lagtime selected for tiC building was 5ns.

3. Clustering was performed on the tICs using k-means clustering algorithms to discretize the space into Markovian states.

4. Lag time for the MSM was calculated by estimating the shortest time at which the timescale of the slowest processes has converged to a particular value (Figures S14A and S14B).

5. To optimize MSM based on the cluster numbers and tIC components on which clustering is performed, we calculated the VAMP-2 score from the MSM, where VAMP stands for Variational Approach for Markov Processes (Figure S15). ¹⁰¹ For a reversible MSM, this score represents the summation of the square of the k slowest eigenvalues, where k is a hyperparameter. Closer the eigenvalue to 1, the corresponding eigenvector captures a slower process. Therefore, we optimize the MSM by maximizing the VAMP-2 score.

6. To validate the Markovian property of our optimized models, Chapman–Kolmogorov test (c-k) test was performed (Figure S16). C-K test states that for a Markov model kth power of P (τ) needs to be equal transition probability matrix determined at kτ time (P (τ)^k ≈ P (kτ)). We showed that differences between the elements of transition probability matrix at higher lag times remain relatively small.

Dimensionality reduction, clustering Markov state model building, and VAMP-2 calculations are performed with the pyEMMA v2.5.6 library. ¹⁰² The optimized MSM for MDMB-FUBINACA unbinding simulations were built with 700 clusters, 7 tiCs and 35 ns of lag time. For HU-210, optimized MSM were build with 800 clusters, 6 tiCs and 35 ns of lag time.

Transition-based reweighting analysis method

Markov State Models have been extensively used to investigate the protein-ligand binding process.^{17,62,103–107} However, these studies were mainly performed for ligands with high off-rates which could be sampled using the unbiased trajectories. For ligand with low off rates, the use of reversible transition matrix would yield incorrect estimates of unbinding kinetics. Therefore, we use The Transition-based reweighting analysis (TRAM)^48,108 method to accurately estimate the unbinding kinetics of new psychoactive substances. TRAM is a thermodynamics and kinetics estimator method, which, unlike MSM, can combine unbiased and biased simulation data to estimate thermodynamics and kinetics. TRAM utilizes the advantages of the local equilibrium approximation of MSM and the benefits of biased simulations to enforce local equilibrium in interstate transitions where it is difficult to attain.

As the simulations are obtained from multiple ensembles (biased and unbiased), it is paramount to classify the MD frames (or the conformations) based on which ensemble it belongs to. Each ensemble represents simulations that are performed with identical Hamiltonian energy functions. Therefore, unbiased simulations are considered as one ensemble, whereas, in umbrella sampling, each biasing window is considered a single ensemble.

Like MSM, in TRAM, the conformational space is also discretized into non-overlapping states. The interstate transitions should follow the following relationship shown in the Equation 7, where is the local free energy of the ith state and kth ensemble. The term is proportional to the stationary density (µ(x)) of state S_i in ensemble k. The µ(x) of each conformation (x) of S_i is weighted with negative exponential of bias energy (b^k(x)) on x in ensemble (Equation 8).

To obtain kinetics and thermodynamics information from TRAM, we have to derive interstate transitions and the stationary density of the entire ensemble (µ(x)), where both the terms follow normalization constraint (Equation 9 and 10). Therefore, there are m²K + X unknown variables. Therefore, to solve these unknown variables, the maximum likelihood approach has been considered, where the likelihood function is defined as Equation 11, which is the combination of the likelihood function of MSM and local equilibrium. This maximum likelihood problem was subjected to the constraints of Equation 7, 9 and 10.

Wu et al. showed that the solution of these maximum-likelihood problem can be turned into system of non-linear algebraic equations (Equation 12, 13 and 14), where is count of the interstate transitions between state S_i and S_j in ensemble k. This system of equations are solved iteratively to estimate and , which provides the prediction of and µ(x) (Equation 15 and 16).

We used the python package pyEMMA v2.5.6 for the practical implementation of TRAM.¹⁰² For calculating transition counts in ensemble , the lagtime of 15 ns was chosen. In the implementation, we need to preprocess each trajectory into three arrays.

1. One of the arrays represents the spatial discretization of each trajectory frame, where each frame belongs to a particular state. Therefore, each element can take values from 0 to m-1 (m is the cluster number). Before the discretization of the space, time-independent component analysis was performed on the biased and unbiased data separately. The number of tIC components for each system was selected based on the number of the tIC components of optimized MSM. Each frame from the unbiased simulation is represented by the unbiased tICs, concatenated with its feature projections on the biased tICs. Similarly, each frame from the biased simulation is represented by its feature projections on the unbiased tic, concatenated with biased tic projection. Therefore, NPS unbinding simulations have 14 tICs, whereas classical cannabinoid un-binding simulations have 12 tICs. The number of clusters is also obtained from the optimized MSMs.

2. Another array represents the corresponding ensemble to which each trajectory frame belongs. There are 300 windows for the umbrella sampling. Therefore, there are 301 ensembles, as the unbiased simulations represent a separate ensemble.

3. Third array represents the corresponding bias potential (b_k(x)) a particular frame feels if it were to be in a particular ensemble. For umbrella sampling, the biased potential is represented as Equation 17, where c_k is selected to be 10 kcal/mol and y_k is the center of each umbrella window.

   

Transition Path Theory

Transition path theory (TPT) analysis is applied to calculate the transition pathway and timescale between different macrostates, representing different configurational spaces in the unbinding process.^109,110 In this work, we define macrostates as a collection of Markovian states present in the area of interest in the unbinding free energy landscape. An essential concept of transition path theory is the committer probability , which is defined as the probability of any Markovian state reaching the final metastable state before it returns to the initial state. Therefore, the Markovian states present in metastable state B has a committer probability of 1. It has been shown that committer probability follows the following system of linear equation as shown in Equation 18, where P_ik is the transition probability between state S_i and S_j as discussed in the previous section.

In this work, the quantity of interest from TPT is the timescale (or rate) between the metastable state transitions as shown in Equation 19, where π_i is the stationary probability of state S_i. TPT calculations were performed by PyEMMA v2.5.6. ¹⁰²

K-L divergence analysis

Kullback–Leibler divergence (K-L divergence) analysis was performed to show the structural differences in protein conformations in different macrostates^4,111. In this study, this technique was used in this work to calculate the difference in the pairwise inverse distance distribution in macrostates. Although K-L divergence is an asymmetric measurement, for this study, we used a symmetric version of the K-L divergence by taking the average of the K-L divergence between two residues. Per residue contribution of K-L divergence is calculated by taking the sum of all the pairwise distances containing that residue. This analysis was performed by in-house Python code.

Trajectory Analysis

Python package GetContacts is used to perform the contact calculation.¹¹² Linear interaction energy analysis was performed to calculate the interaction energy between ligand and receptor using AMBERTools CPPTraj v18.01.^113,114 Trajectory visualization and figure generation are performed with VMD v1.9.3.¹¹⁵

Deep Learning Network for Allosteric Prediction

Neural relational inference (NRI) network was implemented to predict allosteric dependence between the residues in the different parts of the receptors.^54,116 This network is a Variational autoencoder (VAE) comprising encoding and decoding parts.¹¹⁷ The encoder (q_ϕ(z|x)) takes the input Cα coordinates of protein conformations at time t (x_t) and tries to learn the interactions between two residues (z_ij) in the protein as a latent space. The decoder (p_θ(x|z)) network try to regenerate the protein conformation at time t + τ (x_t+τ). Similar to other VAE, the learning process maximizes the evidence lower bound (ELBO) as shown in Equation 20, where p_θ(z) represents the prior distribution for z. Here, the prior distribution is selected as default presented in the original paper, where it is represented as a categorical distribution with K = 4 (P₁ = 0.91, P₂ = 0.03, P₃ = 0.03, P₄ = 0.03).

As shown in Equation, the ELBO consists of two terms. In the first term, further mathematical derivations can show that the first term can be represented as the the reconstruction error (Equation 21), where σ² is variance of the distribution, a user defined parameter.

The second term is also called regularization term which is the the K-L divergence between estimated posterior (q_ϕ(z|x)) and prior distribution (p_θ(z)) (Equation 22). As the prior distribution is a categorical distribution, the K-L divergence becomes entropy of the posterior distribution. We obtained the code for the NRI network from the GitHub implementation and kept most of hyperparameters as default for our training, except for decreasing the hidden layer size to 64.¹¹⁸ From each unbinding simulations, 10 unbiased trajectories were selected where the ligand remain in the bound pose. Each trajectory has a length of 100 ns. Both cases, the τ was selected to be 5 ns. The allosteric weights (posterior probability) were obtained from the validation data (2 trajectories), where training was performed with remaining 8 trajectories (Figure S18). This procedure was repeated three times, where training and validation data were selected randomly.

Standard Binding Free Energy Calculations

To calculate the standard binding free energy from simulation, we adopted a procedure described in Buch et al..⁶² In this procedure, a volumetric correction term is added to the PMF to calculate the final binding free energy (Equation 23). The volumetric correction term is used to predict the free energy at the standard condition (1M) as shown in Equation 24, where V_o corresponds to the volume of a molecule should occupy at the standard condition (1661 ų) and V_u is the volume of the unbounded state in the simulation box. The expression for the PMF contribution of the free energy is shown is in Equation 25, where the denominator of the equation can be represented as .

Therefore, the final derivation of ΔG is shown in Equation 26

In this work, to estimate the free energy (ΔG) x, y, z component of the ligand center of mass is calculated compared to the center of mass of the alpha carbons of binding pocket residues. The three dimensional space was descretized into 25 × 25 × 50 bins and each bin is weighted using TRAM calculated probability density. Depth of the pmf (ΔW) was calculated by averaging the pmf of the 100 bins with highest pmf in the bulk. To evaluate the weighted binding volume , we selected the bins with pmf less than 1 kcal/mol.

Acknowledgements

D.S. acknowledges support from NIGMS MIRA award R35GM-142745 and NSF Early CAREER Award (MCB-1845606). S.D. and D.S. thank folding@home donors for providing computational resources for the study.

Conflict of interest

The authors declare no conflict of interest.

Code and data availability

Unbinding simulation trajectories and topology files that have been used for the analysis can be obtained from https://uofi.box.com/s/gcpafgj5l184c3woqk9m92dk468r58mg. Python scripts and necessary files to generate the figures is provided in the github repository. https://github.com/ShuklaGroup/Dutta_Shukla_Cannabinoid_2023a.git