Mechanistic origin of partial agonism of tetrahydrocannabinol for cannabinoid receptors
Department of Chemical and Biomolecular Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Department of Comparative Biosciences, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Department of Bioengineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Cancer Center at Illinois, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
National Center for Supercomputing Applications, University of Illinois, Urbana, Illinois, USA
Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
NIH Center for Macromolecular Modeling and Bioinformatics, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA
Abstract
Cannabinoid receptor 1 (CB1) is a therapeutically relevant drug target for controlling pain, obesity, and other central nervous system disorders. However, full agonists and antagonists of CB1 have been reported to cause serious side effects in patients. Therefore, partial agonists have emerged as a viable alternative as they can mitigate overstimulation and side effects. One of the key bottlenecks in the design of partial agonists, however, is the lack of understanding of the molecular mechanism of partial agonism itself. In this study, we examine two mechanistic hypotheses for the origin of partial agonism in cannabinoid receptors and predict the mechanistic basis of partial agonism exhibited by Δ9-Tetrahydrocannabinol (THC) against CB1. In particular, we inspect whether partial agonism emerges from the ability of THC to bind in both agonist and antagonist-binding poses or from its ability to only partially activate the receptor. We used extensive molecular dynamics simulations and Markov state modeling to capture the THC binding in both antagonist and agonist-binding poses in the CB1 receptor. Furthermore, we predict that binding of THC in the agonist-binding pose leads to rotation of toggle switch residues and causes partial outward movement of intracellular transmembrane helix 6 (TM6). Our simulations also suggest that the alkyl side chain of THC plays a crucial role in determining partial agonism by stabilizing the ligand in the agonist and antagonist-like poses within the pocket. Taken together, this study provides important insights into the mechanistic origin of the partial agonism of THC.
Untitled section
Keywords: cannabinoid receptors, partial agonist, THC, Markov state model, molecular dynamics, membrane protein
Untitled section
Abbreviations: CB1, cannabinoid receptor 1; GPCR, G-protein-coupled receptor; K-L, Kullback–Leibler; kMC, kinetic Monte Carlo; MD, molecular dynamic; MSM, Markov state model; NAM, negative allosteric modulator; PCA, principal component analysis; TI, thermodynamics integration; THC, tetrahydrocannabinol; TM6, transmembrane helix 6
Article notes
Untitled section
Received 2021 Dec 16; Revised 2022 Feb 14; Collection date 2022 Apr.
Cannabinoid receptor 1(CB1) belongs to the family of class A G-protein-coupled receptors (GPCRs) (1, 2), which modulates diverse cellular signaling processes via intracellular G-proteins (3) and β-arrestins (4). CB1 receptors were first discovered in the last decade of the 20th century as a target of plant cannabinoid molecules (5). Due to its ubiquitous presence in physiological processes, CB1 is an important drug target for the potential treatment of a variety of diseases. In the last 30 years, several synthetic molecules have been designed to target the CB1 for the treatment of pain, obesity, and inflammation (6, 7, 8). CB1 receptor agonists such as MDMB-fubinaca (3) and inverse agonists such as rimonabant (9, 10) have been shown to modulate the receptor activity significantly. However, these designed agonists and antagonists of CB1 also exhibit dangerous side effects. For instance, fubinaca, also known as ”zombie drug,” caused the hospitalization of thousands of patients in New York (11). Rimonabant had to be withdrawn from the market due to its psychological side effects such as depression and anxiety (12). While these designed agonists and antagonists have failed to meet the drug safety guidelines, alternate approaches (e.g., allosteric modulator, partial agonist) can be explored for designing therapeutics. We recently studied the binding of a negative allosteric modulator (NAM), sodium ion (Na+), to cannabinoid receptors (CBs) using molecular dynamics simulation (13). Simulation revealed important differences in binding site and pathway between CB1 and CB2, which can be exploited to design a selective NAM drug. Similarly, a partial agonist dronabinol has been used as an appetite stimulant drug for AIDS patients and an antiemetic drug for chemotherapy (14). Dronabinol is a synthetic form of tetrahydrocannabinol (THC) (Fig. 1), the main psychoactive compound in marijuana, which binds to the CB1 as a partial agonist and affects the endocannabinoid signaling pathway. THC has been shown to demonstrate positive effects for treating Huntington’s disease, Parkinson’s disease, and Alzheimer’s disease (15). Although THC has the potential to become a valuable drug for several diseases, this drug is still banned by the FDA due to its side effects. Therefore, molecular-level understanding of THC and the mechanism by which it partially activates the cannabinoid receptors will inform the design of potential partial agonist drugs targeting CB1 receptor.
Structural studies of CB1 have revealed that toggle switch residue (TRP3566.48 and PHE2003.36) movement by the agonist is crucial for the activation of CB1 (Fig. 2B) (3, 16, 17, 18, 19, 20, 21). Kumar et al. (3) proposed that due to the smaller size of THC as compared with the full agonists, there is less interaction between toggle switch residues and ligand. Furthermore, they proposed that the THC-binding position with downward facing alkyl chain can shift the toggle switch residues and activate the receptor (3). However, the crystal structures and docking studies only provide static interactions of the ligand in the binding pose. These studies do not reveal the mechanism of ligand conformational switching inside the binding pocket. To date, this proposed hypothesis has not been rigorously examined from either the experimental or computational approaches. Therefore, it is difficult to obtain structural understanding of the origin of partial agonism exhibited by THC. A recent study employing metadynamics simulations predicted that a partial agonist GAT228 binds in multiple positions inside the ligand-binding pocket of CB1 (22) due to the large size of the pocket as compared with the ligand volume (18). This prediction is also consistent with the distinct poses of agonist and antagonist in the binding pocket (Fig. 2A). Therefore, it is likely that partial agonist THC might also be stabilized in both the agonist and antagonist-like pose inside the binding pocket (Equations 2 and 3), and only the subset of THC molecules bound to CB1 in the agonist-like pose activates the receptor (Equation 2). This phenomenon would decrease the maximum response of the secondary messenger. The first hypothesis is shown as Equations (1), (2), (3) where R represents the receptor in apo form, RA∗∗ is the agonist bound active state; RPA∗∗ and [RPA] represent partial agonist bound active and inactive state, respectively. In RPA∗∗ and [RPA], partial agonist binds in agonist and antagonist-bound poses, respectively. Agonist and partial agonist are represented as A and PA.
Untitled section
Hypothesis I
Canonical class A GPCR activation is characterized by intracellular TM6 movement, which facilitates the G-protein binding. A partial agonist, salmeterol, binds in the orthosteric pocket in the β2-adrenergic receptor and causes partial movement of TM6 compared with the full movement by an agonist (Fig. 2C) (23, 24, 25). Therefore, we hypothesized that partial activation may happen if the partial agonist stabilizes the receptor in different conformation than the active structure. THC has smaller side chain compared with agonist AM11542 (18). Furthermore, absence of the dimethyl group at the first carbon (C1') of the alkyl chain decreases the interaction with toggle switch residues. Thus, we proposed that THC binding may also cause the partial outward movement of the TM6.
Hypothesis II
The second proposed hypothesis is explained in Equations 4 and 5 where RPA∗ represents partial agonist-bound partially active form of the receptor. Other notations are similar to the reactions 1, 2, 3.
We assessed the validity of these hypotheses by running extensive simulations of THC binding to CB1 (see Experimental procedures section). Using the Markovian property of molecular dynamic (MD) simulations, we built Markov state model (MSM) using the simulation data. From the eigenvalues and vectors of MSM, we predicted the timescale and the thermodynamics of the ligand-binding process. MSM weighted data predict that THC is stabilized in both antagonist and agonist-bound poses in agreement with our first hypothesis. The free energy barrier of ∼1 kcal/mol is estimated for THC to transition between the antagonist and agonist-bound pose, which implies that thermal fluctuations at the human body temperature could easily allow the THC to transition between these poses. In the agonist-bound pose, THC rotates the important toggle switch residue TRP3566.48. The new position of the TRP3566.48 is different compared with active state of CB1. This intermediate position TRP3566.48 leads to partial outward movement of TM6 suggesting that THC can only partially active the receptor according to the second hypothesis. Furthermore, comparison of the ensemble where THC bound in agonist pose with the agonist (AM11542)-bound simulation showed that major differences are in the movement and signaling of the activation microswitches, supporting our results that the partial agonist cannot lead to the full activation of the receptor. This mechanistic study predicts the reason behind the partial agonism behavior of THC compared with other agonists and will aid future drug development targeting cannabinoid receptors.
Results
THC is predicted to be stabilized in both antagonist and agonist-binding poses in the orthosteric pocket of CB1
Active and inactive structures of CB1 reveal that orthosteric-binding site volume undergoes a large change upon activation as compared with other class A GPCRs (18). Comparison of active (PDB ID: 5XRA (17)) and inactive (PDB ID: 5TGZ (18)) structures also reveals that the agonist and antagonist bind in the different regions within the pocket (Fig. 2A). Agonist molecule binds in a region close to the TM5, whereas antagonist binds in the extended pocket formed by TM1 and TM2. In the inactive structure, the downward movement of the N-loop toward the binding pocket separates the agonist and antagonist-binding regions (Fig. 2A). However, partial agonist-bound crystal structure of CB1 or CB2 is not reported in the literature. Therefore, the partial agonist-binding position is not well documented. Preliminary docking studies in both inactive and active structures reveal that THC binds to the agonist-binding region in a similar conformation as other agonists (17, 18). However, docking alone cannot infer the exact binding pose of a ligand, which undergoes a dynamic conformational change within the binding pocket. Therefore, we performed ∼143 μs of MD simulations to predict the THC-binding mechanism in CB1 starting with inactive structure.
To capture the THC-binding process and upward movement of the N-loop, we projected the MD simulation data along the two metrics that characterize these motions (Figs. 3A and S1, A and B). MSM weighted free energy landscape plot (Fig. 3A) shows the predicted movement of THC molecule toward TM5 as indicated by the distance between THC-C1' atom (Fig. 1) and TYR2755.39-Cα atom on TM5. THC is predicted to enter CB1-binding pocket through the space between the space of TM1, TM2, and N-loop (Fig. S2) (16). From MD simulations, THC is further predicted to diffuse inside the pocket and stabilized in antagonist-binding pose where distance between THC(C1') and TYR2755.39(Cα) is 15 to 21 Å (Fig. 3B). We observed two stable local energy wells in the antagonist-binding pose. The free energy well away from agonist-bound pose is named antagonist-like pose 1, and second well is named antagonist-like pose 2. These two minima are separated by the activation barrier of 0.55 ± 0.43 kcal/mol. Superposition of predicted MD structures of THC bound in antagonist-like pose 1 and 2 with inactive structure shows that B and C ring of the tricyclic dibenzopyran group of THC binds in same position as arm 3 of the antagonist, AM6538 (Fig. 3, C and D) (17). The aromatic ring of THC matches with pyrazole ring of the antagonist. However, THC alkyl chain (side chain) orientation varies between the two binding poses. In antagonist-like pose 1, it extends toward the conserved sodium-binding site (13, 26, 27) similar to the arm 1 of antagonist, whereas, in antagonist-like pose 2, it orients itself in the direction of agonist-binding site similar to the arm 2 of the antagonist (Fig. 3, C and D). In both the antagonist-like poses, THC forms stable polar interaction with SER3837.39 and hydrophobic interaction with N-loop (PHE102N−loop, MET103N−loop), TM1 (SER1231.39, ILE1191.35), TM2 (PHE1702.57, PHE1742.61, PHE1772.64), and TM7 (ALA3807.36) (Figs. 4, B and C, S2B and S3A). Previous experiments have shown the importance of PHE1742.61 and PHE1772.64 in the THC binding (28). Relative free energy calculation using alchemical analysis shows that mutations of PHE1742.61 and PHE1772.64 decrease the stability ( kcal/mol and kcal/mol) of the antagonist pose 1 suggesting that these mutations might decrease the THC binding by destabilizing the antagonist-bound pose (29). In these poses, N-loop remains inside the pocket and restricts the movement of THC toward the agonist-binding region. The upward movement of N-loop allows the movement of THC inside the pocket and stabilizes it in the agonist-like pose (Figs. 3E and S16). THC has smaller alkyl side chain than cannabinoid-like full agonists of CB1; thus, first carbon of THC (C1') binds closer to TM5 compared with the agonist, AM11542. In agonist-binding pose, THC(OH) forms polar interaction with SER3837.39 same as other agonists. Furthermore, THC forms extensive hydrophobic interactions with amino acid residues of TM3 (VAL1963.32, LEU1933.29), TM5 (TRP2795.43), TM7 (PHE3797.35), N-loop (MET103N−loop), and ECL2 (PHE268ECL2) as shown in Figure 4D and Fig. S3C.
To capture the timescale of this entire binding process, we ran kinetic Monte Carlo (kMC) simulations on MD data. kMC utilizes MSM transition probability matrix to find the probable pathway for binding (method section). In total, 150 μs long KMC trajectory reveals that entire binding process from unbound to bound poses takes approximately 100 μs (Fig. S4A). From the solution, THC is first stabilized in antagonist-bound pose in ∼50 μs. THC occupies antagonist pose I and II for approximately 30 μs and subsequently moves to the agonist-binding pose (Fig. S4A).
THC chain orientation plays an important role in partial agonism
Alkyl chain of THC plays an important role in binding and activation of CB1. Modification of Alkyl chain leads to change in binding affinity. For example, increasing the chain length of Δ8-THC (structural homolog of Δ9-THC) from five carbons to eight carbons increases the binding affinity from ∼40 nM to ∼8 nM (30). Furthermore, adding a dimethyl group in the first carbon of the alkyl chain is hypothesized to increase the interaction with toggle switch TRP3566.48 (18). To characterize the importance of the alkyl chain of THC, we observed chain dihedral angle (C2-C3-C1'-C2') (Fig. 1) movement during binding. Positive dihedral is crucial to orient the alkyl chain of THC in orthogonal direction of aromatic group as in agonist-like pose (Fig. 3E). The free energy landscape of THC alkyl chain dihedral with respect to binding (Figs. 5A and S5) reveals that THC binds to CB1 in two different chain orientations. If hydrogen atom on C1 faces downward (or positive dihedral angle) during the binding, it moves to agonist-binding pocket with maximum free energy barrier of ∼2 ± 0.4 kcal/mol. This alkyl chain orientation enables THC to pressurize N-loop in upward direction to take orthogonal conformation similar to full agonist. However, if THC enters the receptor with negative dihedral, it is stabilized in the antagonist-like poses. High free energy barrier of ∼4 ± 0.4 kcal/mol is required for THC to move from macrostate 3′ to 4′ (Fig. 5A). In this conformation (macrostate 3′), tricyclic dibenzopyran group of THC is unable shift N-loop upward as the active structure. Hence, alkyl chain of THC shifts the population to macrostate 3, which is more accessible as the free energy barrier is lower than earlier transition. Applying transition path theory (TPT) on MSM states, we calculated effective timescales of these macrostate transitions. TPT provides mean free passage time (MFPT) for the transition between two macrostates by taking into account all possible pathways through the intermediate states. Calculated MFPT shows that transition between state 3′ and 3 (2.4 ± 0.5 μs) is more accessible compared with the 3′ to 4′ (13.2 ± 3.5 μs) (Figs. 5B and S6). The binding of THC to the agonist-like pose shows the orientation of alkyl chain may flip between macrostate 4 and 4′ as they are energetically favorable. However, THC has more tendency to stay in the orthogonal conformation (macrostate 4) compared with macrostate 4′. The conditional probability of THC to be in orthogonal conformation in the agonist-like pose is 0.9 ± 0.0. Our results show that THC is stabilized in agonist and antagonist-like poses that supports our first hypothesis.
THC is predicted to rotate the toggle switch TRP3566.48 in the agonist-binding pocket
Although, we established the prediction that THC can be stabilized in different positions of the orthosteric-binding pocket, it is not clear how THC activates the receptor. Crystal structures of the CB1 receptor in active and inactive position depict that toggle switch residues play an important role in the activation of the receptor. Agonist molecule triggers the movement of toggle switch residue TRP3566.48 movement of the receptor toward the TM5, which consequently leads to outward movement of TM6 (18). It was hypothesized that THC behaves as partial agonist due to the lack of interaction with toggle switch residues (3). However, we observed that binding of THC leads to the rotation of TRP3566.48 (Figs. 6A and S7). The dihedral angle of TRP3566.48 shifts from inactive conformation (χ2 dihedral angle between 60° and 120°) to new intermediate active conformation (−30° to 30°) as well as relatively less stable state (−120° to −60°) (Fig. 6, A and B). We referred energetically favorably accessible states as partially active state 1 and state 2, respectively. The rotation of TRP3566.48 leads to breakage of aromatic interactions with PHE2003.36, and it moves toward TM2 similar to active structure (Fig. 6B). Comparison of representative structures from partially active state 1 with inactive CB2 shows that toggle switch TRP3566.48/TRP2586.48 has similar rotamaric conformation as CB2 inactive structure (Fig. S8). An inverse agonist of CB2, MRI2687, which was predicted to retain the toggle switch residue to similar conformation (31), acts as a partial agonist for CB1. This in turn supports our prediction that THC stabilizes the TRP6.48 in intermediate state to partially activate the receptor.
Although, THC is predicted to rotate toggle switch TRP3566.48 for subset of structures in agonist-binding position, we noticed another stable state around inactive pose of TRP3566.48 (Fig. 6A). In this state, THC is likely unable to shift the toggle switch movement. To explain the reason for the two different orientation of TRP3566.48, we calculated the probability TRP3566.48 rotation with respect to the THC alkyl chain dihedral in agonist-like pose (Fig. 6C). Our calculations show that with negative alkyl chain dihedral, THC can move TRP3566.48 for only 4.0 ± 0.8% structures. With negative dihedral, the first carbon of the alkyl chain (C1') binds far from the toggle switch and cannot induce conformational change (Fig. 5B). However, binding with a positive dihedral leads to rotation of TRP3566.48 for 95.5 ± 0.8% of the structures due to more interactions.
The changes in the toggle switch residue lead to the outward movement of TM6 as shown in the Figure 7A and Fig. S9. When the toggle switch remains in the inactive position, intracellular TM6 can move easily 2 to 3 Å either side inactive structure (low free energy barrier). For this case, TM3-TM6 distance is similar or less than the inactive structure with a probability of 0.74 ± 0.0 (Fig. 7B). However, TRP3566.48 rotation creates a torsion in TM6 and stabilizes the intracellular part of TM6 of the receptor in between inactive and active-like conformation with probability of 0.6 ± 0.0 (Fig. 7B). Comparison of representative structure from partially active state minima with β2-AR inactive structure shows that intracellular TM5, TM6, and TM7 match well (Fig. S10). Therefore, THC in partially activated state forms favorable interaction with TRP3566.48 and stabilizes in intermediate conformation compared with the inactive structure. These results support our second hypothesis that THC can only activate the receptor partially.
Comparison of partial agonist and agonist-bound simulation ensemble using data-driven approach
We have performed 17 μs MD simulations of agonist (AM11542) bound CB1 (Experimental procedures section) and compared the result with the ensemble with the partial agonist bound to the agonist-like pose (32). Principal component analysis (PCA) was implemented on all possible interresidue distances to compare the ligand-bound ensembles of the partial and full agonists. PCA projection onto the first two components shows that the receptor explores different conformational spaces for different ligands (Fig. 8A). As PC1 and PC2 emphasize large amplitude motions, we calculated distance pairs with highest weight contributions for PC1 and PC2. For both agonist and partial agonist-bound ensembles, these distances correspond to the distances between truncated N-terminus and transmembrane domain (Figs. S11 and S12). Binding of agonist and partial agonist leads to outward movement of N-terminus from the orthosteric binding pocket compared with the inactive state (18). As the N-terminus moves out of the binding pocket when CB1 is in the bound state, it fluctuates more because of the reduced interactions with binding pocket residues. Hence, PC1 and PC2 capture these large amplitude motions for both the ensembles. We also performed Kullback–Leibler (K-L) divergence analysis for each interresidue distance to determine the major differences between these two ensembles. K-L divergence analysis shows that the higher divergence occurs in intracellular TM5, TM6, and TM7 (Figs. 8B and S13), depicting the difference in activation microswitches for agonist and antagonist. CB1 has three well-known microswitches: toggle switch movement, intracellular TM6 outward movement (DRY motif salt bridge breaking), and TM7 (NPxxY motif) movement (Fig. S14). Comparison of the agonist bound active-like structure and inactive structure shows that these microswitches rearrange themselves during the protein activation. To capture these motions of the toggle switch, TM6, and TM7, TRP3566.48-ASP1632.50, LYS3436.35-ARG2143.50, and TYR1532.40-TYR3977.53 distances were calculated. Projection of activation microswitch distances onto a 2-D scatter plot shows the difference in activation microswitch movement for partial agonist and agonist-bound ensembles (Fig. 8, C and D). Although the partial agonist affects the movement of these microswitches, it cannot lead to the full activation of the receptor as described earlier. In contrast, an agonist stabilizes the receptor in an active-like pose. Difference in the intracellular microswitch movement might happen due to the difference in allosteric communication pathway between different microswitches. Specifically, we studied how the signaling is transduced from the toggle switch residue to the NPxxY region and compared the pathway for agonist and partial agonist-bound ensemble. An allosteric communication pathway was determined by estimating a mutual information network between residue dihedral angles (33, 34, 35). Using this network, we estimated the shortest pathway between the toggle switch and the NPxxY region using Dijkstra’s algorithm (36). Results show that the shortest paths for signal transduction between toggle switch and NPxxY region are different for partial agonist and agonist (Fig. S15, A and B). For the agonist, we observed an allosteric signaling pathway between the NPxxY region and the toggle switch transduce through the outward movement of intracellular TM6, whereas, due to the lack of intracellular TM6 movement, allosteric signal transduction happens from a different pathway for partial agonist. Therefore, these analyses show how the activation signaling differs between partial agonists and agonists as the ligands affect the toggle switch differently.
Discussion
In this study, we proposed two hypotheses for the partial agonism of THC molecule for Cannabinoid receptor 1. Our first hypothesis is based on different binding position of the THC inside the orthosteric binding pocket. Second hypothesis states that THC may only able to move TM6 partially. To test the hypotheses, we performed unbiased molecular dynamics simulation for THC binding to CB1. Our results support both our hypotheses. Simulations show that during binding THC is stabilized in agonist-like and antagonist-like poses. While binding in the antagonist-like poses, the aromatic group of THC orients itself as arm 2 of the antagonist, whereas alkyl chain of THC can take the conformation of arm 1 and arm 3. During the binding process, alkyl chain orientation is shown to be important factor of determining the binding of THC to agonist pose. When THC enters the receptor with positive side chain dihedral angle (C2-C3-C1'-C2'), it leads to the upward movement of N-loop and favors the THC binding to the agonist-like pose, whereas, with negative dihedral angle, THC increases the free energy barrier for transition from antagonist-like binding pose to agonist-like binding pose. Therefore, chemical modification, which can stabilize the side chain in perpendicular direction to the aromatic ring, may increase the agonistic property of the ligand.
Over the years, various effects on THC alkyl chain modifications on binding and functionality have been proposed. Δ9-THC and Δ8-THC have similar binding affinity and agonistic property. It has been shown that Δ8-THC analog with cis double bond between C1 and C2 (which fixes the side chain in perpendicular direction) with same side chain length increases the agonistic property of the ligand (30). Our simulations also support this experimental observation. With cis bond, THC analog can move to the agonist-like pose with less free energy barrier (Fig. 4A) and therefore has more agonist property, which bolsters our claim in the first hypothesis.
To generalize our first hypothesis that THC can be stabilized at both agonist and antagonist pocket, we performed docking of other partial agonists on representative CB1 structures of THC bound agonist-binding pose and antagonist-binding pose (Fig. S17). Five available partial agonists are selected for GPCRdb database (37): magnolol, AM4089, NMP-4, (S)-Δ3-THC, (R)-Δ3-THC. Docking studies reveal stable docked poses for partial agonists in both the binding pockets. In the antagonist-binding pose, we observed that the partial agonists extend downward toward Na+-binding site similar to THC (13). These partial agonists bind with similar affinity in agonist and antagonist-binding pocket. These docking results show that our hypothesis for partial agonism may be universally valid for other partial agonists for CB1.
Our results also show that the THC can rotate the toggle switch TRP3566.48 conformation when binding to the agonist-like poses. The partially active conformation of CB1 matches well with CB2 inactive structure, which helps to explain yin–yang relation for some ligands between CB1 and CB2. Recently published NMR study on another class A GPCR, AA2AR, has shown that partial agonist rotates TRP2466.48 in a distinct conformation compared with full agonists or antagonists (38). These experimental observations support our results obtained computationally. Furthermore, our simulations reveal that partially active toggle switch movement consequently affects the intracellular side of the receptor and leads to the partial outward movement of TM6. In future study, this partial movement of the TM6 can be experimentally validated by DEER spectroscopy (39, 40, 41, 42).
Overall, our present computational study with MD simulations and MSM is competent with previous experimental results and also provides new predictions for partial agonism of THC for CB1. The findings of the paper are crucial for guiding the chemical modification for cannabinoid ligands for future drug development.
Experimental procedures
System preparation
Inactive structure of CB1 (PDB ID: 5TGZ (17)) is used as a starting structure for MD simulation for the THC binding. This structure does not contain first 98 residues in the N-terminus. However, previous study has shown that most of the N-terminal residues (1–103) can be deleted without abolishing the ligand binding for synthetic cannabinoids (43). Truncated N-terminus and C-terminus and unconnected residues of TM5 and TM6 are capped with neutral terminal residues (acetyl and methylamide groups). Hydrogen atoms are added to protein amino acid residues using reduce command of AMBER package (44). Thermostabilized mutant residues are replaced with original residues using tleap (17). The modified protein structure is embedded in POPC bilayer using CHARMM-GUI (45). Salt concentration of 150 mM (Na+ and Cl−) is to neutralize the system. MD system is solvated using TIP3P water model. 3-D structure of THC is obtained from PubChem in sdf format. Forcefield parameters of THC are obtained using antechamber (46, 47). THC is added to MD system using packmol (48) to generate the starting structure. To perform agonist (AM11542) bound simulation, active-like agonist-bound crystal structure (PDB ID: 5XRA (18)) was selected. Necessary modifications were performed as discussed as for THC-binding simulation. Agonist structure was parameterized using GAFF forcefield. CHARMM-GUI software was used to embed the agonist bound structure in POPC bilayer with NaCl solution (150 mM) in extracellular and intracellular regions.
Simulation details
MD simulations were performed using AMBER18. The MD systems were subjected to minimization with gradient descent and conjugate gradient algorithm for 5000 and 10,000 steps, respectively. Minimized systems were slowly heated from 0 K to 10 K and from 10 K to 300 K in NVT ensemble to increase the temperature at desired level. Each step was done for 1 ns period. To control the pressure of the heated system at 1 bar, NPT ensemble was employed. During modulation of temperature and pressure, protein backbone Cα atoms are restrained with a spring force. Furthermore, total 50 ns of equilibration was performed in NPT ensemble for each system to maintain the temperature and pressure at desired level of 300 K at 1 bar without any restraint force. Production runs were also performed in NPT ensemble. The temperature and pressure of the systems were maintained by Berendsen thermostat and barostat (49). Simulation timestep of 2 fs was used for MD simulation. As hydrogen atom is the lightest atom in the system, vibrational frequency of the hydrogen can be lower than the specified timescale. It can create instability in the system due to the large fluctuation of hydrogen atom. Therefore, SHAKE algorithm (50) was used to put restraint in the movement of hydrogen atom by implementing Lagrangian multiplier. For nonbonded force calculation, 10 Å cutoff distance was used. Periodic boundary condition was applied for all simulations. To consider the nonbonded long-range interactions, particle mesh Ewald method was implemented (51). Simulations were performed with Amber FF14SB (52) and GAFF (46) forcefield parameters.
Adaptive sampling
Ligand binding to GPCR is rare event compared with MD simulations timescale. To capture the entire THC-binding process, adaptive sampling technique is utilized. This approach has been shown to sample the conformational ensemble of a variety of biological systems (27, 53, 54, 55) including the ligand-binding process (13, 56, 57, 58, 59). In this technique, data from one round of simulation are represented by conformational space of protein by using biological relevant features. The features are used to cluster the protein space using k-means clustering algorithm. The structures for next round of simulation are selected randomly from the cluster centers with lowest count of data points. This process is repeated in each round. For this case, intracellular and extracellular helical distances and THC-binding distance from TM5 were used as adaptive sampling feature matrices. Although adaptive sampling helps to parallelize the simulation by sampling from lower probability space, it affects the ensemble distribution and therefore brings sampling error in free energy calculation. To overcome this caveat, MSM (discussed below) is implemented to remove the sampling the bias from the simulation data. A total 143 μs of simulation data were collected using adaptive sampling protocol (Table S2). Agonist-bound simulations were also performed using adaptive sampling protocol using similar matrices as partial agonist THC-binding simulation. A total 17 μs of simulation data were collected for agonist-bound simulation.
Markov state model
The theory of MSM (60, 61, 62) depicts that in a sequence of event, the transition probability of moving from one state to another depends only on the present state and not on the path the system has taken to be there. The trajectory of MD simulation follows the same principle. According to the verlet algorithm, the evolution coordinate and momenta of atoms in system only depend on the present state. Therefore, MD trajectory can be assumed to be Markovian and according to its nature, the probability distribution of the protein conformational ensemble can be calculated from the transition probability matrix between the states. Each element Tij of transition probability matrix calculates the probability of transition using the equation where Cij is the count of jump between the i and j, and Ci is the count of frame in state i. To make the transition probability matrix statistically significant, states with conformation are clustered together into microstates assuming that there is no large energy barrier in the same microstate. However, clustering of conformational space increases the memory of each state, which can invoke “non-markovinity” in our system. This issue can be overcome by increasing the lag time (τ) such that it preserves the Markov property and hence satisfies the equation where and are vectors representing the probability of the microstates. Pyemma (63) package was used to construct MSM. We selected 23 biologically important features to capture the THC binding and protein conformational changes (Table S1 and Fig. S22). Features were transformed into time-lagged independent components (tics) (64, 65) to find the slowest components. Tic components were shown to correlate well with important features for THC binding (Fig. S18, A–C). The transformed data were clustered using k-means clustering algorithm. For MSM building, lag time was chosen by finding out logarithmic convergence of process timescales computed from MSM eigenvalues (process timescale ; λ is the eigenvalue) (Fig. S19). To find out the optimal number of clusters and optimal tic variance, MSM was subjected to VAMP2 (64) scoring to measure the kinetic variance (Fig. S20). Based on transition matrix, MSM predicts the population (stationary probability, πi) of each clustered state, which is needed to calculate thermodynamics (Gi = −kbT logeπi) of the process. Using similar steps, optimal MSM was built for agonist-bound simulation with 65% tic variance and 25 ns of MSM lag time.
Transition path theory
To obtain the kinetics for transition between MSM macrostates, TPT (66, 67) was implemented. TPT calculates the flux between two macrostates (which consists of one or multiply MSM states) using MSM transition matrix (68). The flux (FAB) is given by the equation where is probability that state j will reach macrostate B before A (69). From the flux calculation, we computed MFPT using equation MFPT = τπA/FAB. where πA is the probability the system was in macrostate A. TPT calculation was performed with Pyemma (63) package. Here, the macrostates were defined manually. Five MSM states with highest raw count inside area of interest are defined as macrostates. For example, macrostate 4 in Figure 5A consisted of five MSM states with highest raw count inside the area where THC(C1')-TYR2755.39(Cα) distance is between 9 Å and 13 Å & THC dihedral angle is between 60° and 120°.
Trajectory analysis
All the feature calculation and data processing from the MD simulation were done CPPtraj (70) and MDtraj (71). For the trajectory visualization and analysis, VMD (72) and pymol v2.1 (Schrödinger, LLC) software package were used. To calculate ligand and protein residue contact, GetContacts package (https://getcontacts.github.io/) was implemented.
Alchemical calculation
Alchemical calculations were performed on a reference frame selected from the antagonist-bound pose with thermodynamics integration (TI) technique (73). To prepare the system for TI calculations, alanine mutations were added to the system using tiMerge function of ParmEd software package. Two different systems were prepared for two mutations (PHE1742.61ALA and PHE1772.64ALA). AMBER PMEMD engine was used to carry out TI calculation. Simulation conditions were maintained same as the original simulation condition. For TI, 11 λ windows were prepared between 0 and 1 with a gap of 0.1. Softcore potential was used for van der Waals (vdw) and charge contribution for the dummy atoms. In each window, 2 ns simulation data were generated. Free energy estimation was performed with pyamber package (29).
Kinetic Monte Carlo simulation
kMC simulation is a stochastic method to visualize the evolution of a system based on probability of the transition between different states. In this case, transition probability of different states was obtained from MSM. There are few steps to implement kMC in MSM weighted MD data. First, we generated a random number (R) between 0 and 1. Second, we described the possible transition event from i as discrete cumulation probability distribution. Therefore, the cumulative transition probability between i and j can be written as . If R lies between and , then system transition happens from ith to jth state. Time required for each transition is taken to be the lag time considered for MSM. Same procedure is repeated for desired number of steps. To build a trajectory from kMC state evolution, we picked a random frame at each step from the chosen state. As we were interested to observe the ligand-binding dynamics, kMC simulation is started from a state where THC was 25 Å away from TYR2755.39(Cα).
Data-driven modeling
Data-driven analysis was performed to compare agonist and partial agonist-bound ensembles. For THC, 20,000 frames were selected from clusters where THC is bound to the agonist-bound pose based on cluster probability calculated using MSM. Similarly, 20,000 frames were selected from agonist-bound simulation based on cluster probability calculated using MSM for the agonist-bound simulation. All closest heavy carton distances were calculated using MDtraj (71). PCA was performed using sklearn package of python. Here, symmetrized version of K-L divergence was used to calculate the difference in distance distribution of two ensembles to find out major distinctions. Lastly, to find the activation path, normalized mutual information was calculated for dihedral angles for every residue pair using DihedralMutualInformation class of mdentropy package (74). With this information, we developed a weighted undirected graph where each residue denotes a node. The edges of the graph were created between two residue pairs with the criteria that the residues are within 6 Å from each other and there is nonzero mutual information between two residues. From this graph, we calculated the shortest path between toggle switch (TRP3566.48) and NPxxY region (TYR3977.53) using Dijkstra’s algorithm implementation of the Networkx python package (75, 76).
Docking study
Docking study was performed using Auto dock vina (77) software. 3-D structures of partial agonists are selected from PubChem in sdf format. Antechamber is used to convert the ligands in mol2 format and add partial changes. Then, Auto dock is used to convert the ligands into required pdbqt format. For the docking of the ligands in both agonist and antagonist pose, the box size is calculated from THC-binding structure in respective poses. To dock the ligand in antagonist and agonist-bound poses, MD predicted structures are used.
Error analysis
To determine error in our thermodynamics (free energy, conditional probability) and kinetics (TPT) calculations, we performed bootstrap analysis on MD data (78). In each bootstrap sample, we randomly picked N trajectories, where N is equal to 80% of total number of trajectories. We kept the state labeling same as the original MSM (Fig. S21). For every bootstrap sample, MSM is computed to determine thermodynamics and kinetics. In total, 200 bootstrap samples are generated for error calculations.
Code and data availability
MSM feature matrix, final MSM object, codes, and pdb files used to generate figures can be obtained from https://github.com/ShuklaGroup/THC_binding.git.
Supporting information
This article contains supporting information.
Conflict of interest
The authors declare that they have no conflicts of interest with the contents of this article.
Acknowledgments
This research is performed as a part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) the State of Illinois, and as of December, 2019, the National Geospatial-Intelligence Agency. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.
Funding and additional information
D. S. acknowledges support from NSF Early CAREER Award (MCB-1845606).
Untitled section
Edited by Henrik Dohlman
Supporting information
References
Untitled section
References
- 1.Hilger D., Masureel M., Kobilka B.K. Structure and dynamics of GPCR signaling complexes. Nat. Struct. Mol. Biol. 2018;25:4. doi: 10.1038/s41594-017-0011-7.
- 2.Rosenbaum D.M., Rasmussen S.G., Kobilka B.K. The structure and function of G-protein-coupled receptors. Nature. 2009;459:356–363. doi: 10.1038/nature08144.
- 3.Kumar K.K., Shalev-Benami M., Robertson M.J., Hu H., Banister S.D., Hollingsworth S.A, Latorraca N.R., Kato H.E., Hilger D., Maeda S., Weis W.I., Farrens D.L., Dror R.O., Malhotra S.V., Kobilka B.K., et al. Structure of a signaling cannabinoid receptor 1-G protein complex. Cell. 2019;176:448–458.e12. doi: 10.1016/j.cell.2018.11.040.
- 4.Ibsen M.S., Finlay D.B., Patel M., Javitch J.A., Glass M., Grimsey N.L. Cannabinoid CB1 and CB2 receptor-mediated arrestin translocation: Species, subtype, and agonist-dependence. Front. Pharmacol. 2019;10:350. doi: 10.3389/fphar.2019.00350.
- 5.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.
- 6.Pertwee R.G. Pharmacology of cannabinoid receptor ligands. Curr. Med. Chem. 1999;6:635–664.
- 7.Pertwee R.G. Cannabinoid pharmacology: The first 66 years. Br. J. Pharmacol. 2006;147:S163–S171. doi: 10.1038/sj.bjp.0706406.
- 8.Pertwee R.G., Ross R. Cannabinoid receptors and their ligands. Prostaglandins Leukot. Essent. Fatty Acids. 2002;66:101–121. doi: 10.1054/plef.2001.0341.
- 9.MacLennan S., Reynen P., Kwan J., Bonhaus D. Evidence for inverse agonism of SR141716A at human recombinant cannabinoid CB1 and CB2 receptors. Br. J. Pharmacol. 1998;124:619–622. doi: 10.1038/sj.bjp.0701915.
- 10.Pertwee R.G. Inverse agonism and neutral antagonism at cannabinoid CB1 receptors. Life Sci. 2005;76:1307–1324. doi: 10.1016/j.lfs.2004.10.025.
- 11.Adams A.J., Banister S.D., Irizarry L., Trecki J., Schwartz M., Gerona R. “Zombie” outbreak caused by the synthetic cannabinoid AMB-FUBINACA in New York. N. Engl. J. Med. 2017;376:235–242. doi: 10.1056/NEJMoa1610300.
- 12.Sam A.H., Salem V., Ghatei M.A. Rimonabant: From RIO to ban. J. Obes. 2011;2011:432607. doi: 10.1155/2011/432607.
- 13.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.
- 14.Badowski M.E., Yanful P.K. Dronabinol oral solution in the management of anorexia and weight loss in AIDS and cancer. Ther. Clin. Risk Manag. 2018;14:643. doi: 10.2147/TCRM.S126849.
- 15.Walther S., Halpern M. Cannabinoids and dementia: A review of clinical and preclinical data. Pharmaceuticals. 2010;3:2689–2708. doi: 10.3390/ph3082689.
- 16.Hua T., Li X., Wu L., Iliopoulos-Tsoutsouvas C., Wang Y., Wu M., Shen L., Brust C.A., Nikas S.P., Song F., Song X., Yuan S., Sun Q., Wu Y., Jiang S., et al. Activation and signaling mechanism revealed by cannabinoid receptor-Gi complex structures. Cell. 2020;180:655–665.e18. doi: 10.1016/j.cell.2020.01.008.
- 17.Hua T., Vemuri K., Pu M., Qu L., Han G.W., Wu Y., Zhao S., Shui W., Li S., Korde A., Laprairie R.B., Stahl E.L., Ho J.H., Zvonok N., Zhou H., et al. Crystal structure of the human cannabinoid receptor CB1. Cell. 2016;167:750–762.e14. doi: 10.1016/j.cell.2016.10.004.
- 18.Hua T., Vemuri K., Nikas S.P., Laprairie R.B., Wu Y., Qu L., Pu M., Korde A., Jiang S., Ho J.H., Han G.W., Ding K., Li X., Liu H., Hanson M.A., et al. Crystal structures of agonist-bound human cannabinoid receptor CB1. Nature. 2017;547:468–471. doi: 10.1038/nature23272.
- 19.Shao Z., Yin J., Chapman K., Grzemska M., Clark L., Wang J., Rosenbaum D.M. High-resolution crystal structure of the human CB1 cannabinoid receptor. Nature. 2016;540:602–606. doi: 10.1038/nature20613.
- 20.Shao Z., Yan W., Chapman K., Ramesh K., Ferrell A.J., Yin J., Wang X., Xu Q., Rosenbaum D.M. Structure of an allosteric modulator bound to the CB1 cannabinoid receptor. Nat. Chem. Biol. 2019;15:1199–1205. doi: 10.1038/s41589-019-0387-2.
- 21.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.
- 22.Saleh N., Hucke O., Kramer G., Schmidt E., Montel F., Lipinski R., Ferger B., Clark T., Hildebrand P.W., Tautermann C.S. Multiple binding sites contribute to the mechanism of mixed agonistic and positive allosteric modulators of the cannabinoid CB1 receptor. Angew. Chem. Int. Ed. Engl. 2018;130:2610–2615. doi: 10.1002/anie.201708764.
- 23.Rasmussen S.G.F., DeVree B.T., Zou Y., Kruse A.C., Chung K.Y., Kobilka T.S., Thian F.S., Chae P.S., Pardon E., Calinski D., Mathiesen J.M., Shah S.T.A., Lyons J.A., Caffrey M., Gellman S.H., et al. Crystal structure of the β2 adrenergic receptor–Gs protein complex. Nature. 2011;477:549–555. doi: 10.1038/nature10361.
- 24.Masureel M., Zou Y., Picard L.P., Westhuizen E.V.D., Mahoney J.P., Rodrigues J.P.G.L.M., Mildorf T.J., Dror R.O., Shaw D.E., Bouvier M., Pardon E., Steyaert J., Sunahara R.K., Weis W.I., Zhang C., et al. Structural insights into binding specificity, efficacy and bias of a β2AR partial agonist. Nat. Chem. Biol. 2018;14:1059–1066. doi: 10.1038/s41589-018-0145-x.
- 25.Cherezov V., Rosenbaum D.M., Hanson M.A., Rasmussen S.G.F., Thian F.S., Kobilka T.S., Choi H.-J., Kuhn P., Weis W.I., Kobilka B.K., Stevens R.C. High-resolution crystal structure of an engineered human β 2 -adrenergic G protein–coupled receptor. Science. 2007;318:1258–1265. doi: 10.1126/science.1150577.
- 26.Katritch V., Fenalti G., Abola E.E., Roth B.L., Cherezov V., Stevens R.C. Allosteric sodium in class A GPCR signaling. Trends Biochem. Sci. 2014;39:233–244. doi: 10.1016/j.tibs.2014.03.002.
- 27.Selvam B., Shamsi Z., Shukla D. Universality of the sodium ion binding mechanism in class AG-protein-coupled receptors. Angew. Chem. Int. Ed. Engl. 2018;130:3102–3107. doi: 10.1002/anie.201708889.
- 28.Shim J.-Y., Bertalovitz A.C., Kendall D.A. Identification of essential cannabinoid-binding domains. J. Biol. Chem. 2011;286:33422–33435. doi: 10.1074/jbc.M111.261651.
- 29.Song L.F., Lee T.-S., Zhu C., York D.M., Merz K.M. Using AMBER18 for relative free energy calculations. J. Chem. Inf. Model. 2019;59:3128–3135. doi: 10.1021/acs.jcim.9b00105.
- 30.Bow E.W., Rimoldi J.M. The structure–function relationships of classical cannabinoids: CB1/CB2 modulation. Perspect. Med. Chem. 2016;8:17–39. doi: 10.4137/PMC.S32171.
- 31.Li X., Hua T., Vemuri K., Ho J.H., Wu Y., Wu L., Popov P., Benchama O., Zvonok N., Locke K., Qu L., Han G.W., Iyer M.R., Cinar R., Coffey N.J., et al. Crystal structure of the human cannabinoid receptor CB2. Cell. 2019;176:459–467.e13. doi: 10.1016/j.cell.2018.12.011.
- 32.Fleetwood O., Carlsson J., Delemotte L. Identification of ligand-specific G protein-coupled receptor states and prediction of downstream efficacy via data-driven modeling. Elife. 2020 doi: 10.7554/elife.60715.sa2.
- 33.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.
- 34.Vaidehi N., Bhattacharya S. Allosteric communication pipelines in G-protein-coupled receptors. Curr. Opin. Pharmacol. 2016;30:76–83. doi: 10.1016/j.coph.2016.07.010.
- 35.Nivedha A.K., Tautermann C.S., Bhattacharya S., Lee S., Casarosa P., Kollak I., Kiechle T., Vaidehi N. Identifying functional hotspot residues for biased ligand design in G-protein-coupled receptors. Mol. Pharmacol. 2018;93:288–296. doi: 10.1124/mol.117.110395.
- 36.Kapoor A., Martinez-Rosell G., Provasi D., de Fabritiis G., Filizola M. Dynamic and kinetic elements of μ-opioid receptor functional selectivity. Sci. Rep. 2017;7:11255. doi: 10.1038/s41598-017-11483-8.
- 37.Pándy-Szekeres G., Munk C., Tsonkov T.M., Mordalski S., Harpsøe K., Hauser A.S., Bojarski A.J., Gloriam D.E. GPCRdb in 2018: Adding GPCR structure models and ligands. Nucleic Acids Res. 2018;46:D440–D446. doi: 10.1093/nar/gkx1109.
- 38.Eddy M.T., Martin B.T., Wüthrich K. A2A adenosine receptor partial agonism related to structural rearrangements in an activation microswitch. Structure. 2021;29:170–176. doi: 10.1016/j.str.2020.11.005.
- 39.Jeschke G. DEER distance measurements on proteins. Annu. Rev. Phys. Chem. 2012;63:419–446. doi: 10.1146/annurev-physchem-032511-143716.
- 40.Manglik A., Kim T.H., Masureel M., Altenbach C., Yang Z., Hilger D., Lerch M.T., Kobilka T.S., Thian F.S., Hubbell W.L., Prosser R.S., Kobilka B.K. Structural insights into the dynamic process of β 2 -adrenergic receptor signaling. Cell. 2015;161:1101–1111. doi: 10.1016/j.cell.2015.04.043.
- 41.Mittal S., Shukla D. Maximizing kinetic information gain of Markov state models for optimal design of spectroscopy experiments. J. Phys. Chem. B. 2018;122:10793–10805. doi: 10.1021/acs.jpcb.8b07076.
- 42.Mittal S., Shukla D. Recruiting machine learning methods for molecular simulations of proteins. Mol. Simul. 2018;44:891–904.
- 43.Fay J.F., Farrens D.L. The membrane proximal region of the cannabinoid receptor CB1N-terminus can allosterically modulate ligand affinity. Biochemistry. 2013;52:8286–8294. doi: 10.1021/bi400842k.
- 44.Case D.A., Aktulga H.M., Belfon K., Ben-Shalom I.Y., Brozell S.R., Cerutti D.S., Cheatham T.E., Cisneros G.A., Cruzeiro V.W.D., Darden T.A., Duke R.E., Giambasu G., Gilson M.K., Gohlke H., Goetz A.W., et al. University of California; San Francisco, CA: 2018. AMBER 2018.
- 45.Jo S., Kim T., Iyer V.G., Im W. CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008;29:1859–1865. doi: 10.1002/jcc.20945.
- 46.Wang J., Wolf R.M., Caldwell J.W., Kollman P.A., Case D.A. Development and testing of a general amber force field. J. Comput. Chem. 2004;25:1157–1174. doi: 10.1002/jcc.20035.
- 47.Wang J., Wang W., Kollman P.A., Case D.A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 2006;25:247–260. doi: 10.1016/j.jmgm.2005.12.005.
- 48.Martínez L., Andrade R., Birgin E.G., Martínez J.M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009;30:2157–2164. doi: 10.1002/jcc.21224.
- 49.Berendsen H.J., Postma J.V., van Gunsteren W.F., DiNola A., Haak J.R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984;81:3684–3690.
- 50.Kräutler V., Van Gunsteren W.F., Hünenberger P.H. A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 2001;22:501–508.
- 51.Essmann U., Perera L., Berkowitz M.L., Darden T., Lee H., Pedersen L.G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995;103:8577–8593.
- 52.Maier J.A., Martinez C., Kasavajhala K., Wickstrom L., Hauser K.E., Simmerling C. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015;11:3696–3713. doi: 10.1021/acs.jctc.5b00255.
- 53.Selvam B., Mittal S., Shukla D. Free energy landscape of the complete transport cycle in a key bacterial transporter. ACS Cent. Sci. 2018;4:1146–1154. doi: 10.1021/acscentsci.8b00330.
- 54.Selvam B., Yu Y.-C., Chen L.-Q., Shukla D. Molecular basis of the glucose transport mechanism in plants. ACS Cent. Sci. 2019;5:1085–1096. doi: 10.1021/acscentsci.9b00252.
- 55.Kohlhoff K.J., Shukla D., Lawrenz M., Bowman G.R., Konerding D.E., Belov D., Altman R.B., Pande V.S. Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways. Nat. Chem. 2014;6:15. doi: 10.1038/nchem.1821.
- 56.Shukla S., Zhao C., Shukla D. Dewetting controls plant hormone perception and initiation of drought resistance signaling. Structure. 2019;27:692–702. doi: 10.1016/j.str.2018.12.005.
- 57.Chen J., White A., Nelson D.C., Shukla D. Role of substrate recognition in modulating strigolactone receptor selectivity in witchweed. biorXiv. 2020 doi: 10.1101/2020.07.28.225722. [preprint]
- 58.Aldukhi F., Deb A., Zhao C., Moffett A.S., Shukla D. Molecular mechanism of brassinosteroid perception by the plant growth receptor BRI1. J. Phys. Chem. B. 2019;124:355–365. doi: 10.1021/acs.jpcb.9b09377.
- 59.Lawrenz M., Shukla D., Pande V.S. Cloud computing approaches for prediction of ligand binding poses and pathways. Sci. Rep. 2015;5:7918. doi: 10.1038/srep07918.
- 60.Bowman G.R., Pande V.S., Noé F. Vol. 797. Springer Science & Business Media; Berlin: 2013. (An introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation).
- 61.Prinz J.-H., Wu H., Sarich M., Keller B., Senne M., Held M., Chodera J.D., Schütte C., Noé F. Markov models of molecular kinetics: Generation and validation. J. Chem. Phys. 2011;134:174105. doi: 10.1063/1.3565032.
- 62.Shukla D., Hernández C.X., Weber J.K., Pande V.S. Markov state models provide insights into dynamic modulation of protein function. Acc. Chem. Res. 2015;48:414–422. doi: 10.1021/ar5002999.
- 63.Wehmeyer C., Scherer M.K., Hempel T., Husic B.E., Olsson S., Noé F. Introduction to Markov state modeling with the PyEMMA softwar v1. 0. LiveCoMS. 2018;1:1–12.
- 64.Noé F., Clementi C. Kinetic distance and kinetic maps from molecular dynamics simulation. J. Chem. Theory Comput. 2015;11:5002–5011. doi: 10.1021/acs.jctc.5b00553.
- 65.Pérez-Hernández G., Paul F., Giorgino T., De Fabritiis G., Noé F. Identification of slow molecular order parameters for Markov model construction. J. Chem. Phys. 2013;139 doi: 10.1063/1.4811489.
- 66.Metzner P., Schütte C., Vanden-Eijnden E. Transition path theory for Markov jump processes. Multiscale Model. Simul. 2009;7:1192–1219.
- 67.Weinan E., Vanden-Eijnden E. Towards a theory of transition paths. J. Stat. Phys. 2006;123:503.
- 68.Meng Y., Shukla D., Pande V.S., Roux B. Transition path theory analysis of c-Src kinase activation. Proc. Natl. Acad. Sci. U. S. A. 2016;113:9193–9198. doi: 10.1073/pnas.1602790113.
- 69.Noé F., Schütte C., Vanden-Eijnden E., Reich L., Weikl T.R. Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations. Proc. Natl. Acad. Sci. U. S. A. 2009;106:19011–19016. doi: 10.1073/pnas.0905466106.
- 70.Roe D.R., Cheatham T.E., III 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.
- 71.McGibbon R.T., Beauchamp K.A., Harrigan M.P., Klein C., Swails J.M., Hernández C.X., Schwantes C.R., Wang L.-P., Lane T.J., Pande V.S. MDTraj: A modern open library for the analysis of molecular dynamics trajectories. Biophys. J. 2015;109:1528–1532. doi: 10.1016/j.bpj.2015.08.015.
- 72.Humphrey W., Dalke A., Schulten K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996;14:33–38. doi: 10.1016/0263-7855(96)00018-5.
- 73.Mey A.S., Allen B.K., Macdonald H.E.B., Chodera J.D., Hahn D.F., Kuhn M., Michel J., Mobley D.L., Naden L.N., Prasad S., Rizzi A., Scheen J., Shirts M.R., Tresadern G., Xu H. Best practices for alchemical free energy calculations [article v1.0] Living J. Comput. Mol. Sci. 2020;2:18378. doi: 10.33011/livecoms.2.1.18378.
- 74.Hernández C.X., Pande V.S. MDEntropy: Information-theoretic analyses for molecular dynamics. J. Open Source Softw. 2017;2:427.
- 75.Hagberg A., Swart P., Chult D.S. Exploring Network Structure, Dynamics, and Function Using Networkx. Proceedings of the 7th Python in Science Conference (SciPy2008); Pasadena, CA: 2008.
- 76.Dijkstra E.W. A note on two problems in connexion with graphs. Numer. Math. 1959;1:269–271.
- 77.Trott O., Olson A.J. AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2010;31:455–461. doi: 10.1002/jcc.21334.
- 78.Sultan M.M., Kiss G., Pande V.S. Towards simple kinetic models of functional dynamics for a kinase subfamily. Nat. Chem. 2018;10:903–909. doi: 10.1038/s41557-018-0077-9.
Associated Data
Supplementary Materials
Data Availability Statement
MSM feature matrix, final MSM object, codes, and pdb files used to generate figures can be obtained from https://github.com/ShuklaGroup/THC_binding.git.