Coarse-Grained MD Simulations of Opioid Interactions with the µ -Opioid Receptor and the Surrounding Lipid Membrane

: In our previous studies, a new opioid (NFEPP) was developed to only selectively bind to the µ -opoid receptor (MOR) in inﬂamed tissue and thus avoid the severe side effects of fentanyl. We know that NFEPP has a reduced binding afﬁnity to MOR in healthy tissue. Inspired by the modelling and simulations performed by Sutcliffe et al., we present our own results of coarse-grained molecular dynamics simulations of fentanyl and NFEPP with regards to their interaction with the µ -opioid receptor embedded within the lipid cell membrane. For technical reasons, we have slightly modiﬁed Sutcliffe’s parametrisation of opioids. The pH-dependent opioid simulations are of interest because while fentanyl is protonated at the physiological pH, NFEPP is deprotonated due to its lower p K a value than that of fentanyl. Here, we analyse for the ﬁrst time whether pH changes have an effect on the dynamical behaviour of NFEPP when it is inside the cell membrane. Besides these changes, our analysis shows a possible alternative interaction of NFEPP at pH 7.4 outside the binding region of the MOR. The interaction potential of NFEPP with MOR is also depicted by analysing the provided statistical molecular dynamics simulations with the aid of an eigenvector analysis of a transition rate matrix. In our modelling, we see differences in the XY-diffusion proﬁles of NFEPP compared with fentanyl in the cell membrane.


Introduction
Clinically used opioids, such as morphine or fentanyl, are powerful pain killers, but they also have unwanted and even lethal side effects. The suppression of pain is mediated by these opioid ligands binding to and activating µ-opioid receptors (MOR), which are located on central and peripheral neurons. One possible way to reduce side effects is ensuring that the opioids only activate MOR in the periphery and not centrally (in the brain). Here, one can take advantage of the fact that pain usually goes hand in hand with inflammation in injured tissue. Inflammation lowers the pH value and leads to an acidic chemical environment.
Recent work by our group [1] led to the development of the novel analgesic compound N-(3-fluoro-1-phenethylpiperidin-4-yl)-N-phenylpropionamide (NFEPP), which activates MOR preferentially at acidic extracellular pH levels, as found in injured tissues [2]. Apparently, the low pK a of NFEPP (6.82 [1,3]) precludes its protonation in healthy tissues (pH >7.35; e.g., in the brain) but facilitates its protonation in injured/inflamed microenvironments (low pH, high concentrations of protons). Consistent with the notion that the protonation of opioid ligands is typically required for their binding to opioid receptors [4], NFEPP is active only in injured tissues, while at physiological pH in healthy tissues (e.g., brain), NFEPP deprotonates and loses its ability to bind to MOR. In contrast, the closely related conventional opioid fentanyl (pK a = 8.43 [3,5]) is protonated and able to activate MOR at both low and normal pH, and therefore produces both analgesic effects and central side effects.
Aside from pH, other inflammatory mediators also play important roles. For example, proinflammatory peptides and reactive oxygen species (radicals) can be associated with the formation of disulfide bonds (DSBs), which may modulate the function of G-protein coupled receptors (GPCRs) [6][7][8][9]. In addition, the micropharmacokinetics of ligands penetrating into neuronal membranes should be considered.
A recent study [10] examined the interaction of fentanyl with MOR embedded in a lipid bilayer using coarse-grained (CG) molecular dynamics (MD) simulations. With fully atomistic (FA) simulations, time scales of milliseconds or higher are usually required to observe rare events such as ligand binding. CG force fields combine groups of atoms into interaction sites or beads and thereby diminish the computational costs compared with simulations with FA force fields [11]. Hence, CG MD simulations were employed to deal with such sampling issues by attaining significant computational acceleration relative to FA MD simulations [10,11].
It was earlier noted that fentanyl demonstrates high lipid solubility, and thus, has a tendency to partition into the lipid bilayer without passing through the entire membrane [10]. In the present work, we tweaked the ligand modelling as per the FA atom to CG bead mapping MARTINI 2.0 guidelines [12]. Therefore, the Nd and Qd beads used in the previous study have been substituted with N0 and Q0 beads, respectively. Here, the N (nonpolar) and Q (charged) bead types correspond to particular beads used in the CG representation of deprotonated and protonated ligands, respectively. Furthermore, the suffixes d and 0 denote donor and none hydrogen-bonding abilities, respectively.
Such a modification of the ligand mapping would also alter the intramolecular, as well as intermolecular, interaction of these ligands in the CG system models. Hence, incorporating this modified modelling, the role of pH and the presence of an additional DSB within the MOR were investigated by modelling different scenarios. Afterwards, the interaction of these ligands with the membrane-embedded MOR was studied with the help of CG MD simulations.

