Bias-inducing allosteric binding site in mu-opioid receptor signaling

G-protein-biased agonism of the mu-opioid receptor (μ-OR) is emerging as a promising strategy in analgesia. A deep understanding of how biased agonists modulate and differentiate G-protein-coupled receptors (GPCR) signaling pathways and how this is transferred into the cell are open questions. Here, using extensive all-atom molecular dynamics simulations, we analyzed the binding recognition process and signaling effects of three prototype μ-OR agonists. Our suggested structural mechanism of biased signaling in μ-OR involves an allosteric sodium ion site, water networks, conformational rearrangements in conserved motifs and collective motions of loops and transmembrane helices. These analyses led us to highlight the relevance of a bias-inducing allosteric binding site in the understanding of μ-OR’s G-protein-biased signaling. These results also suggest a competitive equilibrium between the agonists and the allosteric sodium ion, where the bias-inducing allosteric binding site can be modulated by this ion or an agonist such as herkinorin. Notably, herkinorin arises as the archetype modulator of μ-OR and its interactive pattern could be used for screening efforts via protein–ligand interaction fingerprint (PLIF) studies. Agonists and a sodium ion compete for the bias-inducing allosteric binding site that modulates signaling in mu-opioid receptors. Molecular dynamics simulations of the prototype μ-OR agonist suggest a competitive equilibrium involving the agonist and an allosteric sodium ion. Analysis of experimental data from the literature and molecular models provides the structural bases of biased agonism on μ-OR. Agonists and a sodium ion compete for the bias-inducing allosteric binding site that modulates signaling in mu-opioid receptors. Molecular dynamics simulations of the prototype μ-OR agonist suggest a competitive equilibrium involving the agonist and an allosteric sodium ion. Analysis of experimental data from the literature and molecular models provides the structural bases of biased agonism on μ-OR.


