João Marcelo Lamim Ribeiro
Marta Filizola- Department of Pharmacological Sciences, Icahn School of Medicine at Mount Sinai, New York, NY, United States
Effective treatments for pain management remain elusive due to the dangerous side-effects of current gold-standard opioid analgesics, including the respiratory depression that has led to skyrocketing death rates from opioid overdoses over the past decade. In an attempt to address the horrific opioid crisis worldwide, the National Institute on Drug Abuse has recently proposed boosting research on specific pharmacological mechanisms mediated by a number of G protein-coupled receptors (GPCRs). This research is expected to expedite the discovery of medications for opioid overdose and opioid use disorders, leading toward a safer and more effective treatment of pain. Here, we review mechanistic insights from recent all-atom molecular dynamics simulations of a specific subset of GPCRs for which high-resolution experimental structures are available, including opioid, cannabinoid, orexin, metabotropic glutamate, and dopamine receptor subtypes.
Introduction
Pain is a vital, albeit unpleasant, physiological response to tissue damage, but it can become a disease if it strikes in the absence of tissue injury, or continues long after appropriate tissue healing (Ringkamp et al., 2018). As a disease, pain poses an enormous socioeconomic burden on the people who suffer from it, as well as a huge financial strain worldwide. There are several different ways to categorize pain (e.g., chronic, nociceptive, neuropathic, etc.) and treatment decisions depend on the specific type of pain (Chang et al., 2015). For severe and chronic pain, the gold-standard painkillers remain opioid drugs, despite their dangerous side effects and abuse liability.
Overprescription of opioid analgesics in the nineties led to drug misuse, and the consequent “opioid epidemic” or “opioid crisis” in the United States, which has most recently expanded to heroin and other illicit synthetic opioids such as fentanyl and its analogs (Volkow et al., 2019). With an average of 130 Americans dying every day (Centers for Disease Control and Prevention [CDC], 2017), new scientific solutions are desperately needed to effectively manage pain while preventing or treating overdose and opioid use disorder (OUD) manifestations. This recognition recently led the leadership of the National Institutes of Health (NIH) and the National Institutes on Drug Abuse (NIDA) to launch initiatives aimed at accelerating the pace of scientific inquiry that is necessary to address the opioid crisis. One of these initiatives enabled the prioritization of specific mechanisms and pharmacological targets whose study is expected to boost the development of novel drugs that have the highest probability of approval by the Food and Drug Administration (FDA) for the treatment of opioid overdose and OUD (Rasmussen et al., 2019). These “most wanted” mechanisms and targets (Rasmussen et al., 2019), which include several G protein-coupled receptors (GPCRs), were established based on published data and internal studies that the NIDA leadership deemed most promising for the development of improved therapeutics for OUDs.
In the classical view of GPCR-mediated downstream cellular signaling, the receptor transitions into active conformational states which are capable of recruiting and ultimately activating intracellular protein transducers such a G-proteins and β-arrestins. These active states are characterized by specific conformational changes at the intracellular end of the receptor, most notably exemplified by a different extent of outward movement of transmembrane helix 6 (TM6) away from TM3 (e.g., see experimentally determined inactive and active structures of a prototypic GPCR compared to intermediate states in Figure 1). Typically, GPCR activation is mediated by the binding of agonist ligands at the so-called orthosteric binding site, which is the same site where endogenous ligands bind. Antagonist and inverse agonist ligands, on the other hand, shift the conformational equilibrium toward inactive conformations of the receptor while partial agonists are expected to stabilize intermediate conformations between inactive and active states of the receptor. For years, drug design at GPCRs has mostly been focused on optimizing ligands for the receptor orthosteric site. However, by binding non-conserved regions of the receptor and directly affecting the binding and/or efficacy of orthosteric ligands, so-called positive and negative allosteric modulators (PAMs and NAMs, respectively) are receiving more and more attention for the development of improved therapeutics targeting GPCRs. Similarly, so-called biased agonists hold a great potential for drug discovery since they would stabilize receptor conformations that selectively recruit an intracellular protein instead of another, thereby triggering specific biological effects. Figure 2 provides a cartoon depiction of the expected effect of the different types of ligands on the receptor.
Figure 1. A comparison of the representative conformation of the most probable metastable state within an intermediate region of (A) the morphine-bound MOP receptor, where the intermediate state is in light blue, and (B) the TRV-130-bound MOP receptor, where the intermediate state is in light purple, relative to the experimentally determined MOP receptor inactive and active states (light and dark gray, respectively). Note that the most dramatic differences between these conformations stem from the extent of outward movement of TM6 away from TM3, which is one of the most notable conformational changes that has been associated with receptor activation. Images on the right correspond to a 90° rotation of the receptor helical bundle, and represent the view from the intracellular domain.
Figure 2. A cartoon representation of the effects of different ligands on a GPCR. The cooperativity between allosteric and orthosteric ligands can shift their affinities for the receptor, and/or bias the GPCR coupling to a particular intracellular partner.
Indeed, among the NIDA’s ten most wanted medication development priorities in response to the opioid crisis (Rasmussen et al., 2019) are agonists, antagonists, partial agonists, PAMs, and/or NAMs at a number of GPCRs, including orexin-1 or 2, kappa opioid (KOP), GABAB, muscarinic M5, nociceptin opioid peptide (NOP), metabotropic glutamate 2/3, ghrelin, dopamine D3, and cannabinoid CB-1 receptors. Additional NIDA-designated priority medications mediated by GPCRs (Rasmussen et al., 2019) included: (i) serotonin 5HT2C agonists or PAMs, with or without 5HT2A antagonist/NAM activity, (ii) biased μ-opioid (MOP) receptor agonists or PAMs, and (iii) NOP/MOP bifunctional agonists or PAMs.
One of the main obstacles to the development of new therapeutics for pain management or to treat or prevent opioid overdose or OUDs is the limited understanding of the relevant signal transduction mechanisms at the atomic level notwithstanding the recognized role of a number of GPCRs in the regulation of pain transmission and OUD manifestations, as well as the availability of high-resolution experimental structures for several of these GPCRs. Molecular dynamics (MD) simulations can provide a complementary perspective (Ribeiro and Filizola, 2019) on the molecular determinants underlying GPCR-mediated signaling mechanisms involved in pain transmission, respiratory depression, or clinical manifestations of OUD. Availability of more powerful hardware and software has made the use of MD simulations more affordable and available to a larger number of scientists. It is now straightforward for numerous groups to simulate timescales in the microsecond (μs) regime using high-performance computing resources accessible to a large number of academic institutions. Using either standard or enhanced MD simulations, the latter to access even longer timescales or a more extensive sampling, GPCRs are studied in terms of an ensemble of conformations between fully active and inactive states, with a number of factors, such as binding of ligands, lipids, ions, receptors, or intracellular proteins, shifting the equilibrium toward different states. While we refer the reader elsewhere for an overview of strengths and limitations of MD simulations in their application to GPCRs (e.g., Ribeiro and Filizola, 2019), we review here atomically detailed mechanistic insights from MD simulations of high-resolution experimental structures of a number of GPCR subtypes whose study might lead to faster development of medications for the treatment of pain or OUDs (see Table 1 for a summary of all the MD studies reported herein). These GPCRs include opioid, cannabinoid, orexin, metabotropic glutamate, and dopamine receptor subtypes regulating distinctive pharmacological mechanisms. The position of the co-crystallized ligands in the respective high-resolution experimental structures used as a starting point for the MD simulations referenced herein is shown in Figure 3.
Table 1. A compilation of the MD-based studies that have been reported in this review article.
Figure 3. The experimentally determined high-resolution GPCR structures, together with their bound ligands, used in the MD-based studies discussed in this review. Nanobodies and other interacting proteins were removed. PDB 4dkl, The antagonist β-FNA bound to the MOP receptor; PDB 5c1m, The agonist BU72 bound to the MOP receptor; PDB 4djh, The antagonist JDTic bound to the KOP receptor; PDB 6b73, The agonist MP1104 bound to the KOP receptor; PDB 4ea3, The peptide mimetic antagonist compound 24 bound to the NOP receptor; PDB 5tgz, The antagonist AM6538 bound to the CB1 receptor; PDB 5xra, The agonist AM11542 bound to the CB1 receptor; PDB 5u09, The inverse agonist taranabant bound to the CB1 receptor; PDB 4s0v, The antagonist suvorexant bound to the OX2 receptor; PDB 3pbl, The antagonist eticlopride bound to the D3 receptor; PDB 4or2, The negative allosteric modulator FITM bound to the transmembrane domain of mGluR1; PDB 4oo9, The negative allosteric modulator mavoglurant bound to the transmembrane domain of mGluR5; PDB 5cgc, The negative allosteric modulator 3-chloro-4-fluoro-5-[6-(1H-pyrazol-1-yl)pyrimidin-4-yl]benzonitrile bound to the transmembrane domain of mGluR5; PDB 5cgd, The negative allosteric modulator 3-chloro-5-[6-(5-fluoropyridin-2-yl)pyrimidin-4-yl]benzonitrile bound to the transmembrane domain of mGluR5.
Synopsis of MD Simulation Methods Cited Herein
The goal of this section is not to describe in detail the MD simulation methods and tools cited in this review, but rather provide a lay summary for the general audience. Interested readers are referred to the appropriate reviews, cited below, where the methods are described more thoroughly.
The aim of an MD simulation is to provide the time-evolution of a system by solving the appropriate equation(s) of motion. In these equations, the energy interactions between the particles of the system under study must be described. In principle, atomic-level interactions should be handled using quantum mechanics; however, due to the size of typical biological systems, it is often unfeasible to use the fully quantum description, and classical mechanics is used instead (Oren et al., 2001). The typical approach to modeling these interactions is to describe bonded and non-bonded atomic interactions by simple expressions, with different parameters in these expressions representing different atom types (Oren et al., 2001; Ponder and Case, 2003; Lopes et al., 2015; Nerenberg and Head-Gordon, 2018). The determination of an accurate set of parameters for use in these expressions (the so-called force field) is key to properly describe atomic interactions within biological systems, and it is therefore an intensive area of research. In the following sections, we will mention several different force fields that are currently available to the MD practitioner and have been used in the studies reported here, including those corresponding to the names of Assisted Model Building and Energy Refinement (AMBER) (Maier et al., 2015), Chemistry at Harvard Macromolecular Mechanics (CHARMM) (Best et al., 2012), General Amber force field (GAFF) (Wang et al., 2004), CHARMM General Force Field (CGenFF) (Vanommeslaeghe et al., 2010), etc. In addition, tools such as General Automated Atomic Model Parametrization (GAAMP) (Huang and Roux, 2013) have been developed to automatically generate force field parameters for small molecules not accurately described by the aforementioned force fields.
One of the major obstacles in using MD simulations for investigating biological problems is that the timescale for sampling the event of interest is often larger than the times that can be simulated (Bernardi et al., 2015). While the microsecond (μs) regime is nowadays accessible to a large number of MD practitioners, most biological events fall above that threshold (Valsson et al., 2016; Ribeiro et al., 2019). Enhanced sampling methods are designed to provide a faster exploration of the conformational space of the system under study. In this review, we report on studies carried out using two classes of enhanced sampling MD methods. One of these classes, exemplified by metadynamics and Gaussian accelerated molecular dynamics (GaMD) simulations, uses an artificial biasing force to speed up the rate at which the process of interest is sampled, and so long as this is done in a careful manner, the effect of the “bias” can be “reweighed” to recover “unbiased” information. The other class of enhanced sampling MD simulations is exemplified by adaptive sampling protocols, in which successive batches of simulations are started from regions of conformational space that have been undersampled, thus accelerating the sampling of important, but slow, events (Husic and Pande, 2018; Ribeiro et al., 2019). Throughout the remainder of this review we refer to unbiased MD as simulations in which trajectories are propagated without the help of enhanced sampling techniques, although adaptive sampling techniques are, in principle, not biased (Husic and Pande, 2018).
Opioid Receptors
Overdose deaths by prescription, illegal, or synthetic opioids have mostly been attributed to the activation of the MOP receptor, a rhodopsin-like (class A) GPCR located, in part, on brainstem neurons that control respiration. In an attempt to develop improved opioid therapeutics with limited respiratory depression and other unwanted side effects (Janecka et al., 2019), attention has recently shifted to G protein-biased agonists of the MOP receptor. These MOP receptor ligands are believed to produce anti-nociceptive action by stabilizing a receptor conformation that preferentially activates G-protein over β-arrestin, the latter shown to be linked to unwanted side effects (Bohn et al., 1999; Raehal et al., 2005).
Recent MD simulations have been leveraged to reveal the molecular details behind G-protein biased agonism at the MOP receptor (Schneider et al., 2016; Kapoor et al., 2017; Cheng et al., 2018). In particular, oliceridine, also known as TRV-130, a G protein-biased MOP receptor ligand that reached phase III clinical trials for management of moderate to severe pain (Viscusi et al., 2019), has been the subject of a number of MD simulations (Schneider et al., 2016; Kapoor et al., 2017; Cheng et al., 2018). Our group, for instance, investigated the binding of TRV-130 from the bulk solvent to the MOP receptor, as well as its preferred mode of interaction at the crystallographically identified orthosteric binding site, using ∼44 μs of unbiased all-atom MD simulations (Schneider et al., 2016). These MD simulations had the MOP receptor placed in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC)/cholesterol lipid membrane environment and used the CHARMM36 force-field to represent the protein and lipid molecules, and CGenFF for the TRV-130 ligand. The results of these simulations suggested that intermediate binding states of TRV-130 at the so-called vestibule region of the MOP receptor directed ligand access to the orthosteric site, and that two energetically indistinguishable conformations could be adopted in the orthosteric binding pocket. Additional microsecond-scale, unbiased simulations of the MOP receptor bound to TRV-130 or the classical orthosteric opioid drug morphine (Schneider et al., 2016) suggested differences in the allosteric coupling between the MOP receptor orthosteric site and the receptor intracellular region induced by the two different ligands. Notably, we found that residues in direct or water-mediated contact with either ligand did not exhibit a main role in the communication between the orthosteric binding site and the intracellular region of the MOP receptor, notwithstanding their contribution to stable ligand binding at the orthosteric pocket. In addition, unlike the morphine-bound receptor, in which the most contributing residues to the allosteric coupling between the orthosteric binding site and the intracellular region of the MOP receptor resided in both transmembrane (TM) helices TM3 and TM6, the TRV-130 complex did not have strong contributors to the co-information in TM6 (Schneider et al., 2016).
To obtain a more thorough investigation of the molecular details of ligand-induced MOP receptor activation, we recently built a Markov state model (MSM) using over 400 μs of MD simulations of the MOP receptor embedded in a POPC/cholesterol membrane mimetic environment with either morphine or TRV-130 bound at the orthosteric binding site (Kapoor et al., 2017). Here, the CHARMM36 and CGenFF force-fields were also used. The MSM revealed that the conformational landscape of the MOP receptor in complex with either ligand contained several kinetic macrostates (i.e., metastable states) in addition to those corresponding to crystal-like active or inactive conformations of the receptor, defining two different intermediate regions of the conformational space for each ligand-MOP complex. These regions contained different conformational states stabilized by morphine or TRV-130, which may or may not get ever resolved experimentally and yet be useful for the rational design of improved opioids with reduced side effects. Shown as an example in Figure 1 are representative conformations of the most probable metastable states within the intermediate regions available to the simulated morphine-bound and TRV130-bound MOP complexes compared to active or inactive crystallographic states of MOP, which are characterized by a different extent of TM6 outward movement. Another important observation of this MSM analysis was the existence of different activation/deactivation pathways induced by the classical or G protein-biased opioid ligand, which confirmed the substantial difference in the receptor dynamics induced by the two different ligands.
In a recent investigation, unbiased MD simulations were used to study the MOP receptor in a ligand-free form, as well as in complex with TRV-130, the agonist BU72, the antagonist naltrexone (NTX), or the antagonist β-FNA. These simulations were each run for 300 ns with the MOP receptor embedded in a POPC membrane environment, the ligands parameters derived from CGenFF, and both the lipid and receptor molecules described by CHARMM36 force field parameters. Notably, simulations of the TRV-130-MOP receptor complex drew attention to two residues in TM6 and TM7, specifically, Y3267.43 and W2936.48, which had been shown to be important for MOP receptor biased signaling by mutagenesis studies. Superscript residue numbers here and throughout the text refer to the Ballesteros–Weinstein generic numbering scheme (Ballesteros and Weinstein, 1995) wherein the first digit corresponds to the transmembrane helix number and the second digit is a sequence number relative to the most conserved residue in a helix, which is assigned a value of 50. However, corrections to this numbering scheme incorporating structural information have been proposed (Isberg et al., 2015) and will be reported in parenthesis for those residues whose numbering may diverge from Ballesteros-Weinstein’s, such as Y3267.43 (renumbered Y3267.42 by Isberg et al., 2015).
Similar to MOP receptor agonists, centrally acting KOP receptor agonists can be effective in the treatment of pain, but their dysphoric and hallucinogenic side effects have limited their clinical usefulness (Land et al., 2008), shifting focus to the development of peripherally restricted KOP agonists as analgesics with reduced abuse liability (Hasebe et al., 2004) or KOP antagonists for the treatment of substance use disorders (Carlezon and Krystal, 2016). The structural basis of agonism or antagonism at the MOP and KOP receptors has recently been studied using unbiased all-atom MD simulations (Yuan et al., 2015). A total of four ligand-opioid receptor complexes embedded in a POPC membrane environment were simulated, including the KOP receptor in complex with the JDTic antagonist, the MOP receptor complexed with the agonist morphine, and either the MOP or KOP receptor in complex with levallorphan, a morphinan ligand acting as an antagonist at the MOP receptor and an agonist at the KOP receptor (Yuan et al., 2015). The simulations – each 3 μs in length – made use of CGenFF in their description of the ligands, and the CHARMM36 force field for all remaining molecules. In these simulations, the authors found that the amount of water penetration into the interior of the receptors, which is a known characteristic of GPCR activation, was higher when the receptor was complexed with an agonist as opposed to an antagonist (Yuan et al., 2015). In particular, the levallorphan-MOP and JDTic-KOP complexes formed a σ – π stacking interaction with the Y3207.43 (Y3207.42 as per Isberg et al., 2015) residue, which tended to block water penetration into the interior. Solvent accessible surface area calculations on subsequent short MD simulations of several other agonists or antagonists in complex with either the KOP or MOP receptors showed these values were higher for receptors in complex with agonists as opposed to antagonists (Yuan et al., 2015).
The conformational changes induced by 6′-Guanidinonaltrindole (6′-GNTI), a G-protein biased agonist that is selective for the KOP receptor, or the antagonist 5′-Guanidinonaltrindole (5′-GNTI) have recently been studied using unbiased all-atom MD simulations (Cheng et al., 2016). In this work, ∼600 ns MD simulations were performed on the ligand-free KOP receptor, as well as the receptor in complex with either 5′-GNTI or 6′-GNTI, with each system embedded in an explicit POPC membrane environment (Cheng et al., 2016). The MD simulations of the KOP receptor bound to the antagonist 5′-GNTI drew attention to the hydrogen bond between S3247.47 and V691.42 as the basis for the stabilization of the kink angle on TM7 at about 150°, and possibly deriving antagonistic activity. In contrast, the MD simulation of the G-protein biased agonist 6′-GNTI bound to the KOP receptor showed a different value for this kink angle, and highlighted an interaction of the ligand guanidinium group with the E2976.58 residue, together with the steric effect from I2946.55, as key contributors to the rotation of TM6, a known hallmark of GPCR activation (Cheng et al., 2016). The possible absence of guanidinium-E2976.58 interaction in the MOP or the δ-opioid receptor (DOP) receptor due to this residue replacement by a lysine or tryptophan, respectively, was interpreted as the basis for the lack of 6′-GNTI agonism in those opioid receptor subtypes.
More recently, the KOP receptor conformational changes induced by the agonist MP1104 or the antagonist JDTic were investigated using enhanced sampling MD simulations (An et al., 2018). Specifically, using the generalized AMBER force field for the ligands and the AMBER ff14SB force field for the protein, the GaMD method was used to enhance the sampling of long-time, large-scale conformational rearrangements associated with KOP receptor activation by introducing a biasing harmonic potential on certain dihedral angles. The following systems were all simulated in an explicit POPC membrane environment: the ligand-free KOP receptor in an inactive or active conformation, the latter with or without an intracellular protein, the JDTic-inactive KOP receptor complex, and the MP1104-active KOP receptor complex with or without a stabilizing intracellular partner (An et al., 2018). Taken together, the results of these simulations showed that while the agonist stabilized specific functional domains in an active-like conformation, the antagonist shifted the conformational equilibrium toward an inactive conformation. Notably, the inactive ligand-free state of the KOP receptor was the most stable one, in contrast to the ligand-free active form of the receptor, which readily transitioned to an intermediate state characterized by a reduced TM6 outward movement (An et al., 2018). Finally, the results revealed a hydrophobic interaction between Y2465.58 and TM6 in the intermediate metastable state that hindered the transition between the inactive and active conformations of the KOP receptor (An et al., 2018).
The NOP receptor is another opioid target of interest for powerful pain relief with reduced side effects. MD simulations using the AMBER ff99sb force field were recently used to investigate the molecular effect of the novel analgesic cebranopadol (CBP) – which acts as an agonist at both NOP and MOP receptors – on the NOP receptor (Della Longa and Arcovito, 2019). These simulations, run in an explicit POPC membrane environment, used the high-resolution structure of the NOP receptor in complex with the antagonist C24 as a starting point for 1 μs-long simulations of the ligand-free NOP receptor, the C24-NOP receptor complex, and the CBP-NOP receptor complex (Della Longa and Arcovito, 2019). In all cases, the simulations did not sample the large amplitude motions of the transmembrane helices associated with receptor activation even in the presence of the agonist CBP (Della Longa and Arcovito, 2019). The CBP ligand did, however, destabilize the inactive NOP receptor conformation relative to both the ligand-free NOP receptor and the C24-NOP receptor complex such that the NOP receptor bound to CBP could sample a much wider region of the local conformational space (Della Longa and Arcovito, 2019). The authors used these MD simulations to determine some of the earliest microswitches that lead to destabilizing the initial inactive conformation. A histogram of the conformational space of the M1343.36 and W2766.48 residues located in the orthosteric site revealed that a conformational switch to their active-like positions was accessible to the agonist-receptor complex (Della Longa and Arcovito, 2019). In addition, the time evolution of the conformation of these residues showed that in the agonist-receptor complex these active-like states were “locked” in place for the remainder of the simulation (Della Longa and Arcovito, 2019).
Allosteric modulators of opioid receptors (Remesic et al., 2017) constitute another priority area of research with expected higher probability of success in the development of medications in response to the opioid crisis (Rasmussen et al., 2019). NIDA’s “most wanted” allosteric modulators of opioid receptors include MOP PAMs (Rasmussen et al., 2019). One of the reasons why MOP PAMs are of potential interest is that by increasing the potency and/or efficacy of classical opioid drugs, they are expected to produce the same analgesic response achieved by higher doses of opioid drug while simultaneously presenting fewer on-target overdosing risks. Most importantly, these compounds may not be subject to the compensatory mechanisms deriving from chronic MOP activation (e.g., tolerance, dependence, and increased toxicity) because they preserve the temporal and spatial fidelity of signaling in vivo by acting only in the presence of endogenous or other orthosteric ligands (Remesic et al., 2017).
Since experimental high-resolution structures of opioid receptors in complex with allosteric modulators are yet to be published, and automated docking protocols do not yield single binding poses that can be clearly distinguished from the rest, MD simulations can make valuable contributions toward locating allosteric binding sites in opioid receptors, as well as revealing the molecular basis for their binding modes, as we recently demonstrated in an application to the DOP receptor. Specifically, we used metadynamics to simulate the binding of a recently discovered allosteric modulator BMS-986187 of opioid receptors (Burford et al., 2015; Livingston et al., 2018) to the DOP receptor in complex with the orthosteric ligand SNC-80 (Shang et al., 2016). The simulations identified the two most stable binding modes with near-degenerate energies that were discriminated experimentally based on functional studies of normal and mutant receptors (Shang et al., 2016). Figure 4 summarizes the essence of this integrated computational-experimental work, which gave support to the BMS-986187 binding pose in cyan color as the most likely to occur based on the impact of specific mutations (e.g., L/W3007.35) on either the intrinsic binding affinity of the PAM or the affinity/efficacy of the orthosteric ligand.
Figure 4. The two energetically most favorable binding modes of the allosteric modulator BMS-986187 – colored yellow and cyan – at the DOP receptor in complex with the orthosteric ligand SNC-80, which is colored in gray. In sticks are the DOP receptor residues in the putative DOP allosteric site used in mutagenesis experiments to help validate the most likely binding mode. The plots to the right show the effect of mutations on the affinity of the allosteric ligand (top right), the induced changes in affinity of the orthosteric ligand (middle plot), and the changes in orthosteric ligand efficacy when the allosteric ligand is bound (bottom right). The experimental data support the predicted BMS-986187 binding mode in cyan color as the most favorable one. Astars indicate the statistical significance levels as given by Dunnet’s test p values (∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001).
The structural basis for the effect that another allosteric modulator, BMS-986122, has on MOP in complex with either the partial agonist buprenorphine or the agonist methadone, was recently investigated using μs-scale unbiased MD simulations in an explicit membrane mimetic environment (Bartuzi et al., 2019). The results suggested that specific dynamic movements that are characteristic of full receptor activation, such as for instance the bending and rotation of TM7, can be induced by the allosteric modulator even in the presence of a partial agonist at the orthosteric binding site (Bartuzi et al., 2019).
Cannabinoid Receptors
The cannabinoid receptor 1 (CB1) is another important class A GPCR drug target for the development of new analgesics with reduced side effects. For instance, PAMs of this receptor have been shown to suppress pathological pain without producing tolerance or dependence (Slivicki et al., 2018). These properties, alongside their potentially reduced psychoactive side effects due to their lack of intrinsic activity and inherent ceiling efficacy, make CB1 PAMs potentially better therapeutics for inflammation and chronic pain compared to CB1 orthosteric agonists (Wootten et al., 2013). CB1 antagonists have also been reported to exert analgesia in animal models of inflammatory arthritis and hyperalgesia with reduced side effects (Croci and Zarini, 2007; Ueda et al., 2014). Based on these insights, the CB1 receptor is a target of interest for the development of improved therapeutics to combat the opioid crisis (Rasmussen et al., 2019). Luckily, a number of high-resolution experimental structures for either inactive or active CB1 receptors have been made available (Hua et al., 2016, 2017; Shao et al., 2016; Krishna Kumar et al., 2019) and can be used to study CB1-mediated functional mechanisms at a molecular level for the purpose of guiding rational drug discovery.
A recent unbiased MD simulation-based investigation used these structures in order to gain insight into the features that could explain the different efficacy profiles of partial agonists, antagonists, and inverse agonists bound to the CB1 receptor, and could eventually be used for the design of novel therapeutics targeting the CB1 (Jung et al., 2018). Specifically, analysis of these MD simulations made it possible to discriminate the dynamic tendencies of inactive and active CB1 structures in the presence of ligands with different efficacies while the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method was used to assess the contribution of individual ligand-receptor interactions to the binding of the partial agonist Δ9-tetrahydrocannabinol (THC), the antagonist Δ9-tetrahydrocannabivarin (THCV), and the inverse agonist taranabant from ligand binding free energy decompositions of the CB1-ligand complexes. These MD simulations revealed that binding of the inverse agonist to the active CB1 receptor structure made TM1 less rigid, leading to larger root-mean-square deviations of the possible contacts between TM1 and 2, as well as TM1 and TM7, compared to either the partial agonist-bound or the antagonist-bound active CB1 receptor complex. In addition, the simulations drew attention to large conformational changes involving residues Phe2003.36 and Trp3566.48 in the orthosteric binding site. In the inactive CB1 receptor, the inverse agonist taranabant – through a direct interaction with Trp3566.48 – stabilized the conformation of Trp3566.48 and Phe2003.36 with respect to one another while the partial agonist THC and antagonist THCV did not. In contrast, in the CB1 active conformation, taranabant induced a different dynamic behavior for the interaction of Trp3566.48 and Phe2003.36 compared to the partial agonist THC. Furthermore, changes in the binding free energies showed that the partial agonist THC preferred the CB1 active conformation, whereas the simulated inverse agonist taranabant remained more favorably bound to the CB1 inactive conformation via a stable interaction with residue Trp3566.48 (Jung et al., 2018) during the afforded simulation timescale.
In another recent investigation, biased MD simulations were used to probe the binding sites and modes of the CP 55,940 agonist and the GAT228 mixed agonist/PAM (so-called Ago-PAM) to the CB1 receptor (Saleh et al., 2018). Ligand binding events to the CB1 receptor were enhanced using the multiple-walker metadynamics biasing protocol (Saleh et al., 2017) and a funnel-shaped restraint applied to the ligand in the bulk, both for the purpose of aiding convergence (Saleh et al., 2018). The simulation of CP 55,940 binding to the ligand-free CB1 receptor – run for 2 μs to achieve convergence – showed a single, deep minimum along the binding potential of mean force (PMF) (Saleh et al., 2018). The location of this deep minimum corresponded to the orthosteric binding site in the high-resolution structure of the CB1 receptor (Saleh et al., 2018). The binding mode of CP 55,940 was found to reproduce all of the interactions observed in the high-resolution structure of the THC agonist AM11542 bound to the CB1 receptor, except for the interaction with F1742.61, which was replaced by an interaction with residue F102 in the N-terminal region of the receptor (Saleh et al., 2018). In contrast, the binding simulations of the GAT228 Ago-PAM to the ligand-free CB1 receptor, showed two PMF minima corresponding to binding at different sites of the receptor. These two minima had similar affinities (Saleh et al., 2018), which suggests an equilibrium between binding at the two different receptor sites, thus providing structural context to the experimentally observed partial agonistic effect of GAT228 (Saleh et al., 2018). While the PMF global minimum corresponded to GAT228 bound to the orthosteric site via a cluster of hydrophobic interactions with Val1963.32, Leu1933.29, and seven additional Phe residues (Saleh et al., 2018), the other PMF minimum defined a putative CB1 receptor allosteric site (Saleh et al., 2018). Notably, simulations of the binding of GAT228 to the CP 55,940-CB1 receptor complex revealed a 3 Å RMSD displacement of the CP 55,940 binding mode induced by GAT228 binding preferentially at an allosteric site defined by residues W2795.43, Y2755.39, W3566.48, and the N-terminus F268 (Saleh et al., 2018), through a hydrogen bond between the indole hydrogen atom of GAT228 and T1973.33 (Saleh et al., 2018).
Orexin Receptors
The orexin (OX) 1 and 2 receptors, expressed throughout the CNS, are neuropeptide receptors that belong to the β-branch of the rhodopsin-like GPCRs (Yin et al., 2016). Although these receptors are known to be important in regulating mammalian sleep patterns (Yin et al., 2015, 2016; Wacker and Roth, 2016), they have recently received attention in the development of therapeutics to address the opioid crisis (Rasmussen et al., 2019). Although high-resolution experimental structures exist for both the OX1 and OX2 receptors bound to antagonists (Yin et al., 2015, 2016; Suno et al., 2018), recent MD-based studies have focused on the OX2 receptor (Bai et al., 2018; Karhu et al., 2019). The earliest of these studies used 200 ns of unbiased MD simulations of the antagonist suvorexant bound to the OX2 wild-type and N3246.55A mutant receptors embedded in a POPC lipid environment to understand the dynamic interplay between the horseshoe shape pocket of the receptor revealed by crystallography (Yin et al., 2015) and the boat conformation of the ligand at an atomic level of detail (Bai et al., 2018). In line with the notion of a loss of antagonist binding ability and signaling response in the N3246.55A mutant, the results of these simulations showed a distorted horseshoe shape pocket of the N3246.55A mutant of the OX2 receptor compared to wild-type receptor, suggesting that an intact horseshoe shape pocket is required for optimal suvorexant binding and antagonistic activity (Bai et al., 2018).
Molecular determinants of OX2 receptor binding and activation were further investigated in a recent MD-based work (Karhu et al., 2019) focused on comparing the receptor dynamic behavior induced by the agonist Nag26 or the antagonist suvorexant, in addition to predicting the mode of binding of the endogenous ligand orexin-A at the OX2 receptor (Karhu et al., 2019). The microsecond-long unbiased MD simulations of Nag26 or suvorexant bound to the OX2 receptor revealed very different dynamic behaviors between the agonist and antagonist, with the agonist exhibiting much increased flexibility and completely different interaction patterns (Karhu et al., 2019). In particular, while suvorexant induced stabilization of the Q1343.32–Y3547.43 (renumbered Y3547.42 according to Isberg et al., 2015) hydrogen bond, Nag26 promoted counterclockwise rotation of the TM5 extracellular end, influencing the interactions among TM4, 5, and 6 (Karhu et al., 2019).
Dopamine Receptors
The dopamine D3 receptor has also received attention as a drug target for mitigating the opioid crisis (Rasmussen et al., 2019) in large part because of its potential for opioid dependence treatment (Kumar et al., 2016). Recent work on this receptor highlighted the importance of using MD simulations to predict ligand binding at the D3 receptor in agreement with inferences from site-directed mutagenesis (Ferruz et al., 2018). In particular, binding of the antagonist PF-4363467 from the bulk to the dopamine D3 receptor was simulated with ACEMD (Harvey et al., 2009) directed sampling and MSMs generated using High-Throughput Molecular Dynamics (HTMD), using the CHARMM36 force field for the protein and POPC lipid, and ligand force field parameters generated with the GAAMP tool within HTMD. An adaptive sampling protocol was used for these simulations, according to which MSMs were built from successive batches of simulations to identify starting conformations for the next batch, thus affording thorough exploration of the conformational landscape without biasing the potential. A total of over 680 μs of simulation was carried out, resulting in the sampling of two binding paths, which differed in the presence of a second intermediate state in the minor binding path. This intermediate state corresponded to PF-4363467 bound to the D3 receptor at a position between the extracellular vestibule and the orthosteric site. The main structural difference between the predicted PF-4363467-D3 receptor complex and the crystallographically determined eticlopride-D3 receptor complex (Chien et al., 2010) was the formation of an aromatic cryptic pocket between TM5 and TM6 involving residues F3386.44, W3426.48, L3436.49, F3456.51, and F3466.52 deriving from the displacement of residues F1975.47 and F3466.52 (Ferruz et al., 2018).
The SB269652 ligand is a bitopic D2 and D3 receptor ligand with negative allosteric modulation activity (Silvano et al., 2010) that has received research attention for the treatment of drug abuse (Verma et al., 2018). By nature of being bitopic, SB269652 binds to both the orthosteric binding site and an allosteric binding site in the receptor. Previous molecular modeling studies suggested that the tetrahydroisoquinoline (THIQ) moiety of SB269652 binds to the orthosteric binding site via an ionic interaction with D3.32 whereas the indole-2-carboxamide moiety of SB269652 protruded into a putative allosteric site between TM2 and TM7, forming a hydrogen bond with the E2.65 (renumbered E2.64 according to Isberg et al., 2015) residue (Guo et al., 2008; Lane et al., 2014). However, mutagenesis and structure activity relationship studies of SB269652 questioned that this hydrogen bond alone could determine the compound allosteric properties (Lane et al., 2014; Mistry et al., 2015; Shonberg et al., 2015), calling for an in depth dynamics study. Thus, adaptive sampling MD simulations using MSMs were recently used to obtain mechanistic insights into the role of the E2.65 residue in the binding and allosteric properties of SB269652 at both the D2 and D3 receptors (Verma et al., 2018). Specifically, simulations were carried out for the ligand bound to the wild-type D2 receptor, the E2.65A D2 receptor mutant, the wild-type D3 receptor, or the E2.65A D3 receptor – all of which were embedded in a POPC lipid environment – for a total of 76.5 μs, with the simulation time of each complex ranging between 15.9 and 21.3 μs (Verma et al., 2018). The THIQ moiety of SB269652 bound at the orthosteric site was shown to be quite stable in both the D2 and D3 receptors, although subtle differences in its binding poses were observed, due in large part to different interactions between the ligand and the extracellular loop 2 (Verma et al., 2018). This is in agreement with previous chimera mutagenesis results that showed the affinities for the D2 and D3 receptors were different in large part due to the extracellular loop 2 (ECL2) (Silvano et al., 2010). In contrast, the indole-2-carboxamide moiety of SB269652 bound at the allosteric site was shown to undergo significant fluctuations, with the MSM analysis revealing two equiprobable metastable states in the wild-type D2 and D3 receptors, but three different metastable states in both the E2.65A D2 and D3 receptor mutants (Verma et al., 2018). Furthermore, the results suggested that the E2.65 residue mediated the allosteric properties of SB269652 by not only forming a direct hydrogen bond with SB269652, but also by impacting the overall shape and size of the allosteric binding site.
Another recent joint experimental-computational publication also made use of MD simulations to help explain the molecular basis for the binding of bitopic arylamide phenylpiperazine ligands selective for the D3 receptor over D2 (Hayatshahi et al., 2018). Specifically, studies were focused on the prototypic arylamide phenylpiperazine LS-3-134 ligand, which has been found to act as a D3 receptor partial agonist and has also been shown to be 150-fold more selective for the D3 receptor relative to the D2 receptor (Tu et al., 2011). Radioligand binding studies showed that the greatest contribution to the binding energy of the LS-3-134 ligand to the D3 receptor was the phenylpiperazine moiety, but that the arylamide moiety heavily influenced ligand selectivity for the D3 receptor (Hayatshahi et al., 2018). The MD simulations were used to explain the effect that analogs of the piperazine moiety had on the binding affinity. In particular, three different 300 ns unbiased MD simulations were run for LS-3-134 as well as four of its analogs bound to the D3 receptor, revealing that the number of contacts between the protonated nitrogen of piperazine and the D3.32 residue tended to decrease as the size of the piperazine increased, with the exception of only one compound. Calculation of their respective binding energies using MM-PBSA showed that the strength of the electrostatic interaction with the D3.32 residue generally decreased as the size of the piperazine substituent increased. Umbrella sampling simulations were used to generate a PMF aimed at mimicking the unbinding of the ligand protonated nitrogen from the D3.32 residue, and the depth of the bound state in the PMF also agreed with the experimental trend except for one of the analogs (Hayatshahi et al., 2018).
Metabotropic Glutamate Receptors
Metabotropic glutamate receptors (mGluRs) belong to the glutamate (Class C) subfamily of GPCRs given that their endogenous ligand is the neurotransmitter glutamate. Owing to the demonstrated important role of glutamate in pain sensation and transmission, these receptors have been suggested to be promising potential targets for novel pain relieving medications (Pereira and Goudet, 2018). There are eight different subtypes of mGluRs, which are divided into group I (mGluR1 and 5), group II (mGluR2 and 3), and group III (mGluR4, 6, 7, and 8) receptors. While group I mGluRs signal via Gαq, groups II and III signal via Gαi. Several articles suggest opposing effects of the group I vs. group II and III receptors in reference to antinociception, with group I antagonists or group II/III agonists in the spotlight from a drug discovery perspective. In particular, mGluR2/3 agonists or PAMs have been listed among NIDA’s “most wanted” medications in response to the opioid epidemics (Rasmussen et al., 2019).
Like other class C GPCRs, mGluRs are structurally different from rhodopsin-like (class A) GPCRs in that they have a large extracellular domain, also known as Venus Flytrap Domain (VFD), in addition to the 7 transmembrane helical domain (7TMD). Unlike class A GPCRs, the endogenous ligand binding site of mGluRs is located in the extracellular VFD, whereas the transmembrane helical bundle is the primary site for allosteric modulators. Another significant uniqueness is that these receptors are obligate dimers by virtue of a disulfide bond between their VFDs.
Experimental high-resolution structures exist for the 7TMD of mGluR1 (Wu et al., 2014) and the 7TMD and VFT of mGluR5 (Doré et al., 2014; Christopher et al., 2015, 2018; Koehl et al., 2019). There are not, however, published high-resolution structures of the 7TMD of mGluR2 and mGluR3, although their VFT domain has been determined by X-ray crystallography (Muto et al., 2007). Thus, we report here two recent MD-based studies, one for mGluR1 and another for mGluR5, both using their high-resolution experimental structures as starting conformations. To investigate mGluR1 allosteric modulation mechanism at an atomic level of detail appropriate for designing potent NAMs of this receptor, Bai and Yao carried out both biased and unbiased MD simulations on wild-type mGluR1 dimer, as well as its T815M and Y805A mutants, in complex with a NAM known as FITM (4-fluoro-N-(4-(6-(isopropylamino)pyrimidin-4-yl)thiazol-2-yl)-N-methylbenzamide), and embedded in a POPC lipid membrane environment (Bai and Yao, 2016). The simulations used the CHARMM27 force field for the protein and lipid, and the CGenFF force field for FITM. The unbinding PMF, calculated using an adaptive biasing force, revealed an intermediate ligand-binding state along the NAM-mGluR1 dissociation path, stabilized by interactions with residues S735, T748ECL2, C746ECL2, K8117.28 (renumbered K8117.29 according to Isberg et al., 2015) (Bai and Yao, 2016). While hydrogen bonding between FITM and Y805 was identified as a major contributor to stable ligand binding at the crystallographically determined allosteric site, ligand hydrogen bonding to T748 was found to be a crucial factor in stabilizing FITM in an intermediate binding site (Bai and Yao, 2016). Finally, weak interaction analysis of the stabilizing and destabilizing non-covalent interactions using unbiased MD simulations corroborated the importance of van der Waals and hydrogen bond interactions between Y805 and T815 residues in stabilizing the ligand at the binding site.
In another recent unbiased MD investigation, the structural basis for the binding of several NAMs in preclinical or clinical development (mavoglurant, dipraglurant, basimglurant, STX107, and fenobam), as well as three additional NAMs (MPEP, 51D, and 51E) at mGluR5, was elucidated (Fu et al., 2018). The simulations used an explicit POPC membrane in which to embed the ligand-mGluR5 complexes. The force field parameters of the protein atoms were based on the AMBER ff14SB force field, while the force field parameters of the lipid atoms used the Lipid14 force field and ligands were parametrized with the GAFF force field using the Antechamber program. To begin, five NAMs (dipraglurant, basimglurant, STX107, fenobam, and MPEP) were docked to the mGluR5 receptor allosteric site and short 100 ns MD simulations were used to assess the stabilities of the predicted docked poses. Using the final 50 ns of these trajectories, MM/GBSA calculations were carried out to rank the ligands according to their binding free energies. Using the per-residue free energies, eleven residues – I6252.46, I6513.36, S6543.39, P6553.40, L7445.44, W7856.50, F7886.53, M8027.32 (renumbered M8027.33 according to Isberg et al., 2015), V8067.36 (V8067.37 as per Isberg et al., 2015), S8097.39 (S8097.40 as per Isberg et al., 2015), and A8107.40 (A8107.41 as per Isberg et al., 2015) – were identified as the main contributors to the stable binding of all studied NAMs to mGluR5. The apolar nature of most of these eleven residues further suggested that ligands with hydrophobic scaffolds might be better mGluR5 binders.
Concluding Remarks
In this work, we have reviewed recent MD-based investigations of a number of GPCRs that are currently in the spotlight for pain management or to treat or prevent OUD clinical manifestations. With the continued advancements in both computer hardware and MD simulation software, as well as sophisticated tools for analysis of increasingly larger datasets generated by MD simulations, atomic-level insights into the dynamical behavior of GPCRs involved in important pharmacological mechanisms are expected to contribute more and more to the rapid development of therapeutics in response to the opioid crisis.
Author Contributions
JR and MF wrote the manuscript.
Funding
The MF’s lab is funded by National Institutes of Health grants DA045473 and DA038882 for work on opioid receptors. Computations are run on resources available through the Scientific Computing Facility at Mount Sinai and the Extreme Science and Engineering Discovery Environment under MCB080077, which is supported by the National Science Foundation grant number ACI-1053575.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
An, X., Bai, Q., Bing, Z., Zhou, S., Shi, D., Liu, H., et al. (2018). How does agonist and antagonist binding lead to different conformational ensemble equilibria of the κ-Opioid receptor: insight from long-time gaussian accelerated molecular dynamics simulation. ACS Chem. Neurosci. 10, 1575–1584. doi: 10.1021/acschemneuro.8b00535
Bai, Q., Pérez-Sánchez, H., Shi, Z., Li, L., Shi, D., Liu, H., et al. (2018). Computational studies on horseshoe shape pocket of human orexin receptor type 2 and boat conformation of suvorexant by molecular dynamics simulations. Chem. Biol. Drug Des. 92, 1221–1231. doi: 10.1111/cbdd.13181
Bai, Q., and Yao, X. (2016). Investigation of allosteric modulation mechanism of metabotropic glutamate receptor 1 by molecular dynamics simulations, free energy and weak interaction analysis. Sci. Rep. 6:21763. doi: 10.1038/srep21763
Ballesteros, J. A., and Weinstein, H. (1995). “Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors,” in Methods in Neurosciences, eds S. C. Sealfon, and P. M. Conn, (San Diego, CA: Academic Press), 366–428. doi: 10.1016/s1043-9471(05)80049-7
Bartuzi, D., Kaczor, A. A., and Matosiuk, D. (2019). Molecular mechanisms of allosteric probe dependence in μ opioid receptor. J. Biomol. Struct. Dyn. 37, 36–47. doi: 10.1080/07391102.2017.1417914
Bernardi, R. C., Melo, M. C. R., and Schulten, K. (2015). Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta 1850, 872–877. doi: 10.1016/j.bbagen.2014.10.019
Best, R. B., Zhu, X., Shim, J., Lopes, P. E., Mittal, J., Feig, M., et al. (2012). Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles. J. Chem. Theory Comput. 8, 3257–3273. doi: 10.1021/ct300400x
Bohn, L. M., Lefkowitz, R. J., Gainetdinov, R. R., Peppel, K., Caron, M. G., and Lin, F. T. (1999). Enhanced morphine analgesia in mice lacking beta-arrestin 2. Science 286, 2495–2498. doi: 10.1126/science.286.5449.2495
Burford, N. T., Livingston, K. E., Canals, M., Ryan, M. R., Budenholzer, L. M., Han, Y., et al. (2015). Discovery, synthesis, and molecular pharmacology of selective positive allosteric modulators of the delta-opioid receptor. J. Med. Chem. 58, 4220–4229. doi: 10.1021/acs.jmedchem.5b00007
Carlezon, W. A. Jr., and Krystal, A. D. (2016). Kappa-Opioid antagonists for psychiatric disorders: from bench to clinical trials. Depress. Anxiety 33, 895–906. doi: 10.1002/da.22500
Centers for Disease Control and Prevention [CDC], (2017). Wide-Ranging Online Data for Epidemiologic Research (WONDER). Atlanta, GA: CDC.
Chang, D. S., Raghavan, R., Christiansen, S., and Cohen, S. P. (2015). Emerging targets in treating pain. Curr. Opin. Anaesthesiol. 28, 379–397. doi: 10.1097/ACO.0000000000000216
Cheng, J., Sun, X., Li, W., Liu, G., Tu, Y., and Tang, Y. (2016). Molecular switches of the κ opioid receptor triggered by 6′-GNTI and 5′-GNTI. Sci. Rep. 6:18913. doi: 10.1038/srep18913
Cheng, J.-X., Cheng, T., Li, W.-H., Liu, G.-X., Zhu, W.-L., and Tang, Y. (2018). Computational insights into the G-protein-biased activation and inactivation mechanisms of the μ opioid receptor. Acta Pharmacol. Sin. 39:154. doi: 10.1038/aps.2017.158
Chien, E. Y., Liu, W., Zhao, Q., Katritch, V., Han, G. W., Hanson, M. A., et al. (2010). Structure of the human dopamine D3 receptor in complex with a D2/D3 selective antagonist. Science 330, 1091–1095. doi: 10.1126/science.1197410
Christopher, J. A., Aves, S. J., Bennett, K. A., Doreì, A. S., Errey, J. C., Jazayeri, A., et al. (2015). Fragment and structure-based drug discovery for a class C GPCR: discovery of the mGlu5 negative allosteric modulator HTL14242 (3-chloro-5-[6-(5-fluoropyridin-2-yl) pyrimidin-4-yl] benzonitrile). J. Med. Chem. 58, 6653–6664. doi: 10.1021/acs.jmedchem.5b00892
Christopher, J. A., Orgovaìn, Z. N., Congreve, M., Doreì, A. S., Errey, J. C., Marshall, F. H., et al. (2018). Structure-based optimization strategies for G protein-coupled receptor (GPCR) allosteric modulators: a case study from analyses of new metabotropic glutamate receptor 5 (mGlu5) X-ray structures. J. Med. Chem. 62, 207–222. doi: 10.1021/acs.jmedchem.7b01722
Croci, T., and Zarini, E. (2007). Effect of the cannabinoid CB1 receptor antagonist rimonabant on nociceptive responses and adjuvant-induced arthritis in obese and lean rats. Br. J. Pharmacol. 150, 559–566. doi: 10.1038/sj.bjp.0707138
Della Longa, S., and Arcovito, A. (2019). Microswitches for the activation of the nociceptin receptor induced by cebranopadol: hints from microsecond molecular dynamics. J. Chem. Inform Model. 59, 818–831. doi: 10.1021/acs.jcim.8b00759
Doré, A. S., Okrasa, K., Patel, J. C., Serrano-Vega, M., Bennett, K., Cooke, R. M., et al. (2014). Structure of class C GPCR metabotropic glutamate receptor 5 transmembrane domain. Nature 511, 557–562. doi: 10.1038/nature13396
Ferruz, N., Doerr, S., Vanase-Frawley, M. A., Zou, Y., Chen, X., Marr, E. S., et al. (2018). Dopamine D3 receptor antagonist reveals a cryptic pocket in aminergic GPCRs. Sci. Rep. 8:897. doi: 10.1038/s41598-018-19345-7
Fu, T., Zheng, G., Tu, G., Yang, F., Chen, Y., Yao, X., et al. (2018). Exploring the binding mechanism of metabotropic glutamate receptor 5 negative allosteric modulators in clinical trials by molecular dynamics simulations. ACS Chem. Neurosci. 9, 1492–1502. doi: 10.1021/acschemneuro.8b00059
Guo, W., Urizar, E., Kralikova, M., Mobarec, J. C., Shi, L., Filizola, M., et al. (2008). Dopamine D2 receptors form higher order oligomers at physiological expression levels. EMBO J. 27, 2293–2304. doi: 10.1038/emboj.2008.153
Harvey, M. J., Giupponi, G., and Fabritiis, G. D. (2009). ACEMD: accelerating biomolecular dynamics in the microsecond time scale. J. Chem. Theory Comput. 5, 1632–1639. doi: 10.1021/ct9000685
Hasebe, K., Kawai, K., Suzuki, T., Kawamura, K., Tanaka, T., Narita, M., et al. (2004). Possible pharmacotherapy of the opioid kappa receptor agonist for drug dependence. Ann. N. Y. Acad. Sci. 1025, 404–413. doi: 10.1196/annals.1316.050
Hayatshahi, H. S., Xu, K., Griffin, S. A., Taylor, M., Mach, R. H., Liu, J., et al. (2018). Analogues of arylamide phenylpiperazine ligands to investigate the factors influencing D3 dopamine receptor bitropic binding and receptor subtype selectivity. ACS Chem Neurosci. 9, 2972–2983. doi: 10.1021/acschemneuro.8b00142
Hua, T., Vemuri, K., Nikas, S. P., Laprairie, R. B., Wu, Y., Qu, L., et al. (2017). Crystal structures of agonist-bound human cannabinoid receptor CB 1. Nature 547, 468–471. doi: 10.1038/nature23272
Hua, T., Vemuri, K., Pu, M., Qu, L., Han, G. W., Wu, Y., et al. (2016). Crystal structure of the human cannabinoid receptor CB1. Cell 167, 750–762.e14. doi: 10.1016/j.cell.2016.10.004
Huang, L., and Roux, B. (2013). Automated force field parameterization for nonpolarizable and polarizable atomic models based on ab initio target data. J. Chem. Theory Comput. 9, 3543–3556. doi: 10.1021/ct4003477
Husic, B. E., and Pande, V. S. (2018). Markov state models: from an art to a science. J. Am. Chem. Soc. 140, 2386–2396. doi: 10.1021/jacs.7b12191
Isberg, V., de Graaf, C., Bortolato, A., Cherezov, V., Katritch, V., Marshall, F. H., et al. (2015). Generic GPCR residue numbers - aligning topology maps while minding the gaps. Trends Pharmacol. Sci. 36, 22–31. doi: 10.1016/j.tips.2014.11.001
Janecka, A., Piekielna-Ciesielska, J., and Wtorek, K. (2019). Biased agonism as an emerging strategy in the search for better opioid analgesics. Curr. Med. Chem. doi: 10.2174/0929867326666190506103124 [Epub ahead of print].
Jung, S. W., Cho, A. E., and Yu, W. (2018). Exploring the ligand efficacy of cannabinoid receptor 1 (CB1) using molecular dynamics simulations. Sci. Rep. 8:13787. doi: 10.1038/s41598-018-31749-z
Kapoor, A., Martinez-Rosell, G., Provasi, D., De Fabritiis, G., and Filizola, M. (2017). Dynamic and kinetic elements of μ-Opioid receptor functional selectivity. Sci. Rep. 7:11255. doi: 10.1038/s41598-017-11483-8
Karhu, L. V., Magarkar, A., Bunker, A., and Xhaard, H. G. (2019). Determinants of orexin receptor binding and activation–a molecular dynamics study. J. Phys. Chem. B 123, 2609–2622. doi: 10.1021/acs.jpcb.8b10220
Koehl, A., Hu, H., Feng, D., Sun, B., Zhang, Y., Robertson, M. J., et al. (2019). Structural insights into the activation of metabotropic glutamate receptors. Nature 566, 79–84. doi: 10.1038/s41586-019-0881-4
Krishna Kumar, K., Shalev-Benami, M., Robertson, M. J., Hu, H., Banister, S. D., Hollingsworth, S. A., et al. (2019). Structure of a signaling cannabinoid receptor 1-G protein complex. Cell 176, 448–458.e12. doi: 10.1016/j.cell.2018.11.040
Kumar, V., Bonifazi, A., Ellenberger, M. P., Keck, T. M., Pommier, E., Rais, R., et al. (2016). Highly selective dopamine D3 receptor (D3R) antagonists and partial agonists based on eticlopride and the D3R crystal structure: new leads for opioid dependence treatment. J. Med. Chem. 59, 7634–7650. doi: 10.1021/acs.jmedchem.6b00860
Land, B. B., Bruchas, M. R., Lemos, J. C., Xu, M., Melief, E. J., and Chavkin, C. (2008). The dysphoric component of stress is encoded by activation of the dynorphin kappa-opioid system. J. Neurosci. 28, 407–414. doi: 10.1523/JNEUROSCI.4458-07.2008
Lane, J. R., Donthamsetti, P., Shonberg, J., Draper-Joyce, C. J., Dentry, S., Michino, M., et al. (2014). A new mechanism of allostery in a G protein–coupled receptor dimer. Nat. Chem. Biol. 10, 745–752. doi: 10.1038/nchembio.1593
Livingston, K. E., Stanczyk, M. A., Burford, N. T., Alt, A., Canals, M., and Traynor, J. R. (2018). Pharmacologic evidence for a putative conserved allosteric site on opioid receptors. Mol. Pharmacol. 93, 157–167. doi: 10.1124/mol.117.109561
Lopes, P. E., Guvench, O., and MacKerell, A. D. Jr. (2015). Current status of protein force fields for molecular dynamics simulations. Methods Mol. Biol. 1215, 47–71. doi: 10.1007/978-1-4939-1465-4_3
Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255
Mistry, S. N., Shonberg, J., Draper-Joyce, C. J., Klein Herenbrink, C., Michino, M., Shi, L., et al. (2015). Discovery of a novel class of negative allosteric modulator of the dopamine D2 receptor through fragmentation of a bitopic ligand. J. Med. Chem. 58, 6819–6843. doi: 10.1021/acs.jmedchem.5b00585
Muto, T., Tsuchiya, D., Morikawa, K., and Jingami, H. (2007). Structures of the extracellular regions of the group II/III metabotropic glutamate receptors. Proc. Natl. Acad. Sci. U.S.A. 104, 3759–3764. doi: 10.1073/pnas.0611577104
Nerenberg, P. S., and Head-Gordon, T. (2018). New developments in force fields for biomolecular simulations. Curr. Opin. Struct. Biol. 49, 129–138. doi: 10.1016/j.sbi.2018.02.002
Oren, M. B., Alexander, D. M. Jr., Benoit, R., and Masakatsu, W. (2001). Computational Biochemistry and Biophysics. Boca Raton, FL: CRC Press.
Pereira, V., and Goudet, C. (2018). Emerging trends in pain modulation by metabotropic glutamate receptors. Front. Mol. Neurosci. 11:464. doi: 10.3389/fnmol.2018.00464
Ponder, J. W., and Case, D. A. (2003). Force fields for protein simulations. Adv. Protein Chem. 66, 27–85. doi: 10.1016/s0065-3233(03)66002-x
Raehal, K. M., Walker, J. K., and Bohn, L. M. (2005). Morphine side effects in beta-arrestin 2 knockout mice. J. Pharmacol. Exp. Ther. 314, 1195–1201. doi: 10.1124/jpet.105.087254
Rasmussen, K., White, D. A., and Acri, J. B. (2019). NIDA’s medication development priorities in response to the Opioid Crisis: ten most wanted. Neuropsychopharmacology 44, 657–659. doi: 10.1038/s41386-018-0292-5
Remesic, M., Hruby, V. J., Porreca, F., and Lee, Y. S. (2017). Recent advances in the realm of allosteric modulators for opioid receptors for future therapeutics. ACS Chem. Neurosci. 8, 1147–1158. doi: 10.1021/acschemneuro.7b00090
Ribeiro, J. M. L., and Filizola, M. (2019). Allostery in G protein-coupled receptors investigated by molecular dynamics simulations. Curr. Opin. Struct. Biol. 55, 121–128. doi: 10.1016/j.sbi.2019.03.016
Ribeiro, J. M. L., Tsai, S. T., Pramanik, D., Wang, Y., and Tiwary, P. (2019). Kinetics of Ligand-Protein dissociation from all-atom simulations: are we there yet? Biochemistry 58, 156–165. doi: 10.1021/acs.biochem.8b00977
Ringkamp, M., Dougherty, P. M., and Raja, S. N. (2018). “Anatomy and physiology of the pain signaling process,” in Essentials of Pain Medicine, eds H. T. Benzon, S. N. Raja, S. M. Fishman, and S. S. Liu, (Amsterdam: Elsevier), 3.–10.
Saleh, N., Hucke, O., Kramer, G., Schmidt, E., Montel, F., Lipinski, R., et al. (2018). Multiple binding sites contribute to the mechanism of mixed agonistic and positive allosteric modulators of the cannabinoid CB1 receptor. Angewandte Chem. 130, 2610–2615. doi: 10.1002/anie.201708764
Saleh, N., Ibrahim, P., Saladino, G., Gervasio, F. L., and Clark, T. (2017). An efficient metadynamics-based protocol to model the binding affinity and the transition state ensemble of G-protein-coupled receptor ligands. J. Chem. Inform. Model. 57, 1210–1217. doi: 10.1021/acs.jcim.6b00772
Schneider, S., Provasi, D., and Filizola, M. (2016). How oliceridine (TRV-130) binds and stabilizes a μ-opioid receptor conformational state that selectively triggers G protein signaling pathways. Biochemistry 55, 6456–6466. doi: 10.1021/acs.biochem.6b00948
Shang, Y., Yeatman, H. R., Provasi, D., Alt, A., Christopoulos, A., Canals, M., et al. (2016). Proposed mode of binding and action of positive allosteric modulators at opioid receptors. ACS Chem. Biol. 11, 1220–1229. doi: 10.1021/acschembio.5b00712
Shao, Z., Yin, J., Chapman, K., Grzemska, M., Clark, L., Wang, J., et al. (2016). High-resolution crystal structure of the human CB1 cannabinoid receptor. Nature 540, 602–606. doi: 10.1038/nature20613
Shonberg, J., Draper-Joyce, C., Mistry, S. N., Christopoulos, A., Scammells, P. J., Lane, J. R., et al. (2015). Structure-activity study of N-((trans)-4-(2-(7-cyano-3,4-dihydroisoquinolin-2(1H)-yl)ethyl)cyclohexyl)-1H-ind ole-2-carboxamide (SB269652), a bitopic ligand that acts as a negative allosteric modulator of the dopamine D2 receptor. J. Med. Chem. 58, 5287–5307. doi: 10.1021/acs.jmedchem.5b00581
Silvano, E., Millan, M. J., la Cour, C. M., Han, Y., Duan, L., Griffin, S. A., et al. (2010). The tetrahydroisoquinoline derivative SB269, 652 is an allosteric antagonist at dopamine D3 and D2 receptors. Mol. Pharmacol. 78, 925–934. doi: 10.1124/mol.110.065755
Slivicki, R. A., Xu, Z., Kulkarni, P. M., Pertwee, R. G., Mackie, K., Thakur, G. A., et al. (2018). Positive allosteric modulation of cannabinoid receptor Type 1 suppresses pathological pain without producing tolerance or dependence. Biol. Psychiatry 84, 722–733. doi: 10.1016/j.biopsych.2017.06.032
Suno, R., Kimura, K. T., Nakane, T., Yamashita, K., Wang, J., Fujiwara, T., et al. (2018). Crystal structures of human orexin 2 receptor bound to the subtype-selective antagonist EMPA. Structure 26, 7–19.e5. doi: 10.1016/j.str.2017.11.005
Tu, Z., Li, S., Cui, J., Xu, J., Taylor, M., Ho, D., et al. (2011). Synthesis and pharmacological evaluation of fluorine-containing D3 dopamine receptor ligands. J. Med. Chem. 54, 1555–1564. doi: 10.1021/jm101323b
Ueda, M., Iwasaki, H., Wang, S., Murata, E., Poon, K. Y., Mao, J., et al. (2014). Cannabinoid receptor type 1 antagonist, AM251, attenuates mechanical allodynia and thermal hyperalgesia after burn injury. Anesthesiology 121, 1311–1319. doi: 10.1097/ALN.0000000000000422
Valsson, O., Tiwary, P., and Parrinello, M. (2016). Enhancing important fluctuations: rare events and metadynamics from a conceptual viewpoint. Annu. Rev. Phys. Chem. 67, 159–184. doi: 10.1146/annurev-physchem-040215-112229
Vanommeslaeghe, K., Hatcher, E., Acharya, C., Kundu, S., Zhong, S., Shim, J., et al. (2010). CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 31, 671–690. doi: 10.1002/jcc.21367
Verma, R. K., Abramyan, A. M., Michino, M., Free, R. B., Sibley, D. R., Javitch, J. A., et al. (2018). The E2. 65A mutation disrupts dynamic binding poses of SB269652 at the dopamine D2 and D3 receptors. PLoS Comput. Biol. 14:e1005948. doi: 10.1371/journal.pcbi.1005948
Viscusi, E. R., Skobieranda, F., Soergel, D. G., Cook, E., Burt, D. A., and Singla, N. (2019). APOLLO-1: a randomized placebo and active-controlled phase III study investigating oliceridine (TRV130), a G protein-biased ligand at the micro-opioid receptor, for management of moderate-to-severe acute pain following bunionectomy. J. Pain Res. 12, 927–943. doi: 10.2147/JPR.S171013
Volkow, N. D., Jones, E. B., Einstein, E. B., and Wargo, E. M. (2019). Prevention and treatment of opioid misuse and addiction: a review. JAMA Psychiatry 76, 208–216. doi: 10.1001/jamapsychiatry.2018.3126
Wacker, D., and Roth, B. L. (2016). An alerting structure: human orexin receptor 1. Nat. Struct. Mol. Biol. 23, 265–266. doi: 10.1038/nsmb.3198
Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035
Wootten, D., Christopoulos, A., and Sexton, P. M. (2013). Emerging paradigms in GPCR allostery: implications for drug discovery. Nat. Rev. Drug Discov. 12, 630–644. doi: 10.1038/nrd4052
Wu, H., Wang, C., Gregory, K. J., Han, G. W., Cho, H. P., Xia, Y., et al. (2014). Structure of a class C GPCR metabotropic glutamate receptor 1 bound to an allosteric modulator. Science 344, 58–64. doi: 10.1126/science.1249489
Yin, J., Babaoglu, K., Brautigam, C. A., Clark, L., Shao, Z., Scheuermann, T. H., et al. (2016). Structure and ligand-binding mechanism of the human OX 1 and OX 2 orexin receptors. Nat. Struct. Mol. Biol. 23, 293–299. doi: 10.1038/nsmb.3183
Yin, J., Mobarec, J. C., Kolb, P., and Rosenbaum, D. M. (2015). Crystal structure of the human OX 2 orexin receptor bound to the insomnia drug suvorexant. Nature 519, 247–250. doi: 10.1038/nature14035
Keywords: GPCRs, opioid crisis, molecular dynamics, pain, opioid use disorder
Citation: Ribeiro JML and Filizola M (2019) Insights From Molecular Dynamics Simulations of a Number of G-Protein Coupled Receptor Targets for the Treatment of Pain and Opioid Use Disorders. Front. Mol. Neurosci. 12:207. doi: 10.3389/fnmol.2019.00207
Received: 31 May 2019; Accepted: 07 August 2019;
Published: 23 August 2019.
Edited by:
Meritxell Canals, University of Nottingham, United KingdomReviewed by:
Leonardo Pardo, Autonomous University of Barcelona, SpainHugo Gutiérrez De Teran, Uppsala University, Sweden
Copyright © 2019 Ribeiro and Filizola. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Marta Filizola, marta.filizola@mssm.edu