Fully Atomistic Molecular Dynamics Simulations
In the the FA and CG MD simulations, the difference between inflamed and healthy tissues was modelled by changes in pH and with the formation of a DSB [6]. The protonation states of the different amino acid residues of the MOR were determined using the PROPKA predictor tool [13,14] with respect to the system acidity values of pH 6.5 and pH 7.4, corresponding to the inflamed and healthy system states, respectively.
For molecular modelling, the mouse MOR structure was procured from the RCSB database (Protein Data Bank (PDB) [15]: 6DDF [16]). Interestingly, in the G i -stabilised active states of MOR, the TM6 helix of the MOR is displaced towards the TM7 helix [16]. Furthermore, a recent in silico study has identified the preferential TM6 and TM7 conformational changes in the course of fentanyl-induced MOR activation [17]. Hence, in the mouse MOR, CYS 292 6.47 of the TM6 helix and CYS 321 7.38 of the neighbouring TM7 helix were selected for the introduction of an additional DSB at pH 6.5 to model the effect of reactive oxygen species on the MOR within inflamed tissues [6][7][8][9]. Sulfur atoms of these two cysteine residues are at a distance of 0.987 nm in the native mouse MOR crystal structure [16]. Thus, the dynamics of the MOR were examined with and without an additional CYS 292 6.47 − CYS 321 7.38 DSB (for terminology, please refer to [18]) in the receptor at pH 6.5 and pH 7.4, respectively.
The receptors were inserted into the POPC (1-palmitoyl-2-oleoyl-sn-glycero-3phosphocholine) bilayer models using the CHARMM-GUI Membrane Builder [19]. Similarly to [3], MD simulations were performed with GROMACS 2019.6 [20], using the CHARMM36m force field for the ligands [21], receptor [22] and lipids [23]. As an explicit solvent, the CHARMM TIP3P water model [24] was used. Sodium and chloride counterions were added to neutralise the excess charge and obtain a salt concentration of 0.15 M. All simulation systems went through consecutive minimisation, equilibration and production runs, using the GROMACS scripts generated by CHARMM-GUI [19]. First, the systems were energy minimised with steepest descent algorithms, followed by six-step equilibration runs. The first two runs were performed in the NVT (constant particle number, volume, and temperature) ensemble and the remaining runs in the NPT (constant particle number, pressure, and temperature) ensemble. Restraint forces were applied to the receptor, lipids, and water molecules, and Z-axis positional restraints were placed on lipid atoms to restrict their motion along the X-Y plane. These restraints were progressively reduced during the equilibration process. Ultimately, unrestrained NPT production runs of 1 µs were performed, with periodic boundary conditions along all three orthonormal directions.
The particle mesh Ewald (PME) method [25] was employed to calculate long-range Coulombic interactions, with a 1.2 nm cut-off for real-space interactions. A force-switch function was implemented for the Lennard-Jones interactions, with a smooth cut-off from 1.0 to 1.2 nm. The temperature was maintained at 310 K using the Nosé-Hoover thermostat [26,27]. System pressure was kept at 1 bar with the Parrinello-Rahman barostat [28], using a semi-isotropic scheme where pressures along the X-Y-directions and the Z-direction were coupled separately. The coupling constant and compressibility of the barostat were set to 5 ps and 4.5 × 10 −5 bar, respectively. The LINCS algorithm [29] was used to constrain the covalent bonds between hydrogen and other heavy atoms, allowing a simulation time step of 2 fs. Production run trajectories were saved every 1 ns, and processed with GRO-MACS analysis tools to generate the required information. VMD software [30] was used for visualisation.
The protonated form of fentanyl, and the protonated and deprotonated forms of NFEPP [1], were sketched and parameterised using the CHARMM-GUI Ligand Reader & Modeler [31]. To parameterise these ligands in MARTINI, initially, 1 µs FA MD simulations of a single ligand in TIP3P water and 0.15 M NaCl were conducted using the CHARMM36m force field [10]. System pressure was maintained at 1 bar with the Parrinello-Rahman barostat using an isotropic scheme, where pressures along all three orthonormal directions were coupled together. The rest of the simulation parameters were kept similar to those described earlier for the MD simulations of MOR embedded in POPC bilayer.

