Molecular mechanisms of fentanyl mediated β-arrestin biased signaling.

The development of novel analgesics with improved safety profiles to combat the opioid epidemic represents a central question to G protein coupled receptor structural biology and pharmacology: What chemical features dictate G protein or β-arrestin signaling? Here we use adaptively biased molecular dynamics simulations to determine how fentanyl, a potent β-arrestin biased agonist, binds the μ-opioid receptor (μOR). The resulting fentanyl-bound pose provides rational insight into a wealth of historical structure-activity-relationship on its chemical scaffold. Following an in-silico derived hypothesis we found that fentanyl and the synthetic opioid peptide DAMGO require M153 to induce β-arrestin coupling, while M153 was dispensable for G protein coupling. We propose and validate an activation mechanism where the n-aniline ring of fentanyl mediates μOR β-arrestin through a novel M153 "microswitch" by synthesizing fentanyl-based derivatives that exhibit complete, clinically desirable, G protein biased coupling. Together, these results provide molecular insight into fentanyl mediated β-arrestin biased signaling and a rational framework for further optimization of fentanyl-based analgesics with improved safety profiles.


Introduction
Activation of the μ-opioid receptor (μOR), a class A G protein coupled receptor (GPCR), by opiates leads to clinically desired antinociceptive properties [1], but also presents many unwanted on-target effects including addiction and potentially fatal respiratory depression. Opiates represent the largest class of prescribed drugs for the management of severe pain [2], and their subsequent abuse has led to the fastest growing public health epidemic in United States history with the Centers for Disease Control reporting a ten-fold increase in lethal overdoses from 2013 to 2017 [3]. Driven primarily by fentanyl, a highly potent and easily synthesized synthetic opiate [4], the unprecedented rate of opioid overdoses has prompted intensive research on the development of novel analgesics without the risks of standard opiates [5].
Fentanyl, and more commonly prescribed opiate alkaloids such as oxycodone, stimulate desirable antinociceptive activity by interacting with the μOR's orthosteric ligand binding site. Ligand bound μOR, in turn, mediates heterotrimeric Gi/o complex signaling to inhibit cyclic adenosine monophosphate (cAMP) production [6]. Activated μOR can also couple to β-arrestins through phosphorylation events on their Cterminal tail to promote both receptor internalization and desensitization [7]. This β-arrestin signaling axis has been hypothesized as the source for many unwanted opiate-induced side-effects as β-arrestin knockout mice exhibit prolonged morphine induced analgesia, protection against respiratory depression, opiate dependence and tolerance, and constipation [8][9][10][11].These studies suggest that decoupling of G protein from arrestin signaling through biased agonism may offer a clinically viable approach to developing safer opioid analgesics. While this hypothesis is based on knock-out models, it remains an open question whether this effect can be replicated by a small molecule.
The development of biased agonists to selectively modulate GPCR-stimulated G protein or arrestin signaling pathways for the μOR has been met with varied success [12][13][14]. Specifically, oliceridine (TRV130) provides similar antinociceptive activity to morphine while retaining a high potential for dependence and abuse with no respiratory safety advantages over morphine [15]. Another novel opiate, PZM21, was similarly shown to be less potent than morphine while eliciting undesirable respiratory depression and dose-dependent tolerance [16]. This limited success of designing G protein biased μOR agonists is ultimately hindered by a poor mechanistic understanding of how a ligand's chemical structure influences the receptors conformational state to mediate specific signaling pathways.
To gain structural insight into biased signaling, here we sought to understand the mechanism of μOR recognition and activation by fentanyl, a β-arrestin bias agonist. As ligand binding and unbinding events often occur on the timescales unreachable by standard simulations, we introduced a computationally efficient mollified adaptive biasing potential (mABP) [17] in to the graphics processing unit (GPU) enabled AMBER molecular dynamics (MD) engine to accelerate sampling of rare binding events. Starting from solvent, we simulated binding events for fentanyl, and two of its more potent derivatives carfentanil and lofentanil, to the μOR. These simulations uncovered a common pose within the receptors orthosteric site consistent with a wealth a historical selectivity and structure-activity relationship studies on the fentanyl scaffold. Within this shared pose, fentanyl induced a unique rotomeric conformation of M153 3.36 (Ballesteros-Weinstein numbering), which subsequently displaces W295 6.48 , a conserved "microswitch" critically important for receptor activation in all class A GPCRs [18]. Mutational studies of M153 3.36 reveal its role as a specific "microswitch" for μOR β-arrestin signaling as both fentanyl and DAMGO, a synthetic opioid peptide, require M153 3.36 for β-arrestin but not Gi complex coupling. We further synthesized fentanyl derivatives designed to be unable to module the M153 3.36 β-arrestin microswitch. Consistent with our computationally derived hypothesis, these compounds exhibit no detectable levels of β-arrestin coupling while retaining partial agonist-like Gi coupling. Our study provides a structural basis of fentanyl mediated β-arrestin bias signaling and provides chemical insight for the design of safer fentanyl-based analgesics.