Introduction
G-protein-coupled receptors (GPCRs) are seven-transmembrane (TM) helical proteins that trigger signaling pathways, mediated by interactions with heterotrimeric G-proteins and/or with β-arrestins. Preferential activation to either pathway is possible, leading to different physiological effects. Agonists that have this ability are called biased agonists, and the effect is named functional selectivity. Other modes of GPCR activation are intracellular activation, dimerization activation, transactivation, or biphasic activation [1].
Biased activation of the mu-opioid receptor (μ-OR) towards G-protein is relevant for the design of painkillers with reduced detrimental side effects. For example, activation of β-arrestin signaling pathways is suggested to lead to tolerance, desensitization, addiction, constipation, and respiratory suppression [2]. Albeit not fully understood, these effects seem to be related to β-arrestin recruitment or receptor internalization. Moreover, the coexistence of diverse and complex GPCR activation types is an active research topic. This complexity is further evidenced by the ability of GPCRs to adopt different conformations. Experimental evidence shows that different agonists can stabilize these conformations, leading to new conformational landscapes [3][4][5]. This behavior is counterintuitive, based on the paradigm that all agonists would conduct to only one active conformation [6].
Besides the biased ligands, experimental evidence shows that a conserved allosteric binding site plays a crucial role in the activation of class A GPCRs that include μ-OR [7]. Modulation of the GPCRs via this allosteric site is related to downstream signaling and physiological effects [8,9]. The high-resolution (1.8 Å) X-ray crystal structure of the adenosine A2A receptor revealed a highly ordered cluster of atoms, composed of one sodium ion, water molecules, and conserved amino acids [8] located on TM2 and TM3, such as D 2.50 and S 3.39 . The importance of an allosteric sodium binding site nearby D 2.50 has been reviewed [7][8][9], particularly on the receptor activation process and G-protein/arrestin signaling bias. Also, the structure of the delta-opioid receptor (δ-OR), in the inactive state, showed an electron density distribution located in an allosteric site in the vicinity of D 2.50 , corresponding to a sodium ion and also unveiled N 3. 35 as an important amino acid in this binding site [9]. In agreement with these findings, mutations N 3.35 V and N 3. 35 A decreased and abolished δ-OR agonist sodium effect, respectively. Interestingly, in both cases, the δ-OR arrestin-mediated signaling pathway is increased. However, both mutations did not affect the canonical signaling, which is mediated by G-proteins. This suggests that while the allosteric sodium ion can modulate the recruitment of β-arrestin, the G-protein-receptor binding only depends on the agonist [9]. Initial experiments showed a dependence on sodium concentration for the affinity of agonists but not for antagonists [10,11]. Antagonists/ weak partial agonists are also affected by D 2.50 mutations, which become β-arrestin-biased agonist [7]. Other small molecules such as amiloride can compete with the allosteric sodium site [9,12,13].
Mutations at D 2.50 and N 3.35 in δ-OR at varied concentration of Na + affect recruitment of β-arrestin, most likely by affecting the TM7 helical bundle [9]. Similarly, the relevance of TM7 movement is suggested in argininevasopressin type 2 receptor (V2R) [14] using fluorescence spectroscopy; and in β2-adrenergic receptors (β2AR) [15] characterized by 19 F-NMR experiments. These studies revealed that movements at TM6 helix are required for G-protein signaling. In contrast, β-arrestin-biased ligands impact the TM7-H8 region. This information could be transferable across class A GPCRs [7,8,16,17] where the allosteric sodium site is preserved. To note, experimental information regarding the allosteric sodium site in μ-OR is scarce [18,19], but many computational studies have shown the collapse of the allosteric sodium site in the activation and the egress pathways of the sodium ion in the presence of agonists and modulators [19][20][21]. These studies also suggest the lack of spontaneous release of the sodium ion from the bundle which required random accelerated molecular dynamics simulations to explore possible outlets from the allosteric site [19].
This important behavior has motivated the search for allosteric modulators in μ-OR. For example, ignavine increases agonist response for DAMGO and has positive modulatory activity in morphine and DAMGO assays. Furthermore, ignavine produces an analgesic effect in vivo through opioid response, proved by tail-flick and tail pressure tests [22]. However, the mechanistic modulation transmission in opioid receptors is still unclear. The positive allosteric modulator (PAM) BMS-986122 of μ-OR, which suppresses the β-arrestin recruitment in conjunction with agonists [23], can decrease the ability of NaCl to modulate the binding affinity of agonists and has been proposed to act as an indirect competitor of the allosteric sodium ion [18]. Additionally, molecular dynamics (MD) simulations of BMS-986122 suggested that residues on TM7 (7.53 and 7.35) are involved in the positive allosteric signal transmission [24]. Notably, the structure of BMS-986122 contains motifs related to C2-benzyloxy of herkinorin; this is in agreement with our suggested binding mode of herkinorin involving modulation of TM7 [25].
Comparing the X-ray structures of μ-OR in the active [26] and inactive [27] states allows analyzing helical rearrangements related to the activation process. These rearrangements are common to other members of the class A GPCR superfamily (β-2AR, rhodopsin, M2R) [7,8,16], such as helical movements at TM6 and TM7. Moreover, changes in water-mediated interactions involving amino acids inside the receptor's cavity are observed when the activation is established [26]. The active state of the receptor is characterized by an ordered network involving water molecules. Notably, careful analysis of the water network shows that it is functionally relevant, transferring information from the extracellular side (ligand binding site) through the receptor [26], and reaching the interior of the cell. This last step is characterized by collective movements of intracellular regions of the TMs and loops. This information relates to the exquisite modulation of the receptor that results from ligand binding.
Interestingly, herkinorin, a non-nitrogenous μ-OR agonist, does not produce β-arrestin2 recruitment and μ-OR cell internalization, even when the G-proteincoupled receptor kinase 2 (GRK2) is overexpressed. This is remarkable because GRK2 is known to recognize and phosphorylate activated GPCRs [28,29], promoting β-arrestin recruitment and internalization [30]. In turn, morphine-activated μ-OR generates limited β-arrestin2 recruitment and undergoes low μ-OR internalization under normal conditions. However, if GRK2 is overexpressed, internalization is readily promoted. In contrast, DAMGO (H-Tyr-D-Ala-Gly-N(Me)Phe-Gly-ol), a selective μ-OR peptide agonist [31], recruits β-arrestin2 and undergoes μ-OR internalization under normal conditions [29]. We previously proposed a binding mode of herkinorin into μ-OR. Our model suggests that herkinorin's C2-benzyloxy group locates in the allosteric sodium binding site. Binding free energy calculations of the herkinorin's proposed binding mode are in full agreement with experimental data [25]. To note, small structural modifications affect binding selectivity across opioid receptor subtypes as well as biased signaling preference. For instance, salvinorin A, a precursor of herkinorin, is a selective agonist of the kappa opioid receptor (κ-OR) and significantly recruits arrestin [32,33].
In this work, we present the structural basis of biased agonism in μ-OR with particular emphasis on the allosteric sodium site, what we call bias-inducing allosteric binding site. Taking advantage of the signaling differences for herkinorin, morphine and DAMGO, molecular dynamic simulations were performed with these agonists, while conditions in the allosteric sodium site were studied. We observed specific conformations for the NPxxY conserved motif [9], associated with the modulation at the bias-inducing allosteric binding site. Also, we concluded that this allosteric site can modify the water distribution in the intracellular side of μ-OR, in the presence of the allosteric sodium ion. We found that herkinorin can modulate this allosteric binding site preserving one conformation of NPxxY. In contrast, morphine and DAMGO systems in the absence of allosteric sodium display two extra different conformations that can be responsible for the non-canonical signaling due to the lack of modulation in this allosteric site. Finally, the conjunction between agonist and the bias-inducing allosteric binding site can regulate intracellular loops movement and the final part of transmembrane helices that are responsible for coupling macromolecular structures.

