pubs.acs.org/jcim Article Structural Assessment of Agonist Efficacy in the μ‑Opioid Receptor: Morphine and Fentanyl Elicit Different Activation Patterns Adrián Ricarte, James A. R. Dalton,* and Jesús Giraldo* Cite This: J. Chem. Inf. Model. 2021, 61, 1251−1274 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. ABSTRACT: Over the past two decades, the opioid epidemic in the United Downloaded via UNIV AUTONOMA DE BARCELONA on March 22, 2021 at 08:43:07 (UTC). States and Canada has evidenced the need for a better understanding of the molecular mechanisms of medications used to fight pain. Morphine and fentanyl are widely used in opiate-mediated analgesia for the treatment of chronic pain. These compounds target the μ-opioid receptor (MOR), a class A G protein-coupled receptor (GPCR). In light of described higher efficacy of fentanyl with respect to morphine, we have performed independent μs- length unbiased molecular dynamics (MD) simulations of MOR complexes with each of these ligands, including the MOR antagonist naltrexone as a negative control. Consequently, MD simulations totaling 58 μs have been conducted to elucidate at the atomic level ligand-specific receptor activity and signal transmission in the MOR. In particular, we have identified stable binding poses of morphine and fentanyl, which interact differently with the MOR. Different ligand−receptor interaction landscapes directly induce sidechain conformational changes of orthosteric pocket residues: Asp1493.32, Tyr1503.33, Gln1262.60, and Lys2355.39. The induced conformations determine Asp1493.32−Tyr3287.43 sidechain−sidechain interactions and Trp2956.48−Ala2425.46 sidechain−backbone H-bond formations, as well as Met1533.36 conformational changes. In addition to differences in ligand binding, different intracellular receptor conformational changes are observed as morphine preferentially activates transmembrane (TM) helices: TM3 and TM5, while fentanyl preferentially activates TM6 and TM7. As conformational changes in TM6 and TM7 are widely described as being the most crucial aspect in GPCR activation, this may contribute to the greater efficacy of fentanyl over morphine. These computationally observed functional differences between fentanyl and morphine may provide new avenues for the design of safer but not weaker opioid drugs because it is desirable to increase the safety of medicines without sacrificing their efficacy. ■ INTRODUCTION Chronic pain is a major public health problem with a high MOR Gi protein-biased agonist TRV130 (oliceridine).9,10 Oliceridine has been recently approved (August 2020) by the prevalence and impact on quality of life.1,2 Opiate-mediated US Food and Drug Administration (FDA) under the trade analgesia has shown to be the most efficacious treatment for name Olinvyk for short-term intravenous use in hospitals and chronic pain,3 which primarily targets the μ-opioid receptor other controlled settings.11 Moreover, in the prescribing (MOR), a class A G protein-coupled receptor (GPCR) information leaflet, it is stated that oliceridine is a full opioid encoded by the Oprm1 gene.4 Upon binding of opioid agonist and is relatively selective for the MOR. It is also analgesics, ligand-specific conformational changes occur in the recognized that, depending on the dose used, adverse MOR,5 which lead to activation of the receptor and coupling reactions, including respiratory, and CNS depression may with intracellular heterotrimeric Gi protein or β-arrestin,5 thus appear. Importantly, it is also stated that the precise mediating the propagation of the receptor message through mechanism of the analgesic action is unknown although different signaling pathways. It remains elusive precisely how specific CNS opioid receptors are thought to play a role in the ligand−receptor interactions at the binding site influence analgesic effects of the drug.12,13 Thus, although drug discovery MOR conformational dynamics and hence different signal programs based on biased agonism seem to be a promising transmissions. While the therapeutic effect of opioid analgesics strategy, in particular for the MOR, they still need further is mainly attributed to MOR activation through Gi protein signaling,5−8 their side effects, opioid-induced hyperalgesia, constipation, respiratory depression, and analgesic tolerance, Received: August 5, 2020 have mostly been linked to β-arrestin signaling.5−8 Thus, Published: January 15, 2021 agonists biased toward the Gi protein pathway would be ideal to be considered for drug development as potential future medications. Following this research direction, high expect- ations were generated by preclinical in vivo studies on the © 2021 American Chemical Society https://dx.doi.org/10.1021/acs.jcim.0c00890 1251 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article knowledge and experimentation both at the molecular and mMOR crystal structure presents, compared to the inactive quantitative pharmacology levels.14−18 In this regard, following crystal structure, (i) large outward movement of TM6 relative the strategy of MOR biased agonism, several studies have to TM3, (ii) smaller inward movements of TM5 and TM7, examined how agonists targeting the MOR interact with (iii) breakage of the H-bond between Arg1653.50 and specific residues to selectively regulate coupling of Gi protein Thr2796.34 (equivalent of the classical Arg3.50−Asp/Glu6.30 or β-arrestin,19−23 which in some cases have led to the GPCR ionic lock, which involves an acidic amino acid at the discovery of new analgesics with reduced negative side position 6.30 that is lacking in the MOR)43 and formation of a effects.19,24 Yet, although some structural studies of varying H-bond between Arg1653.50 and Tyr2525.58, (iv) inward MOR agonist efficacy have been addressed,25,26 there is still a movement of the N7.49PxxY7.53 motif on TM7 toward lack of complete understanding of the precise nature of TM5,5,7,19,30 and (v) upward axial movement of TM3.44 In ligand−MOR interactions responsible for observed pharmaco- general, previous MD studies were performed by removing logical profiles. In this regard, further studies at the atomic non-native Nb39 and T4 lysozyme, as well as bound ligands level including agonists with different efficacies can shed light BU-72 and β-FNA, from activated and inactive states of the on MOR activation mechanisms. This was the reason for mMOR, respectively.5,29,30 In particular, a model was provided selecting morphine and fentanyl for the present study. Fentanyl to show how efficacy depends on small chemical differences in is a high-efficacy MOR agonist, whereas morphine is a lower- structurally similar morphine-like ligands.5 A common feature efficacy MOR agonist, which, depending on the assay of bound opioid ligands revealed in these studies5,7,19,29,30,45 is performed, can behave as a partial agonist.27,28 It is expected the protonation of their tertiary amine, which is necessary to that the generated knowledge can be later incorporated into form a conserved salt bridge with Asp1473.32 on TM3 in the drug discovery programs addressed to avoid ligands with orthosteric pocket. This residue is also located within the H- unwanted effects yet conserving the analgesic capacity of the bond distance with Tyr3267.43 on TM7 in the mMOR,46 thus most powerful drugs. favoring interaction between TM3 and TM7 and activation of Biophysical studies of the MOR have led to the the receptor.44,47−49 This can be observed in activated and determination of three different crystal structures of the inactive crystal structures of the mMOR,5,29 where the distance receptor: inactive,29 activated,5 and fully active30 structures of between these residues is closer in the activated crystal the murine MOR (mMOR). In 2012, the inactive crystal structure than in the inactive one. In addition, computational structure of the mMOR29 was obtained by X-ray crystallog- studies of the mMOR have reported specific protein−ligand raphy, co-crystallized with the morphine-like antagonist β-FNA interactions such as the plausible existence of two possible covalently bound to Lys2335.39 (superscript numbers refer to tautomers of His2976.52, which can have relevance in ligand the Ballesteros and Weinstein generic numbering scheme31). function,19 as well as Tyr3267.43 and Gln1242.60 H-bonding This structure provided the first high-resolution insight of the interactions, which may also determine ligand activity.6,46 MOR,29 which enabled the application of structure-based However, computational techniques are yet to satisfactorily computational approaches: rational targeted drug design,32−34 explain varying MOR agonist efficacy and potency. receptor oligomerization studies with other GPCRs,35,36 and Opioid analgesics, the prototypical pain killers, differ from study of the mechanism by which sodium ions prevent large- each other in their structural scaffold and potency. Of all of scale movement of transmembrane (TM) helix 6 away from them, the naturally occurring compound morphine has been TM3, which inhibits receptor activation.37,38 In 2015, an considered the reference MOR agonist to which other opioid activated crystal structure of the mMOR bound to the analgesics are compared.50 In the last decade, fentanyl, a morphine-like agonist BU-72,5 with intracellular binding of a synthetic opioid 50−100 times more potent than mor- G protein mimetic camelid antibody fragment (Nb39), was phine,51−53 and fentanyl analogs have been related to a rapid resolved by X-ray crystallography.5 G protein mimetic increase in the number of opioid overdose deaths because of nanobodies have been shown to be of great utility for legal (such as transdermal patches, which led to an increase in stabilization of active-like states of GPCRs.39−41 The crystal deaths up to 201054,55) or illicit manufacture.51−53,56,57 This structure of the activated mMOR shed light on the structural has contributed to what is called the opioid epidemic, which features of MOR activation that are responsible for the efficacy has affected mainly the United States and Canada and, in a of most therapeutic analgesics. In addition, it allowed the lower degree, Europe.58,59 Fentanyl differs from morphine description of an extensive reorganization of the protein−water primarily in its pharmacokinetic properties,50 constituting a polar network, required in the full activation process,5 which is highly lipid-soluble drug compared to low lipid solubility of associated with an efficient allosteric coupling between the morphine, allowing it to penetrate the blood−brain barrier receptor orthosteric pocket and Gi protein-coupling inter- faster.50,53 Therefore, fentanyl has a faster onset but shorter face.5,42 However, it was not until 2018 that a structure of the duration of action than morphine.50,53 Following receptor fully active state of the mMOR was achieved by cryo-electron activation, differences in morphine and fentanyl structural microscopy.30 This receptor structure with bound peptide scaffolds determine specific ligand-dependent receptor inter- agonist DAMGO and human Gi protein has provided new nalization processes.60 It has been described60 that high- insights into the mechanism of intracellular binding to the efficacy opiates like fentanyl induce rapid MOR phosphor- MOR.30 ylation and internalization compared to morphine. However, Computational techniques such as molecular dynamics 50−100 times potency difference51−53 between fentanyl and (MD) simulations have shown promise in further under- morphine has been widely discussed by in vitro studies23,61−64 standing the complexities of MOR signaling mechanisms. For and its potential relation with fentanyl β-arrestin-biased example, several MD simulation studies of the behavior has not been clearly established. In this regard, a MOR,5,7,19,25,26,30 initialized from inactive29 and activated5 recent literature65 study supports that both fentanyl and mMOR crystal structures, have been worthy in assessing some morphine work as unbiased agonists66−68 toward Gi protein relevant features of the activation process. Briefly, the activated and β-arrestin signaling pathways.61−63 Apparently, contra- 1252 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 1. Crystal and docked ligands. Morphine and fentanyl binding poses. (a) 2D binding poses of respective inactive, activated, and fully active co-crystallized ligands: β-FNA, BU-72, and DAMGO, and docked ligands of interest: morphine, fentanyl, and naltrexone. (b) 3D conformation of fentanyl and morphine representative stable binding poses observed during MD simulations (purple and pink, left and right, respectively). Residues of the hMOR in close contact with fentanyl or morphine (<3.5 Å) are colored in dark green or dark red, respectively. Residues colored in light green or light red, respectively, are not in close contact but displayed for comparison purposes. dictory conclusions on biased signaling may depend not only mathematical models and parameters employed. Thus, while on the variety of biological assays used but also on the the Emax/EC50 sigmoidal-fitting parameter was chosen in 1253 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling ■ pubs.acs.org/jcim Article Vasudevan et al.,65 a comparison between operational τ/KA RESULTS and τ parameters was made in Burgueño et al.64 In the latter Homology Models of the hMOR Remain Conforma- study, while fentanyl was found to be G protein-biased under tionally Stable in Control MD Simulations. Native the τ/KA parameter, an opposite bias toward the β-arrestin activated and inactive hMOR homology models were signaling pathway was found when using the τ parameter.64 generated from respective activated and inactive mMOR Moreover, it has been shown, by measuring the antinociceptive crystals5,29 by renumbering sequences, mutating non-con- and respiratory depressant effects of some MOR agonists, that served residues, and adding non-crystallized segments of the the low intrinsic efficacy of some opioid ligands can explain inactive hMOR N-terminus, inactive hMOR intracellular loop, their improved side effect profile.69 Additionally, a strong and activated hMOR C-terminus (see Methods). In order to correlation between measures of efficacy for receptor validate the accuracy of activated and inactive hMOR models, activation, G protein coupling, and β-arrestin recruitment we performed two control replicas of 2 μs-length unbiased MD was found for some MOR agonists, including those that had simulations with co-crystallized agonist BU-725 and three previously been described as biased.69 A detailed review on the relationship between agonist efficacy, biased agonism, and control replicas of 3 μs-length with docked antagonist therapeutic window can be bound in Gillis et al.18 naltrexone, respectively. As co-crystallized antagonist β-FNA From the above data, it can be concluded that the of the inactive mMOR29 is covalently bound to Lys2335.39, it is mechanisms of MOR activation and signal transduction are not an ideal control antagonist for our purposes here, but not yet fully understood. Although downstream signaling because of high similarity between the molecular structures of events occur after a number of steps following the initial drug− naltrexone and β-FNA (Figure 1), β-FNA can be easily receptor interactions at the receptor binding site, a chemical substituted for naltrexone. Therefore, in order to test docking connection is expected to be present between initial and final accuracy, BU-72 was re-docked into our hMOR active model, signaling states. Thus, understanding of the distinct receptor while naltrexone was docked into our inactive hMOR model conformational changes induced by high- and low-efficacy and compared to co-crystallized β-FNA. Docking results reveal MOR agonists, and their differential involvement in the a top-ranked BU-72 binding pose with a root mean square activation of the receptor, can be informative to identify the deviation (RMSD) of 0.4 Å from its co-crystallized position structural features that cause different MOR functionalities. and a top-ranked naltrexone binding pose with an RMSD of This knowledge is crucial for the rational design of safer but 0.7 Å compared to the common heavy atoms of co-crystallized not weaker drugs. β-FNA. Trajectories of the hMOR with bound BU-72 and The experimentally observed functional MOR response to naltrexone show stable receptor conformations (average morphine and fentanyl is a consequence of molecular events RMSDs of 1.8 Å (0.4 SD) and 2.8 Å (0.2 SD), respectively that occur first in the receptor as a result of the interactions (Figures S1 and S4)) and ligand binding poses (average established by these ligands with residues of the receptor RMSDs of 1.7 Å (0.6 SD) and 2.2 Å (0.6 SD), respectively binding pocket, which are then transmitted from the (Figures S1 and S2)). These results indicate that both extracellular to the intracellular side of the cell membrane. In activated and inactive hMOR models remain near their initial this context, in the present study, we performed unbiased μs- conformation throughout long-timescale MD trajectories, length MD simulations to complete a detailed atomic-level which supports their use as reliable starting points to evaluate study of the receptor conformational space induced or selected the binding effect of ligands of interest: morphine and fentanyl. in the MOR by morphine and fentanyl in order to explain the Binding Poses of Morphine and Fentanyl Differ. In variations seen in opiate efficacy.27 Receptor activation by order to computationally examine at the atomic level the MOR either induction or selection of receptor active states has been efficacy differences between fentanyl and morphine, we first recurrently under debate70 though they may coexist and be docked each ligand into our models of activated and inactive part of the same thermodynamic cycle.71 Here, both models hMOR models. Starting from these protein−ligand systems, are considered, which will allow noteworthy structural we performed three replicas of 3 μs-length unbiased MD comparisons. In order to understand how opioid analgesics simulations to allow for ligand-dependent receptor conforma- transmit their molecular signal and generate effects in humans, tional change. From the respective morphine- and fentanyl- we consider it appropriate to perform our study in the native bound MD trajectories, we identified one stable binding pose human MOR (hMOR) and across multiple independent MD for each ligand (Figure 1), consistent across all replicas and trajectories. Our results show that fentanyl displays a predominantly selected over other alternative temporary poses completely different binding pose from morphine, which (Figure S2). According to RMSD across MD simulations leads to divergent effects on the orthosteric pocket conforma- (Figure S2), and considering a threshold of ± 2.5 Å to define tional arrangement. The fact that both morphine and fentanyl ligand stability, we observe that both morphine and fentanyl could activate both Gi protein and β-arrestin signaling are stable after the first microsecond of simulation time and pathways raises the question of which signaling pathway until the end of each trajectory, independently of the initial corresponds to the receptor conformational changes found receptor state (average ligand RMSDs of 1.3 Å (0.5 SD) and computationally. This is a question for which there is not a 2.2 Å (0.5 SD), respectively). Notable initial ligand instability, definite answer with the present knowledge. Nevertheless, to limited to the first microsecond of replicas starting from the avoid speculations, structural comparison of computational receptor activated state with bound fentanyl and replicas #1 structures was made with the current MOR crystal and #2 starting from the receptor inactive state with bound structures,5,29,30 which are related to the Gi signaling pathway. morphine, can be explained by receptor−ligand relaxation and Our results show that extracellular conformational changes of adaptation (Figure S2). Therefore, both ligands are bound in a the receptor compromise the ability of its intracellular side to stable fashion, and any differences in ligand stability in the consistently reach (induce) or maintain (select) the Gi protein- activated state are probably a result of better fit of the mMOR signaling fully active state.30 crystal structure for morphine-like ligands, as well as 1254 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article differences in the molecular scaffold between morphine and establishes its most stable interactions with the N-terminus, fentanyl (Figure 1). On the other hand, difference of size extracellular regions of TM helix 5, TM6, and TM7, and a between co-crystallized morphine-like ligands of activated5 and section of TM3 facing TM5 (Figure 1 and Figure S3). On the inactive29 mMOR models, BU-72 (smaller) and β-FNA other hand, fentanyl establishes more frequent contact with (larger), respectively, creates a wider orthosteric pocket in extracellular loop (ECL) 2, extracellular regions of TM2 and the inactive hMOR model. This may allow ligands larger than TM7, and a section of TM3 facing TM2 (Figure 1 and Figure morphine, like fentanyl, to obtain its preferred binding pose S3). Interestingly, fentanyl comes into contact with faster in this receptor state. Nevertheless, irrespective of initial Cys219ECL2, which is directly involved in an important receptor conformation, both morphine and fentanyl reach structural disulfide bond with Cys1423.25, which may have consistent stable poses in the orthosteric pocket over replicated implications for conformational change in TM3,72 especially as trajectories, which indicate satisfactory docking accuracy this helix has previously been described to undergo an upward axial movement in receptor activation of class A GPCRs.44,73 (Figure 1). In addition to fentanyl and morphine, in the present study, we Differences are not only observed in the conformation of have performed three additional replicas of 3 μs-length morphine and fentanyl but also in orthosteric pocket residues unbiased MD simulations with the MOR-bound antagonist that are in close contact with them (ligand−residue distance naltrexone starting from the activated hMOR state, which <3.5 Å). The distribution of residues most frequently making supports the comparison between morphine and fentanyl. contact, summarized in Table 1, indicates that morphine Naltrexone, similar to morphine and fentanyl, overall maintains a stable disposition in the binding pocket across replicas Table 1. Crystal Ligands, Morphine, Fentanyl, and (Figure S2). However, when compared to control replicas Naltrexone Protein−Ligand interactionsa where naltrexone is bound in the inactive hMOR state, it ligand unique interactions common interactions presents a slightly worse receptor accommodation (average RMSDs of 2.6 Å compared to 2.2 Å, Figure S2) as might be crystal ligands expected for its antagonist condition. Similarities between beta- mTyr1483.33 mAsp1473.32 funaltrexamine naltrexone and morphine scaffolds lead to similar interaction (inactive crystal) mGlu2295.35 environments (Figure 1 and Table 1). Thus, contact between mLys2335.39 naltrexone and the hMOR is similarly and frequently formed mVal2365.42 with TM5, TM6, TM7, and a section of TM3 facing TM5 mTrp2936.48 (Table 1 and Figure S3). However, differently from morphine, BU-72 (activated mPro59N‑ter mAsp1473.32 naltrexone comes into close contact with Tyr3287.43, a residue crystal) mLeu1292.65 mTyr3267.43 also involved in the interaction between fentanyl and the DAMGO (active mGln1242.60 mAsp1473.32 receptor in our MD simulations. As these interactions with crystal) mCys217ECL2 mTyr3267.43 Tyr3287.43 are common to both a potent agonist and an mLys3036.58 antagonist, it could be speculated that they are apparently mTrp3187.35 associated with affinity rather than with efficacy. However, a mIle3227.39 detailed analysis of the simulations reveals a fine structural docked ligands tuning of this residue, which allows the functional distinction morphine hAsp56N‑term hSer57N‑term between agonists and antagonists (see below). hLys2355.39 (i) hAsp1493.32 (i) (ii) (iii) Fentanyl and Morphine Induce Different Receptor hIle2986.51 hTyr1503.33 (i) Activation Pathways. According to RMSD of the trans- hVal3026.55 membrane domain (TMD), the hMOR experiences varying hTrp3207.35 (iii) degrees of backbone conformational change across our MD hIle3247.39 simulations (Table S1 and Figure S4). This is an expected fentanyl hGln1262.60 (iii) hSer57N‑term result as the MOR is known to be a highly dynamic hAsn1292.63 hAsp1493.32 (i) (ii) (iii) protein,5,7,29,30 as shown when comparing both activated and hVal1453.28 hTyr1503.33 (i) inactive crystals5,29 with the fully active crystal structure30 hCys219ECL2 (iii) (TMD RMSDs of 1.4 and 3.0 Å, respectively). Making use of a hTyr3287.43 (ii) (iii) threshold of ± 3.5 Å, calculated from TMD RMSD of crystal naltrexone hMet1533.36 hAsp1493.32 (i) (ii) (iii) structures, we quantified the conformational changes between hLys2355.39 (i) hTyr1503.33 (i) active and inactive states with respect to either its initial state hVal2385.42 (i) or fully active crystal. In general terms, receptor conformation hHis2996.52 is relatively stable (≤ 3.5 Å) in all replicas starting from the hVal3026.55 activated state with the bound agonist (morphine or fentanyl) hTyr3287.43 (ii) (iii) with an average TMD RMSD of 3.2 Å (0.4 SD) from its initial a state and 2.2 Å (0.4 SD) from the fully active crystal30 (Table Residues in close contact (<3.5 Å) listed with respective murine (m) S1 and Figure S4). On the contrary, in MD simulations or human (h) numbering following the Ballesteros and Weinstein starting from the inactive state with the bound agonist generic scheme,31 for crystal structures and MD systems, respectively. (morphine or fentanyl), TMD conformations remain mostly Residues interacting with morphine, fentanyl, and naltrexone are included when frequency of contact is >20% of trajectory time during inactive-like with an average TMD RMSD of 2.8 Å (0.3 SD) at least two MD simulation replicas, independent of the initial state of from its initial state and distant (> 3.5 Å) from the fully active the receptor. Similarity between residues in contact across crystals and crystal30 (TMD RMSD of 3.6 Å (0.4 SD), Table S1 and Figure MD simulations are labeled as common to the (i) inactive crystal, (ii) S4). However, the receptor in replicas #2 and #3 with bound activated crystal, and (iii) fully active crystal. fentanyl or replica #2 with bound morphine starting from the 1255 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article inactive state presents occasional TMD RMSDs of 2.8 and 2.9 Å or 2.6 Å from the fully active crystal,30 respectively, which indicate transient receptor activation toward active-like states (≤3.5 Å, Figure S4). Larger conformational changes in simulations starting from the activated state than those starting from the inactive state are explained by the lack of the intracellular partner in our models. In more detail, computa- tional results also suggest that both activated and inactive crystal structures5,29 (and the models generated from them) are better adapted to morphine than fentanyl because of the similarity between morphine and the co-crystallized ligand scaffolds, which implies that the receptor has to undergo bigger conformational changes to properly bind fentanyl. TMD conformational changes can be compared with those observed when the antagonist naltrexone is bound, where a clear tendency toward inactivation is observed when starting from the activated state (average RMSD of 3.9 Å (0.3 SD) from the initial state, Table S1) and a preservation of the inactive state is found in control replicas starting from the inactive state Figure 2. Conformational change of the TM domain in MD (average RMSD of 2.8 Å (0.2 SD) from the initial state, Table simulations of the inactive hMOR. Extracellular (top) or intracellular S1). Surprisingly, when comparing naltrexone-bound receptor (bottom) views of representative receptor conformations from replicas #3 and #1 starting from the inactive state with bound simulations with the fully active crystal,30 the average RMSD fentanyl (left, in green) or morphine (right, in red), respectively, achieved is only 3.0 Å (0.4 SD) or 3.7 Å (0.3 SD) in compared with the fully active crystal structure (blue). trajectories starting from activated or inactive states, respectively (Table S1). In order to better describe receptor intracellular conformational changes involved in receptor replica #1, the receptor maintains an average active-like activation, four metrics were employed relative to the protein position of +0.3 Å (0.2 SD, Figure 3, Table S1, and Figure center, considering only the two last microseconds of MD S5). In comparison, all replicas starting from the activated state simulations: (i) TM3, which moves upward, (ii) intracellular with bound morphine, despite transient active-like TM3 tips of TM5 and (iii) TM6, which move inward and outward, conformations, show an average less active-like TM3 position respectively; and (iv) NPxxY motif of TM7,43 which rotates of only +0.2 Å (0.2 SD, Figure 3, Table S1, and Figure S5). toward the core of the helical bundle. Meanwhile, a fifth metric Surprisingly, the tendency observed in replicas starting from was used for assessing intracellular separation between TM3 the activated state, where fentanyl selects for more upward and TM6: the distance between residues Arg1673.50 and conformations of TM3 than morphine, is not conserved in Thr2816.34 (Arg−Thr), which is equivalent to the ionic lock in replicas starting from the inactive hMOR. Under this class A GPCRs.43,74−76 Thresholds for these metrics were condition, 2/3 trajectories with bound morphine (replicas #1 based on fully active,30 activated,5 and inactive29 crystal and #3) induce an active-like average shift of +0.3 Å (0.2 SD), structures (normalized according to the inactive crystal while replica #2 shows an even larger average of +0.6 Å (0.2 structure), with a representative sign, which represents the SD), in accordance with the fully active crystal30 (Figure 3, movement involved with respect to the receptor center Table S1, and Figure S5). On the contrary, none of the replicas (inward or outward movement, or downward or upward starting from the inactive state with bound fentanyl presents movement, “−” or “+” sign, respectively): +0.3 Å for TM3, average TM3 movements ≥ +0.3 Å, resulting in mostly −0.8 Å for TM5 and NPxxY, +1.7 Å for TM6, and + 10.0 Å for selecting an inactive-like TM3 conformation (Figure 3, Table Arg−Thr distance as utilized in a previous study.7 These S1, and Figure S5). Interestingly, under these conditions, thresholds were fine-tuned with naltrexone-bound MD replica #3 shows a downward TM3 movement of −0.4 Å (0.2 simulation results to classify active- or inactive-like con- SD, Figure 3, Table S1, and Figure S5), which suggests that the formations of the five intracellular metrics. specific contact fentanyl establishes with Cys219ECL2 does not Upward axial movement of TM3 (Figures 2 and 3, Table S1, reliably invoke TM3 upward axial movement, at least when and Figure S5) is a key indicator of receptor activation, as starting from the inactive state. Rather, direct interactions shown previously in other GPCRs.77,78 Comparison between between the ligand and TM3 appear more influential, as mMOR crystal structures shows that TM3 adopts a more evidenced by the more numerous interactions morphine makes upward position when the receptor is activated: by distances of with TM3 relative to fentanyl. In terms of other TM helices, +0.2 Å and +0.6 Å in activated5 and fully active30 crystals with the distance between the protein center and intracellular tip of respect to the inactive one,29 respectively. Based on these TM5 (see Methods) decreases in activated5 or fully active differences and antagonist-bound MD simulations (average crystal structures30 with respect to the inactive structure29 by movement of +0.1 Å (0.2 SD) and + 0.3 Å (0.2 SD) in replicas −1.5 Å or −2.1 Å, respectively. Considering these crystal starting from activated or inactive states, respectively), we movements and naltrexone-bound MD simulations (average considered a threshold value of +0.3 Å. Our trajectories show movement of +0.5 Å (0.6 SD) and +0.2 Å (0.5 SD) in replicas that replicas #2 and #3 starting from the activated hMOR with starting from activated or inactive states, respectively), a bound fentanyl present an average TM3 upward axial threshold of −0.8 Å can therefore be established for movement of +0.6 Å (0.2 SD), more upward than the initial determining active- or inactive-like conformations of TM5. receptor activated state and adopting a conformation similar to In MD simulations starting from the activated hMOR with the fully active crystal30 (Figure 3, Table S1, and Figure S5). In either fentanyl or morphine bound, as expected, TM5 stabilizes 1256 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 3. Intracellular conformation in MD replicas defined by the Arg1673.50−Thr2816.34 distance, vertical TM3 movement, and lateral NPxxY motif, TM5, and TM6 movements. (a) Structural comparison between the fully active crystal30 (blue) and either activated5 or inactive29 crystals (top, pink or purple, left or right, respectively), or representative conformation of MD-generated replicas #3 and #1 starting from the inactive state with bound fentanyl or morphine, respectively (bottom, green or red, left or right, respectively). As a function of time, (b) active- or inactive-like states of the five different activation-related intracellular features of each replica starting from the activated or inactive hMOR with bound fentanyl, morphine, or naltrexone (active-like conformations of: Arg1673.50−Thr2816.34, NPxxY, intracellular TM6, intracellular TM5, and TM3, in purple, green, yellow, orange, or blue, respectively). White gaps represent inactive-like conformations of each respective metric. See the Methods section for thresholds used to define active- or inactive-like states of intracellular features and Figure S5,S6 for respective metric individual time-dependent analysis. (selects) the active-like conformation with normalized (Figure S5). Likewise, in trajectories starting from the inactive distances mostly between −0.8 Å and − 2.1 Å (Figure 3, state with bound fentanyl or morphine, TM5 does not obtain Table S1, and Figure S5). An exception is replica #3 with an active-like conformation and remains inactive (Figure 3, bound morphine, where TM5 fails to reliably select an active- Table S1, and Figure S5). Moreover, in replica #3, TM5 like conformation and instead permits a change toward an undergoes an average change of +1.7 Å with a maximum of inactive-like conformation with a maximum of +0.7 Å at 2.8 μs +2.9 Å (Figure S5), which places TM5 in a more outward 1257 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article conformation with respect to the receptor core (Figures 2 and active30 crystals present more inward positions of the NPxxY 3, Table S1, and Figure S5). This may be the result of an motif than the inactive crystal 29 (−1.7 and −0.8 Å, outward movement of TM6 (see below), which temporarily respectively). Consequently, in addition to antagonist-bound influences TM5 through their mutual interhelical contact, an MD simulations (average movement of −0.8 Å (0.6 SD) and effect observed in active-like crystal structures.5,30 Exception- +0.1 Å (0.3 SD) in replicas starting from activated or inactive ally, replicas #1 and #2 when morphine is bound in MD states, respectively), a threshold of −0.8 Å can be defined and simulations starting from the inactive state present transient or applied for determining active- or inactive-like conformations sustained inward TM5 movement induction with maximum of this intracellular region (Figure 3, Table S1, and Figure S6). and average displacements of −1.7 Å and −1.5 Å (0.8 SD), Because of computational limitations, our MD simulations do respectively (Figure S5). Altogether, these results indicate that not include bound Gi protein; therefore, it might be expected TM5 generally remains in an active-like conformation if the that activated-like NPxxY conformations could be selected or receptor begins in that state, or if initially inactive, morphine is induced over fully active conformations. During MD more effective than fentanyl at inducing active-like con- simulations starting from the activated state, NPxxY is broadly formations of TM5 (Figures 2 and 3, Table S1, and Figure S5). selected in its original conformation with an average distance Similarly, activated5 and fully active30 crystal structures present of −2.5 Å (0.5 SD), independent of the agonist bound (Figure a more outward TM6 intracellular tip conformation (see 3, Table S1, and Figure S6). On the other hand, all replicas Methods) relative to the inactive crystal,29 with TM6 starting from the inactive state with bound fentanyl induce a respectively moving +2.1 Å and +2.0 Å away from the protein conformational change so that NPxxY achieves a stable active- center. Based on these differences and naltrexone-bound MD like position (average of −1.4 Å, 0.6 SD), which is the greatest simulations (average movement of +0.8 Å (0.5 SD) and +0.6 Å in replica #3 (average of −1.8 Å, 0.7 SD, Figure 3, Table S1, (0.4 SD) in replicas starting from activated or inactive states, and Figure S6). Conversely, when morphine is bound, only respectively), a distance threshold of +1.7 Å can therefore be replica #1 starting from the inactive state induces an active-like used for defining active- or inactive-like conformations of position that crosses the −0.8 Å threshold (average of −1.5 Å, TM6. This helix selects an active-like outward conformation 0.3 SD, Figure S6), whereas other replicas select an inactive- with an average distance of +3.5 Å (0.6 SD) in trajectories like conformation. These results suggest that fentanyl starting from the activated hMOR, independent of the agonist consistently induces active-like conformations of both TM7 bound (Figure 3, Table S1, and Figure S6). However, a and the NPxxY motif stronger than morphine (Figure 3). The transient deactivation of TM6 occurs when it falls into an intracellular distance between Arg1673.50 and Thr2816.34, inactive-like conformation in replica #1 with bound morphine equivalent to the ionic lock in class A GPCRs,43,74−76 can be at 1.8 μs, recording a distance of +1.4 Å (Figure S6). Likewise, classified as being in an active- or inactive-like conformation by in all replicas starting from the inactive hMOR with bound applying a threshold value of +10.0 Å, as used previously.7 The morphine, TM6 does not demonstrate any outward move- Arg−Thr distance is broadly selected in its original state in ment, maintaining only an average distance of +0.4 Å (0.4 SD) respective MD simulations starting from activated or inactive with respect to the inactive crystal and below the +1.7 states, independent of the agonist bound, resulting in threshold (Figure 3, Table S1, and Figure S6). However, with bound fentanyl, active-like conformations of TM6 are induced respective averages of 13.0 Å (0.9 SD) and 6.5 Å (0.5 SD, in replicas #2 and #3 with maximum values of +2.3 Å and +3.2 Figure 3, Table S1, and Figure S7). However, in replica #3 Å, respectively (Figures 2 and 3, Table S1, and Figure S6). when starting from the inactive state with bound fentanyl, this While in replica #2 these conformational changes are observed interaction is seen to transiently break at 1.8 μs when the transiently, a more sustained conformational change occurs in distance reaches 10.1 Å (Figure S7). In the replicas where replica #3 with an average TM6 movement of +1.9 Å (0.6 SD) naltrexone is bound, Arg1673.50 and Thr2816.34 reach and over the last microsecond. Altogether, these results indicate stabilize distances <10 Å over the entire two last microseconds that while both fentanyl and morphine are able to select an of MD simulations, started from either control inactive or active-like conformation of TM6 in simulations starting from activated states, as expected by the antagonist nature of the activated state, uniquely fentanyl is able to induce TM6 naltrexone, contrary to the overall active-like states selected by activation when starting from the inactive state, which leads agonists when starting from the activated state (Figure 3, Table TM6 to adopt conformations more similar to activated and S1, and Figure S7). Similarly, naltrexone induces an inward fully active crystals5,30 (Figures 2 and 3, Table S1, and Figure movement of TM6, in parallel to an outward movement of S6). This may represent an indication of the greater efficacy of TM5 and the NPxxY region, in simulations starting from the fentanyl relative to morphine. activated state, thus adopting, as expected, inactive-like Rearrangement of TM helices during receptor activation conformations of these metrics (Figure 3, Table S1, and leads to specific intracellular conformational changes such as Figure S7). This result indicates that the threshold used clearly the NPxxY motif on TM7 and Arg167 3.50 −Thr281 6.34 differentiates between active- or inactive-like states of this interaction (Figure 3). The fully active crystal structure of interaction. Altogether our results strongly indicate that the mMOR30 allows a better understanding of the conforma- receptor activation by means of TM6 outward movement tional changes adopted by these regions when the intracellular (Figures 2 and 3, Table S1, and Figure S5-S7) is more Gi protein partner is bound to the receptor. Interestingly, in frequently observed with bound fentanyl, which is uniquely the fully active crystal,30 the NPxxY motif adopts an alternative able to initiate in a μs time period, than with morphine. conformation closer to the inactive crystal structure29 than to Therefore, in addition to conformational changes of helical the activated crystal structure.5 This difference may be due to regions mentioned above, the different frequency of activation the non-native binding of Nb39 to the receptor. In these of these intracellular metrics is consistent with the crystal structures, the NPxxY position can be described experimentally observed difference in efficacy between according to the protein center. As such, activated5 and fully morphine and fentanyl. 1258 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 4. Density maps of the outward movement of TM6 and conformation of the NPxxY motif with respect to the protein center. Density maps of fentanyl or morphine (top or bottom row, respectively) of the three replicas started from activated or inactive states (left or right column, respectively). Inactive, activated, and fully active crystals are represented for comparison by dark blue, light blue, or green, respectively. Density maps show conformations more or less frequently induced represented with a hot to cold color gradient. All these TMs and intracellular metric movements (average of the absence of TM3 upward axial movement, which, in tendencies summarized in Table S1) observed during receptor addition to TM6 outward movement, is required to properly (in)activation are indicative of conformational changes at the break intracellular TM3−TM6 interactions. Altogether, these intracellular side of the receptor in system dynamics starting results are indicative of different receptor activation pathways either from activated or inactive states. As shown in Figure 3, induced by morphine (TM3 and TM5) and fentanyl (TM6 some of these activation features co-exist in time. Therefore, and TM7). their simultaneous presence represents a highly active It has been widely discussed if fentanyl acts as a β-arrestin conformation of the receptor. However, the activation of one biased agonist or behaves as a balanced agonist.18,23,61−68 As feature prior to another can be an indication of the exploration done for other GPCRs,77,78 we performed a density map of a certain process or pathway of receptor activation. Our (Figure 4) of the receptor intracellular state achieved results indicate that all active-like conformational intracellular considering two of the main representative intracellular metrics are mostly preserved in MD simulations starting from metrics, which conformationally differ more between morphine the activated agonist-bound receptor state, which is biologically and fentanyl: (i) outward movement of the intracellular side of reasonable because morphine and fentanyl are both agonists. TM6 and (ii) conformation of NPxxY with respect to the However, agonist-specific differences are consistently observed protein center. This representation (Figure 4) suggests that in during receptor activation from the inactive state. While with our MD simulations, we reach a common active-like state with bound morphine, the hMOR largely presents an upward axial either fentanyl or morphine bound. It should be considered movement of TM3 and exclusively induced TM5 active-like that our models starting from the activated hMOR were conformations, with bound fentanyl, the hMOR shows specific modeled from its respective crystal,5 which contains a G activation of TM6 and stronger activation of the NPxxY motif protein mimetic camelid nanobody subsequently removed in on TM7 (Figure 3). Interestingly, replica #3 starting from the our model. Therefore, the conformations achieved in our inactive state with bound fentanyl, which most frequently trajectories can be conditioned by the initial crystal used, presents active-like conformations of TM6 and the NPxxY which may limit the observations of other alternative motif, is the only replica that yields transient active-like conformations. However, in a timescale of three microseconds, Arg1673.50−Thr2816.34 distances (Figure 3). The non-preser- our simulations have not identified receptor conformational vation of this Arg1673.50−Thr2816.34 separation may be a result regions consistent with biased agonism, in agreement with a 1259 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 5. Orthosteric pocket conformational changes of the hMOR in MD simulations. (a) and (b) Predominant trans (∼180.0°) or cis conformation (∼300.0°) of (a) Gln1262.60 χ2 dihedral angle or (b) Lys2355.39 χ2 dihedral angle, in representative MD conformations with bound fentanyl or morphine (purple or pink, hMOR residues in green or red, left and right, respectively). Structural images also show different Asp1493.32 χ1 dihedral angle cis and trans conformations selected, and Asp1493.32−ligand, Tyr1503.33−ligand, or Asp1493.32−Tyr3287.43 H-bond formation. (c) Conformation frequency of the Asp1493.32 trans χ1 dihedral angle or Gln1262.60 and Lys2355.39 cis χ2 dihedral angle selected with bound fentanyl, morphine, or naltrexone, independent of the initial state of the receptor (green, red, or purple, respectively). (d) H-bond occupancy of Asp1493.32− Tyr3287.43, Lys2355.39−ligand, or Tyr1503.33−ligand interaction of simulations with bound fentanyl, morphine, or naltrexone, independent of the initial state of the receptor (green, red, or purple, respectively). (e) Structural comparison of Met1533.36 conformation between the fully active crystal30 (blue) and either replica #3 starting from the inactive state with bound morphine or fentanyl (red or green, left or right, respectively). (f) Frequency of Met1533.36 gauche− conformation in trajectories with bound fentanyl, morphine, or naltrexone, starting from activated or inactive states (green, red, or purple, respectively). recent experimental study.65 It is worth noting that conclusions more prone to β-arrestin binding is a speculation to be on biased agonism from experimental functional studies may considered in further studies. depend on the parameter used for its quantification, where the Morphine and Fentanyl Mediate Different Orthos- proposed ΔΔlog(τ/KA), ΔΔlog(τ), or ΔΔlog(Emax/EC50) bias teric Pocket Conformational Changes. To ascertain which factors have led to different biased profiles when applied to conformational changes in the orthosteric pocket, induced or fentanyl, with Gi protein bias, and β-arrestin with and without selected by morphine and fentanyl (see Figure 1), are most bias, respectively.23,65 Interestingly, density maps of simu- responsible for observed intracellular conformational differ- ences, we performed an energetic analysis of residues in close lations starting from the inactive state show that fentanyl contact with each ligand (Figure S8). Accordingly, from reaches more intermediate intracellular states than morphine, residues in close contact with both ligands (summarized in the latter presenting a larger density closer to the inactive Table 1), fentanyl makes stronger energetic interaction with crystal conformation29 (Figure 4). This result may suppose an Asp56N‑term, Gln1262.60, Cys219ECL2, and Tyr3287.43, whereas additional indication of the differences in efficacy observed morphine makes stronger interactions with Tyr1503.33 and between both ligands.51−53 On the other hand, the possibility Lys2355.39 (Figure S8). Both ligands establish the highly that the higher trend of fentanyl of generating intermediate energetic and conserved MOR-opiate salt bridge with states may facilitate the generation of receptor conformations Asp1493.32, as previously identified in other studies.5,7,19,29,30,45 1260 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article From a conformational perspective, we observe (i) H-bonds In addition to Asp1493.32−ligand interaction, greater (either direct or water-mediated) between the ligand and energetic interactions can also be observed between Asp1493.32 or Tyr1503.33 and between residues Asp1493.32 and Tyr1503.33 and morphine with respect to bound fentanyl Tyr3287.43, (ii) specific χ1 dihedral angles of Asp1493.32 and (average potential energy (P.E.) of −5.6 kcal/mol (2.4 SD) or Met1533.36 (directly determined by conformational changes of − 1.9 kcal/mol (2.7 SD), respectively, Figure S8). Differences Tyr1503.33), and (iii) specific χ2 dihedral angle selections of in the energetic landscape between fentanyl and morphine Gln1262.60 and Lys2355.39 sidechains (Figure 5 and Figures relate to their different scaffold, which presents a single oxygen S9−S12). Taking these metrics one by one, morphine and group in the former (carbonyl group) and three oxygen groups fentanyl have both been described to make a highly conserved in the latter (one ether and two hydroxyl groups, Figure 1). MOR-opiate electrostatic interaction with Asp1493.32, a residue This distinction translates to a different H-bond occupancy located on TM3.5,7,19,29,30,45 In our MD simulations, this (either direct or water-mediated) between Tyr1503.33 and interaction between the ligand and Asp1493.32 is highly fentanyl or morphine with average occupancies of 32.2% (7.2 energetically favorable with bound fentanyl but even greater SD) or 89.2% (9.4 SD), respectively (Figure 5). In the case of with bound morphine, independent of the initial receptor state the antagonist naltrexone, which presents a scaffold similar to (average interaction energy of −78.1 or −94.9 kcal/mol, morphine (Figure 1), it adopts a Tyr1503.33−naltrexone H- respectively, Figure S8). Despite this energetic difference bond occupancy of 56.1% (9.6 SD, Figure 5). Tyr1503.33 and between morphine and fentanyl, both ligands maintain a stable morphine or naltrexone hydroxyl/carbonyl groups are interaction (<4.5 Å) over three microseconds between their considered either donors or acceptors when evaluating amine group and Asp1493.32 (Figure 5 and Figure S9). In Tyr1503.33−ligand H-bond presence/absence. The higher H- addition, Asp1493.32 presents different χ1 dihedral angle bond occupancy between morphine and Tyr1503.33 enhances conformations when morphine or fentanyl is bound in the upward axial movement of TM3 (Figure 3, Table S1, and hMOR (Figure 5 and Figure S9). Two different Asp1493.32 Figure S5) and exerts a stabilizing effect in the bottom of the conformations are observable in our MD simulations: (i) a receptor orthosteric pocket, as observed with the conformation trans conformation (χ1 dihedral angle >115.0°), predom- of Met1533.36 (Figure 5 and Figure S10). In MD simulations inantly induced in 5/6 trajectories with bound fentanyl starting from the activated state, independent of the agonist (average frequency of 89.2%, Figure 5 and Figure S9), which bound, the Met1533.36 χ1 dihedral angle generally adopts a positions the sidechain closer to TM2 and (ii) a gauche− gauche− conformation (∼60.0°). Likewise, in MD simulations conformation (χ1 dihedral angle <115.0°) induced by starting from the inactive state with bound morphine, morphine in all trajectories (average frequency of 93.7%, Met1533.36 gauche− conformation is selected with a frequency Figure 5 and Figure S9), which positions the sidechain away of 93.2% (8.8 SD, Figure 5 and Figure S10). However, when from TM2. Similar to morphine-bound replicas, when the fentanyl is bound, alternative conformations of Met1533.36 (χ1 antagonist naltrexone is bound to the hMOR, an overall > 90.0°) are more frequently induced with a frequency of Asp1493.32 gauche− conformation is induced, independent of 65.3% (4.9 SD, Figure 5 and Figure S10). Interestingly, in the initial state of the receptor (Figure 5 and Figure S9). Under naltrexone-bound conditions, Met1533.36 gauche− conforma- this condition, conformational change from gauche− to trans tion is overall selected only with a frequency of 58.3% (24.5 conformation is only observed in replica #3 starting from the SD, Figure 5 and Figure S10), which indicates that this residue inactive state at 2.2 μs. As a result of the trans conformation conformationally fluctuates more with bound naltrexone than predominantly observed with bound fentanyl, Asp1493.32 is when fentanyl or morphine is bound in the hMOR, probably as within closer H-bonding distance with Tyr3287.43 located on a result of the direct contact naltrexone makes with this TM7.46 Activated and fully active crystal structures5,30 show sidechain (Table 1 and Figure S3). Selection of χ1 gauche− closer Asp1493.32−Tyr3287.43 interaction distances than the appears to be associated with aforementioned active-like inactive crystal structure29 (3.2 and 3.4 Å, respectively), and conformations of TM3 (Figure 3, Table S1, and Figure S5). this interaction also differs in our MD simulations with either This process can be linked with fentanyl-mediated TM6 bound morphine or fentanyl (Figure 5 and Figure S9). activation where the receptor core undergoes reorganization Throughout 4/6 trajectories with bound fentanyl, independent and becomes more flexible. Once this transition ends, the of the initial receptor state, Asp1493.32−Tyr3287.43 maintains a stability of the core can be recovered, as indicated by closer interaction (<4.5 Å) with an H-bond occupancy of trajectories starting from the activated state (Figure 5 and 89.5% (36.1 SD, Figure 5 and Figure S9). Conversely, across Figure S10). Conversely to the low H-bond occupancy all trajectories with bound morphine, only transient between fentanyl and Tyr1503.33 in TM3, which conditions interactions between Asp1493.32 and Tyr3287.43 are observed the conformation of Met1533.36, an energetically favorable with a H-bond occupancy of 24.1% (22.4 SD, Figure 5 and interaction with Gln1262.60 on TM2 is established with an Figure S9). Similarly, the gauche− conformation overall average P.E. of −3.0 kcal/mol (4.4 SD). This contrasts with an selected when naltrexone is bound leads to an Asp1493.32− unfavorable average P.E. of +2.7 kcal/mol (2.2 SD) displayed Tyr3287.43 H-bond occupancy of 41.1%, independent of the by morphine (Figure S8). Subsequently, Gln1262.60 shows starting receptor state (19.5 SD, Figure 5 and Figure S9). conformational differences when morphine or fentanyl is These results are supported by the low H-bond occupancy bound in the hMOR (Figure 5 and Figure S11). In general, observed when the antagonist naltrexone is bound into the two different Gln1262.60 χ2 dihedral angles are observed in our orthosteric pocket of the receptor (Asp1493.32−Tyr3287.43 H- MD simulations: (i) a cis and (ii) trans (∼300.0° and ∼ 180.0°, bond occupancy of 41.6% across all replicas, Figure 5), which, respectively). By using a threshold of ± 240°, Gln1262.60 can despite being in close contact with Tyr3287.43 like fentanyl, be categorized according to these two conformations (Figure 5 presents larger Asp1493.32−Tyr3287.43 distances. Therefore, and Figure S11). Although both conformations are observable fentanyl mediates a stronger Asp1493.32−Tyr3287.43 H-bond with each ligand, its frequency is noticeably different from one than morphine in the orthosteric pocket. to another. Systems with bound fentanyl predominantly induce 1261 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article the trans conformation 86.0% of the time (10.5 SD, Figure 5 and Figure S11), which positions Gln1262.60 toward the orthosteric pocket, allowing a three-way sidechain interaction with Tyr3287.43 and Asp1493.32 and supporting closer Asp−Tyr interaction (Figure 5 and Figure S13). By contrast, morphine selects the Gln1262.60 cis conformation with a frequency of 80.5% (13.5 SD, Figure 5 and Figure S11), which positions it toward the membrane and negatively affects the Asp−Tyr interaction (Figure 5 and Figure S13). Similar to morphine, naltrexone selects the Gln1262.60 cis conformation with a frequency of 87.3% (10.0 SD, Figure 5 and Figure S11). Finally, Lys2355.39−agonist interaction is the only one established between morphine and TM5 and is present at least 20.0% of the time in 5/6 trajectories (Table 1 and Figure S2). This interaction shows energetic differences between morphine and fentanyl with an average P.E. of +6.1 kcal/mol (6.9 SD) or +16.5 kcal/mol (3.8 SD), respectively (Figure S7). Similar to Tyr1503.33, the Lys2355.39 amino group establishes an interaction with any of morphine’s oxygen groups (ether, hydroxyl groups #1 and #2) with a H-bond occupancy of 63.5% (13.2 SD) but generally avoids interaction with fentanyl’s carbonyl group (H-bond occupancy of 7.2% (6.9 SD), Figure 5 and Figure S12). When naltrexone is bound, which shares a similar scaffold with morphine (Figure 1), Lys2355.39 H-bonds naltrexone with an occupancy of 43.9% Figure 6. Trp2956.48 conformational change in MD simulations with respect to MOR crystals. Comparison between different conforma- (18.6 SD). In this specific interaction, ligand oxygen groups are tions of Trp2956.48 (χ2 dihedral angle represented by spheres) of fully only considered as acceptors. Accordingly, Lys2355.39 displays active, activated, and inactive crystals (salmon, blue and purple, two different χ2 dihedral angles, either adopting a cis or trans respectively) and (i) representative MD-generated hMOR conforma- conformation (∼300.0° or ∼ 180.0°, respectively) where a tions with bound fentanyl (green, middle row) of replica #2 starting threshold of ± 240° can be used to categorize these two from the activated state, replicas #3 and #2 starting from the inactive sidechain conformations (Figure 5 and Figure S12). state (from left to right, respectively) or, ii) with bound morphine Independent of the bound ligand, the trans conformation is (red, bottom row), replica #3 starting from the activated state, or the most commonly selected with the sidechain interacting replicas #1 and #2 starting from the inactive state (from left to right, with the extracellular solvent rather than with the ligand. respectively). However, different rates of cis conformation are induced between fentanyl, morphine, and naltrexone, with average Trp2956.48 trans conformation is stabilized within 1250 ns percentages of 11.5% (5.1 SD), 38.0% (13.9 SD), and 40.0% (Figures 6 and 7 and Figure S14). This enables a novel H- (14.9 SD), respectively (Figure 4 and Figure S12). In this case, bond between Trp2956.48 and Ala2425.46 (sidechain amino the cis conformation permits interaction with the bound ligand. group and backbone carbonyl group, respectively) with an This indicates that morphine mediates different conforma- average occupancy of 97.7% (1.1 SD), which is beneficial for tional dynamics of Lys2355.39 (Figure S13), which, accordingly, maintaining active conformations of TM5 and TM6 (Figure helps to stabilize the binding pose of morphine inside the 7). Conversely, all replicas starting from the inactive state with orthosteric pocket of the receptor. bound fentanyl, as well as replicas #2 and #3 where morphine Conformational Changes of Trp2956.48 Inducing is bound, select an inactive-like Trp2956.48 gauche− con- Different Bottom-Orthosteric Pocket Receptor Inter- formation (Figure 6 and Figure S14). This conformation actions Are Relevant in TM5−TM6 Packing. During enables H-bonding with the adjacent residue on TM7, activation of the MOR, TM6 moves outward, which, together Asn3307.45, which enhances TM6−TM7 interaction and avoids with conformational rearrangement of TM3, TM5, and TM7, the interaction of this residue with water molecules, thus creates solvation of the intracellular cavity5 and conformational inhibiting TM6 activation (Figure 7). Interestingly, in replica changes in specific residues on TM6, such as Trp2956.48.5 By #1 starting from the inactive state with bound morphine, the classifying sidechain conformation by dihedral angles, gauche− activated-like trans conformation is induced at 400 ns, which is (0−120°), trans (120−240°), or gauche+ (240−360°) in the only replica in this category where the NPxxY motif yields mMOR crystals, Trp2956.48 χ2 changes from gauche− to trans an active-like conformation, observed from 200 ns onward to gauche+ in inactive, activated, and fully active crystals, (Figure S6). This enables sidechain−backbone Trp2956.48− respectively (values of 78.6°, 121.3°, and 342.9°, respectively). Ala2425.46 H-bonding between TM5−TM6 with an occupancy In our MD simulations, we observe two stable Trp2956.48 of 68.0%, which assists TM5 activation after 500 ns (Figure conformations (Figure 6 and Figure S14): (i) an activated-like S5), even though TM6 activation is not observed. Despite trans conformation at 225.6° (6.9 SD) and (ii) an inactive-like active-like conformations of the NPxxY motif being present in gauche− conformation at 97.1° (4.1° SD). We do not observe all replicas starting from the inactive state with bound fentanyl, the gauche+ conformation probably because, due to computa- only replica #3 transiently induces the trans conformation of tional limitations, we have not included an interacting Gi Trp2956.48 during the first 500 ns before returning to gauche− protein in our simulations. In all trajectories starting from the (Figure S14). This conformational change results in only 6% activated hMOR state, independent of the agonist bound, the Trp−Ala H-bonding occupancy, but repeated sidechain 1262 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 7. TM5−TM6 interaction in MD simulations: Trp2956.48 and Ala2425.46. MD-generated receptor structures showing the interaction of Trp2956.48 with Ala2425.46 or Asn3307.45, in (a) replica #3 or #2 starting from inactive or activated states with bound fentanyl or morphine, respectively (left or right, green or red, respectively) or (b) replica #2 starting from the inactive state with bound fentanyl (green). (c−e) Per replica average (c) Trp2956.48 χ2 dihedral angle, (d) Trp2956.48−Ala2425.46 distance, and (e) Trp2956.48−Ala2425.46 H-bond occupancy (shades of green, red, or purple for replicas #1−3 with bound fentanyl, morphine, or naltrexone, respectively). switching is enough to trigger the gradual activation of TM6, these differences, closer interactions between Trp2956.48 and which is observed more strongly later in this trajectory (Figure Ala2425.46 are observed with active-like conformations of TM5, S6). Indeed, Trp2956.48 and Ala2425.46 mostly interact via TM6, and the NPxxY motif (Figure 7), which demonstrates hydrophobic contact for the rest of the trajectory (closest the importance of rearrangements in TM5−TM6 packing. For distance: 4.7 Å, Figure S14) rather than H-bonding because comparison purposes, naltrexone-bound trajectories induced Trp2956.48 remains in its gauche− conformation. This suggests or selected gauche− conformations from respective activated or that although bound morphine is capable of eliciting similar inactive states with percentages of 74.6 and 88.8%, respectively, conformational changes in Trp2956.48, switching from gauche− with the closest Trp2956.48−Ala2425.46 distance of 8.4 Å to trans, its effect can be different depending on other residues, achieved in control simulations starting from the inactive state. such as Met1533.36. In particular, it is noticeable in replica #3 This distance, larger than that observed in fentanyl-bound with bound fentanyl that Met1533.36 χ2 conformational trajectories, suggests a lack of TM5−TM6 packing, which switching occurs at the same time as Trp2956.48 (Figures S10 favors the inactive state of the receptor. In addition, we observe and S13) with these switches in Met1533.36 conformation different distributions between Asp1493.32−Tyr3287.43 H- becoming more frequent as TM6 moves outward. Despite bonding in the orthosteric pocket and Trp2956.48−Ala2425.46 1263 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article interaction, where average Trp−Ala distances of 3.7 Å (1.8 Thr2816.34 distance state (p < 0.001 in all cases, Figure S16), SD) or 6.8 Å (3.0 SD) are reported when the Asp−Tyr H- showing overall Trp2956.48−Ala2425.46 closer distances in bond is formed or broken, respectively (Figure S13). This active- than inactive-like states. Close Trp2956.48−Ala2425.46 demonstrates the interconnectivity of residue−residue inter- distances (see Methods), which are mainly attributed to helical actions in the orthosteric pocket, which may have different rearrangement, enhance the establishment of hydrophobic implications for intracellular receptor conformations. contact between both residues. These results highlight the Intracellular Conformational Changes Correlate with importance of these two residues becoming physically close Ligand-Specific Conformational Changes in the Re- irrespective of their H-bonding status. ceptor Orthosteric Pocket. In order to elucidate how Altogether, these results indicate that, independent of the ligand-specific orthosteric pocket conformational changes ligand bound and the starting conformation of the receptor for determine intracellular receptor conformations, we performed MD simulations, mostly all orthosteric metrics studied in this a number of statistical analyses (Figures S15 and S16) with study statistically reveal significant differences between active- respect to (in)activation of specific intracellular receptor and inactive-like states of the intracellular metrics used, thus metrics: TM3 upward-, TM5 inward-, TM6 outward-, and establishing an (in)activation communication connection NPxxY motif inward-movement and Arg1673.50−Thr2816.34 mechanism through the entire receptor. The differences distance. The chi-square (chi2) test was chosen to examine found between agonists in their propensity to select or induce the association between these intracellular metrics and specific receptor conformations through these mechanistic orthosteric pocket metrics: sidechain−sidechain Asp1493.32− structural features may explain the differences in efficacy Tyr3287.43, Lys2355.39− or Tyr1503.33−ligand, sidechain− observed experimentally.27,51−53 backbone Trp2956.48−Ala2425.46 H-bond formation and Internal Receptor Polar-Water Network Differs Met1533.36 gauche− conformation. The chi2 analyses involve between Morphine and Fentanyl. A polar network two-way frequency tables of binary variables (presence/ mediated by water molecules has been described to be absence of H-bonds or Met1533.36 gauche− conformation involved in signal transmission from the extracellular to the and active/inactive state of each intracellular metric, Figures intracellular side of the mMOR.5 An average water density S15 and S16). In addition, the potential association between map of the hMOR reveals differences in this network between the Trp2956.48−Ala2425.46 distance and active- or inactive-like activated/inactive states in morphine/fentanyl-bound MD state of the selected intracellular metrics was measured using simulations (Figure 8). As might be expected, the different Student’s t-tests (Figure S16). All these statistical analyses were molecular size and binding pose of fentanyl compared to performed independently of the initial state of the receptor and morphine affects the number of water molecules that can enter the ligand bound. Our objective was not to find differences between ligands but to detect general mechanistic trends of receptor activation. Because we are combining data from different ligands and MD simulations with different starting states, we should take these tests as exploratory rather than confirmatory of the tested hypotheses. Nevertheless, the consistency found in the results suggests these analyses are valuable to realize the apparent relationships between the selected structural features. Chi2 tests show that differences exist between active- and inactive-like intracellular receptor conformations and Asp1493.32−Tyr3287.43, Tyr1503.33−ligand, Lys2355.39−ligand, Trp2956.48−Ala2425.46 H-bond formation and Met1533.36 gauche− conformation (p < 0.05 in all cases, Figures S15 and S16). The absence of statistical significance was found for the association between Tyr1503.33−ligand H- bond formation and the intracellular TM6 state (p = 0.051, Figure S15) and between Lys2355.39−ligand H-bond formation and either TM3 or TM5 movements (p = 0.9 or 0.7, respectively, Figure S15). Interestingly, the fact that this last H- bond formation is not associated with movements of TM5 (Figure S15) suggests that Lys2355.39 affects the position of the ligand in the orthosteric pocket, which in turn determines the conformation of other residues (such as Asp1493.32), rather than directly transmitting the signal through TM5. Large chi2 values observed in the association between the sidechain− backbone Trp2956.48−Ala2425.46 H-bond and intracellular metrics (values between 252.7 and 582.9, which are extremely significant because 2 × 2 contingency tables contain 1 degree Figure 8. Polar water network in MD simulations. Comparison of the water network connecting extracellular and intracellular sides of the of freedom and significant P < 0.05 values are reached if chi2 > hMOR with bound fentanyl or morphine (green or orange, 3.841, Figure S16) highlight the relevance of this interaction in respectively). Receptor conformations from replicas #1 and #1 the (in)activation of the hMOR. Similarly, Student’s t-tests starting from the activated hMOR (snapshots at 2.4 μs and 3.0 μs, show that the Trp2956.48−Ala2425.46 distance statistically respectively) with (a) fentanyl and (b) morphine and replicas #3 and differentiates between active and inactive TM3, TM5, TM6, #1 starting from the inactive hMOR (pictures at 2.2 μs and 1.8 μs, and NPxxY motif conformational states and the Arg1673.50− respectively) with (c) fentanyl and (d) naltrexone. 1264 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article the orthosteric pocket. Therefore, we observe more waters signaling to either one effector or another. This bar-coded when morphine is bound than with fentanyl, independent of phosphorylation signaling affects indirectly G protein coupling the initial receptor state (Figure 8). In terms of the receptor because arrestin- and G protein-coupled receptor populations core or intracellular regions, two different water clusters can be coexist and the increase of one population decreases the identified: (i) between Trp2956.48 and Asn3307.45 on TM6 and other.82 However, phosphorylation patterns may vary between TM7, respectively, and (ii) adjacent to the intracellular GPCRs, in particular for the MOR. For this receptor, it was N7.49PxxY7.53 motif on TM7. In MD simulations with bound found that a single threonine (T180) in ICL2 is fundamental fentanyl, greater solvation is observed at both these regions, in for agonist-dependent receptor phosphorylation and subse- particular near the NPxxY motif, which has previously been quent arrestin binding, activation, and signalling.83 In the described as important for active conformations of the mMOR present study, the examination of MD simulations of the or class A GPCRs, in general.5,79−81 This is consistent with hMOR with either fentanyl or morphine bound did not find stabilization of an active-like receptor conformation or receptor differences with respect to T180. It is worth noting that with activation from the inactive state, which is observed in replica present computational means, it is not possible to assess which #3 with bound fentanyl (Figure 8). On the contrary, with active receptor conformations are more prone to be bound morphine, these same water clusters are less phosphorylated, which ultimately lead to β-arrestin signaling. pronounced, even with greater solvation in the orthosteric However, at the timescales and starting models used, our MD pocket, which suggests comparatively reduced receptor simulations suggest the presence of a single common active- activation or increased destabilization of the active state. like receptor state stabilized by both morphine and fentanyl, in This is the case even when an active-like trans conformation of comparison to multiple conformations observed in other Trp2956.48 is induced, which allows sidechain−backbone H- GPCRs such as the adenosine A2A receptor.78 This leads us to bonding with Ala2425.46 (Figure 8), as observed in replica #1 consider an unbiased behavior of both morphine and fentanyl, when starting from the inactive state. Furthermore, in this where both ligands achieve a common active-like receptor state same replica, despite corresponding activation of the NPxxY by differently activating ligand-specific conformations of the motif (Figure 3), the lack of TM6 conformational change as a receptor. Nevertheless, it cannot be completely discarded that whole allows fewer water molecules to access this region the higher propensity of fentanyl to induce receptor compared to simulations with bound fentanyl. In the same intermediate states could facilitate the generation of con- direction, when the antagonist naltrexone is bound, solvation formations predisposed to β-arrestin binding. In addition, it is of the two aforementioned regions is reduced even in worth noting that reliable quantification of biased agonism simulations starting from the activated state, which facilitates needs the use of proper parameters, and in this regard, the hydrophobic interactions, TMD rearrangements, and inactiva- following has been proposed: the use of operational efficacy tion of the receptor (Figure 8). This indicates that proper (τ) versus the transduction coefficient (τ/KA) in the case of conformational movements of both TM6, most importantly comparing ligands producing different maximum responses or, through Trp2956.48, and the NPxxY motif are necessary for in other words, when partial agonism is present64 or, moreover, formation of a water column in the receptor core, which is the inclusion of constitutive receptor activity in the thought to be necessary for activation of the MOR.5,79−81 mathematical modeling in those cases where basal receptor ■ DISCUSSION The recent spate of deaths from fentanyl and its derived response is observed.17,84 On the contrary, other authors have applied the Emax/EC50 parameter for bias calculation of μ- opioid agonists.65 Thus, it may happen that the mathematical compounds51−53,56,57 raises special interest in how these high- models used for agonist bias quantification, in particular in efficacy agonists interact with the hMOR compared to lower- those systems where pharmacological complexity either at the efficacy agonists such as morphine. Questions about potential ligand (partial agonism) or receptor (constitutive activity) differences between morphine- and fentanyl-based signaling levels is present, may have permitted the occurrence of lead to the need for better understanding of how these two conflicting results in some cases. Nevertheless, a connection ligands modulate MOR conformation and function and how between functional experimental results yielded by ligands this could instigate stronger responses from the receptor51−53 from particular receptors and the molecular interactions of or even stimulate different downstream pathways.25,64 There these ligands with the residues in the receptor binding pocket are conflicting experimental data about whether fentanyl is a is expected. The observed functional response is a consequence biased agonist or not. It was proposed through in vitro studies of molecular events that occur first within the receptor and are that fentanyl is more biased toward β-arrestin activation than later propagated to the effector system. morphine,23,25,64 but the recent literature65 did not find The results obtained in this study greatly depend on the significant bias of neither fentanyl nor a collection of fentanyl interactions established by morphine and fentanyl in the MOR analogues toward Gi protein or β-arrestin signaling pathways. orthosteric pocket, which to be properly determined requires a Interestingly, opposite bias profiles were obtained for fentanyl sufficiently accurate initial docking pose, as has been previously depending on the model parameter used for bias quantifica- demonstrated in MD simulations of other GPCR-ligand tion: Gi protein bias or β-arrestin bias when τ/KA or τ complexes.86 In this regard, our accurate docking of co- parameters were respectively chosen.64 MD simulations cannot crystallized ligands (or their close analogues) gives us reason to give a definite answer to this problem. GPCR-dependent β- have confidence in our docking protocol. We consider that, in arrestin signaling involves the phosphorylation of particular order to properly characterize the efficacy of fentanyl and serine or threonine residues at the intracellular regions of the morphine in humans, it is important to study the hMOR receptor by GPCR kinases (GRKs). In a recent study on the instead of the mMOR, thus avoiding protein−ligand dopamine D1 receptor,82 it was found that intracellular loop 3 interactions with non-conserved residues. At the time our (ICL3) phosphorylation affects arrestin binding and activation models were generated, only crystals of the activated5 and through various phosphorylation patterns, which direct the inactive29 mMOR were available. Models of activated or 1265 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article inactive hMOR states were successfully validated by high morphine-like scaffold of agonist BU-72 and antagonist β-FNA conformational stability in MD simulations of the receptor and in the activated or inactive mMOR,5,29 respectively. Con- bound control ligands: BU-72 and naltrexone. However, a versely, the native fentanyl binding pose is still unclear, with third mMOR crystal structure30 was later released, which previously proposed binding poses differing in their details the fully active receptor conformation coupled to orientation,26,90−96 and no consensus having been reached. human Gi protein. Despite not providing data of the β-arrestin Mutagenesis-based studies46,97 first characterized the binding signaling pathway, it supposes a novel point of reference for landscape of fentanyl and fentanyl derivatives involving comparison with our MD simulations in terms of Gi protein Asp1493.32, Tyr1503.33, Asn1523.35, Trp3207.35, His3217.36, and signaling. Therefore, this crystal structure was used as an Tyr3287.43. During the last two decades, the use of computa- independent positive control, which, in addition to the use of tional techniques has revealed new details in fentanyl binding, co-crystallized BU-72-bound activated hMOR and naltrexone- identifying different residues that might potentially inter- bound inactive hMOR, allows us to quantify different act.25,26,91,96 Despite fentanyl being more prone to change in conformational changes induced or selected in the hMOR by its initial conformation during our MD simulations compared fentanyl or morphine. The question arises which is the proper to morphine, the same stable binding mode of fentanyl is reference state for the system. Both the apo receptor87 and an reached in the first microsecond of all replicas, especially in the antagonist-bound receptor may, in principle, be suitable for inactive receptor state. This is interesting because a recent this purpose. We chose the latter condition because of the high report investigating fentanyl binding in MD simulations of the stability of an antagonist-bound inactive receptor state, as well mMOR failed to obtain a stable binding pose in the inactive as the potential to deactivate an active receptor state to the state unless sodium ions (a known negative allosteric inactive. However, utilizing the apo receptor to investigate the modulator) were co-bound.26 Such ions are not required in apparently low reported basal activity of the hMOR85 is an our study. As fentanyl and morphine are both agonists, it might interesting area for future study, in particular regarding be expected that they should favor binding of the activated potential allosteric modulation by anionic phospholipids as state over the inactive; however, tight induced-fit of the has been reported for other homologous class A GPCRs.77,78,88 activated crystal structure for its co-crystallized morphine-like Under the simple framework of the two-state model of agonist BU-725 appears to enable precise docking of morphine L but partially hinders fentanyl. Conversely, the larger size of co- receptor activation: R ↔ R *, with R and R*, the inactive and crystallized antagonist β-FNA in the inactive crystal29 may active receptor states and L, the interconversion equilibrium favor a faster stabilization of fentanyl in this state compared to constant L = [R*]/[R], the agonist intrinsic efficacy α can be the activated one. Nevertheless, once stable ligand binding is αL seen either through the induction branch AR ↔ AR* or reached within our MD simulations starting from activated or K K/ α inactive states, morphine and fentanyl make different ligand− through the selection branches: A + R ↔ AR and A + R* ←→ receptor interactions. In particular, morphine interacts more AR* of the thermodynamic cycle, where K is the dissociation frequently with the N-terminus, TM5, TM6, and TM7, coming equilibrium constant. Thus, induction and selection ap- into contact with residues: Asp56N‑ter, Lys2355.39, Ile2986.51, proaches are equivalent in terms of intrinsic efficacy α within Val3026.55, Trp3207.35, and Ile3247.39, while fentanyl interacts the context of four receptor species in equilibrium. However, with ECL2, TM2, and TM7, coming into contact with the situation can be more complex when receptor states Gln1262.60, Asn1292.63, Val1453.28, Cys219ECL2, and Tyr3287.43. include ensembles of protein conformations. Conversion of AR However, there are similarities as both morphine and fentanyl into AR* may involve different intermediates with different interact with Ser57N‑term, Tyr1503.33, and Asp1493.32. kinetics. The kinetic component can be a limiting factor in The different molecular scaffolds of morphine and fentanyl some cases making more complex the correspondence between lead to different binding modes, which are associated with induction and selection approaches. Considering that the different conformational arrangements in the orthosteric inactive receptor R is, normally, the major species in the pocket and which may be related to their different absence of an agonist, it is expected that an agonist making pharmacological efficacies. The widely described5,7,17,29,30,45 “easier” the conversion of AR into AR* will reflect this electrostatic interaction with Asp1493.32 has been observed as molecular feature onto its pharmacological profile. Our essential for binding of μ-opioid ligands, and the Asp1493.32 simulations suggest that this could be the case in the conformation resulting from it plays a direct role in receptor comparison between morphine and fentanyl. Therefore, in (in)activation as shown by mutagenesis experiments.26,48,49 the present study, receptor activation by either induction or Likewise, H-bonding between Asp1493.32 and Tyr3287.43 in the selection of receptor active states has been considered. This orthosteric pocket has been shown to be relevant for MOR has been possible because both inactive and active MOR activation.26 Further experimental studies26,46,95 have shown structures were available. Therefore, when starting from a that Tyr3287.43 has a significant effect on ligand potency and specific simulation initial state (activated or inactive hMOR) if agonist-induced receptor activation, specifically for fentanyl, the receptor achieves the opposite state (inactive- or active-like and is conserved and functional in δ- and κ-opioid receptors as conformations, respectively) the induced-fit approach is demonstrated by mutagenesis.98 Other mutagenesis studies followed, otherwise, if the initial state is preserved, we can have identified Tyr1503.33 to be relevant for agonist binding state that a selection approach occurs. This makes the dynamic affinity,95,97 including fentanyl. In addition, mutation of the structural analysis more robust and reliable. It should be taken conserved residue Tyr1393.33 in the κ-opioid receptor alters into account that because we are using receptor states obtained ligand potency.98 Likewise, a recent mutagenesis study99 has for particular ligands, there is a bias toward these states for shown that Cys219ECL2 likely comes into contact with fentanyl those ligands structurally resembling the crystallized ones and is relevant in its binding but does not affect receptor (BU72 and β-FNA). In this regard, the morphine binding pose activation. This is an interesting result because it implies a is widely known42,89,90 because of its similarity with the molecular separation of the affinity and efficacy concepts, 1266 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article which is difficult if not impossible to obtain by parameter Asp1493.32−Tyr3287.43 H-bonding through a three-way side- estimation from functional studies as the estimated binding chain interaction thus favoring receptor activation.44,47 In constants in operational models include the receptor activation comparison, the Gln1262.60 cis sidechain conformation, component.100 Therefore, our morphine binding pose is in predominantly selected by morphine, shows lower accordance with previously described binding modes,42,89,90 Asp1493.32−Tyr3287.43 H-bonding and less favorable receptor co-crystallized BU-725 and β-FNA,29 and side-directed muta- activation.44,47 Differences in Gln1262.60 conformation with genesis, which determined Val3026.55 or Trp3207.35 to be fentanyl and morphine and their influence on Asp1493.32− crucial for the decrease101 or increase29 of morphine affinity, Tyr3287.43 H-bonding may partly explain their difference in respectively. Despite this, to our knowledge, no mutagenesis efficacy.51−53,56,57 The sidechain conformation of Asp1493.32, studies have been performed for residues: Asp56 N‑ter, which is dictated by the different binding poses of morphine Lys2355.39, and Ile3247.39, even though they have been and fentanyl, is further dictated by different interactions the implicated in morphine binding before.42,90 On the other ligands make with Lys2355.39. As both ligands contain an hand, the fentanyl binding pose presented in our study, which oxygen acceptor group(s), the sidechain amino group of presents protein−ligand interactions similar to those shown Lys2355.39 can theoretically interact with either. In relevant previously,25,91 is in agreement with the first point mutation crystals, Lys2355.39 covalently binds antagonist β-FNA,29 does experiments45,46,97 involving residues Asp1493.32, Tyr1503.33, not interact with agonist BU-72,5 and makes an H-bond with and Tyr3287.43 but not Asn1523.35, Trp3207.35, and His3217.36. agonist DAMGO.30 Differences in the molecular scaffold As far as we are aware, no experimental mutagenesis has been between ligands lead to higher Lys2355.39−ligand H-bond performed on Gln1262.60 and Asn1292.63, which we predict to occupancy with bound morphine, which results in different interact with fentanyl. The relative stability of our observed conformations of Lys2355.39. Accordingly, Lys2355.39−ligand binding poses of morphine and fentanyl, irrespective of H-bonding is related to enhanced stabilization of morphine receptor conformation, leads us to conclude that they are relative to fentanyl in the orthosteric pocket, thus leading sufficiently accurate. As an extra validation of the fentanyl Asp1493.32 to adopt a gauche− sidechain conformation when binding pose, we docked fentanyl to the new DAMGO-bound morphine is bound. Therefore, our results indicate that fully active crystal structure of the mMOR,30 which might be Asp1493.32−Tyr3287.43 H-bonding, previously described to be expected to yield greater accuracy. This results in a similar essential for MOR activation,44,47 is determined by ligand docking pose to that obtained with our hMOR models and interactions with Asp1493.32 and Gln1262.60 and by its position MD simulations (Figure S17), thus adding an extra layer of in the orthosteric pocket influenced by its interaction with confidence. Lys2355.39. Conformational dispositions in the orthosteric pocket lead Oxygen H-bond acceptor groups on morphine and fentanyl to alteration of other receptor residues such as Met1533.36 and not only establish interactions with TM5 but also make H- Trp2956.48, which assist receptor activation.26,102,103 More bonds with Tyr1503.33 on TM3, either directly or through a specifically, different conformations of Met1533.36 have been water molecule. Similar to Lys2355.39, morphine has a higher associated with a specific ligand-dependent microswitch for Tyr1503.33−ligand H-bond occupancy than fentanyl, which MOR β-arresting signaling.25 Although in our results we have may enhance TM3 activation and exert a stabilizing effect on not found multiple active-like states, we observe different the receptor orthosteric pocket. This includes the crystal conformations of Met1533.36 more frequently induced when gauche sidechain conformation of Met1533.36, which is largely either fentanyl or morphine is bound into the receptor induced in simulations starting from the activated state and orthosteric pocket. In addition, Trp2956.48 has previously been when morphine is bound. This conformation stabilizes observed to rotate during activation of the mMOR26,103 and hydrophobic interactions that occur at the bottom of the class A GPCRs, in general,74 and has been commonly named a orthosteric pocket restricting the transition from one state to “toggle switch”. Such “toggling” can result in different H- the other. In this respect, signal transmission along TM3 from bonding between TM5, TM6, and TM7 and solvation of the Tyr1503.33 to Met1533.36 of the hMOR is similar to that of receptor core.26,103 In particular, it has been described that other GPCRs such as the nociceptin receptor.102 On the fentanyl stabilizes different rotameric states of Trp2956.48 contrary, selection of the Met1533.36 trans conformation by compared to morphine or the apo receptor.26 In our study, fentanyl when the receptor is in an inactive-like state suggests a we identify ligand-specific orthosteric pocket changes in (i) destabilization of the hydrophobic interactions, which supports ligand interactions with Gln1262.60, Asp1493.32, Tyr1503.33, and transition from an inactive receptor state to a more active one Lys2355.39 including concomitant sidechain rotameric changes, by disconnecting TM3 from TM6. These results indicate that (ii) interhelical Asp1493.32−Tyr3287.43 H-bonding, (iii) side- fentanyl has a facility to change Met1533.36 from an inactive to chain rotameric changes in Trp295 6.48 and Met153 3.36 active-like conformation more frequently than morphine. In (influenced by conformational change of Tyr1503.33), and addition to Asp1493.32−Tyr3287.43 and Tyr1503.33−ligand H- (iv) proximity and interhelical H-bonding between bonding and associated conformational change of Met1533.36, Trp2956.48−Ala2425.46 (influenced by conformational change the conformational toggling of Trp2956.48 has previously been in Trp2956.48). In terms of how these conformational changes described in the activation of the MOR26 and δ-opioid are connected, we observe that morphine and fentanyl engage receptor.104 The Trp2956.48 sidechain also adopts different with Asp1493.32 via different gauche− or trans sidechain conformations in inactive,29 activated,5 and fully active crystals conformations, respectively. As trans Asp1493.32 is within the of the mMOR.30 In our MD simulations with either bound H-bonding distance with Tyr3287.43, as expected, we find morphine or fentanyl, Trp2956.48 adopts an inactive-like gauche higher Asp1493.32−Tyr3287.43 H-bond occupancy with bound or an activated-like trans conformation depending on the fentanyl. Likewise, a trans sidechain conformation of receptor state and ligand bound. The trans position resembles Gln1262.60, which is induced predominantly when fentanyl is the conformation observed in the fully active crystal bound, interacts with Tyr328 7.43 , further stabilizing structure,30 adopting a perpendicular orientation to TM6 but 1267 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article with the indole NH pointing toward TM5 instead of TM7. and receptor starting state. This communication has been This conformation is in general agreement with a recent study proposed to operate, at least in part, through a polar network implementing MD simulations of the mMOR with docked mediated by water molecules, which changes depending on the fentanyl but differs from those observed with docked receptor state.5,79−81 Here, we have identified two water morphine.26 The reason for this discrepancy over morphine clusters, which differ between active- and inactive-like states or is unclear but could be due to differences between mouse and between morphine/fentanyl bound systems: (i) intracellular human receptor models or co-bound sodium ions, which might solvation near the NPxxY region on TM7 and (ii) solvation in negate agonist activity.26 In our simulations starting from the the receptor core between Asn3307.45 and Trp2956.48. In inactive hMOR, the trans conformation of Trp2956.48 is only particular, the former has been described to occur during obtained twice, once with bound morphine and the other more activation of the hMOR and class A GPCRs in general.5,53,81 transiently with fentanyl. This shows that activation of the This is partly because TM6 conformational change involved in hMOR, in particular TM6 conformational change, is a the transition from the inactive to near-active receptor state, as relatively difficult process to capture in a microsecond period. observed more with bound fentanyl, leads to greater solvation This contrasts with the relative ease of receptor activation in intracellular receptor regions, including the NPxxY motif, as observed in MD simulations of other class A GPCRs, such as well as the receptor core near Trp2956.48. On the contrary, the CB177 or adenosine A2A.78 In some respects, this fits with the smaller size of morphine allows greater solvation of the low constitutive activity displayed by the hMOR in orthosteric pocket, which may negatively affect its binding pharmacological experiments,85 which suggests a highly stable stability in the inactive receptor state. Furthermore, the inactive state. When Trp2956.48 adopts the trans conformation, energetic state of these waters and how fentanyl potentially it allows interhelical H-bonding with Ala2425.46, a position displaces more of them in the orthosteric pocket relative to which has been shown to be relevant for TM5 flexibility in morphine can be the topic of future study. class A GPCRs.81 Conversely, the gauche conformation of In summary, our dynamics of the hMOR are determined by Trp2956.48 H-bonds Asn3307.45 when whole TM6 movement is ligand binding to activated5 and inactive29 receptor models, not present and the distance between Trp2956.48 and Ala2425.46 based on mMOR crystal structures, including extracellular/ is large. As a result, conformational change and activation of intracellular loops and modeled N-termini. The recently TM5 and TM6 are highly related to close Trp2956.48− published fully active mMOR cryo-EM structure30 allows for Ala2425.46 interaction, either by H-bond formation or helical an independent positive control of MD-generated receptor rearrangement. Together with Asp1493.32−Tyr3287.43 and structures. In addition, MD simulations with the antagonist Tyr1503.33−ligand H-bonding, as well as Met1533.36 conforma- naltrexone bound represent a negative control that allows us to tional change in the receptor orthosteric pocket, these are the better describe differences between fentanyl and morphine. factors that most strongly govern intracellular receptor Likewise, the constitution of the geometry of our hMOR conformation and where differences between morphine and models is validated by the mutual ligand and receptor fentanyl are most clear. conformational stability in control systems with BU-72 or In our MD simulations, conformational changes of TM naltrexone bound, respectively, to activated or inactive states. helices differently occur with bound morphine or fentanyl. Two different physically stable binding poses of morphine and While morphine largely induces or selects conformational fentanyl in the hMOR mediate different sets of protein−ligand activation changes of TM3 and TM5, fentanyl more frequently interactions (in particular involving Gln1262.60, Lys2355.39, achieved active-like conformations of TM6 and TM7. Asp1493.32, Lys2355.39, and Tyr1503.33), which differentially Conformational rearrangements of these TM helices are change orthosteric pocket conformation (in particular reflected in changes in specific intracellular regions such as Asp1493.32−Tyr3287.43 H-bonding, Met1533.36 orientation, Arg1673.50−Thr2816.34 interaction between TM3 and TM6 and Trp2956.48−Ala2425.46 interaction). These result in and the NPxxY motif on TM7. However, in the most part, we different intracellular TM3, TM5, TM6, and TM7 conforma- observe that the distance between Arg1673.50 and Thr2816.34 tional changes, including movement of the NPxxY motif and and conformation of the NPxxY motif remain constant when Arg1673.50−Thr2816.34 separation. Fentanyl has a stronger the receptor is already in the activated state with bound effect on TM6 and TM7 conformation, while morphine fentanyl or morphine, which means that both agonists preferentially affects TM3 and TM5. As conformational change generally sustain the active receptor state, which is consistent in TM6 is critical for GPCR activation and G protein binding with their agonist character. However, activation of the in particular, this likely explains the enhanced receptor activity receptor from the inactive state through breakage of elicited by fentanyl in vivo, which has 50−100-fold higher Arg1673.50−Thr2816.34 interaction and conformational change potency than morphine.51−53,56,57 The structural insights of the NPxxY motif is more frequently observed with bound gained in the present study can be used in future work with fentanyl than with morphine. This suggests that fentanyl more a spread of μ-opioid targeting drugs, to identify those strongly induces conformational changes in the NPxxY motif interactions and conformational changes associated with through enhanced Asp1493.32−Tyr3287.43 H-bonding and particular ligand efficacy and the disjunction between signaling enables TM3−TM6 separation through disruption of networks through the receptor, which probably determine the Met1533.36 conformation and Trp2956.48−Ala2425.46 transient proficiency in modulation of the receptor intracellular pocket and Gi protein or β-arrestin recognition. ■ H-bonding and further close contact. This may be indicative of its higher efficacy observed at the experimental level51−53 and also reflects the allosteric communication that operates METHODS between the orthosteric pocket and intracellular regions of Homology Modeling. The crystal structures of activated the hMOR. Nearly all orthosteric metrics analyzed in this study and inactive states of the mMOR (PDB ids: 5C1M5 and have shown statistically significant association with the state of 4DKL29) were used to model respective hMOR states using intracellular metrics studied, independent of the ligand bound Chimera v1.11105 and Modeller v9.16.106 Non-native camelid 1268 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article G protein mimetic nanobody Nb39 and T4L lysozyme fusion MD Simulations. MD simulations of each hMOR system protein were removed from respective activated5 and inactive29 starting from either activated or inactive states with bound mMOR crystal structures, as well as co-crystallized opioids: morphine, fentanyl, or naltrexone were performed in triplicate BU-725 and β-FNA,29 respectively. In order to homology using the CHARMM36 force-field114 with ACEMD116 on model human wt-sequence (Uniprot ID: P35372)107 hMOR specialized GPU-computer hardware. Each replica was energy- models, respective activated5 and inactive29 mMOR crystal minimized for 2000 steps followed by 28 ns of equilibration at structures were used as templates and with renumbered 300 K and 1 atm, with positional harmonic restraints on sequences, so Gly52 was renumbered to Gly54, and so. Using protein and ligand heavy atoms progressively released over the Chimera105 and the Dunbrack 2010 rotamer library,108 first 8 ns of equilibration and thereafter continued without sidechain modifications were performed on S55R, H56D, restraints. This was followed with an unbiased continuous Q61P, V68I, N139T, V189I, and I308V, which correspond to production run under the same conditions for 3 μs for each nonconserved residues between species. The most probable replica. In addition, the control hMOR system with bound BU- rotamer with the fewest steric clashes was chosen in each 72 was run in duplicate for 2 μs after the same energy- instance. Using Modeller,106 noncrystallized segments of the minimization and equilibration process as systems of interest, inactive hMOR N-terminus (residues Gly54 to Ile73), inactive making a total additive simulation time of 58 μs. hMOR intracellular loop 3 (ICL3; residues Val264 to Lys271), MD Simulation Analysis. In order to evaluate activated and activated hMOR C-terminus (residues Arg347 to Ile354) and inactive hMOR model accuracy, we respectively analyzed were homology modeled using the other receptor structure as a the conformational stability of the receptor helical domain and template. Subsequently, activated and inactive hMOR models co-crystallized ligand BU72 or naltrexone (substituted for co- were energy-minimized in the AMBER14SB force-field109 crystallized β-FNA) making use of RMSD measurements with using Chimera in vacuum conditions105 to optimize internal VMD software v1.9.2.117 Analysis of the three microseconds of interactions. each trajectory with bound morphine, fentanyl, or naltrexone Ligand Docking. The structures of three ligands of interest (including control replicas of the naltrexone-bound inactive for the present study, morphine and fentanyl, and negative state) was performed using VMD117 with the following control naltrexone, were extracted from Pubchem110 and protocol: (i) ligand stability in inactive and activated hMOR docked in activated and inactive hMOR models using models and (ii) analysis of respective stable binding poses; (iii) Autodock4.2.6.111 Specifically, morphine, fentanyl, and nal- comparison between ligand-mediated hMOR orthosteric trexone were prepared by respectively selecting their bioactive pocket conformational changes and their correlation with (iv) TM helical movement; (v) intracellular receptor isomer, protonating their amino group (consistent with conformational state and (vi) conformational changes of the physiological pH), and assigning a net charge of +1 in order rest of the receptor; (vii) water-network differences associated to correctly interact with Asp1493.32. This particular protein− with different receptor states. Evaluation of conformational ligand interaction is in accordance with co-crystallized μ-opioid stability was measured by the analysis of the two last agonists (BU-72; DAMGO),5,30 antagonist (β-FNA),29 and microseconds of each replica of the systems generated, related structural data.91 In addition to the naltrexone-bound whereas residue-specific conformational changes were studied inactive state of the hMOR, for control purposes, co- considering the entire simulation time-period. In detail, RMSD crystallized BU-72 was re-docked into the activated hMOR measurements of fentanyl or morphine in their respective MD state with its amino group protonated as other ligands. simulations were used to monitor ligand stability. In addition, Docking grid points were extended to cover total orthosteric all fentanyl, morphine, and naltrexone conformations achieved pocket volumes, respectively. The selected docked conforma- during respective trajectories were classified into five different tion of each ligand in each receptor state represents the top hit clusters using a cutoff of 2.0 Å using VMD Clustering Tool.117 identified by the best predicted affinity, which also makes A list of residues that came into contact with either morphine correct interaction with Asp1493.32. Subsequent docked or fentanyl (<3.5 Å) in each replica of each system and their ligand−receptor complexes were energy-minimized in the frequency were extracted making use of a smoothing/averaging AMBER14SB force-field109 using Chimera105 in vacuum TCL script118 executed in VMD.117 Electrostatic, van der conditions to optimize protein−ligand interactions. Waals, and hydrophobic interactions between the ligand and MD System Setup. Two control systems, BU-72 bound to residues with frequencies ≥20.0% in two or more replicas were the activated hMOR and naltrexone bound to the inactive analyzed for energetics using NAMD Energy plugin119 v1.4 hMOR, and five different systems of interest were generated: within VMD.117 RMSD analysis of the hMOR TM domain was morphine-bound activated and inactive hMOR, fentanyl- performed with respect to respective receptor starting bound activated and inactive hMOR, and naltrexone-bound conformation or with the fully active crystal30 to assess activated hMOR; each embedded separately into a POPC receptor stability or (in)activation. A threshold of ± 3.5 Å from membrane using the CHARMM-GUI web-based interface112 the initial receptor state as well as activated and inactive crystal and solvated with TIP3P water molecules. All models were conformational divergence from fully active crystal30 was used oriented in the membrane according to the OPM database113 to determine significant receptor conformational changes. entry of the activated mMOR (id: 5c1m).5 The disulfide Intracellular receptor conformational changes associated with bridge between Cys1423.25 and Cys219ECL2 was maintained in activation/deactivation were analyzed by making use of all systems, and charge neutralizing ions (0.15 M KCl) were VMD117 in terms of movement with respect to the protein automatically introduced in order to make a net system charge center of (i) TM3, (ii) the intracellular tip of TM5, (iii) the of 0. CHARMM-GUI112 automatically generated membrane, intracellular tip of TM6, and (iv) the NPxxY motif on TM7, in water, and protein parameters according to the CHARMM36 addition to (v) distance measurement between alpha carbons force-field114 and ligand parameters in accordance with of residues Arg1673.50 and Thr2816.34. TM3 upward axial CHARMM General Force Field (CGenFF) v.1.0.0.115 movement was measured by z axis offset difference between 1269 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article TM3 and the protein center of mass (COM). Helix tips were initial state of the receptor or the ligand bound, thus making a defined according to the second and third-last intracellular total sample of 900 receptor−ligand conformations. This helical turns (TM5: residues 251−261, TM6: 274−284) and sample was considered sufficient to construct two-way distances measured relative to the protein COM. Taking in contingency tables to evaluate the Chi2 score and to perform consideration differences between crystal structures5,29,30 in Student’s t-test analysis for exploring general mechanistic their positions of TM3, TM5, and TM6 as well as MD trends. P < 0.05 was considered statistically significant. IBM simulations with the antagonist naltrexone bound, relevant SPSS Statistics 20.0.0 was used for statistical analyses. internal distances measured were normalized with respect to the inactive crystal,29 with respective distance thresholds of +0.3, −0.8, and + 1.7 Å used to classify active or inactive helix ■ * ASSOCIATED CONTENT sı Supporting Information conformations. Active- or inactive-like conformations of The Supporting Information is available free of charge at intracellular metrics such as the NPxxY motif and https://pubs.acs.org/doi/10.1021/acs.jcim.0c00890. Arg1673.50−Thr2816.34 interaction were defined by distance from the protein COM (using −0.8 Å threshold) or by +10.0 (Figures S1−S17) Stability of BU-72 bound in the Å for Arg1673.50−Thr2816.34 inter-residue distance. Residues activated hMOR model during control MD simulations; showing differences in their P.E. of interaction between bound morphine, fentanyl, and naltrexone conformational fentanyl and morphine were analyzed according to (i) clustering and RMSD in MD simulations; proportion sidechain χ1 dihedral angle; (ii) χ2 dihedral angle; (iii) P.E. of residues in close contact with morphine, fentanyl, and of interaction with nearby residues; (iv) distribution with naltrexone in MD simulations; entire TMD receptor respect to conformational changes of other residues; and (v) conformational change; TM helix movement with correlation with intracellular receptor conformational changes. respect to the protein center in MD simulations; Sidechain dihedral angles were classified according to a intracellular Arg1673.50−Thr2816.34 distance in MD threshold between predominantly observed χ1 and χ2 dihedral simulations; average potential energy per replica of angle conformations. Thus, a ± 115.0° value was used as a protein−ligand and sidechain−sidechain interactions in threshold for classifying Asp1493.32 χ1 dihedral angle MD simulations; fentanyl docking in active MOR crystal conformation (gauche and trans conformations), ± 90.0° structures; statistical analyses; distribution of orthosteric threshold for the Met1533.36 χ1 dihedral angle (gauche or pocket conformational changes; and (Table S1) average alternative conformations), and ± 240.0° threshold for and SD of metrics mentioned in the text (PDF) Lys2355.39 and Gln1262.60 χ2 dihedral angles (predominant cis and trans conformations of 300.0° and 180.0°, respectively). Electrostatic intersidechain or ligand−sidechain interactions ■ AUTHOR INFORMATION Corresponding Authors were specified according to ± 4.5 Å distance threshold. H- James A. R. Dalton − Laboratory of Molecular bonding between morphine, fentanyl, or naltrexone with Neuropharmacology and Bioinformatics, Unitat de Asp1493.32 was selected according to a threshold of ± 2.5 Å. Bioestadística and Institut de Neurociències, Universitat Distances/H-bonding was measured between (i) fentanyl, Autònoma de Barcelona, 08193 Bellaterra, Spain; Unitat de morphine, or naltrexone amino groups and Asp1493.32 gamma Neurociència Traslacional, Parc Taulí Hospital Universitari, carbon (considered as the center of the Asp1493.32 carboxylate Institut d’Investigació i Innovació Parc Taulí (I3PT), Institut group), (ii) fentanyl, morphine, or naltrexone oxygen accept- de Neurociències, Universitat Autònoma de Barcelona, 08193 or/donor atoms (a single carbonyl group in the former or one Bellaterra, Spain; Instituto de Salud Carlos III, Centro de ether and two hydroxyl groups in morphine or one ether, one Investigación Biomédica en Red de Salud Mental, hydroxyl group and one carbonyl group in naltrexone) and the CIBERSAM, 28029 Madrid, Spain; Email: James.Dalton@ Tyr1503.33 hydroxyl group or the Lys2355.39 amino group, (iii) uab.es Asp1493.32 gamma carbon and the Tyr3287.43 hydroxyl group, Jesús Giraldo − Laboratory of Molecular Neuropharmacology and (iv) the Trp2956.48 indole NH group and the Ala2425.46 and Bioinformatics, Unitat de Bioestadística and Institut de backbone carbonyl group. Distance between residues Neurociències, Universitat Autònoma de Barcelona, 08193 Trp2956.48 and Ala2425.46 was measured between their Bellaterra, Spain; Unitat de Neurociència Traslacional, Parc respective sidechain COMs. Water density maps were created Taulí Hospital Universitari, Institut d’Investigació i Innovació using VolMap within VMD117 by calculating the weighted Parc Taulí (I3PT), Institut de Neurociències, Universitat atomic density of water molecules averaged over the last two Autònoma de Barcelona, 08193 Bellaterra, Spain; Instituto microseconds of each MD replica. de Salud Carlos III, Centro de Investigación Biomédica en Statistical Analysis. Chi-square (Chi2) and Student’s t- Red de Salud Mental, CIBERSAM, 28029 Madrid, Spain; test statistical analyses were performed to assess the statistical orcid.org/0000-0001-7082-4695; Email: Jesus.Giraldo@ significance of the potential association between orthosteric uab.es pocket metrics (sidechain−sidechain Asp1493.32−Tyr3287.43, Tyr1503.33−ligand, Lys2355.39−ligand, and sidechain−back- Author bone Trp2956.48−Ala2425.46 H-bond presence or absence, Adrián Ricarte − Laboratory of Molecular Met1533.36 gauche conformation selection, or Trp2956.48− Neuropharmacology and Bioinformatics, Unitat de Ala2425.46 distance, respectively) with the intracellular Bioestadística and Institut de Neurociències, Universitat conformations achieved (active- or inactive-like conformation Autònoma de Barcelona, 08193 Bellaterra, Spain; Unitat de of TM3, TM5, TM6, NPxxY, and Arg1673.50−Thr2816.34 Neurociència Traslacional, Parc Taulí Hospital Universitari, distance). To perform these analyses, we used 50 representa- Institut d’Investigació i Innovació Parc Taulí (I3PT), Institut tive snapshots of each replica, extracted every 40 ns of the last de Neurociències, Universitat Autònoma de Barcelona, 08193 two microseconds of respective trajectories, independent of the Bellaterra, Spain; Instituto de Salud Carlos III, Centro de 1270 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Investigación Biomédica en Red de Salud Mental, (12) OLINVYK (oliceridine) injection, for intravenous use, CII. CIBERSAM, 28029 Madrid, Spain Olinvyk web site. https://olinvyk.com/docs/OLINVYK_Final_ Label_Ver%20002_Nov2020_1103.pdf (accessed Dec 3, 2020). Complete contact information is available at: (13) Lambert, D.; Calo, G. Approval of oliceridine (TRV130) for https://pubs.acs.org/10.1021/acs.jcim.0c00890 intravenous use in moderate to severe pain in adults. Br J Anaesth. 2020, 125, e473−e474. Funding (14) Gurevich, V. V.; Gurevich, E. V. Biased GPCR signaling: This project has received funding from the European Union’s Possible mechanisms and inherent limitations. Pharmacol. Ther. 2020, Horizon2020 research and innovation programme under grant 211, 107540. agreement No 848068 and from Ministerio de Ciencia, (15) Wingler, L. M.; Skiba, M. A.; McMahon, C.; Staus, D. P.; Kleinhenz, A. L. W.; Suomivuori, C. M.; Latorraca, N. R.; Dror, R. O.; Innovación y Universidades (Spain) under grant agreement Lefkowitz, R. J.; Kruse, A. C. Angiotensin and biased analogs induce SAF2017-87199-R. This publication reflects only the authors’ structurally distinct active conformations within a GPCR. Science view, and the European Commission is not responsible for any 2020, 367, 888−892. use that may be made of the information it contains. (16) Kenakin, T.; Christopoulos, A. Signalling bias in new drug Notes discovery: detection, quantification and therapeutic impact. Nat. Rev. The authors declare no competing financial interest. Drug Discovery 2013, 12, 205−216. ■ (17) Hall, D. A.; Giraldo, J. A method for the quantification of biased signalling at constitutively active receptors. Br. J. Pharmacol. REFERENCES 2018, 175, 2046−2062. (1) Steglitz, J.; Buscemi, J.; Ferguson, M. J. The future of pain (18) Gillis, A.; Sreenivasan, V.; Christie, M. J. Intrinsic efficacy of research, education, and treatment: a summary of the IOM report opioid ligands and its importance for apparent bias, operational “Relieving pain in America: a blueprint for transforming prevention, analysis and therapeutic window. Mol. Pharmacol. 2020, 119, care, education, and research”. Transl. Behav. Med. 2012, 2, 6−8. No. 119214. (2) Breivik, H.; Collett, B.; Ventafridda, V.; Cohen, R.; Gallacher, D. (19) Manglik, A.; Lin, H.; Aryal, D. K.; McCorvy, J. D.; Dengler, D.; Survey of chronic pain in Europe: prevalence, impact on daily life, and Corder, G.; Levit, A.; Kling, R. C.; Bernat, V.; Hübner, H.; Huang, X. treatment. Eur. J. Pain. 2006, 10, 287−333. P.; Sassano, M. F.; Giguère, P. M.; Löber, S.; Da, D.; Scherrer, G.; (3) O’Brien, T.; Christrup, L. L.; Drewes, A. M.; Fallon, M. T.; Kobilka, B. K.; Gmeiner, P.; Roth, B. L.; Shoichet, B. K. Structure- Kress, H. G.; McQuay, H. J.; Mikus, G.; Morlion, B. J.; Perez- based discovery of opioid analgesics with reduced side effects. Nature Cajaraville, J.; Pogatzki-Zahn, E.; Varrassi, G.; Wells, J. C. European 2016, 537, 185−190. Pain Federation position paper on appropriate opioid use in chronic (20) Azzam, A. A. H.; McDonald, J.; Lambert, D. G. Hot topics in pain management. Eur. J. Pain. 2017, 21, 3−19. opioid pharmacology: mixed and biased opioids. Br. J. Anaesth. 2019, (4) Wang, J. B.; Johnson, P. S.; Persico, A. M.; Hawkins, A. L.; 122, e136−e145. Griffin, C. A.; Uhl, G. R. Human mu opiate receptor. cDNA and (21) Piekielna-Ciesielska, J.; Ferrari, F.; Calo, G.; Janecka, A. genomic clones, pharmacologic characterization and chromosomal Cyclopeptide Dmt-[D-Lys-p-CF3-Phe-Phe-Asp]NH2, a novel G assignment. FEBS Lett. 1994, 338, 217−222. protein-biased agonist of the mu opioid receptor. Peptides 2018, (5) Huang, W.; Manglik, A.; Venkatakrishnan, A. J.; Laeremans, T.; 101, 227−233. Feinberg, E. N.; Sanborn, A. L.; Kato, H. E.; Livingston, K. E.; (22) Pasquinucci, L.; Turnaturi, R.; Calò, G.; Pappalardo, F.; Ferrari, Thorsen, T. S.; Kling, R. C.; Granier, S.; Gmeiner, P.; Husbands, S. F.; Russo, G.; Arena, E.; Montenegro, L.; Chiechio, S.; Prezzavento, M.; Traynor, J. R.; Weis, W. I.; Steyaert, J.; Dror, R. O.; Kobilka, B. K. O.; Parenti, C. (2S)-N-2-methoxy-2-phenylethyl-6,7-benzomorphan Structural insights into μ-opioid receptor activation. Nature 2015, compound (2S-LP2): Discovery of a biased mu/delta opioid receptor 524, 315−321. agonist. Eur. J. Med. Chem. 2019, 168, 189−198. (6) Roeckel, L. A.; Utard, V.; Reiss, D.; Mouheiche, J.; Maurin, H.; (23) Schmid, C. L.; Kennedy, N. M.; Ross, N. C.; Lovell, K. M.; Yue, Robé, A.; Audouard, E.; Wood, J. N.; Goumon, Y.; Simonin, F.; Z.; Morgenweck, J.; Cameron, M. D.; Bannister, T. D.; Bohn, L. M. Gaveriaux-Ruff, C. Morphine-induced hyperalgesia involves mu Bias factor and therapeutic window correlate to predict safer opioid opioid receptors and the metabolite morphine-3-glucuronide. Sci. analgesics. Cell 2017, 171, 1165−1175.e13. Rep. 2017, 7, 10406. (24) Grim, T. W.; Acevedo-Canabal, A.; Bohn, L. M. Toward (7) Kapoor, A.; Martinez-Rosell, G.; Provasi, D.; de Fabritiis, G.; directing opioid receptor signaling to refine opioid therapeutics. Biol. Filizola, M. Dynamic and kinetic elements of μ-opioid receptor Psychiatry. 2020, 87, 15−21. functional selectivity. Sci. Rep. 2017, 7, 11255. (25) de Waal, P. W.; Shi, J.; You, E.; Wang, X.; Melcher, K.; Jiang, (8) Raehal, K. M.; Walker, J. K.; Bohn, L. M. Morphine side effects Y.; Xu, H. E.; Dickson, B. M. Molecular mechanisms of fentanyl in β-arrestin 2 knockout mice. J. Pharmacol. Exp. Ther. 2005, 314, mediated β-arrestin biased signaling. PLoS Comput. Biol. 2020, 16, 1195−1201. No. e1007394. (9) DeWire, S. M.; Yamashita, D. S.; Rominger, D. H.; Liu, G.; (26) Lipiński, P. F. J.; Jarończyk, M.; Dobrowolski, J. C.; Sadlej, J. Cowan, C. L.; Graczyk, T. M.; Chen, X. T.; Pitis, P. M.; Gotchev, D.; Molecular dynamics of fentanyl bound to μ-opioid receptor. J. Mol. Yuan, C.; Koblish, M.; Lark, M. W.; Violin, J. D. A G protein-biased Model. 2019, 25, 144. ligand at the μ-opioid receptor is potently analgesic with reduced (27) Morgan, M. M.; Christie, M. J. Analysis of opioid efficacy, gastrointestinal and respiratory dysfunction compared with morphine. tolerance, addiction and dependence from cell culture to human. Br. J. J. Pharmacol. Exp. Ther. 2013, 344, 708−717. Pharmacol. 2011, 164, 1322−1334. (10) Schneider, S.; Provasi, D.; Filizola, M. How oliceridine (TRV- (28) Kelly, E. Efficacy and ligand bias at the μ-opioid receptor. Br. J. 130) binds and stabilizes a μ-opioid receptor conformational state Pharmacol. 2013, 169, 1430−1446. that selectively triggers G protein signaling pathways. Biochemistry (29) Manglik, A.; Kruse, A. C.; Kobilka, T. S.; Thian, F. S.; 2016, 55, 6456−6466. Mathiesen, J. M.; Sunahara, R. K.; Pardo, L.; Weis, W. I.; Kobilka, B. (11) FDA Approves New Opioid for Intravenous Use in Hospitals, K.; Granier, S. Crystal structure of the μ-opioid receptor bound to a Other Controlled Clinical Settings. FDA news release web site. morphinan antagonist. Nature 2012, 485, 321−326. https://www.fda.gov/news-events/press-announcements/fda- (30) Koehl, A.; Hu, H.; Maeda, S.; Zhang, Y.; Qu, Q.; Paggi, J. M.; approves-new-opioid-intravenous-use-hospitals-other-controlled- Latorraca, N. R.; Hilger, D.; Dawson, R.; Matile, H.; Schertler, G. F. clinical-settings (accessed Dec 3, 2020). X.; Granier, S.; Weis, W. I.; Dror, R. O.; Manglik, A.; Skiniotis, G.; 1271 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Kobilka, B. K. Structure of the μ-opioid receptor-Gi protein complex. multiple transmembrane domains. Biochemistry 2008, 47, 10576− Nature 2018, 558, 547−552. 10586. (31) Ballesteros, J. A.; Weinstein, H. [19] Integrated methods for the (48) Li, J. G.; Chen, C.; Yin, J.; Rice, K.; Zhang, Y.; Matecka, D.; de construction of three-dimensional models and computational probing Riel, J. K.; DesJarlais, R. L.; Liu-Chen, L. Y. Asp147 in the third of structure-function relations in G protein-coupled receptors. Rec. transmembrane helix of the rat μ opioid receptor forms ion-pairing Molar. Bio. 1995, 25, 366−428. with morphine and naltrexone. Life Sci. 1999, 65, 175−185. (32) Filizola, M.; Devi, L. A. Structural biology: How opioid drugs (49) Befort, K.; Tabbara, L.; Bausch, S.; Chavkin, C.; Evans, C.; bind to receptors. Nature 2012, 485, 314−317. Kieffer, B. The conserved aspartate residue in the third putative (33) Carroll, F. I.; Dolle, R. E. The discovery and development of transmembrane domain of the delta-opioid receptor is not the anionic the N-substituted trans-3,4-dimethyl-4-(3′-hydroxyphenyl)piperidine counterpart for cationic opiate binding but is a constituent of the class of pure opioid receptor antagonists. ChemMedChem 2014, 9, receptor binding site. Mol. Pharmacol. 1996, 49, 216−223. 1638−1654. (50) James, A.; Williams, J. Basic Opioid Pharmacology - An Update. (34) Yuan, Y.; Zaidi, S. A.; Elbegdorj, O.; Aschenbach, L. C.; Li, G.; Br J Pain. 2020, 14, 115−121. Stevens, D. L.; Scoggins, K. L.; Dewey, W. L.; Selley, D. E.; Zhang, Y. (51) O’Donnell, J. K.; Halpin, J.; Mattson, C. L.; Goldberger, B. A.; Design, synthesis, and biological evaluation of 14-heteroaromatic- Gladden, R. M. Deaths involving fentanyl, fentanyl analogs, and, U- substituted naltrexone derivatives: pharmacological profile switch 47700 - 10 States, July-December 2016. MMWR Morb. Mortal. Wkly. from mu opioid receptor selectivity to mu/kappa opioid receptor dual Rep. 2017, 16, 1197−1202. selectivity. J. Med. Chem. 2013, 56, 9156−9169. (52) Mounteney, J.; Giraudon, I.; Denissov, G.; Griffiths, P. (35) Fujita, W.; Gomes, I.; Devi, L. A. Revolution in GPCR Fentanyls: Are we missing the signs? Highly potent and on the rise signalling: opioid receptor heteromers as novel therapeutic targets: in Europe. Int. J. Drug Policy. 2015, 26, 626−631. IUPHAR review 10. Br. J. Pharmacol. 2014, 171, 4155−4176. (53) Vardanyan, R. S.; Hruby, V. J. Fentanyl-related compounds and (36) Lee, C. W.; Ho, I. K. Pharmacological profiles of oligomerized derivatives: current status and future prospects for pharmaceutical μ-opioid receptors. Cell 2013, 2, 689−714. applications. Future Med. Chem. 2014, 6, 385−412. (37) Shang, Y.; LeRouzic, V.; Schneider, S.; Bisignano, P.; Pasternak, (54) Grissinger, M. Inappropriate prescribing of fentanyl patches is G. W.; Filizola, M. Mechanistic insights into the allosteric modulation still causing alarming safety problems. Pharm. Ther. 2010, 35, 653− of opioid receptors by sodium ions. Biochemistry 2014, 53, 5140− 654. 5149. (55) Schifano, F.; Chiappini, S.; Corkery, J. M.; Guirguis, A. (38) Livingston, K. E.; Traynor, J. R. Disruption of the Na+ ion Assessing the 2004-2018 fentanyl misusing issues reported to an binding site as a mechanism for positive allosteric modulation of the international range of adverse reporting systems. Front. Pharmacol. mu-opioid receptor. Proc. Nat. Acad. Sci. U. S. A. 2014, 111, 18369− 2019, 10, 46. 18374. (56) Lyons, P. J.; Rivosecchi, R. M.; Nery, J. P.; Kane-Gill, S. L. (39) Rasmussen, S. G.; Choi, H. J.; Fung, J. J.; Pardon, E.; Casarosa, Fentanyl-induced hyperalgesia in acute pain management. J. Pain P.; Chae, P. S.; Devree, B. T.; Rosenbaum, D. M.; Thian, F. S.; Palliat. Care Pharmacother. 2015, 29, 153−160. Kobilka, T. S.; Schnapp, A.; Konetzki, I.; Sunahara, R. K.; Gellman, S. (57) Warner, M.; Trinidad, J. P.; Bastian, B. A.; Minino, A. M.; H.; Pautsch, A.; Steyaert, J.; Weis, W. I.; Kobilka, B. K. Structure of a Hedegaard, H. Drugs most frequently involved in drug overdose nanobody-stabilized active state of the β(2) adrenoceptor. Nature deaths: United States, 2010-2014. Natl. Vital Stat. Rep. 2016, 65, 1− 2011, 469, 175−180. 15. (40) Ring, A. M.; Manglik, A.; Kruse, A. C.; Enos, M. D.; Weis, W. (58) Skolnick, P. The Opioid Epidemic: Crisis and Solutions. Annu. I.; Garcia, K. C.; Kobilka, B. K. Adrenaline-activated structure of β2- Rev. Pharmacol. Toxicol. 2018, 58, 143−159. adrenoceptor stabilized by an engineered nanobody. Nature 2013, (59) Verhamme, K. M. C.; Bohnen, A. M. Are We Facing an Opioid 502, 575−579. (41) Kruse, A. C.; Ring, A. M.; Manglik, A.; Hu, J.; Hu, K.; Eitel, K.; Crisis in Europe? Lancet Public Health. 2019, 4, e483−e484. Hübner, H.; Pardon, E.; Valant, C.; Sexton, P. M.; Christopoulos, A.; (60) Anselmi, L.; Jaramillo, I.; Palacios, M.; Huynh, J.; Sternini, C. Felder, C. C.; Gmeiner, P.; Steyaert, J.; Weis, W. I.; Garcia, K. C.; Ligand-induced μ opioid receptor internalization in enteric neurons Wess, J.; Kobilka, B. K. Activation and allosteric modulation of a following chronic treatment with the opiate fentanyl. J. Neurosci. Res. muscarinic acetylcholine receptor. Nature 2013, 504, 101−106. 2013, 91, 854−860. (42) Sounier, R.; Mas, C.; Steyaert, J.; Laeremans, T.; Manglik, A.; (61) Zheng, H.; Loh, H. H.; Law, P. Y. Beta-arrestin-dependent mu- Huang, W.; Kobilka, B. K.; Déméné, H.; Granier, S. Propagation of opioid receptor-activated extracellular signal-regulated kinases (ERKs) conformational changes during μ-opioid receptor activation. Nature translocate to nucleus in contrast to G protein-dependent ERK 2015, 524, 375−378. activation. Mol. Pharmacol. 2007, 73, 178−190. (43) Huang, P.; Visiers, I.; Weinstein, H.; Liu-Chen, L. Y. The local (62) Mori, T.; Kuzumaki, N.; Arima, T.; Narita, M.; Tateishi, R.; environment at the cytoplasmic end of TM6 of the μ opioid receptor Kondo, T.; Hamada, Y.; Kuwata, H.; Kawata, M.; Yamazaki, M.; differs from those of rhodopsin and monoamine receptors: Sugita, K.; Matsuzawa, A.; Baba, K.; Yamauchi, T.; Higashiyama, K.; introduction of an ionic lock between the cytoplasmic ends of helices Nonaka, M.; Miyano, K.; Uezono, Y.; Narita, M.. Usefulness for the 3 and 6 by a L6.30(275)E mutation inactivates the μ opioid receptor combination of G protein- and β-arrestin-biased ligands of μ-opioid and reduces the constitutive activity of its T6.34(279)K mutant. receptors: Prevention of antinociceptive tolerance. Mol. Pain. 2017. Biochemistry 2002, 41, 11972−11980. (63) Kovoor, A.; Celver, J. P.; Wu, A.; Chavkin, C. Agonist induced (44) Dalton, J. A.; Lans, I.; Giraldo, J. Quantifying conformational homologous desensitization of μ-opioid receptors mediated by G changes in GPCRs: glimpse of a common functional mechanism. protein-coupled receptor kinases is dependent on agonist efficacy. BMC Bioinf. 2015, 16, 124. Mol. Pharmacol. 1998, 54, 704−711. (45) Kaserer, T.; Lantero, A.; Schmidhammer, H.; Spetea, M.; (64) Burgueño, J.; Pujol, M.; Monroy, X.; Roche, D.; Varela, M. J.; Schuster, D. μ opioid receptor: novel antagonists and structural Merlos, M.; Giraldo, J. A complementary scale of biased agonism for modeling. Sci. Rep. 2016, 6, No. 21548. agonists with differing maximal responses. Sci. Rep. 2017, 7, (46) Mansour, A.; Taylor, L. P.; Fine, J. L.; Thompson, R. C.; No. 15389. Hoversten, M. T.; Mosberg, H. I.; Watson, S. J.; Akil, H. Key residues (65) Vasudevan, L.; Vandeputte, M.; Deventer, M.; Wouters, E.; defining the mu-opioid receptor binding pocket: a site-directed Cannaert, A.; Stove, C. P. Assessment of structure-activity relation- mutagenesis study. J. Neurochem. 1997, 68, 344−353. ships and biased agonism at the Mu opioid receptor of novel synthetic (47) Xu, W.; Sanz, A.; Pardo, L.; Liu-Chen, L. Y. Activation of the opioids using a novel, stable bio-assay platform. Biochem. Pharmacol. mu opioid receptor involves conformational rearrangements of 2020, 177, No. 113910. 1272 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article (66) Kalvass, J. C.; Olson, E. R.; Cassidy, M. P.; Selley, D. E.; (85) Seifert, R.; Wenzel-Seifert, K. Constitutive activity of G protein- Pollack, G. M. Pharmacokinetics and pharmacodynamics of seven coupled receptors: cause of disease and common property of wild- opioids in P-glycoprotein-competent mice: assessment of unbound type receptors. Naunyn-Schmiedeberg’s Arch. Pharmacol. 2002, 366, brain EC50,u and correlation of in vitro, preclinical, and clinical data. 381−416. J. Pharmacol. Exp. Ther. 2007, 323, 346−355. (86) Söldner, C. A.; Horn, A. H. C.; Sticht, H. A Metadynamics- (67) Bobeck, E. N.; Haseman, R. A.; Hong, D.; Ingram, S. L.; Based Protocol for the Determination of GPCR-Ligand Binding Morgan, M. M. Differential development of antinociceptive tolerance Modes. Int. J. Mol. Sci. 2019, 20, 1970. to morphine and fentanyl is not linked to efficacy in the ventrolateral (87) Bruzzese, A.; Dalton, J. A. R.; Giraldo, J. Statistics for the periaqueductal gray of the rat. J. Pain. 2012, 13, 799−807. analysis of molecular dynamics simulations: providing P values for (68) Trescot, A. M.; Datta, S.; Lee, M.; Hansen, H. Opioid agonist-dependent GPCR activation. Sci. Rep. 2020, 10, 19942. pharmacology. Pain Physician. 2008, 11, 133−153. (88) Bruzzese, A.; Gil, C.; Dalton, J. A. R.; Giraldo, J. Structural (69) Gillis, A.; Gondin, A. B.; Kliewer, A.; Sanchez, J.; Lim, H. D.; insights into positive and negative allosteric regulation of a G protein- Alamein, C.; Manandhar, P.; Santiago, M.; Fritzwanker, S.; Schmiedel, coupled receptor through protein-lipid interactions. Sci. Rep. 2018, 8, F.; Katte, T. A.; Reekie, T.; Grimsey, N. L.; Kassiou, M.; Kellam, B.; 4456. Krasel, C.; Halls, M. L.; Connor, M.; Lane, J. R.; Schulz, S.; Christie, (89) Serohijos, A. W.; Yin, S.; Ding, F.; Gauthier, J.; Gibson, D. G.; M. J.; Canals, M. Low intrinsic efficacy for, G protein activation can Maixner, W.; Dokholyan, N. V.; Diatchenko, L. Structural basis for explain the improved side effect profiles of new opioid agonists. Sci. mu-opioid receptor binding and activation. Structure 2011, 19, 1683− Signaling 2020, 13, No. eaaz3140. 1690. (70) Changeux, J. P.; Edelstein, S. Conformational selection or (90) Berríos-Cárcamo, P.; Quintanilla, M. E.; Herrera-Marschitz, M.; induced fit? 50 years of debate resolved. F1000 Biol. Rep. 2011, 3, 19. Vasiliou, V.; Zapata-Torres, G.; Rivera-Meza, M. Racemic Salsolinol (71) Giraldo, J. Agonist induction, conformational selection, and and its enantiomers act as agonists of the μ-Opioid receptor by mutant receptors. FEBS Lett. 2004, 556, 13−18. activating the Gi protein-adenylate cyclase pathway. Front. Behav. (72) Lans, I.; Dalton, J. A. R.; Giraldo, J. Helix 3 acts as a Neurosci. 2017, 10, 253. conformational hinge in Class A GPCR activation: An analysis of (91) Spahn, V.; Del Vecchio, G.; Labuz, D.; Rodriguez- interhelical interaction energies in crystal structures. J. Struct. Biol. Gaztelumendi, A.; Massaly, N.; Temp, J.; Durmaz, V.; Sabri, P.; 2015, 192, 545−553. Reidelbach, M.; Machelska, H.; Weber, M.; Stein, C. A nontoxic pain (73) Tehan, B. G.; Bortolato, A.; Blaney, F. E.; Weir, M. P.; Mason, killer designed by modeling of pathological receptor conformations. J. S. Unifying family A GPCR theories of activation. Pharmacol. Ther. Science 2017, 355, 966−969. 2014, 143, 51−60. (92) Subramanian, G.; Paterlini, M. G.; Portoghese, P. S.; Ferguson, (74) Rosenbaum, D. M.; Rasmussen, S. G.; Kobilka, B. K. The D. M. Molecular docking reveals a novel binding site model for fentanyl at the μ-opioid receptor. J. Med. Chem. 2000, 43, 381−391. structure and function of G protein-coupled receptors. Nature 2009, (93) Jarończyk, M.; Lipiński, P. F. J.; Dobrowolski, J. C.; Sadlej, J. 459, 356−363. The FMO analysis of the molecular interaction of fentanyl derivatives (75) Trzaskowski, B.; Latek, D.; Yuan, S.; Ghoshdastider, U.; with the μ-opioid receptor. Chem. Pap. 2017, 71, 1429. Debinski, A.; Filipek, S. Action of molecular switches in GPCRs– (94) Ellis, C. R.; Kruhlak, N. L.; Kim, M. T.; Hawkins, E. G.; theoretical and experimental studies. Curr. Med. Chem. 2012, 19, Stavitskaya, L. Predicting opioid receptor binding affinity of 1090−1109. pharmacologically unclassified designer substances using molecular (76) Vogel, R.; Mahalingam, M.; Lüdeke, S.; Huber, T.; Siebert, F.; docking. PLoS One 2018, 13, No. e0197734. Sakmar, T. P. Functional role of the “ionic lock” - an interhelical (95) Dosen-Micovic, L.; Ivanovic, M.; Micovic, V. Steric interactions hydrogen-bond network in family A heptahelical receptors. J. Mol. and the activity of fentanyl analogs at the μ-opioid receptor. Bioorg. Biol. 2008, 380, 648−655. Med. Chem. 2006, 14, 2887−2895. (77) Díaz, Ó .; Dalton, J. A. R.; Giraldo, J. Revealing the mechanism (96) Gentilucci, L.; Tolomelli, A.; De Marco, R.; Artali, R. Molecular of agonist-mediated cannabinoid receptor 1 (CB1) activation and docking of opiates and opioid peptides, a tool for the design of phospholipid-mediated allosteric modulation. J. Med. Chem. 2019, 62, selective agonists and antagonists, and for the investigation of atypical 5638−5654. ligand-receptor interactions. Curr. Med. Chem. 2012, 19, 1587−1601. (78) Bruzzese, A.; Dalton, J. A. R.; Giraldo, J. Insights into adenosine (97) Xu, H.; Lu, Y. F.; Partilla, J. S.; Zheng, Q. X.; Wang, J. B.; Brine, A2A receptor activation through cooperative modulation of agonist G. A.; Carroll, F. I.; Rice, K. C.; Chen, K. X.; Chi, Z. Q.; Rothman, R. and allosteric lipid interactions. PLoS Comput. Biol. 2020, 16, B. Opioid peptide receptor studies, 11: involvement of Tyr148, No. e1007818. Trp318 and His319 of the rat mu-opioid receptor in binding of mu- (79) Yuan, S.; Filipek, S.; Palczewski, K.; Vogel, H. Activation of G selective ligands. Synapse 1999, 32, 23−28. protein-coupled receptors correlates with the formation of a (98) Yan, F.; Mosier, P. D.; Westkaemper, R. B.; Stewart, J.; continuous internal water pathway. Nat. Commun. 2014, 5, 1−10. Zjawiony, J. K.; Vortherms, T. A.; Sheffler, D. J.; Roth, B. L. (80) Tomobe, K.; Yamamoto, E.; Kholmurodov, K.; Yasuoka, K. Identification of the molecular mechanisms by which the diterpenoid Water permeation through the internal water pathway in activated salvinorin A binds to kappa-opioid receptors. Biochemistry 2005, 44, GPCR rhodopsin. PLoS One 2017, 12, 1−14. 8643−8651. (81) Venkatakrishnan, A. J.; Ma, A. K.; Fonseca, R.; Latorraca, N. R.; (99) Zhang, P.; Johnson, P. S.; Zöllner, C.; Wang, W.; Wang, Z.; Kelly, B.; Betz, R. M.; Asawa, C.; Kobilka, B. K.; Dror, R. O. Diverse Montes, A. E.; Seidleck, B. K.; Blaschak, C. J.; Surratt, C. K. Mutation GPCRs exhibit conserved water networks for stabilization and of human mu opioid receptor extracellular ″disulfide cysteine″ activation. Proc. Natl. Acad. Sci. U. S. A. 2019, 116, 3288−3293. residues alters ligand binding but does not prevent receptor targeting (82) Kaya, A. I.; Perry, N. A.; Gurevich, V. V.; Iverson, T. M. to the cell plasma membrane. Brain Res. Mol. Brain Res. 1999, 72, Phosphorylation barcode-dependent signal bias of the dopamine D1 195−204. receptor. Proc. Natl. Acad. Sci. U. S. A. 2020, 117, 14139−14149. (100) Roche, D.; van der Graaf, P. H.; Giraldo, J. Have many (83) Celver, J. P.; Lowe, J.; Kovoor, A.; Gurevich, V. V.; Chavkin, C. estimates of efficacy and affinity been misled? Revisiting the Threonine 180 is required for G protein-coupled receptor kinase 3- operational model of agonism. Drug Discovery Today 2016, 21, and beta-arrestin 2-mediated desensitization of the mu-opioid 1735−1739. receptor in Xenopus oocytes. J. Biol. Chem. 2001, 276, 4894−4900. (101) Sader, S.; Anant, K.; Wu, C. To probe interaction of morphine (84) Zhou, B.; Hall, D. A.; Giraldo, J. Can Adding Constitutive and IBNtxA with 7TM and 6TM variants of the human μ-opioid Receptor Activity Redefine Biased Signaling Quantification? Trends receptor using all-atom molecular dynamics simulations with an Pharmacol. Sci. 2019, 40, 156−160. explicit membrane. Phys. Chem. Chem. Phys. 2018, 20, 1724−1741. 1273 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article (102) Della Longa, S.; Arcovito, A. ″In silico″ study of the binding of two novel antagonists to the nociceptin receptor. J. Comput.-Aided Mol. Des. 2018, 32, 385−400. (103) Fowler, C. B.; Pogozheva, I. D.; Lomize, A. L.; LeVine, H., 3rd; Mosberg, H. I. Complex of an active mu-opioid receptor with a cyclic peptide agonist modeled from experimental constraints. Biochemistry 2004, 43, 15796−15810. (104) Tryoen-Tóth, P.; Décaillot, F. M.; Filliol, D.; Befort, K.; Lazarus, L. H.; Schiller, P. W.; Schmidhammer, H.; Kieffer, B. L. Inverse agonism and neutral antagonism at wild-type and con- stitutively active mutant delta opioid receptors. J. Pharmacol. Exp. Ther. 2005, 313, 410−421. (105) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera–a visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605−1612. (106) Webb, B.; Sali, A. Protein structure modeling with MODELLER. Methods Mol. Biol. 2014, 1137, 1−15. (107) The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2018, 46, 2699. (108) Shapovalov, M. V.; Dunbrack, R. L., Jr. A smoothed backbone- dependent rotamer library for proteins derived from adaptive kernel density estimates and regressions. Structure 2011, 19, 844−858. (109) 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. (110) Kim, S.; Thiessen, P. A.; Bolton, E. E.; Chen, J.; Fu, G.; Gindulyte, A.; Han, L.; He, J.; He, S.; Shoemaker, B. A.; Wang, J.; Yu, B.; Zhang, J.; Bryant, S. H. PubChem substance and compound databases. Nucleic Acids Res. 2016, 44, 202−213. (111) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 30, 2785−2791. (112) 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. (113) Lomize, M. A.; Lomize, A. L.; Pogozheva, I. D.; Mosberg, H. I. OPM: orientations of proteins in membranes database. Bioinformatics 2006, 22, 623−625. (114) Huang, J.; MacKerell, A. D. J. CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data. J. Comput. Chem. 2013, 34, 2135−2145. (115) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D., Jr. CHARMM General Force Field (CGenFF): A force field for drug-like molecules compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671−690. (116) Harvey, M. J.; Giupponi, G.; De Fabritiis, G. ACEMD: Accelerating biomolecular dynamics in the microsecond time scale. J. Chem. Theory Comput. 2009, 5, 1632−1639. (117) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol Graphics. 1996, 14, 33−38. (118) Saam, J. Trajectory Smooth 1.1 TCL-script n.d. available at http://www.ks.uiuc.edu/ Research/vmd/script_library/scripts/trajec- tory_smooth/. (119) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. ■ NOTE ADDED AFTER ASAP PUBLICATION This paper was originally published ASAP on January 15, 2021. The TOC graphic and Figure 1 were replaced, and the paper was reposted on February 25, 2021. 1274 https://dx.doi.org/10.1021/acs.jcim.0c00890 J. Chem. Inf. Model. 2021, 61, 1251−1274