Coarse-Grained Molecular Dynamics Simulations
The protein structure coordinates of the mouse MOR structure (PDB: 6DDF) were converted to CG MARTINI 2.2 representation [32] using Martinize2 [12,[32][33][34][35]. An additional CYS 292 6.47 − CYS 321 7.38 DSB was incorporated into the MOR at pH 6.5, as mentioned earlier for the FA simulations. The ElNeDyn elastic network [35] with a spring force constant of 500 kJ · mol −1 · nm −2 and a cutoff of 0.9 nm was applied to the CG protein structure [36]. The FA ligand trajectories were extracted from the single ligand simulations. Atom-to-bead mapping of the ligands was performed with the CG Builder tool [37], and as shown in Figure 1a,b, each atom in the FA representation of the ligand was assigned to an appropriate CG bead [12]. The CG ligand parametrisation was performed using the Bartender program [38][39][40][41], with the FA ligand trajectories supplied as reference. A temperature value of 310 K and a system dielectric value of 74.21 [42] were provided during the parametrisation.
The CG MOR was then embedded in a membrane bilayer with a lipid composition, as listed in Table 1 [12,36,[43][44][45][46][47][48][49][50], using the insane script [44]. The average areas per phospholipid were set to 0.59 nm 2 and 0.55 nm 2 for the inner and the outer leaflets, respectively [51]. The initial size of the system box was maintained as 10.0 × 10.0 × 12.5 nm 3 . The MORbilayer complex was afterwards solvated in polarisable water [34] with 0.15 M NaCl. For the simulations with opioid ligands, 4 molecules of the ligand were then added to the system using the PACKMOL package [52]. In Figure 1c, a snapshot of the initial state of a CG system setup before the MD simulation is presented, consisting of four deprotonated NFEPP molecules along with the MOR at pH 7.4. The CG MD simulations were performed using the MARTINI 2.3P force field [53]. The systems were first energy minimised over 50,000 steps using the steepest descent algorithm, then equilibrated with NVT ensemble and then NPT ensembles for a total duration of 10 ns. Constraints were applied to the MOR and were progressively reduced during the equilibration process. Subsequently, production simulations with a 20 fs timestep were performed in NPT ensembles for 5 µs. The FA ligands were transformed into CG ligands using MARTINI 2.0 mapping scheme [12], as depicted for the deprotonated NFEPP. The B4 CG bead is Q0 for protonated fentanyl and protonated NFEPP molecules. The B5 CG bead is C2 for protonated fentanyl molecules. Other CG beads in protonated fentanyl and protonated NFEPP molecules are the same as those for deprotonated NFEPP molecules. (c) Initial CG simulation setup of a membrane-embedded MOR at pH 7.4, surrounded by solvent and four deprotonated NFEPP molecules shown in black, red, grey and violet.
The temperature was set to 310 K and maintained using the V-rescale thermostat [54]. The pressure was set to 1 bar and maintained using the Parrinello-Rahman barostat. Semi-isotropic pressure coupling was applied for these systems, with pressure along the X-Y-directions and the Z-direction coupled separately. The barostat coupling constant and compressibility were set to 12 ps and 3 × 10 −4 bar, respectively. For calculation of nonbonded interactions, a pair list was generated using the Verlet scheme with a buffer tolerance of 0.005. The Coulombic terms were calculated using PME and a real-space cutoff of 1.1 nm [55]. The relative dielectric constant was set to 2.5, as recommended for the polarizable water model [34]. For the VDW terms, a cutoff scheme with a cutoff value of 1.1 nm was employed along with the Verlet cutoff scheme for the potential-shift. The LINCS algorithm was implemented to constrain the bonds to their equilibration values. The simulation trajectories were analysed using the GROMACS analysis tools and MATLAB software [56]. The trajectory visualisation was performed in VMD.

Square Root Approximation
The trajectories generated a sequence of states. Due to the large amount of data and the high-dimensional structure of the state space, besides a direct analysis of the states, we also aim to carry out an analysis of the distribution of opioids in the membrane. Therefore, we aim to describe the ligand transitions inside the membrane by a rate matrix Q, giving the kinetics between collections of dynamically similar states. The states generated by the trajectory are Boltzmann distributed, and thus, dynamically similar states are adjacent. For assembling the matrix Q, we decompose the state space into Voronoi cells (ω i ) i=1,...,n such that each of the Voronoi cells has at least one of the states in it. In the next step, we count the number of states w i in each Voronoi cell ω i . For neighbouring cells, the entries Q ij of the rate matrix Q are then computed by This method, which is based on a cascade of approximations of the actual transition rates, is called square-root-approximation (SqrA) [57][58][59].

Results
The distance between the TM6 and TM7 helices of the membrane-bound MOR was tracked across the FA and CG simulations without any ligand molecule. In Figure 2c,d, it could be observed that the distances between the helices fluctuate by up to 0.2 nm over the simulation period. However, the pH 6.5 systems with an additional CYS 292 6.47 − CYS 321 7.38 DSB (Figure 2b) had the helices in a greater proximity compared with pH 7.4 systems without an additional DSB (Figure 2a). Significantly, this trend was found to be similar in both FA and CG simulations. Thus, the CG models were able to exhibit similar characteristics, as noted in the corresponding FA simulations. Hence, to access faster time scales, CG models of membraneembedded MOR were analysed with MD simulations. Furthermore, the interaction of ligand molecules with the surrounding environment was also studied using CG modelling.
The spatial distributions of ligand molecules with regard to the receptor are depicted in Figure 3. Each cross corresponds to a ligand position in consecutive time steps of the MD simulation. In Figure 3a,b, we map the XY-distribution at pH 6.5 compared with pH 7.4. The Z direction histogram is shown in Figure 3c. In contrast to the opioids which are activating the MOR, the non-active NFEPP at pH 7.4 enters the membrane more deeply (blue line). This results in less availability of NFEPP at the membrane surface at pH 7.4. After examining the data presented in Figure 3, we decided to explore the interaction between the NFEPP molecules and the MOR at pH 7 in greater detail. Several distances were calculated between these two groups to elicit further information, as depicted in Figure 4. From the distances across the XY plane between the NFEPP molecules and the MOR in Figure 4a, no significant trends could be observed. Therefore, we decided to analyse the movement of ligands along the Z axis, to check their movement across the simulation system. For this purpose, the Z axis values of the centre of mass coordinates were subtracted from the commensurate values of the NFEPP molecules. Although not much could be inferred from these calculations for three out of the four ligands shown in Figure 4b, the remaining ligand, henceforth referred to as "Ligand 1", was able to transport to a region below the centre of mass of the MOR. As the central portion of the MOR is embedded within the lipid bilayer, this indicated the entry of Ligand 1 into the bilayer region. Encouraged by this observation, we decided to calculate the minimum interatomic distances between the NFEPP molecules and the MOR. Based on the data from Figure 4c, we decided to calculate the radial distribution functions between Ligand 1 and the different residue groups of the MOR for the last 250 ns of the corresponding simulation trajectory. During this time frame, not only is Ligand 1 in close proximity to the MOR, but it is also present within the membrane bilayer region of the system.  From the trajectory snapshot depicting the interaction of NFEPP molecules with the MOR at pH 7 in Figure 5a, it was observed that, apart from Ligand 1, other ligands were in the vicinity of the upper leaflet of the membrane bilayer. Only Ligand 1, which was able to cross the solvent-membrane interface as previously noted in Figure 4b, appears to be closer to the lower leaflet than the upper leaflet of the bilayer. Afterwards, the same trajectory frame is displayed in Figure 5b, with the amino acid residues on the MOR surface categorised into four types based on their steric and Coulombic properties: acidic (negatively charged: Asp, Glu), basic (positively charged: Arg, Lys), nonpolar (neutral: Ala, Gly, Ile, Leu, Met, Phe, Pro, Trp and Val) and polar (neutral: Asn, Cys, Gln, Ser, Thr, Tyr, deprotonated His and protonated Asp) [30,60]. From this trajectory snapshot, it could be observed that deprotonated NFEPP prefers interaction with both the nonpolar and polar amino acid residues on the membrane-embedded region of the MOR surface. It shows relatively diminished interaction with the acidic and basic surface residues, primarily located at the solvent-accessible extremities of the MOR. Similar trends could be observed in Figure 5c, where radial distribution functions were calculated between Ligand 1 and these four residue groups. Significantly, once inside the membrane bilayer, the deprotonated NFEPP was found to demonstrate a stable interaction with the amphiphilic MOR surface, as opposed to a lipophilic interaction with the lipid bilayer tails. Furthermore, upon analyzing the interaction between Ligand 1 and the MOR in greater detail, it was observed that interaction with the residue PHE 156 3.41 in the TM3 helix of membrane-embedded MOR surface played an important role in stabilizing the interaction between Ligand 1 and MOR, as depicted in Figure 5d. It is not easy to quantify the differences between the distributions of the ligands from Figure 3a,b. In order to make these differences more visible, we divided the XYdirection of the membrane into equally sized cells and used SqrA in order to derive the diffusion profile of the ligands by generating a rate matrix Q of the process; see Figure 6. The term "diffusion profile" refers to the fact that the SqrA is basically a discretised Fokker-Planck diffusion operator [61]. For the SqrA, the correspondence between eigenvalues and time scales has been discussed [62]. In this sense, the graphs Figure 6a-d represent the χ-function of the slowest process, i.e., the gradient of this function points in the direction of the slowest processes. Compared with the diffusion profile of a ligand, which does not interact with the receptor Figure 6e, the χ-function of the NFEPP molecule at pH 6.5 and at pH 7.4 has a very different profile. It has two plateaus (metastable states) instead of four. The fentanyl molecule shows a kind of χ-function that is more like a free diffusion. This too shows that NFEPP has a higher tendency to interact with the MOR when moving inside the membrane. The second highest eigenvalue of the diffusion rate matrix corresponds to the characteristic time scale of the movement [62]. Again in Figure 6f, the biggest difference is between NFEPP at pH 6.5 and NFEPP at pH 7.4. Thus, the protonation state of NFEPP seems to have a substantial impact on its dynamical behaviour within the membrane.

Discussion
As coarse-graining combines three to four, or sometimes more, heavy atoms and their associated hydrogen atoms into a single bead [12], the fentanyl CG models employed in this study could also be applicable to opioid molecules with a chemical structure similar to fentanyl [10]. The same could be said for NFEPP and the halogenated fentanyl derivatives with a structure close to NFEPP. Furthermore, an improved CG modelling scheme and force field MARTINI 3.0 is currently available, with reportedly better representation of intrapeptide and interpeptide interactions along with protein-lipid interactions [63]. Additionally, an enhanced scheme for modelling small molecules, such as fentanyl and NFEPP, has been recently developed for this force field [64]. However, at present, the MARTINI 3.0 parameters for Chol, DPG3, DPSM and POP2 are not available [10,36,65]. Therefore, it was decided to use the MARTINI 2.2 and MARTINI 2.0 modelling schematics for the MOR and ligands, respectively, along with the MARTINI 2.3P force field.
The present study aimed to shed some light on the generic interaction of fentanyl and NFEPP opioid molecules with MOR and the lipid bilayer as a function of the surrounding physicochemical microenvironment. To access greater time scales, CG simulations were performed. It was observed that charged ligand molecules were unable to penetrate into the lipid bilayer. However, the deprotonated or neutral NFEPP at pH 7 was able to enter the membrane after crossing the membrane head group-lipid tail interface. However, upon forming a stable contact with the amphiphilic region of the MOR surface embedded within the bilayer, the ligand molecule preferred interaction with the MOR surface, instead of returning back to the hydrophobic environment in the midst of the lipid tail segments. Such amphiphilic behaviour has been also observed for deprotonated fentanyl in a recent FA MD simulation study, where a fentanyl molecule, initially placed between the tail segments of lipid molecules of the bilayer, was able to reach the membrane-solvent surface and interact with the solvent water molecules by reorienting itself within the membrane bilayer [66]. Furthermore, similar association of a BPTU ligand from the solvent region to the membrane-embedded region of a G-protein coupled receptor via the membrane bilayer has also been observed [67].
The pH-dependent membrane absorption of ligand molecules has been reported earlier [68,69], where ligands demonstrated increased surface uptake with increasing extracellular pH values, especially at pH greater than the ligand's pK a values. However, more recent in vivo studies are required to further corroborate these findings. Significantly, fluorinated analogues of fentanyl with a similar pK a were found to exhibit equivalent behaviour in vivo [70,71]. However, NFEPP and other fluorinated analogues, with pK a values lower than the physiological pH of 7.4, were found to exhibit highly pH-dependent biological activity near their pK a values, similar to the behaviour demonstrated by NFEPP in the current study. Therefore, the localisation of deprotonated NFEPP molecules near the MOR surface after entering via the lipid bilayer could diminish its bio-availability in the extracellular environment, which might also reduce its likelihood to interact with the MOR binding region and trigger its subsequent activation [3].
It is already known that NFEPP and fentanyl at pH 7.4 behave differently with regard to receptor activation. It has also been shown that NFEPP at pH 7.4 and NFEPP at pH 6.5 have different activation profiles [1]. However, the differences between fentanyl and NFEPP at pH 6.5 have rarely been analysed using molecular simulation methods in the literature. Our novel observations suggest that, besides the binding affinities, there might be more general differences between fentanyl and NFEPP with regard to interaction with the cell membrane. The fluorine atom of NFEPP increases the amphiphilic character of fentanyl. The process of approaching and binding to the MOR could also be different for fentanyl and NFEPP at pH 6.5, which has not been taken into account in pre-clinical studies so far, but is visible by comparing Figure 6a   Data Availability Statement: Code and data are available at https://github.com/user8490/opioid.

Conflicts of Interest:
The authors declare no conflict of interest.