Results
The μ-OR biased signaling mechanism, proposed in this work, initiated at the orthosteric binding site continues via an allosteric binding site [26], then reaches a hydrophobic barrier (framed by DRY and NPxxY conserved motifs) [34,35] to finalize in the C-terminal helix H8, as depicted in Fig. 1. To explore different activation mechanisms, we analyzed agonists with a variety of downstream signaling effects, such as β-arrestin recruitment and receptor internalization, as well as the impact of an allosteric sodium ion. The three agonists studied here, shown in Fig. 1, are herkinorin, morphine, and DAMGO; we incorporated naloxone as the antagonist reference. Based on βarr2-GFP and confocal microscopy assays, herkinorin does not recruit β-arrestins2 and does not internalize μ-OR; DAMGO recruits β-arrestins2 and internalizes μ-OR, and morphine undergoes these processes when GRK2 is overexpressed [2,29]. To note, based on non-imaging-based β-arrestin recruitment assays, the degree of biased agonism of herkinorin has shown dependence on experimental conditions. The presence of fusion tags directly attached to the C terminus of GPCRs may explain the differences observed in these experiments [36,37]. For clarity, the systems proposed in this work as predominant are highlighted in italics throughout the text.
Binding poses were selected based on docking studies; see description in the Supplementary Material. MD simulations were carried out in triplicates. The allosteric sodium ion effect was studied for all systems except for herkinorin. The herkinorin C2-benzyloxy group occupies the allosteric sodium site, near N150 3.35 [25]. In the absence of the allosteric sodium ion, the D 2.50 amino acid was set as protonated, as suggested in previous pKa studies of GPCRs [38]. The simulation stability was assessed following the RMSD of the protein backbone and the ligand (Figures S1 and S2) clearly, all the proteins reached equilibrium. Interestingly, the conformational stability of DAMGO (+Na + ) is affected by the presence of allosteric sodium.
In what follows, we present analyses of the different systems regarding conformational changes in conserved motifs, differences in water distributions, and the dependence of the hydrophobic barrier and helical motions on the ligand studied. This information is confronted with experimental information available.

