This chapter contains the description of our MD simulation results: (i) the fentanyl and morphine binding modes supplemented by mutagenesis data, (ii) the description of the differences in the dynamic behavior of the binding site residues, and (iii) the principle components analysis. The discussion of presented results finishes this chapter.

Fentanyl binding mode

We undertook several approaches for the determination of the fentanyl binding mode. Our idea was to find the binding mode by bridging molecular docking and molecular dynamics (following the workflow similar to the one proposed in ref. 40). The ligand was docked either manually or automatically to the binding site of: i) the inactive receptor without a sodium cation present in the allosteric site (4DKL, FENT-IN), ii) the inactive receptor with a sodium cation present in the allosteric site (4DKL with sodium at Asp1142.50, FENT-IN-NA), and iii) the active receptor (FENT-ACT), see “Binding mode determination”.

Molecular dynamics runs of length up to at least 500 ns were performed in all cases. The simulations with unstable behavior of fentanyl in the binding site were discarded. For example, no stable fentanyl positioning was found in FENT-IN simulations, and therefore this configuration was omitted from further investigations. Finally, a single binding mode to which several FENT-IN-NA simulations and FENT-ACT simulations had converged was chosen for the extended MD runs. Therefore, similar interactions in the binding pocket for both FENT-IN-NA (Fig. 4) and FENT-ACT (Fig. SM-1) were observed.

In this binding mode, the protonated amine of the piperidine ring interacts stably with Asp1473.32, as expected of high-affinity μOR agonists. The N-phenethyl chain is oriented toward the intracellular side, while the 4-N-phenylpropanamide chain toward the extracellular side. Both these chains retain some residual mobility and in some simulations a sort of rotational alternation can be observed without the overall change of the binding mode (Fig. 3a and b).

Fig. 3 The fentanyl side-chains retain some residual mobility in the binding pocket: a 4-N-phenylpropanamide mobility, b N-phenethyl mobility (conformers are represented with different colors). Helices of the receptor are colored from N-terminal (blue) to C-terminal (red) Full size image

Figure 4 presents the intermolecular contacts of the fentanyl-μOR complex (3D and 2D projections). Apart from the mentioned charge-assisted hydrogen bonding pipNH+…Asp1473.32, the piperidine ring also forms hydrophobic interactions with Gln1242.60 and for some time with Ile3227.39 or Tyr3267.43. The N-phenethyl chain in the most frequent arrangement is localized so that the aromatic ring of fentanyl forms hydrophobic contacts with Tyr3267.43 for some of the simulation periods with Ile2966.51 and Ile3227.39 (dependent on the rotamer of the N-phenethyl chain) and/or with Trp2936.48 (dependent on the rotameric state of the Trp2936.48 side chain). As to the 4-N-phenylpropanamide terminus, two groups of residues interact either with the phenyl ring or with the propionyl chain of fentanyl. The first grouping involves: Gln1242.60, Trp133ECL1, Val1433.28, Cys217ECL2, Asp216ECL2 and the second one: Ile1443.27, Tyr1483.33, Thr218ECL2, Leu219ECL2, Leu2325.38.

Fig. 4 Intermolecular contacts of fentanyl in the receptor cavity presented for FENT-IN-Na+: a in 3D view and (all helices˗gray, Leu and Ile˗green, Asp˗red, Trp˗brown, Tyr, Gln˗orange, His˗blue) and b 2D view. Yellow point on left part of ligand is the sodium ion Full size image

Additionally, we observed that the 4-N-phenylpropanamide may interact with His54 of the receptor N-terminus in the active structure (FENT-ACT) where the flexible N-terminus of the receptor is present. One can take into account that experimentally, there is some evidence that N-terminus is important for fentanyl’s binding (lowering the K D by ten times upon N-term deletion) but not for that of morphine [41].

Morphine binding mode