Implementation of fABAMBER
Expanding upon our prior work to introduce a computationally efficient, scalable mollified ABP (mABP) scheme in the GROAMCS MD engine [19], here we introduced mABP with overfill protection [20] to the graphics processing unit (GPU) accelerated version of pmemd.cuda [21] for AMBER16 [22] which is available for download at Github (https://github.com/ParkerdeWaal/fABAMBER). By minimizing communication between GPU and CPU routines, our mABP implementation, termed as fast-Adaptive-Biasing-AMBER (fABAMBER), introduces nearly no overhead to the native MD engine for small and large simulation systems and provides an approximately 3-fold increase in simulation throughput compared to AMBER's Nonequilibrium Free Energy (NFE) toolkit ( Figure 1a). Beyond application to ligand binding presented here, fABAMBER can perform free energy calculations using dihedral, distance, and rootmean-square deviation (RMSD) collective variables (CVs) similar to the AMBER's NFE toolkit and the PLUMED 2 [23] plugin for GROMACS. fABAMBER thus allows users to invoke mABP methodologies while leveraging the full computational efficiency of the AMBER's GPU accelerated pmemd.cuda engine.
Additionally, both fABAMBER and fABMACS utilize a recent advancement, called overfill protection, that prevents the adaptive bias from driving the system over unrealistic barriers [20]. We discuss the significance of this feature for receptor-ligand binding simulations below.

Ligand binding to β2AR and μOR
To assess the accuracy of ligand binding pose prediction for GPCRs using fABAMBER, we first simulated binding of carazolol, an inverse agonist, to the β2 adrenergic receptor (β2AR) ( Figure S1) [24]. Ligand binding to β2AR has been previously achieved through brute force computation using ANTON, a specialized computer designed to maximize throughput of unbiased MD steps, and thus serves as a reasonable test-case [25]. Carazolol was initially placed in solvent above β2AR. The ligand-receptor conformation space was projected onto a 2D collective variable (CV) space comprising two root-meansquared deviations (RMSDs) where the ligand was decomposed into its flexible isopropylamino (CV1) and rigid carbazole moieties (CV2; Figure S1). All simulations were run for approximately 2 µs to provide ample sampling of ligand binding and unbinding events. Importantly, the underlying mABP is blinded to any exogenous information including the crystallographic ligand position and simulations were not manually stopped once the correct pose was found. This system setup was uniformly applied to all simulations performed here.
Initial simulations without bias overfill protection resulted in severe structural deformation as ligand pushed through the receptor's secondary structure (data not shown). In a series of overfill protection calibration experiments ( Figure S2), simulations performed at low (11 and 14 kcal/mol) and high bias (20 kcal/mol) potential fill-limits failed to identify prominent binding events. At a fill-limits of 17 kcal/mol, two of four simulations identifying a common low energy conformation within the receptors orthosteric site which overlapped with the expected crystallographic position with a RMSD less than 0.5 Å ( Figure S1).
Importantly, this pose successfully recaptured both crystallographic hydrogen bonding and salt-bridge interactions.
We next simulated the binding of BU72, a potent morphine-derived agonist [26], to the crystal structure of inactive μOR (Figures 1b & S3) [27]. After a series of calibration experiments of fill-limits were performed similar to that of β2AR (data not shown), we selected a 15 kcal/mol fill-limit for all subsequent μOR simulations. Prominent energetic minima were identified within the orthosteric site in two of four simulation replicates ( Figure S4). The representative ligand conformation captured near-atomic level crystallographic interactions between BU72's morphine core and the μOR ( Figure 1d). Notably, a salt bridge between the ligand's protonated amine and D149 3.32 and a water mediated hydrogen bonding network between the distal hydroxyl group to Y150 3.33 and H299 6.52 . Slight positional variance of the pendant phenyl ring was observed as we simulated the connecting carbon in its sp 3 tetrahedral configuration while the crystal structure unexpectedly contained a near-planar sp 3 species [28].

Fentanyl binding to the μOR
Encouraged by successful pose identification for carazolol and BU72, we next sought to determine how fentanyl (ED 50 = 0.011 mg/kg), as well as two more potent derivatives, carfentanil (ED 50 = 0.00032 mg/kg) and lofentanil (ED 50 = 0.0006 mg/kg), engage and activate the human μOR ( Figure 2a) [29]. Carfentanil and lofentanil both differ from fentanyl through a shared carbomethoxy group at the 4-position of the core piperidine ring while lofentanil has an additional cis-methyl group at the 3-position. From six 2-μs mABP simulation replicates for each compound, a common low-energy conformation emerged across all three ligands (Figures 2b,c & S4-S7). All fentanyl derivatives were oriented in a same latitudinal manner in the orthosteric binding site, forming extensive hydrophobic contacts between their n-alkyl phenyl and propylamide groups with transmembrane helices (TM) 2/3 and 3/5/6, respectfully. The n-aniline ring was oriented downwards towards the receptor's hydrophobic core above M153 3.36 . The piperidine ring is angled towards TM3, allowing for a salt bridge to form between its protonated amine and D149 3 Figure S8) [28].

Internal dynamics, affinity, and selectivity
Chemical modifications that both minimize ligand conformational entropy and stabilize productive ligandreceptor contacts have been hypothesized to increase ligand potency [19]. To assess how the carbomethoxy and 3-cis methyl groups of carfentanil and lofentanil affect ligand rigidity, we computed the free energy of rotation around a dihedral defined between the n-aniline and the piperidine ring. Fentanyl's largely flat landscape comprised three energetic minima with no barriers greater than 1.1 kcal/mol ( Figure   3a). Two global minima were identified within a shared super basin spanning approximately -2 and 0 radians, separated by a negligible saddle, and a secondary local minimum, S3, located at 2 radians.
Addition of a 4-carbomethoxy group (carfentanil) shifted the positions of S1 and S2 closer to -2 and 0 radians respectfully while increasing the barrier separating these two states from 0.1 to 0.9 kcal/mol. The depth of S3 shrank in favor of the primary super basin containing S1 and S2 which also saw a slight increase in its bounding barriers from 1.1 to 1.25 kcal/mol. Secondary addition of a 3-cis methyl group (lofentanil) drastically altered the dihedral landscape in favor of a single global minimum, S2, that is bounded by steep asymmetrical barriers of approximately 2.5 and 3.0 kcal/mol. Simulation of 3-cis methyl fentanyl similarly shifted the dihedral landscape in favor of the S2 conformation ( Figure S9). Across our simulations of fentanyl, carfentanil, and lofentanil, all ligands adopted the energetically favored S2 dihedral conformation state within the receptors orthosteric site (Figure 3b).
We next sought to determine the structural basis of fentanyl μOR selectivity over OR and κOR ( Figure   3c) [30]. Only two residues differ between the μOR and OR orthosteric sites: N129 2.63 is a lysine and   Figure 3d). These selectivity profiles correlate well with both ligand propensity to adopt the productive S2 dihedral conformation and the degree of steric hindrance to each receptor [30].