Agonist and allosteric modulation affect conformations around N332
NPxxY is a highly conserved motif for class A GPCRs; it is located in the intercellular region of TM7 and is involved in the activation processes and modifications of a hydrophobic barrier [39]. Experimental data show that mutations in the NPxxY motif of β2AR or bradykinin B2 (B2R) affect coupling with G-proteins, exhibiting phosphorylation, and subsequent sequestration of these receptors [38,[40][41][42]. The conformation of NPxxY in μ-OR active crystal structure shows that N332 7.49 and Y336 7.53 residues interact with the backbone carbonyl group of V285 6.40 through a specific water-mediated polar network ( Figure  S3). This interaction has also been observed in X-ray crystal structures of rhodopsin [26]. Notably, in the inactive state of μ-OR, the water-mediated interaction between Y336 7.53 and V285 6.40 is broken. We analyzed the interactions involving the NPxxY motif and the side-chain dihedral angles χ 1 and χ 2 of N332 7.49 , which were sampled during the simulations. In our systems, the residues Y 7.53 , N 7.49 (of NPxxY) and V 6.40 interact similarly to those observed in crystal structures of class A GPCRs ( Fig. 2 and Figure S3). Our data are also in agreement with the number and location of specific water molecules. X-ray crystal structures show water-mediated interactions between N332 7.49 and V285 6.40 , involving two water molecules for active states and one water molecule for inactive states. Crystal structures also show different water-mediated interactions for the residue Y336 7.53 . While in the inactive state this residue interacts with N 1.50 ; in the active form, it interacts with V 6.40 .
For the dihedral angle χ 1 , two distributions centered at −72 and −162 degrees were observed in the systems that involve agonists ( Figure S4). In the naloxone(+Na + ) system only the distribution at −72 degrees is favored. Remarkably, the dihedral angle χ 1 at −162 degrees was more populated in the morphine(−Na + ) and DAMGO(−Na + ) systems with respect to morphine(+Na + ) and DAMGO(+Na + ) systems which populated distributions at −72 degrees. In turn, the herkinorin(−Na + ) system retained the shown to the right. As reference, the trend of β-arrestin2 recruitment is shown in grayscale underneath the chemical structures. Darker color represents higher recruitment based on βarr2-GFP assays [29] conformation at −72 degrees as the system in the presence of the allosteric sodium ion. The experimental information discussed previously shows that χ 1 angle is found at −163 degrees for the cryo-EM structure 6DDE and −76 degrees for the X-ray crystal 5C1M of μ-OR.
All in all, this information suggests that, in the absence of sodium ion, herkinorin is modulating the NPxxY motif (TM7), thus affecting downstream signaling. In contrast, morphine and DAMGO are dependent on allosteric sodium ion modulation, which is lost in the activation process. Thereby herkinorin seems to be a sodium-independent modulator of G-protein signaling. It is expected that the μ-OR in the active state, obtained bound to a G-protein mimetic nanobody (Nb39), retains NPxxY configurations, just as it does when it interacts with G-proteins. The fact that herkinorin retains the conformation observed by NPxxY in the crystal structure agrees with our model.
Let us describe the second and third conformations observed for the χ 2 dihedral angle. Our simulations show that morphine(−Na + ) and DAMGO(−Na + ) populated χ 2 at 85 or −120 degrees. Additionally, water molecules mediating N332 7.49 −V285 6.40 interaction are relocated, and they no longer interact with Y336 7.53 . Furthermore, the interaction between D114 2.50 and N332 7.49 nitrogen was mediated by one water molecule (Fig. 2). This suggests that different conformations of NPxxY are involved in different downstream signaling pathways, such as phosphorylation in Y336 7.53 when it is not participating in the water network interaction. Non-canonical signaling has been tested experimentally through Y336 7.53 phosphorylation of μ-OR [44].
For the observed changes in NPxxY, we propose a dynamic equilibrium in morphine(+Na + )/morphine(−Na + ) and DAMGO(+Na + )/DAMGO(−Na + ) systems. The equilibrium seems to be favored towards morphine(−Na + ) or DAMGO(−Na + ). Moreover, experimentally DAMGO promotes internalization of the receptor more easily than morphine. The conformation that seems to be related to this signaling pathway is the one with χ 2 dihedral angle at either 85 or −120 degrees. Similar conclusions can be drawn from the analysis of χ 1 . These configurations are favored in morphine(−Na) and DAMGO(−Na) systems. We then explored if water networks were also distinctive between agonists.

Water network is involved in activation but not in biased agonism
Crystal structures of GPCRs have shown functionally relevant water networks. Inspection of our simulations provided relevant water-mediated intercommunications in the orthostatic site. In all the systems studied, the most remarkable water-mediated interaction was primarily observed between H297 6.52 and the ligand. The exception is for DAMGO(-Na + / +Na + ), which could make this contact in the presence or absence of water molecules. In the active state, the μ-OR crystal structure shows a water-mediated interaction between H297 6.52 and its co-crystal agonist BU72 [26]. This interaction is stabilized by the backbone carbonyl group of K233 5.39 , similarly to our systems of morphine(+Na ± Na + ) and herkinorin. However, herkinorin(−Na + ) retained this interaction through the carbonyl group of C4-methyl ester ( Figure S5). Additionally, H297 6.52 interconnects with D114 2.50 by four water molecules on average; this interaction was observed for all the systems studied except DAMGO(−Na + ), where W293 6.48 disrupts the network. If the sodium ion was present, it contributed to mediate this interaction. Other different networks in the crystal structures are the interaction mediated by one or two water molecules, observed between ligands and amine group of K233 5.39 or the water-mediated interactions between D114 2.50 -D147 3.32 amino acids that were stabilized by neighbors amino acid. Thus, water networks are relevant for the activation process, but it is not clear if they are involved in the transmission of biased signaling information.

Distribution of water molecules and hydrophobic barriers through the receptors
Differences in the distribution of water molecules across the different agonists were two-fold, from one side the probability of finding voxels (cavities) inside of the receptor and on the other side on the topology of the water distribution itself, as shown in Fig. 3. Voxels are colored by the probability of finding water molecules across the simulations. Red and blue colors represent high and low probability, respectively. The topology of the volume occupied by water molecules was characterized by a narrow region (Fig. 3). This narrow region involves interactions between  Figure S6). Interestingly, R165 3.50 and Y336 7.53 are part of the conserved DRY and NPxxY motifs, respectively. These two structural motifs are the limits of a large hydrophobic barrier formed by approximately six amino acids on the inactive state of class A GPCRs where in the activation process, once the hydrophobic barrier opens, hydrogen bonding interactions emerge between DRY and NPxxY motifs, and the hydrophobic gap rearranges, as shown experimentally by Schertler's group [34,35]. In our work, for the naloxone(+Na + ) system, the probability of finding water molecules on the intracellular part of μ-OR is low (no green and red regions, Fig. 3). The distance between R165 3.50 -Y336 7.53 is larger than the one observed in the agonist systems (Table 1), supporting a hydrophobic barrier kept in the inactive state. Interestingly, the distribution of water molecules in the intracellular part of μ-OR was denser for morphine(−Na + ), DAMGO(−Na + ), and herkinorin(−Na + ) complexes than the system in the inactive state or morphine(+Na + ) and DAMGO(+Na + ). Correspondingly, our data show comparable conformations to those observed at M(L) 3.46 , R 3.50 , and Y 7.53 in μ-OR and rhodopsin crystal structures in the active state ( Figure S6).