In the case of morphine, the determination of the binding mode was straightforward, after only a very minor readjustment from the pose obtained by manual docking. In both the active and the inactive receptor (MOR-IN-NA, MOR-ACT), morphine keeps a very stable positioning throughout all simulations, and therefore similar interactions in the binding pocket for both MORPH-IN-NA (Fig. 5) and MORPH-ACT (Fig. SM-2) were observed.

Fig. 5 Intermolecular contacts of morphine in the receptor cavity presented for MORPH-IN-Na+: a in 3D view (all helices˗gray, Val and Ile˗green, Asp˗red, Trp˗brown, Tyr˗orange, His˗blue) and b 2D view. Yellow point on left part is the sodium ion Full size image

The intermolecular contacts (see Fig. 5) between the ligand and μOR involves first and foremost the interaction of the protonated amine with Asp1473.32. The oxygens of the phenol, hydroxy, and ether groups are oriented toward the extracellular side and exposed to interactions with the solvent molecules hydrating the binding site. The water molecules in the direct vicinity of these oxygens form a water network contacting for example Tyr1483.33. Morphine is also involved in a number of stable hydrophobic interactions with Met1513.36, Val2365.42, Trp2936.48, His2976.52, Ile2966.51, Val3006.55, Trp3187.35, Ile3227.39, Asn1503.35, and Lys2335.39. Such a binding mode is generally similar to the ones found in crystallography for a 4,5-α-epoxy-morphinan antagonist β-funaltrexamine [9] or for a more distant agonist derivative BU72 [10]. It is also consistent with other modeling studies that dealt with morphine so far [42, 43].

The data emerging from the μOR mutagenesis

The description of the binding modes of both ligands is supported by mutagenesis data:

(i) First and foremost, both studied ligands exhibit an ionic interaction with Asp1473.32. The mutation of this residue for alanine or asparagine (able to form H-bonds but not ionic interactions) lowers the μOR affinity of morphine more than 30 times [44]. Such a drastic drop is not observed with the Asp147Glu mutation [44], where the side-chain is negatively charged (as in Asp), but longer by one –CH 2 – unit which is apparently suboptimal, still demonstrating, however, the necessity of the N+H…COO− interaction in μOR binding. No data on the Asp1473.32 mutations are available for fentanyl; however, it is reasonable to suppose that the importance of the interaction also holds true for 4-anilidopiperidines (as a general rule in aminergic GPCRs). “Carba”-derivative of fentanyl was reported to bind with much weaker affinity [45]. (ii) Both ligands interact either directly (H-bonds) or indirectly (H-bonds through a bridging water molecule) with Tyr1483.33. This corresponds to the effect of Tyr148Phe mutation [44], upon which such an interaction is not possible (there are no atoms in the side-chain able to form an H-bond), and which was reported to bring about a threefold drop in morphine affinity and a 2- to 7-fold drop in affinity of ohmefentanyl stereoisomers. (iii) An identical mutation on Tyr3267.43, namely Tyr326Phe, was reported to cause a significant deterioration in fentanyl affinity (more than 72 times), while a not so dramatic one for morphine (about six times) [46]. In our simulations, for neither of the ligands there exists a direct H-bonding interaction with Tyr3267.43. On the other hand, Tyr3267.43 is involved in the interaction with Asp1473.32 (previously called “the 3-7 lock”). The interaction is either direct or mediated by a water molecule (or dissolved). Our simulations show (these results are described in “Differences in the dynamic behavior of the binding site residues”) that this contact is tighter for fentanyl. Hypothetically, the presence of the H-bonding ability in this residue is important for adequate positioning of the side chains for μOR interactions with fentanyl, and therefore a mutation depriving such a function impacts fentanyls’ affinity to a much greater extent than that of morphine. (iv) A mutation where the effect on affinity and potency is clear for morphine (5 x ↓ affinity) but not for fentanyl (2 x ↓ affinity) is Trp318Leu [47]. In our simulations, fentanyl did not form stable contacts with Trp3187.35. Contrarily, morphine did form hydrophobic interactions, so an exchange for a smaller side-chain resulting in a smaller contact surface should consequently bring about some lowering of the affinity.