Structure activity relationship
To date, several hundred fentanyl analogs have been synthesized, producing a wealth of historical structure-activity-relationship (SAR) data [29]. We next sought to contextualize this SAR data by docking compounds into a rigid receptor configuration obtained from the simulated fentanyl bound state. The This n-alkyl tetrazol group is shifted sligntly compared to fentanyl and sufentanil's conjugated rings and makes largely unfavorable contacts within the hydrophobic TM2/3 binding cleft ( Figure 4a). Remifentanil, like carfentanil (ED 50 = 0.00032 mg/kg), contains a 4-carbomethoxy-substituted piperidine ring; however the n-alkyl phenyl ring is replaced with a secondary polar carbomethoxy group. The polar n-alkyl carbomethoxy group is buried down towards the receptors core's where the hydroxyl group can form a hydrogen bond Y324 7.43 ( Figure 4b).
Addition of a methyl group to the piperidine scaffold has yielded one of the more potent fentanyl derivatives, 3-cis-methylfentanyl (ED 50 = 0.00058 mg/kg) [31]. Within the lofentanil bound pose, the 3-cis methyl group occupies a hydrophobic space between I146 3.29 , D149 3.32 , and found between the S1 and productive S2 dihedral conformations without incurring steric clash with the receptor ( Figure S9). Addition of a 2-cis methyl group (ED 50 = 0.665 mg/kg) introduces significant steric clash against I146 3.29 to preclude interaction between the ligand's protonated amine and D149 3.32 .
Recent attempts to synthesize bivalent fentanyl derivatives commonly employ extension or functionalization of the amide or n-alkyl groups at the cost of decreased analgesic activity [29,32,33].
Our bound pose accommodates these large extensions without altering position of the n-aniline or piperidine rings (Figure 4f). Attempts to model conformationally restricted fentanyl derivatives introduced significant steric clash and prevented the ligand from occupying the identified binding pose (Figure 4g), consistent with the fact that these molecules exhibit insignificant analgesic activity [34,35]. Together, our bound pose of fentanyl and its derivatives provide a rational, structure based explanation of SAR for the 4-anilidopiperdine series compounds.