The hydrophobic barrier in the herkinorin/μ-OR system
Notably, the herkinorin(−Na + ) system has a short and modified hydrophobic region (enclosed by DRY and NPxxY motifs) under the parameters implemented in the trj_cavity tool. This could be explained by the close interaction between R165 3.50 and Y336 7.53 in the herkinorin(−Na + ) system (Table 1). MD simulations of oliceridine in μ-OR (which mitigates the β-arrestin recruitment) highlight the relevance of the interaction with R165 3.50 in the transmission of information initiated by biased ligands [45].
The changes in the N332 7.49 conformation and the water-mediated interactions at NPxxY (involving Y336 7.53 ) did not seem to affect the short hydrophobic barrier in all systems with agonist due to no significant differences for R165 3.50 -Y336 7.53 distances. It is consistent with the activation processes in GPCR. Although water distribution differences were observed in the intracellular part of μ-OR with agonists, this can be attributed to the presence of the allosteric sodium ion. These observations suggest that the changes observed in the conformation of N332 7.49 in the absence of the allosteric sodium ion do not affect the activation mechanism but could be responsible for alternative downstream signaling.

Structural helical and loop motions: functional mode analysis (FMA)
We performed an extended analysis of the μ-OR backbone chain's structural collective motions and preservation of helicity over TM7 information. Typically, principal component analysis (PCA) provides compacted information about projections or principal components (PCs), which means separated motions are contributing to the mean square fluctuation (MSF) of atomic coordinates. To get a large fraction (80-90%) of atomic MSF, one should analyze the first 10-20 PCs in protein simulation [46,47]; however, analysis of 10 PCs is demanding. In contrast, FMA builds a vector formed by PCs (up to 50), representing a collective motion (just one projection) that is directly related to a time-dependent biological function. In our work, the collective motion was calculated as ewMCM (ensembleweighted maximal correlation motion) related to the RMSD of the backbone chain. ewMCM is more dependent on PCs and protein energetic maps than MCM (maximally correlation motion) [46]. Table 2 shows the number of selected PCs to obtain the best model (R m ) and the best cross-validation (R c ) related to the RMSD of the protein backbone of μ-OR. Figure 4 and Figure S7 show the ewMCM oscillation between extremes p min a and p max a collective motions. As mentioned earlier, our systems reached equilibrium, as shown by the steady changes in backbone-RMSD values. Thus, the collective motion is thought of as harmonic oscillation.

ICL2, ICL3, ECL2, and H8 motion
The collective motions in morphine(−Na + ) and DAMGO(−Na + ) presented remarkable oscillations in the intracellular loops (ICL) ICL2 and ICL3 (Fig. 4). These two intracellular loops, along with TM3, TM5, and TM6, are described as the primary interactive sites in GPCR-G-protein complexes [43]. In contrast, herkinorin(−Na + ) stabilized ICL2 diminishing the oscillations in the three simulations that might help G-protein signaling. As morphine(+Na + ) and DAMGO(+Na + ) also restrained the movements in ICL2, the modulation in the allosteric sodium site may affect the ICL2 oscillations. Morphine(+Na + ) and DAMGO(+Na + ) showed oscillations at TM1 (Figure S7), which was only seen on the naloxone(+Na + ) system, indicating instability of the active state of μ-OR. This information supports the idea that the equilibrium morphine(+Na + )/morphine(−Na + ) favors morphine(−Na + ) system and that the DAMGO(+Na + )/DAMGO(−Na + ) equilibrium favors DAMGO(−Na + ). The extracellular loop (ECL) ECL2 in herkinorin(−Na + ) is modified concerning the other systems, leading to more oscillations. ECL2 has been linked to the positive allosteric modulation of ligand binding in μ-OR to G-protein signaling as the binding mode of the biased agonist oliceridine [45,48].
Interestingly, the DAMGO(−Na + ) system has a preponderant collective backbone motion over the H8 helix compared to morphine(−Na + ) and herkinorin(−Na + ) systems. Under normal experimental conditions, DAMGO induces μ-OR internalization and phosphorylation at T 370 , and S 375 amino acid of H8 could explain recruitment by arrestins [29]. Moreover, the computational work of oliceridine suggests that F347(H8) is one of the most relevant amino acids in the transmission of information in μ-OR [45]. This means sizable and specific H8 motion would facilitate the sampling of conformations that produce phosphorylation, which can, in turn, induce receptor recruitment.
DAMGO(−Na + ) system shows functionally relevant motion at H8, and the conformation at N332 7.49 (in NPxxY motif) was different from DAMGO(+Na + ), morphine(+Na + ), and herkinorin(−Na + ), suggesting that the conformation of N332 7.49 would have influenced TM7 downstream signaling to affect H8 motion. All in all, systems started signaling in the orthosteric-position of agonist and the modulation of the allosteric sodium site, passing by the modification of NPxxY, which altered H8 motion. Although the herkinorin(−Na + ) system had no allosteric sodium ion, it did not produce larger motion at H8, particularly noticeable in two of the simulation replicas, explained by the retention of N332 7.49 conformation. This is similar to the one observed in the active state crystal structure because the C2-benzyloxy group of herkinorin modulated the allosteric sodium site preventing H8 and NPxxY changes. Interestingly, in the herkinorin(−Na + ) simulation where H8 displacement was observed (simulation 2), the ECL2 did not suffer any change.