A generalized glance at both ligands

In a generalized view (Fig. 6), both morphine and fentanyl interact with TM3, TM5, TM6, and TM7. The positioning of fentanyl and morphine overlap only to a minor extent (Fig. 7), even though some of the contacts are present in both cases (Figs. 4 and 5). Notably however, the extent of these contacts does differ in both cases.

Fig. 6 Overlap of fentanyl and morphine positioning in the μOR binding pocket. Helices of the receptor are colored from N-terminal (red) to C-terminal (blue) Full size image

Fig. 7 Both ligands contact the helices at different residues: a fentanyl (blue) and b morphine (yellow) in the μOR. The residues that are contacted by the ligand are marked in red (simplified depiction) Full size image

Fentanyl’s atoms are in the proximity (closer than 7.0 Å) of as many as seven TM3 residues, while for morphine it is only four residues. The situation is somehow reverse if TM6 is considered. Here, morphine is able to contact five TM6 residues, and fentanyl only three. Regarding TM7, both ligands are in contact to a similar extent, but its place is different. Morphine interacts closer to the extracellular side (Trp3187.35), and fentanyl at the bottom of the binding pocket (Tyr3267.43). Fentanyl, as it stretches toward the extracellular side with the 4-N-phenylpropanamide chain, is also able to reach ECL1, ECL2, and the N-terminus (present in 5C1M structure).

Differences in the dynamic behavior of the binding site residues

The comparison of the experimental structures in the inactive and active crystal forms has revealed main “micro-switches” to explain the inhibition/activation mechanism of GPCR. There are putative “Trp rotamer toggle switch” and breaking/creating of “3–7 lock” [48, 49]. Apart from contacting different residues of the μOR binding pocket than morphine, fentanyl contacts the elements common to both ligands but induces some different effects. These include the conformation of: Trp2936.48, Asp1473.32, and Tyr3207.43.

The first of the mentioned residues, highly conserved amino-acid Trp2936.48, seems particularly important since it has been thought to be directly involved in the activation process in the so-called “Trp rotamer toggle switch” (called “transmission switch” too [49]). This “micro-switch” hypothesis was based on the MC simulation of the TM6 in cannabinoid CB 1 receptor [50, 51]. More recent MD simulations performed on inactive and activated μOR crystal structures seem to support it [52,53,54]. This residue links the agonist binding site with the movement of TM5 and TM6 through the rearrangement of the TM3–5-6 plane. In our simulations, the differences in χ 1 and χ 2 angles in Trp2936.48 are clearly seen for morphine and fentanyl interactions with μOR. Therefore, we monitor the process through the Trp2936.48 movement, which can be evidenced by the time course of the χ 1 and χ 2 angles.