M153 mediates β-arrestin signaling
Confident in our computationally derived poses of fentanyl, carfentanil, and lofentanil, we next sought to understand the mechanism of ligand induced biased activation of the receptor. Focusing only on the orthosteric pocket, comparison of the DAMGO-bound active μOR-Gi protein complex cryo-EM structure to our model of the receptor bound to the β-arrestin biased agonist fentanyl revealed largely similar side chain arrangements with exception of M153 3.36 and W295 6.48 ( Figure 5a) [36]. Notably M153 3.36 is pushed downward by fentanyl's n-aniline ring, but not in the DAMGO structure, to adopt a rotomeric conformation which directly displaces W295 6.48 . This fentanyl induced M153 3.36 rotomer is also not observed in the antagonist bound or morphine derived agonist bound structures of the μOR ( Figure S10) [27,28]. Thus, we hypothesized that fentanyl's unique displacement of M153 3.36 may mediate its β-arrestin biased signaling.
To validate the role of M153 3.36 in β-arrestin biased signaling, we constructed a series of hydrophobic mutations at the M153 3.36 position of decreasing sidechain size and assessed ligand-induced assembly of the receptor with Gi protein complex or with arrestin using a direct interaction NanoBiT assay [37]. All five mutations substantially reduced the EC 50 of Gi coupling for both DAMGO and fentanyl, which displayed a correlation with the size of side chain of the mutated residue ( Figure 5b). Interestingly, while fentanyl requires M153 3.36 to achieve Gi coupling above 55%, DAMGO was less sensitive to side-chain replacement and retained near 100% coupling for larger amino acid side chains M153F/L/I 3.36 .
Unexpectedly, β-arrestin coupling for both fentanyl and DAMGO exhibit a shared dependence upon M153 3.36 . For the mutations assessed, all but larger amino acid side chains M153F 3.36 and M153F/L 3.36 for DAMGO and fentanyl, respectfully, failed to elicit coupling, and those that retained coupling saw severely reduced maximal levels ( Figure 5b). Together, these results suggest that M153 3.36 is required for βarrestin, but not Gi coupling.