Helicity
During the MD simulations, the degree of helicity was modified ( Figure S8). For example, the herkinorin(−Na + ) system had lower degrees of helicity, such as morphine(−Na + ) and DAMGO(−Na + ), at W318 7.35 (Fig. 4). Also, the herkinorin(−Na + ) system showed loss of helicity at C336 7.53 and L335 7.52 , which are more stable compared to morphine(−Na + ) and DAMGO(−Na + ) systems. Because of N332 7.49 proximity, the unstable helicity in this region could be affected by the χ 2 angle conformations of N332 7.49 in morphine(−Na + ) and DAMGO(−Na + ). Interestingly, it has been proposed that W318 7.35 and C336 7.53 are involved in the TM7 transmission of mechanistic signaling, creating flexions at TM7 in the presence of a PAM [24]. By this means, herkinorin could induce positive downstream signaling over TM7.

The relation of binding mode with functional selectivity in μ-OR. Protein-ligand interactions
In previous studies, we [25] and others [49,50] have explored protein-ligand interaction fingerprints (PLIF) to classify and identify biased vs. non-biased ligands. Here, we generated a PLIF based on the contacts obtained from the MD simulations and compared biased vs. non-biased μ-OR agonists. Protein-ligand contacts within a radius of 4.0 Å were taken into account. Analysis of the contacts for herkinorin(−Na + ), morphine(−Na + ), and DAMGO(-Na + ) systems (Fig. 5), showed amino acids that make contacts during the entire simulation time, for example, D147 3.32 , I296 6.51 , V300 6.55 , W318 7.35 , and I322 7.39 . We also observed that these contacts were present for the naloxone(+Na + ) system ( Figure S9) during more than 90 percent of the simulation time. Interactions involving these amino acids are found in active and inactive crystal BU72-μ-OR and could be a reasonable representation of the orthosteric binding site. Prevalent and unique contacts of the herkinorin(−Na + ) system corresponded to two allosteric binding sites. The first one encompasses the amino acids A117 2.53 , V143 3.28 , and C217 ECL2 (TM2, TM3 and ECL2), which have been reported for the biased agonist oliceridine (TRV-130) (Q124 2.60 , V143 3.28 , and C217 ECL2 ) [45]. Oliceridine promotes G-protein signaling [51], similar bias-agonists mitragynine and NAP also interact with these allosteric amino acids [36,52,53]. The second binding site for herkinorin included contacts with N150 3.35 and S329 7.46 , corresponding to the allosteric sodium region. Consistently, a representative pose of oliceridine from the MD simulations makes contact with S329 7.46 through its methoxy-thiophene moiety [45].

Binding mode is affected by the allosteric sodium ion
Morphine has the same PLIF profile in the presence or absence of the allosteric sodium ion with a small variation in L233 5.39 , but DAMGO did not ( Figure S9). The lack of allosteric sodium ion DAMGO researched extracellular regions on ECL2 and TM5 making contacts with T218 ECL2 , L219 ECL2 , and L232 5.38 , whereas the DAMGO(+Na + ) system makes contacts with N127 2.63 , W133 ECL1 , and V143 3.28 , interacting with the extracellular part of TM2-TM3 and ECL1. It is important to note that these contact differences were mostly due to placements of Phe and Gly-ol group in DAMGO. If DAMGO(−Na + ) system predominates in the equilibrium DAMGO(+Na + )/DAMGO(−Na + ), β − arrestin signaling will be favored because, in contrast to biased agonists, the extracellular region of TM2 and TM3 cannot be reached all time.