In the inactive 4DKL crystal structure of receptor, Trp2936.48 has dihedral angles χ 1 = −77° and χ 2 = 79°, while in the active crystal structure 5C1M, these values read −87° and 121°, respectively. The simulations of the empty receptor started from either the active (3 out of 3 APO-ACT replications) or the inactive (2 out of 3 APO-IN-NA replications) converge to values of about −70° and 110°, respectively (Fig. 8a and b). However, one APO-IN-NA simulation results in another rotamer of Trp2936.48 (χ 1 ~ 180° and χ 2 ~ 110°, presented in Fig. 8b). The first option is unexceptionally present in MORPH-IN-NA and MORPH-ACT simulations. In the case of FENT, both possibilities are sampled and additionally, a third rotamer appears (χ 1 ~ 180° and χ 2 ~ −110°, see Fig. 8c). The −70°/110° rotamer is associated with a particular, stabilized conformation of the N-phenethyl chain [t, g−, g−], and with the chain flexibility at ϑ 2 and ϑ 4 ([t/g−/g+; t,g−/g+] conformers (for ϑ definitions see Scheme 1).

Fig. 8 The Trp2936.48 rotameric states found in the simulations: a MORPH and the majority of APO simulations; b and c FENT simulations. Note that rotameric states of Phe2896.44 are very similar in all cases Full size image

To better understand the structural changes for Trp2936.48, the distributions of two dihedral angles are presented in Fig. 9a, b, and the time evolution of these angles is illustrated in Fig. 9c, d.

Fig. 9 Distributions of Trp2936.48 dihedral angles: a χ 1, b χ 2 . The histograms are made based on the last 300 ns for three replications (900 ns in total); The plots c and d show time dependences of χ 1 and χ 2 for Trp2936.48 during 1.2 μs simulations. Each component picture demonstrates a separate simulation. The top, middle, and bottom row correspond to APO, FENT, and MORPH, respectively. Observe that only for FENT are different values of the torsion angles stabilized Full size image

The positions of the χ 1 and χ 2 torsional angles in the APO, MORPH, and FENT simulations can be easily read from the angle distributions (Fig. 9a and b). Observe that in the APO and FENT simulations based on the inactive receptor, two positions of the χ 1 angle are found, while in the other simulations only one value is populated (Fig. 9a). It seems to be important that MORPH simulations indicate the same value of χ 1 for both inactive and active receptor structures. Moreover, the χ 1 angle in APO simulations of the active receptor structure is practically the same. In contrast, in FENT simulations based on the active structure of receptor the χ 1 angle is very clearly different and equal to ca. −160° (Fig. 9a). On the other hand, the χ 2 torsional angle value is split only for the FENT simulations based on the inactive receptor (Fig. 9b). Interestingly, the χ 2 values in APO and MORPH simulations are basically identical in the inactive or active receptor. Again, the value of the χ 2 angle is very clearly different than that in APO and MORPH simulations, and it equals ca. −110° (Fig. 9a). The exceptional values of the χ 1 and χ 2 torsional angles in the FENT simulations evidently distinguish the fentanyl interactions with μOR from those of morphine.

The time dependences of the χ 1 and χ 2 torsional angles in the APO, MORPH, and FENT simulations (Fig. 9c and d) confirm the image given by the angle distributions (Fig. 9a and b). Indeed, there are two χ 1 angle values for the APO and FENT simulations based on the inactive receptor and one for all the other simulations (Fig. 9c). As before, the χ 1 angle indicated by the time dependence in the FENT simulations of the active receptor structure is different than that in the APO and MORPH ones. Similar conclusions can be drawn based on the time dependences of the χ 2 torsional angle (Fig. 9d).

Additionally, we determined the number of water molecules within 5 Å from Asp1142.50, as the hydration of Asp1142.50 is connected with the rotation of Trp2936.48, which opens the space toward the cell interior for water influx. The distribution of the number of water molecules is presented in Fig. 10. It appeared that the number of water molecules in all APO, FENT, and MORPH simulations in the inactive receptor structure is similar (from 4 to 11, within the statistical fluctuations) and does not differentiate the ligands. Also, the number of water molecules and distribution shape in the APO, FENT, and MORPH simulations in the active receptor structure are practically identical (from 4 to 9) and do not differentiate fentanyl from morphine.

Fig. 10 The distribution of the number of water molecules within 5 Å from any atom of Asp1142.50. The histograms are made based on data from the last 300 ns of three replications (in total 900 ns) Full size image

The χ 1 and χ 2 dihedral angles in Phe2896.44 and Trp2936.48 are known to be correlated [44]. The time dependence of the torsional angles in Phe2896.44 is presented in Fig. 11a. Two positions of the χ 1 (ca. −180° −100°) are observed in the inactive receptor. However, for APO and MORPH the χ 1 equal to ca. −100° occurs rather occasionally, whereas for FENT χ 1 = −100° is stable for quite a significant part of the simulations (Fig. 11a). On the other hand, in the active receptor for both: APO, MORPH, and FENT simulations only χ 1 = −180° is present (Fig. 11a). In the inactive receptor, two positions of the χ 2 (ca. −100° and 90°) are observed for APO and MORPH simulations but for FENT ones, in the earlier stage of simulations, up to 0.9 μs, from 2 to 5 values of χ 2 can be observed, but eventually the same two angles seem to be populated ca. −100° and 90 ° (Fig. 11b). In the active receptor, the same two positions of the χ 2 (ca. −100° and 90°) are observed for all simulations, but for APO and MORPH, flips between the two angles are very frequent, while for FENT the flips are occasional. The frequency of changes of the χ 1 and χ 2 dihedral angles in Phe2896.44 in FENT with respect to APO and MORPH seem to be the main difference of interaction between fentanyl and morphine with μOR.

Fig. 11 The Phe2896.44 a χ 1 and b χ 2 dihedral angles time dependences of the χ 1 and χ 2 during 1.2 μs simulations. Each component picture demonstrates a separate simulation. The top, middle, and bottom rows correspond to APO, FENT, and MORPH, respectively Full size image

We also observe changes of the hydrogen bond between Asp1473.32 and Tyr3207.43. It has been speculated to be involved in “the 3–7 lock”, i.e., the breaking of the link between TM3 and TM7 during the activation, which is realized by the breaking of a hydrogen bond between Asp1473.32 and Tyr3207.43. It has been suggested to play a key role in the activation of the μOR. Opening of “the 3–7 lock” was the first switch proposed to activate the rhodopsin. Here, the lock behavior is demonstrated in the form of distributions (Fig.12a) and time evolutions of the distance between Asp1473.32 and the Tyr3267.43 (Fig. 12b). When fentanyl is in the binding pocket, the Tyr3267.43 hydroxyl group remains closer to Asp1473.32 than when morphine is bound. It is particularly visible in FENT-IN-NA simulations, where distribution of the Tyr3267.43 O(H)…CG Asp1473.32 distance is focused around 5.0 Å, while for MORPH-IN-NA there is a shift toward higher values and for APO-IN-NA it is usually distributed around 7.0 Å (Fig. 12). The time evolutions support the image yielded by distributions (Fig. 12b).

Fig. 12 a The distributions of Tyr3267.43 O…CG Asp1473.32 distance. The histograms are made based on data from the last 300 ns of three replications (in total 900 ns). b The plots showing time evolution of the distance (in nm) during 1.2 μs simulations. Each component picture demonstrates a separate simulation. The top, middle, and bottom row correspond to APO, FENT, and MORPH, respectively. It seems that time dependencies of the Tyr3267.43 O…CG Asp1473.32 distance do not differentiate the APO, FENT, and MORPH systems Full size image

Additionally, except for the distribution of the Tyr3267.43 O…CG Asp1473.32 distance, the observation of the Asp1473.32 χ 1 dihedral angle changes are very important. These parameters, the distribution of the distance and the time dependence of Asp1473.32 χ 1 dihedral angle for Asp1473.32 during simulations, are presented in Fig. 13a, b.

Fig. 13 a The distributions of Asp1473.32 χ 1 dihedral angle. The histograms are made based on data from the last 300 ns of three replications (in total 900 ns). b The plots showing time dependences of the χ 1 dihedral angle for Asp1473.32 during 1.2 μs simulations. Each component picture demonstrates a separate simulation. The top, middle, and bottom row correspond to APO, FENT, and MORPH, respectively Full size image

In this case, the distributions for the inactive receptor almost do not change in the active one (Fig. 13a) for the two ligands, morphine and fentanyl. Notice, however, that only one angle is populated for MORPH, while there are two for FENT and APO simulations. The time evolution supports the image yielded by distributions (Fig. 13b).

Principal component analysis

The structural variability observed in the simulations was also investigated by the means of principal component analysis [55]. For the sake of analysis, only the last 300 ns of the simulations were taken into account. The data from different runs were pooled so to obtain a common subspace.

The first principal component (PC1) accounts for as much as 46% of the observed variance. The further four PCs are associated with 6% (PC2 and PC3), 4% (PC4), and 3% (PC5) variability. The plot of PC2 vs PC1 reveals clear clustering (Figs. 14 and 15). The points corresponding to the IN-NA runs gather at the low extremum of PC1, while those of the ACT runs at the high extremum thereof. Thus, it seems that the movement associated with PC1 corresponds to the large helical movements on the activation pathway. Indeed, it is what is revealed by the visual inspection of the trajectory projections on the average structure from PCA along the PC1 (Fig. 16).

Fig. 14 The plot of the PC2 vs PC1 values. APO-IN-NA (red), FENT-IN-NA (blue), MORPH-IN-NA (green), APO-ACT (black), FENT-ACT (orange), MORPH-ACT (pink) Full size image

Fig. 15 Superposition of the distributions of: a PC1 values, b PC2 values. The histograms are made based on data from the last 300 ns of three replications (in total 900 ns) Full size image

Fig. 16 The trajectory projections on the average structure from PCA along: a the PC1, b the PC2. Helices of the receptor are colored from N-terminal (blue) to C-terminal (red) Full size image

Notably, there is also some less pronounced grouping within the ‘active’/‘inactive’ clusters. The APO-IN-NA generally has the lowest values. FENT-IN-NA groups in the middle between APO-IN-NA and MORPH-IN-NA. This could be interpreted that APO-IN-NA runs do not go toward the activation (no G protein and no agonist present), while in the runs with the agonist it occurs at least partially. Regarding the ACT simulations, PC1 values are rather concentrated, with only some left-tailing of the APO-ACT. An interpretation for that could be that with the lack of agonist in the binding site, APO-ACT deactivates faster than the structures with fentanyl or morphine.

The yet clearer clustering within the IN-NA runs is obtained if the second PC is taken into account. PC2 values for FENT-IN-NA group to the higher extremum, while APO-IN-NA and MORPH-IN-NA are evenly distributed. The clustering in PC2 values is not so pronounced when ACT simulations are considered. The extraction of the motion associated with PC2 is not easy. However, the clustering could suggest that some receptor elements on the activation pathway inactive–active might be different in the presence of morphine or fentanyl.

To finish, we would like to present a link between the differences in dynamic behavior and the differences in signaling pathway efficacy. The activation mechanism of GPCRs involves the rearrangement of hydrophobic residues below the orthosteric binding pocket (“hydrophobic core”). This rearrangement enables large helical movements that lead to the formation of the water channel from the extracellular to the intracellular side. A residue that connects the ligand binding pocket and the hydrophobic core is Trp6.48 (highly conserved, Trp or Phe present in over 80% of GPCR A family proteins). The changes of the rotameric states of this residue have been suggested to form a “microswitch” that participates in the early activation stages by influencing the packing of the hydrophobic core residues and opening a gate for the water influx from the bulk toward the center of the receptor [56]. Presumably, if some agonists stabilize different sets of rotamers of Trp6.48, as is the case in our simulations, it can result in differences in the course of the activation process, yielding different efficacies with regard to binding the intracellular partners.

A similar suggestion could be elucidated for the observed differences in the interaction between Tyr3267.43 and Asp1473.32. Tyr326Phe (human receptor numbering) mutations have also been shown to influence opiate receptor bias with yet a different pattern. Tyr → Phe exchange disables the ability to form the hydrogen bond between the Asp3.32 and mutated amino acid in TM7, which had been formerly designated as a “3–7 lock”. In FENT-IN-NA simulations the lock seems ‘tighter’, with either direct H-bond interaction or one mediated by a water molecule happening more often than in APO-IN-NA or MORPH-IN-NA. This seems due to the different arrangement of fentanyl’s piperidine with respect to Asp1473.32, as well as a direct interaction of its N-phenethyl with the Tyr3267.43 aromatic ring.

The question of if and how these differences could lead to structural differences at the intracellular binding site cannot be answered on the basis of the simulations of the timescale presented here. An early hallmark of such possible differences may be different hydration of the allosteric sodium site, which is larger in FENT-IN-NA than MORPH-IN-NA (Fig. 10). More experimental work is necessary to elucidate this hypothesis.