Fentanyl's n-aniline ring mediates β-arrestin signaling
Given the fact that fentanyl's n-aniline ring is directly packed against M153 3.36 , a key residue for β-arrestin signaling, we next sought to determine if fentanyl's n-aniline ring is required for β-arrestin recruitment ( Figure 6a). A series of fentanyl derivatives (FD-1, -2, -3) was synthesized where the n-aniline ring was replaced with smaller aliphatic functional groups, which were predicted to have less packing effects on the rotomeric conformation of M153 3.36 (Figure 6b).
All three derivatives exhibited significantly reduced EC 50 and maximal Gi coupling compared to fentanyl, where ligand potency increased as a function of aliphatic group size (Figure 6c). One compound in particular, propyl substituted FD-3, retained partial-agonist properties comprising a mid-nanomolar EC 50 with a maximal coupling efficiency that is approximately 67% that of morphine [16] and comparable to the M153I 3.36 . All of our derivatives were completely Gi biased and failed to elicit any detectable level of βarrestin coupling up to 4 μM ( Figure 6c). Taken together with our prior mutational analysis, these experiments suggest that the n-aniline ring is required for fentanyl mediated β-arrestin, but not Gi coupling, and this effect is mediated through an interaction with M153 3.36 .

Discussion
Enhanced sampling methods for molecular dynamics (MD) are often invoked to accelerate simulations of rare events for biological systems [38], such as ligand binding and release. In this study, we implemented a previously derived efficient mABP scheme [17] with minimal communication to the GPU enabled AMBER MD engine, which significantly reduced the computational time to sample rare events while providing high-density simulation throughput.
During our initial benchmarking and system testing simulations of carazolol binding to β2AR, we found that overfill protection [20] was essential to our success. Overfill protection prevents trajectory "spoiling", arising from the bias potential pushing the ligand through the protein causing irreversible structural deformation events. Interestingly, a relationship between bias potential fill-limit and successful ligand binding events was observed. Simulations performed at lower fill-limit resulted in no binding events, whereas those performed at higher fill-limit resulted in many ligand entry events into the orthosteric pocket, but failed to identify any dominant bound conformation. For μOR binding simulations, we arrived at a fill-limit lower than that used for β2AR simulations. Unlike the orthosteric site of β2AR which is partially occluded by the receptor's extracellular loop 2, the μOR has no such barrier above the orthosteric site. These experiments suggest a careful fill-limit calibration is required to ensure that the bias potential is high enough to sample rare events without spoiling the trajectory. Likely, this value varies between GPCRs and potentially between different ligands.
Simulations of fentanyl, carfentanil, and lofentanil identified a common shared pose within the receptor's orthosteric site that makes chemically favorable interactions with the μOR. As the ligand's scaffold comprising the n-aniline, piperidine, and n-alkyl phenyl rings remains unchanged, the degree of withinand between-ligand convergence across independent mABP simulations largely improves the statistical confidence in the bound conformation. Importantly, each ligand's protonated amine was found to form a salt-bridge with D149 3.32 . As with other aminergic GPCRs, interaction with this fully conserved D 3.32 provides a critical anchoring point for ligand orientation and is an essential feature for receptor activation [18]. We note that recent attempts to "dock" fentanyl to the μOR did not observed this key interaction conserved for class A biological amine receptors [39]. For studies that correctly identify the D 3.32 salt bridge, the resulting pose is sterically incompatible with fentanyl derivatives bearing n-alkyl or amide extensions [40,41]. The mABP approach employed here offers a major advantage over virtual screening by employing a fully dynamic, unrestrained, explicitly solvated membrane environment [42]. For these reasons, this technique has also proved valued in systems with induced-fit mechanics [43].
Beyond their described roles of increasing the ligand-receptor interaction surface area, we also sought to understand how modifications to the piperidine ring affect conformational heterogeneity around the naniline and piperidine rings. All tested piperidine additions reduced conformational heterogeneity, with the 3-cis methyl group contributing a majority of this effect over 4-carboxymethyl. As lofentanil exhibited only one energetically favorable dihedral conformation, we hypothesize the S2 conformation is productive for both ligand binding and receptor activation (Figures 3a,b). Indeed, all of our μOR bound poses for fentanyl, carfentanil and lofentanil adopted the S2 dihedral state within the receptors orthosteric site.
These simulations not only inform on the structural influence of piperidine scaffold modifications, but also provide strong support in favor of the common binding conformation for fentanyl and its derivatives in the orthosteric site.
Fentanyl and its derivatives display robust affinity and selectivity for the μOR over the OR and κOR despite limited perturbations to the orthosteric site amino acid composition [30]. By modelling the OR orthosteric site onto our fentanyl bound pose, we proposed N129K 2.63 to be the largest contributor to the decreasing binding affinity. The larger side chain of the lysine perturbs the shape of the orthosteric binding pocket. Unlike the OR, the only possible W320Y 7.35 rotomer in κOR introduces direct steric clash with the ligand's rigid core. Based on this structural analysis, an ordered ranking affinity for fentanyl derivatives was predicted as μOR > OR > κOR, consistent with observed in-vivo potencies [30].
We next sought to provide a rational structure-based explanation for historical SAR of the fentanyl scaffold. Modifications to the 4-position on the piperidine ring were found to increase binding affinity by both increasing the ligand-receptor interface and propensity to adopt the productive S2 dihedral conformation, consistent with the increased efficacy seen in sufentanil over fentanyl. This increased potency can be offset by replacing the n-alkyl group with a polar moiety, which introduces repulsive interactions to the hydrophobic TM2/3 cleft. This change accounts for the affinity decrease from sufentanil to alfentanil and from carfentanil to remifentanil. These polar groups may also interfere with formation of a hydrogen bond between D149 3.32 and Y324 7.43 observed in both crystal and cryo-EM structures of agonist bound μOR [28,36]. As agonists stabilize this conserved hydrogen bond [18,44], we speculate that its disruption by polar n-alkyl groups will negatively affect receptor activation. In addition, fentanyl's bound conformation provides a rational explanation for alkyl substations around the piperidine scaffold. We further explored compatibility of fentanyl analogs with amide and n-alkyl extensions. While these functionalizations greatly diminishes μOR activity [32,33], we successfully docked examples of both ligand types in conformations that do not alter the core-ligand interaction interface comprising the naniline and piperidine rings. As we previously noted, prior studies have been unable to identify bound conformation-compatible amide and n-alkyl extensions [40,41]. Lastly, we found that rigid n-aniline and fused-alkyl fentanyl derivatives prevent ligand binding through the introduction of severe steric clash, consistent with their inability to induce antinociceptive properties [34,35]. Together, our bound pose of fentanyl and its derivatives provides a rational, structure-based explanation of SAR for the 4anilidopiperdine series compounds which further support the bound conformation of fentanyl as observed in our simulations.
Comparison of fentanyl and its derivatives to structures of agonist-and antagonist-bound μORs revealed a unique M153 3.36 rotomeric conformation previously unobserved [27,28,36]. Mediated by the common n-aniline ring, the rotomers of M153 3.36 directly influences the conserved W295 7.35 microswitch, critical for class A GPCR receptor activation. Mutational studies revealed a unique sensitivity to M153 3.36 perturbations for fentanyl-mediated but not DAMGO-mediated, maximal Gi coupling levels. As fentanyl's n-aniline ring is buried deeply within the orthosteric site to displace M153 3.36 unlike DAMGO, this sensitivity is consistent with each ligand's spatially distinct binding footprint. We thus hypothesized that fentanyl's ability to displace M153 3.36 is required for β-arrestin, but not Gi, coupling. We reasoned that replacing the bulky n-aniline ring of fentanyl with smaller aliphatic groups would diminish the ligand's βarrestin biased signaling profile. Indeed, all of the fentanyl analogs failed to induce β-arrestin coupling while one compound in particular, propyl-substituted FD-3, retained partial Gi coupling levels with nanomolar efficacy, comparable to M153 3.36 mutations. We propose a mechanism where the n-aniline ring of fentanyl mediates μOR β-arrestin through a novel M153 3.36 microswitch. Together, our study provides structural insight into fentanyl-mediated β-arrestin biased signaling and provides chemical insights to the design of novel fentanyl-based analgesics that exhibit complete bias toward Gi coupling.
In summary, through adaptively biased molecular dynamics simulations we have derived a common binding pose for the fentanyl derivatives within the µOR orthosteric ligand binding pocket. This binding pose is consistent with the vast SAR data from large numbers of fentanyl-related compounds. Moreover, this fentanyl binding mode has allowed us to determine the critical role of a single residue, M153 3.36 , for µOR β-arrestin biased signaling, which is further validated by fentanyl derivatives that are capable of mediating G protein signaling but specifically devoid of β-arrestin biased signaling. These results demonstrate the level of our understanding for the biased signaling in the µOR-fentanyl system. Beyond application to the μOR-fentanyl system, our findings also have broader impacts for our understanding on biased signaling in other class A GPCRs. Similarly, subtle changes in the chemical structures of β 2 AR ligands [45] or single residue mutation in the serotonin receptor 5-HT 2B [46] can also alter the biased signaling specificity by these receptors, demonstrating that GPCR signaling specificity can be altered by very minor changes either in the receptor itself or in the chemical structures of the bound ligands.