Discussion
In summary, analysis of water distribution, conformational _ in NPxxY, and water networks suggested that DRY, NPxxY motifs and intracellular part of μ-OR are affected by agonist presence and the modulation of the allosteric sodium binding site. The main alteration in the NPxxY motif (in the hydrophobic limit) was on the conformation around the χ 2 angle of N332 7.49 amino acid and the water-mediated network. The reason for this alteration seems to be related to the modulation of the allosteric sodium pocket but not in the orthosteric site. The binding modes observed for morphine and DAMGO did not affect the conformation at the NPxxY motif. The conformation involving N 7.49 and water molecules in DAMGO(−Na + ) and morphine(−Na + ) systems gave a conformation at the NPxxY motif that may explain β-arrestin signaling pathways. Notably, herkinorin(−Na + ) retained a conformation at the χ 2 dihedral angle as the one observed in the crystallographic structure, in the active state. In the presence of herkinorin, there is no need for an allosteric sodium ion to control χ 1 and χ 2 angles, granted that the C2-benzyloxy group occupies this pocket. This site could be called, G-protein bias-inducing allosteric binding site since it can be modulated by Na + or biased ligands reaching that site. The impact of information transmission on the μ-OR's intracellular backbone was mainly observed in H8 helix and ICL3 loop movements by FMA. In this analysis of collective harmonic motion in TMs and loops, morphine(−Na + ) and DAMGO(−Na + ) produce changes in the ICL2 while herkinorin(−Na + ) restrains the loop movement. H8 helix in DAMGO(−Na + ) has larger displacement, particularly noticeable in two of the simulations. Although H8 did not have a complete sequence in our system, it can be suggested that this motion could affect the rest of the helix. We propose an equilibrium between morphine(+Na + )/morphine(−Na + ) and DAMGO(+Na + )/DAMGO(−Na + ). This equilibrium may alter the transduction of the signaling pathway. The main outcome is the understanding of the modulation of the allosteric sodium pocket, exemplified by herkinorin(−Na + ) through the C2-benzyloxy group. New μ-OR biased ligands can be designed by maintaining the interacting pattern obtained from the PLIF studies in the allosteric sodium region. This work highlights the relevance of the allosteric modulation in the extracellular part of class A GPCRs, in agreement with previous studies [54]. It has been proposed that, in contrast to non-biased agonists, biased agonists with specific and extended binding modes convey a closure of the binding pocket. Biased agonists could also promote an equivalent effect for the closure of the biased-inducing allosteric binding site. Our data show that the sodium ion or herkinorin induce the retention of conformations in conserved residues of class A GPCRs. The C2-benzyloxy group of herkinorin or the sodium ion participate in the closure of the allosteric site. As a result, specific water network interactions that connect the extracellular and intracellular regions of μ-OR are also induced by the allosteric modulation.

Molecular alignments
The flexible alignment was carried out through the software molecular operating environment (MOE) developed by Chemical Computing Group Inc. This process is useful to find structural similarity features through property density functions [55,56]. Co-crystal ligands BU72 and β-FNA (from PDBIDs 5C1M and 4DKL) [27] were used as structural ligand references in the orthosteric binding site of active and inactive μ-OR, respectively. On the other hand, agonist morphine and antagonist naloxone were placed into μ-OR with molecular alignment since they are structurally related to ligand references. The calculation parameters were adjusted by previous work [25].

Docking
Automated docking was performed for ligand agonists DAMGO and herkinorin. Processed with the software Internal Coordinate System, (ICM) 3.8-3 44 . All non-receptor molecules were removed, and the binding site was defined using the ICM PocketFinder algorithm 45

Molecular dynamics simulations setup
Active and inactive states of μ-OR were obtained from their corresponding X-ray coordinates (PDBIDs 5C1M and 4DKL, respectively) [9,27]. Amino acid sequence range was from M65 to F347 in the active state, whereas the range in the inactive state was from M65 to I351. First, all non-receptor molecules were removed. Then, acetyl and methylamide capping groups were placed on N-terminal and C-terminal amino acids, and the disulfide bond between the amino acids C140-C217 was inserted. The missing amino acids on ICL3 (MLSGSK sequence) in the inactive state were included with ProMod Version 3.70 implemented in SWISS-MODEL [60][61][62]. Histidine residues were simulated as the neutral tautomers. Protonation states of all amino acids were assigned at pH 7 except for D114 2.50 , protonated only in allosteric sodium ion absence.

