How μ-Opioid Receptor Recognizes Fentanyl

The opioid crisis has escalated during the COVID-19 pandemic. More than half of the overdose-related deaths are related to synthetic opioids represented by fentanyl which is a potent agonist of mu-opioid receptor (mOR). In recent years, crystal structures of mOR complexed with morphine derivatives have been determined; however, structural basis of mOR activation by fentanyl-like synthetic opioids remains lacking. Exploiting the X-ray structure of mOR bound to a morphinan ligand and several state-of-the-art simulation techniques, including weighted ensemble and continuous constant pH molecular dynamics, we elucidated the detailed binding mechanism of fentanyl with mOR. Surprisingly, in addition to the orthosteric site common to morphinan opiates, fentanyl can move deeper and bind mOR through hydrogen bonding with a conserved histidine H297, which has been shown to modulate mOR’s ligand affinity and pH dependence in mutagenesis experiments, but its precise role remains unclear. Intriguingly, the secondary binding mode is only accessible when H297 adopts a neutral HID tautomer. Alternative binding modes and involvement of tautomer states may represent general mechanisms in G protein-coupled receptor (GPCR)-ligand recognition. Our work provides a starting point for understanding mOR activation by fentanyl analogs that are emerging at a rapid pace and assisting the design of safer analgesics to combat the opioid crisis. Current protein simulation studies employ standard protonation and tautomer states; our work demonstrates the need to move beyond the practice to advance our understanding of protein-ligand recognition.

database. S7 The CHARMM-GUI web server S8 was then used to construct the system of mOR embedded in a POPC lipid bilayer. The disulfide bond observed in the crystal structure was imposed between Cys140 and Cys217. All titratable residues were fixed in the standard protonation states. All histidine residues were set to the neutral N δ tautomer (HID) as in the default setting of CHARMM. S9 The system was first energy minimized using the steepest descent followed by conjugate gradient algorithm. The system was then gradually heated from 0 to 310 K over 200 ps with harmonic restraints on the protein heavy atoms, lipids, and bound water molecules (same as Step 1 in Table S1). Following heating, the system was equilibrated for 117 ns, during which time the various harmonic restrains were gradually reduced to zero ( Table S1). The final dimension of the system was ∼87×87×113 Å 3 . The final snapshot was used for constructing the fentanyl-bound mOR model and as the starting structure for the replica-exchange CpHMD simulations of apo mOR (CpH-apo).
Relaxation of the docked fentanyl-mOR complex in a lipid bilayer. The fentanylbound mOR model was prepared by superimposing a top fentanyl binding pose from a previous docking study S10 onto the final snapshot from the aforementioned equilibration simulation of apo active mOR. The docked pose showed a salt bridge between the piperidine nitrogen and Asp147 and an aromatic stacking between the phenethyl ring and His297. The docked model was equilibrated for 115 ns, using harmonic restraints imposed on the protein and fentanyl heavy atoms in the first 5 ns (details see Table S2).
During the unrestrained part of the simulation, the root-mean-square deviation (RMSD) of the fentanyl heavy atom positions with respect to the docked structure steadily increases in the first 20 ns and stabilizes at about 6 Å in the remainder of the 110 ns simulation (Fig.   S1A). Concomitant with the fentanyl RMSD increase, the minimum distance between fentanyl's piperidine nitrogen and His297's imidazole nitrogen decreases from 10 to about 7 Å (Fig. S1D); however, the salt bridge between the positively charged piperidine amine S-2 and the negatively charged Asp147 remains largely stable except for occasional excursions ( Fig. S1C), consistent with a previous µs simulation study. S11 The final snapshot was used as the starting configuration for the WE simulations.
Weighted ensemble MD simulations. Weighted ensemble (WE) is a path sampling protocol that uses splitting and merging trajectories to enhance sampling of rare events.
Briefly, the configuration space is divided into bins based on a predetermined progress coordinate and a fixed number of walkers (trajectories) per bin is targeted. At the beginning of the simulation, walkers are initiated from a single bin and after a specified time interval, resampling is performed by evaluating the number of walkers per bin, and for bins with less than desired number of walkers, the walker is replicated (or split), and for bins with more than desired number of walkers, the walkers are pruned. Thus, over time, more bins are sampled and the simulation progresses along the progress coordinate. Details theory and algorithm can be found elsewhere. S12-S14 Two WE MD simulations were carried out starting from the equilibrated fentanyl-mOR complex structure, with His297 fixed in either HIE (WE-HIE) or HID (WE-HID) state. The Python-based tool WESTPA S13 was used to control the WE protocol and data storage.