Molecular dynamics simulation setup and equilibration
All-atom atomospheric simulations of the human β2AR (PDB Code: 2RH1) and μOR (PDB Code: 4DKL) were performed in the NPT ensemble with periodic boundary conditions and the CHARMM36m forcefield [47]. From the inactive crystal structures, both receptors were prepared by removing all non-GPCR protein chains, fusions partners, and heteroatoms including the crystallographic ligand. Mouse μOR was humanized and subjected to 5,000 rounds of refinement using Modeller9.18 [48].ICL3 was left unbuilt for both receptors. Protonation states of titratable residues were set in agreement with prior simulations of inactive GPCRs [49] and histidines were modelled with an explicit proton on their epsilon nitrogen. The allosteric sodium interacting with D 2.50 was modeled to further mimic the receptors inactive conformation [50]. Each receptor was then capped with neutral acetyl and methylamine groups and embedded into a

Ligand binding simulations and pose extraction
To accelerate sampling of ligand binding, we implemented our previously described mABP scheme [17] with overfill protection [20] into the GPU-enabled pmemd.cuda molecular dynamics [21] engine for both AMBER16 [22]. Two RMSD CVs were defined between each ligand and the receptor (Figures S1&S3).
The receptor reference points of CV1 and CV2 are two dynamically updating centers-of-geometry comprising Cα atoms found below the orthosteric site in TM2/3/7 and TM3/5/6 respectfully. From the geometric center of the CV1 and CV2 references, a receptor center is dynamically calculated and serves as the origin for a cylindrical restraint (radius = 18 Å). This restraint is combined with a maximal allowable CV value of 22 Å to prevent unwanted ligand diffusion away from the receptor and into periodic images.
The collective variable space was discretized into a 480 x 480 grid and ranged from 0 to 60 Å for a bin width of 0.125 Å. Biasing parameters were b = 0.9, c = 0.005/ δt, α = 20. mABP production simulations of in replicates of 2 μs were performed in the NTP ensemble. All simulation parameters and trajectories are available upon request. fABAMBER patches for AMBER16 are can be found on Github (https://github.com/ParkerdeWaal/fABAMBER).
Representative pose extraction for each ligand was performed as follows: The free-energy landscape of simulations that identified energetic minima within the orthosteric site were combined and averaged.
Simulations without ligand binding events were discarded. A square region comprising the averaged minima was defined and all trajectory frames within this CV space were extracted into a combined subtrajectory comprising only the ligand. The most representative ligand configuration was determined using the DBSCAN clustering algorithm [52] within CPPTRAJ [53] where each cluster required a minimum of 25 points with a distance cutoff between points of 0.7 Å.

Fentanyl dihderal simulations
Fentanyl and its derivatives were individually solvated in a box of TIP3P waters with a minimum of 8 Å padding and neutralized with the addition of a single sodium ion. Each system was equilibrated in the NVT and NPT ensembles for 10 and 10 ns respectfully using largely the same temperature and pressure coupling schemes used for GPCR simulations where only the pressure-coupling scheme was changed to isotropic. A single common dihedral CV for all ligands was defined by atoms connecting the n-aniline and piperidine scaffold and discretized into a 300 x 300 grid ranging from 0 to 2π radians. Biasing parameters were b = 0.8, c = 0.1/ δt, α = 5. mABP simulations lasting 500 ns were performed in the NTP ensemble for each ligand.

Fentanyl derivative docking
Flexible ligand docking of fentanyl derivatives into a rigid fentanyl bound μOR receptor snapshot was performed using DOCK6.9 [54]. Prior to docking, the W320 7.35 rotomer was adjusted to accommodate 4position substations to the piperidine scaffold, each ligand was energy minimized, and partial charges were derived using CHIMERA [55].
For β-arrestin recruitment assay, the cDNAs of human μOR (wild type and M153 mutants), as well as pre- For Gi activation assay, the cDNAs of receptors were subcloned into a modified pfastbac1 vector (Invitrogen), which contained an expression cassette with a haemagglutinin (HA) signal sequence at the N terminus and LgBiT at the C terminus. The human cDNA of Gαi1, Gβ1 and Gγ2 were separately inserted into pfastbac1 vector, while the SmBiT was fused to C terminus of Gβ1 subunit.

Sf9 cell membrane preparation and Gi activation assay
High-titer recombinant baculovirus (>10 9 virus particles per ml) of μOR (wild type and mutants), Gαi1, Gβ1 and Gγ2 was obtained using Bac-to-Bac Baculovirus Expression System (Invitrogen) as previously described [56,57]. Cell suspensions were cultured for 4 days while shaking at 27 ℃ to generate P1 virus. A modified NanoBiT assay was performed to evaluate Gi activation [37]. Briefly, sf9 membrane containing μOR or its mutants were added into 384-well white plates (Greiner) at 20 μl/well. NanoLuc substrates (Promega, after 1:20 dilution) were added at 10 μl/well followed by incubation by 10 μl of 4×compound solutions for 15 min at room temperature. The luminescence was measured at 700 nm using a Perkin-Elmer EnVision plater reader. Results were plotted and analyzed in GraphPad Prism 7.0.

AD-293 cell transfection and β-arrestin recruitment assay
The human embryonic kidney AD-293 cells (Agilent) were cultured in growth medium that consisted of Dulbecco's modified Eagle's medium (DMEM) with 10% fetal bovine serum (Invitrogen) at 37 ℃ in the presence of 95% O 2 and 5% CO 2 . Cells were plated at a density of 5×10 5 cells/well of a 6-well plate for overnight. Using FuGENE ® HD (Promega), cells were transfected with 1.5 μg receptor plasmid and 1.5 μg β-arrestin2 (3A) plasmid per well and then cultured for 24 hours.
A modified NanoBiT assay was performed to evaluate β-arrestin recruitment [37]. On the day of assay, transfected cells were harvested and resuspended using Opti-MEM (Invitrogen). The cell density was adjusted to 4×10 5 cells/ml before plating cells into 384-well white plates (Greiner) per 20 μl/well. 10 μl/well of NanoLuc substrate (Promega, after 1:20 dilution) was added followed by addition of 10 μl/well indicated concentration of 4×compounds solutions. After 15 min incubation, luminescence was recorded in Perkin-Elmer EnVision plate reader, and curves were generated and analyzed using Graphpad Prism 7.0.

Fentanyl derivative synthesis
Methods for synthesis of compound FD-1: The appropriate anhydride (1.2 equiv) was added slowly to a stirred 0.5 M solution of compound 1a (1 mmol) and Et 3 N (2.0equiv) in CH 2 Cl 2 at RT. After having been stirred for 5 h, the reaction mixture was quenched with water. The organic layer was washed with saturated NaHCO 3 (aq), 1 M HCl, and brine.

Methods for synthesis of compound FD-2:
The appropriate anhydride (1.2 equiv) was added slowly to a stirred 0.5 M solution of compound 1b (1 mmol) and Et 3 N (2.0equiv) in CH 2 Cl 2 at RT. After having been stirred for 5 h, the reaction mixture was quenched with water. The organic layer was washed with saturated NaHCO 3 (aq), 1 M HCl, and brine.
The extract was dried over Na 2 SO 4 and concentrated to provide amides 2b which was used directly.MS Methods for synthesis of compound FD-3: The appropriate anhydride (1.2 equiv) was added slowly to a stirred 0.5 M solution of compound 1c (1 mmol) and Et 3 N (2.0equiv) in CH 2 Cl 2 at RT. After having been stirred for 5 h, the reaction mixture was quenched with water. The organic layer was washed with saturated NaHCO 3 (aq), 1 M HCl, and brine.
The extract was dried over Na 2 SO 4 and concentrated to provide amides 2c which was used directly.MS

Data availability statement:
Parameters, initial coordinates, and simulation restart files for all systems are available in the fABAMBER github page. All trajectories are available upon request .