Molecular dynamics simulations
Simulations were carried out using GROMACS 5.0.4, and three simulations were performed for each system [77]. Energy minimizations were conducted using the steepest descent algorithm. Equilibration and MD runs were performed with the leap-frog integration algorithm at 310 K and 1 bar. The systems were equilibrated through seven steps decreasing harmonic restrictions on protein, lipid, and ligands (Table S1). Runs corresponding to production step consisted of 500 ns with 2 fs time step applying the Nose-Hoover thermostat (relaxation time 2 ps) and the semi-isotropic Parrinello-Raman barostat (relaxation time 5 ps). Bond lengths to hydrogen atoms were constrained with LINCS. Receptor-ligand, lipids, and ion-solvent sets were considered as three separated temperature-coupling groups. Non-bonded interactions used the Verlet scheme with a 12.0 Å cut-off. Long-range electrostatic interactions were computed using the particle mesh Ewald (PME) method, and the Lennard-Jones (LJ) interaction smoothly decayed to zero crossing 10 Å distances. The initial conformations of morphine and naloxone were based on the molecular alignment to the crystallized ligand coordinates in the active and inactive states, respectively. The conformation and binding mode of DAMGO and herkinorin were selected from docking studies, preceding binding models, and experimental information [25,57,58,[78][79][80].

Water molecules distribution in μ-OR
The analysis was performed with the tool trj_cavity [81] implemented within the GROMACS package. This tool is handy to identify and detect cavities in MD trajectories with voxels (volume pixel), including solvent analysis. The number of surrounding protein directions was set to six, ideal to find internal protein cavities. Water molecules density distribution was extracted according to occupation frequency of water in a cavity voxel through 500 ns MD.

Amino acid contact
The amino acid contact analysis was undertaken with the MDAnalysis library in the python language [82,83]. Protein-ligand interaction contacts were calculated on 500 ns MD. First, we selected amino acids for each 0.5 ns time step with distances less than 4.0 Å from ligand atoms. Afterward, frequency contact was extracted for each amino acid over all the simulation time. Finally, contacts with less than 5% of frequency were not considered.

Functional mode analysis
We analyzed the collective atomic motions of μ-OR using the methodology known as functional mode analysis (FMA) [46]. This technique allows detecting collective motion related to a function such as a root mean square deviation (RMSD). The maximally correlated motion (MCM) extracted from FMA represents a functional state in terms of a single collective coordinate. Then, the motions that display the largest possibility to induce a change in the functional quantity are estimated from the given protein ensemble.
Hence, we diagonalized the covariance matrix C to obtain eigenvector e i and eigenvalues 2 i to build principal components p i .
FMA uses e i to establish an MCM that depends on the time by the projection p a along the normalized collective vector a correlated to the change in the RMSD or some function f through Pearson's coefficient R.
The ensemble-weighted maximally correlated motion (ewMCM) is related as the most probable collective motion v * generating a functional state p * a that is an arbitrary and fixed p a value between extremes p min a and p max a .
p * a = v * ⋅ a For each value of pa*, a set of pi* is given and related by the expression.
For each value of p * a a set of p * i is given and related by the expression We can visualize the ewMCM for the given ensemble according to the expression Since p a (t) is a collective motion dependent on time, p min a and p max a are motions in the minimum or maximum of RMSD around the average. If the system reaches equilibrium (constant RMSD), an oscillatory motion will be observed. We visualized the μ-OR backbone (−N-C α -CO-) ewMCM correlated to RMSD in all systems p * i e i since this weighted collective motion is more dependent on PCs and protein energetic maps than just the MCM [46]. We used 500 ns MD to perform FMA and then selected the best PCs set based on correlation coefficients of models and cross-validation. Models were built from the first 400 ns, while the last 100 ns completed the cross-validation. Table 2 gives the PCs number used in the FMA.

Helicity
Data were extracted by GROMACS 5.0.4 inbuilt tools gmx helix, based on an ideal helix that is formed by at least three consecutive amino acids with ideal Φ and Ψ dihedral angles (axis of rotation on N-C α and C α -C C=O , respectively) [84,85]. For each amino acid, both time percentage and α-helix structure on 500 ns were calculated from MD.

Conclusion
The molecular dynamic simulations for μ-OR performed here show that the modulation of the bias-inducing allosteric binding site influences conformational changes in the conserved NPxxY motif of class A GPCRs. Also, the water distribution in the intracellular part of μ-OR can be attributed to this allosteric site. In the absence of the allosteric sodium ion, the agonist herkinorin can modulate the bias-inducing allosteric binding site and can preserve the conformation of NPxxY (TM7). This can explain the biased agonism of herkinorin for G-protein signaling in μ-OR. In contrast, morphine and DAMGO systems in the absence of the allosteric sodium ion display two different NPxxY conformations that can be responsible for the lack of G-protein canonical signaling. Lastly, the conjunction between agonist and the modulation of the bias-inducing allosteric binding site in μ-OR appears essential for the regulation of movements at intracellular loops and transmembrane helices. These regions directly interact and recognize G-proteins or arrestins.
AMM performed analysis, investigation, visualization, and writingreview and editing, and KMM contributed to conceptualization, project administration, resources, funding acquisition, supervision, and writing-review and editing.

Declarations
Conflict of interest There are no conflicts to declare.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.