Continuous constant pH molecular dynamics (CpHMD)
We applied the membraneenabled hybrid-solvent continuous constant pH molecular dynamics (CpHMD) method S15,S16 to determine the protonation states of all titratable residues in the mOR systems. In this method, conformational dynamics is propagated in explicit solvent and lipids, while the solvation forces for propagating titration coordinates are calculated using the membrane generalized-Born GBSW model S17 based on the conformations sampled in explicit solvent. To accelerate convergence of the coupled conformational and protonationstate sampling, a replica-exchange protocol in the pH space is used. S15 The membraneenabled hybrid-solvent CpHMD method S15,S16 has been validated for pK a calculations and pH-dependent simulations of transmembrane proteins. S16,S18,S19 The detailed protocols can be found here S20 Three sets of replica-exchange CpHMD simulations were performed, starting from the equilibrated apo active mOR structure, equilibrated fentanyl-mOR complex in the D147binding mode, and the fentanyl-mOR complex in the H297-binding mode obtained from the WE-HID simulation in which the FEN-H297 distance is less than 3.5 Å. The pH replica-exchange protocol included 16 replicas in the pH range 2.5-9.5 with an increment of 0.5 unit. A GB calculation was invoked every 10 MD steps to update the titration coordinates. In the GB calculation, the default settings were used, consistent with our S-4 previous work. S15 Each pH replica underwent molecular dynamics in the NPT ensemble with an aggregate sampling time of 320 ns. All Asp, Glu, and His sidechains as well as the piperidine amine were allowed to titrate. The model pK a 's of Asp, Glu, and His are 3.8, 4.2, and 6.5, respectively, S21 while that of the piperidine amine in fentanyl is 8.9. S22 Molecular dynamics protocol. The temperature and pressure were maintained at 310 K and 1 atm by the Langevin thermostat and Monte Carlo barostat, respectively in the simulations with the Amber program, S1 while the modified Hoover thermostat and Langevin piston coupling method were used in the simulations with the CHARMM program. S9 Longrange electrostatics was treated by the particle-mesh Ewald (PME) method S23 with a realspace cut-off of 12 Å and a sixth-order interpolation with a 1.6-Å −1 grid spacing. The van der Waals interactions were smoothly switched to zero between 10 and 12 Å. Bonds involving hydrogen atoms were constrained using the SHAKE algorithm to enable a 2-fs timestep. Analysis was performed using CPPTRAJ program. S24 For analysis of the equilibrium simulations starting from the D147-and H297-binding modes, the last 200 or 400 ns data was used, respectively.
Analysis. Tanimoto coefficient between profiles A and B is calculated using the following equation: where x iA or x iB denotes the contact value between fentanyl and residue i of mOR from profile A or B, respectively. x i has a value of 1 if a contact exists and 0 otherwise.
S-5 NPT 50.0 0 0 0 a Force constant in the harmonic restraint potential for the backbone heavy atoms in kcal/mol/Å 2 . b Force constant in the position/dihedral restraint potential. The unit is kcal/mol/Å 2 for the position restraints or kcal/mol for dihedral angle restraints. mOR Wat referes to the crystal waters and seven additional water added using the program DOWSER. S6 Table S2: Steps in the equilibration simulation to relax the docked fentanyl-mOR complex structure Step Ensemble Time (   Only a lower bound to the pK a of the piperidine amine of fentanyl is given, as the amine remains charged in the entire simulation pH range 2.5-9.5. S-7   , and His297 (magenta) from the simulations with HIE297 (left) and HIP297 (right). C. Time series of the fentanyl rotation angle (blue) and vertical orientation angle (gold) from the simulations with HIE297 (left) and HIP297 (right). Fentanyl rotation angle is measured by the angle between N2-H vector and the y-axis. Fentanyl vertical orientation angle is measured by the angle between the N2-C4 vector and the membrane normal (z-axis). S-19