Structural Insights into σ1 Receptor Interactions with Opioid Ligands by Molecular Dynamics Simulations

Despite considerable advances over the past years in understanding the mechanisms of action and the role of the σ1 receptor, several questions regarding this receptor remain unanswered. This receptor has been identified as a useful target for the treatment of a diverse range of diseases, from various central nervous system disorders to cancer. The recently solved issue of the crystal structure of the σ1 receptor has made elucidating the structure–activity relationship feasible. The interaction of seven representative opioid ligands with the crystal structure of the σ1 receptor (PDB ID: 5HK1) was simulated for the first time using molecular dynamics (MD). Analysis of the MD trajectories has provided the receptor–ligand interaction fingerprints, combining information on the crucial receptor residues and frequency of the residue–ligand contacts. The contact frequencies and the contact maps suggest that for all studied ligands, the hydrophilic (hydrogen bonding) interactions with Glu172 are an important factor for the ligands’ affinities toward the σ1 receptor. However, the hydrophobic interactions with Tyr120, Val162, Leu105, and Ile124 also significantly contribute to the ligand–receptor interplay and, in particular, differentiate the action of the agonistic morphine from the antagonistic haloperidol.


Introduction
The growing number of crystal structures of the receptors has shed a strong light on their structural features, the arrangement of the transmembrane helices, ligand-binding interactions, and specificities. Very recently, two high-resolution crystal structures of the human transmembrane σ 1 receptor in complex with two ligands (a high-affinity antagonist PD144418 and either agonist or inverse agonist 4-IBP; PDB ID: 5HK1 and 5HK2, respectively) have been published for the first time by Schmidt et al. [1]. The amino acid sequence of the σ 1 receptor does not resemble any of the other mammalian proteins. To the best of our knowledge, up until now, no molecular dynamics (MD) simulations using the crystal σ 1 receptor structure has been published.
The crystal data show a triangular structure constituted of three associated units with a single transmembrane domain for each protomer (Figure 1) located at the N-terminus [1]. The three associated protomers, consisting of four α-helices each, form a flat triangle with a transmembrane domain at each corner. The receptor binding pocket is hydrophobic but the receptor-ligand interactions known structural data are available for the first time. The question that needs to be addressed is how is the σ 1 receptor binding site able to accommodate chemically diverse ligands?
The design and development of new potent and selective σ 1 ligands are important for medicinal chemistry. The pharmacophoric requirement for the selective and optimal ligand binding to the σ 1 receptor is the presence of the following: (i) the piperidine nitrogen atom, which matches the positive ionizable residue in the binding site, (ii) the amide oxygen atom, which is able to accept a hydrogen bond from a donor residue of the receptor, and (iii) the two aromatic rings, which can be engaged in hydrophobic interactions. The family of fentanyl derivatives satisfies these requirements and hence they could possibly have a significant binding affinity towards the σ 1 receptor.
The site-directed mutagenesis assays and the crystal structure of σ 1 show that Glu172 and Asp126 are essential for the ligand binding [1]. Glu172 is hydrogen bonded by the protonated ligand's amine group, while Asp126, instead of interacting with the ligand, forms an intramolecular hydrogen bond with Glu172. To figure out the structural details of ligand binding to the σ 1 receptor, we chose seven ligands. Historically, the σ receptors were classified as opioid receptors; thus, some of the chosen ligands are simultaneously typical µ-opioid receptor ligands: fentanyl (FENT), morphine (MORPH), haloperidol (HALO), (R,S)-pentazocine (PTZ), and (R,S)-phenazocine (PHZ). There are two reasons for such a choice of ligands: (i) Firstly, most evidence to date has dealt with the therapeutic possibility of the σ 1 receptor antagonists (i.e., haloperidol) for the treatment of neuropathic pain since opioid enhances the effectiveness of medical treatment [2,7]; (ii) Secondly, the σ 1 receptor belongs now to a set of known experimental structures [1]; therefore, it is poised for a reevaluation of the interaction of these "old" ligands toward a putative binding partner.
In this paper, we present and discuss the results of the molecular dynamics simulations of the σ 1 receptor-ligand system based on the σ 1 receptor crystal structure [1]. The aim is to sketch the mechanism of inducing distinct responses of the σ 1 receptor to the chosen seven ligands. The study includes the following: (i) Carrying out independent MD simulations for the top-ranked poses of each ligand; (ii) Examining the frequency of contacts and patterns between different ligand moieties and the receptor; ligand and different receptor residues; maps of different ligand moieties vs. different receptor residues; (ii) Determining the binding energies (∆G bind ) and decomposing them according to the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) approach to identify the most important contributions from polar and nonpolar contacts. Decomposition of ∆G bind refers to terms such as ionic bond (or salt bridge) and hydrogen bond where the latter contains also attractive van der Waals and hydrophobic interactions essential for stabilizing the structure; (iv) Searching for a relationship/trend between ∆G bind and hydrophobicity/hydrophilicity of the residues through the linear regressions based on the studied ligands being treated as a training set.

Method
To investigate the interactions between the σ 1 receptor and the ligands, an automatic, multi-stage protocol was designed.

Ligand Preparation
The ligands (see Scheme 1) were prepared in the following fashion. First, sketch structures were drawn in the JSME applet. JSME is a 2D sketcher; however, it allows the drawing of bonds extending above and below the drawing plane. 2D sketches were automatically converted to naive 3D models, which took into consideration only atomic hybridization, and submitted to the Automated Topology Builder (ATB) [35]. The ATB server returned the optimized 3D structures and corresponding GROMACS [36] topologies. The retrieved files were used as input for the MD simulations.

Receptor Preparation
The receptor structure was taken from the 5HK1 entry of the PDB database [1]. This deposition contains the σ 1 receptor homotrimer bound with three PD144418 ligand molecules, two lipid molecules, and sulfate ions. It is the higher resolution antagonist-bound crystal structure, one of two published by Schmidt et al. [1]. Chain B of the receptor was extracted only from the 5HK1 structure for further modeling. Additionally, the N-terminal helix (residues 1-36), which anchors the receptor in the membrane, was removed as it does not interact with either the ligand or its surroundings. The structure was refined using Modeller 9.18 [37] as a few residues were missing several atoms. The protein backbone was kept intact during this procedure. Figure 1 presents the σ 1 receptor ligand complex.

Molecular Docking of Ligands to the σ 1 Receptor
The ligands were docked in silico to the crystal structure of the σ 1 receptor using the AutoDock program [38]. The binding pocket, which was manually selected with the 'Grid Box' tool from AutoDockTools, ran on the crystallographic structure 5HK1 with its native ligand PD144418. The procedure consisted of two steps: (1) setting up a minimal box to completely enclose the ligand and (2) extending the box size by approximately 50% in all directions. Docking was performed with the default AutoDock settings as its main purpose was to prepare the initial configuration for the MD simulations. Three docking runs yielded 8-10 docking poses each out of which up to three were selected for further processing [39]. The selection criterion was the presence of the hydrogen bond between the ligand's protonated NH + moiety and the C=O group of the Glu172 receptor's residue. Such a bonding was observed in the crystallographic structure and is thought to be the key to the interaction with the σ 1 receptor ligands.

Molecular Dynamics Simulation of the Receptor-Ligand Complex
In order to refine complex structures obtained from the docking procedure for each pose, we ran a fully flexible 80 ns MD simulation using the GROMACS program [36]. Production simulation stages were preceded by the energy minimization and the NVT and NPT equilibrations. Simulations were conducted with explicit solvent molecules at T = 300 K and δt = 2 fs time step.

The Formal Analysis of the MD Simulations and Presentation
The intention of the MD analysis was to find key similarities and differences in the ligand binding patterns. To this aim, 3000 structural snapshots, separated by 20 ps intervals, were extracted from the terminal with 60 ns on each trajectory. The interaction "fingerprint" was obtained from IChem analysis [40], providing both information on the occurrence and type of the interaction. This analysis returned to us the averaged contact maps depicting "fingerprint" patterns of the ligands binding to the receptor. We kept the default IChem parameters [40]. The receptor response to the ligand was monitored by the root mean square deviations (RMSD) and root mean square fluctuations (RMSF), which were calculated based on the binding site residues ( Figure S1).

MM/PBSA Calculations for Binding Free Energy
To evaluate the strength of the interactions between the ligand and the σ 1 receptor, the binding energy (∆G bind ) was calculated within the framework of the molecular mechanics/ Poisson-Boltzmann surface area (MM/PBSA) calculation approach [41]. It is a validated methodology combining the molecular mechanics energy and implicit solvation models. The enthalpic contribution to the free energy reflects the strength and specificity of the interactions between both partners, and includes ionic, electrostatic, hydrogen and halogen bonds, weak C-H . . . π hydrogen bond and the interactions with the prevailing dispersion components: π-π stacking, van der Waals and hydrophobic interactions [42,43]. Our main aim has been to compare the differences in the ligands' binding to the σ 1 receptor; therefore, we did not include the calculation of the entropic terms ∆S and refer only to the relative binding energies rather than to the absolute binding energies [44]. Moreover, as the binding energy is highly sensitive to the input parameters, the MM/PBSA model is more relevant in only calculating the relative differences. In order to base only on statistically converged results, independent short MD simulations starting with initial different velocities rather than long single simulations were run.
The components of the electrostatic ∆E elect , the van der Waals ∆E vdW , the polar ∆E polarsol and the nonpolar ∆E nonpolarsol calculated as SASA [45] will be discussed in Section 3.3. The decomposition of the interaction energies on a per-atom basis was also calculated. All these components were averaged over a series of snapshots from the corresponding trajectories. From the terminal 60 ns MD trajectory, we extracted 100 snapshots every 0.6 ns from the production trajectories and ran the MM/PBSA protocol on them using the GROMACS package [45].

Simulation of the PD144418 Ligand
In order to confront our methodology with experimental data, we performed the docking procedure on the PD144418 ligand that we had extracted previously from the 5HK1 PDB crystallographic structure. Our main goal with this task was to check if the ligand known to form a stable complex with the σ 1 receptor would reside within the pocket throughout the simulation while maintaining crucial interactions. The pre-simulation procedure was exactly the same as for the previously examined seven ligands. The simulation was performed under the same conditions as before but for a longer time (100 ns). The RMSD between the best-docked pose and the crystallographic structure was 0.59 Å and the crucial interaction of Glu172 with the protonated nitrogen atom in the ligand was retained in 99% of the trajectory snapshots. These all lead us to the conclusion that our procedure is valid for providing atom-level insight into the interaction of the σ 1 receptor with fentanyl ligands.

Dynamical Properties of the σ 1 Receptor-Ligand Interactions
The binding pocket includes polar (p) and hydrophobic (h) residues: Asp126(p), Trp89(p), Met93(h), Leu95(h), Leu105(h), Leu182(h), Phe107(h), Ile124(h), Trp164(p), Tyr103(p), Ser117(p), His116(p). Actually, the binding pocket is quite hydrophobic and its interior is occluded from the solvent [1]. The path of the ligand entry, as well as its exit, remains unknown [1]. It seems most likely that the ligands enter into the binding pocket in dynamic fluctuations of the system. Throughout this paragraph, we focus on the direct contacts between the ligand functional groups and the residues. Figure 2 illustrates the final snapshots of the MD simulations obtained for the seven ligands. Key interactions in the final snapshots are described below. First of all, the intramolecular salt bridge with the COO − oxygen atom of Glu172 is present in all systems close to Asp126. In general, one can distinguish three groups of residues in the binding pocket: (i) The first group is localized around Asp126 and forms the intermolecular salt bridge between the protonated nitrogen atom (N-H + ) in the ligands and the carboxyl (COO − ) group of the receptor Glu172 residue. We include also the stable hydrogen bonds between the carboxyl C=O group of Thr181 with the OH groups of the ligands in this region. (ii) The second region is placed around the Tyr103 residue and plays a determining role in the π-π stacking interactions between Tyr103 and the aromatic/heteroaromatic rings of FENT, HALO, MORPH, PTZ, and PHZ. (iii) It is noticeable that the residues (Phe107, Trp164, Ile178, Val162, Leu105, Leu182, and Ala185) from the third group belong to the hydrophobic pocket. They yield the required van der Waals and hydrophobic interactions, arranging appropriately the hydrophobic moieties of the ligands. A more detailed description of the ligand-receptor interaction can be found in Table S1. The conformational variations of the receptor upon the complex formation are presented in terms of the RMSD of the atoms' positions with respect to their initial positions ( Figure S1). The plots show that after about 20 ns of an initial increase of fluctuations, the receptor reaches equilibrium. The binding modes of the σ 1 receptor-ligand are presented in Figure S2.

Ligand-Receptor Binding in Terms of the Frequency of Ligand-Receptor Contacts
The ligand-receptor interactions can be presented in a more quantitative way using the patterns of the frequency of contacts between: (i) the ligand moieties with the receptor; (ii) the receptor residues with the whole ligand; (iii) the ligand-residue contact maps averaged over the MD trajectory.
The initial step of analysis of the contacts consists of the choice of moieties in the ligands participating in the interaction with the receptor (for abbreviations of moieties in the ligands, see Figure 3).
The search for the most frequently formed contacts led us to analyze the ligand-protein contact maps (Figures 3 and 4). Figure 4A presents the normalized frequency of the receptor-ligand interactions, while Figure 4B on the right illustrates the frequency normalized by taking into account the hydrophobicity scale, published in [46]. The positive and negative signals correspond to the hydrophobic and hydrophilic residues, respectively. The first three ligands represent σ 1 agonist (FENT), σ 1 antagonist (HALO), and the ligand of reduced affinity to the σ 1 receptor (MORPH). These ligands, apart from polar interactions (permanent salt bridge with the carboxylic group of the Asp126 and hydrogen bonds) are characterized by the π-π and hydrophobic contacts.
In the case of FENT, the most frequent contacts are with hydrophobic moieties (Ar2, Et, Ar1, chair, Link) while the contacts with hydrophilic moieties (NH + and C=O) are rare (Figure 3a). In Figure 4Aa, the highest peaks are observed for five receptor residues: Leu105, Ile 124, Ile178, Leu186, and Glu172. Out of them, only Glu172 is hydrophilic. The other residues interact through the alkyl groups and significantly contribute to hydrophobic interactions. Several other residues that interact less frequently, such as Met93, Val162, and Leu199, additionally strengthen the hydrophobic interactions. Moreover, common interactions with Phe107, Phe184, and Phe200 introduce the π-π stacking interactions with the two fentanyl phenyl rings and increase the overall hydrophobic reception of this ligand. This interaction picture in which the dominant role is played by hydrophobic interactions that are strengthened by π-π stacking interactions and where hydrogen bonding plays a minor role would support the classification of fentanyl as a hydrophobic σ 1 ligand. This is clearly supported by the relevant interaction spectra shown in Figures 3 and 4A,B. Therefore, it may be attributed to the simplistic correlations between hydrophobicity and binding affinity. This allows us to perceive the importance of the hydrophobic character of amino acid chains seen in Figure 4.
In the case of HALO, the most frequent contacts are with the hydrophobic moieties (Ar1, Ar2, Link, F, chair, Cl) while the contacts with the hydrophilic moieties (NH + , C=O and OH) are less frequent (Figures 3b and 4A,B). The highest peaks are observed in Figure 4A,B for the following receptor residues: Val162, Leu105, Tyr103, Trp89, Ser117, Val190, and Thr181. Notice that this time interactions with the hydrophobic (Leu, Val, and Leu) and hydrophilic (Trp, Thr, Ser, and Tyr) residues are found, although the hydrophilic ones are not as frequent and are of lower intensity than the hydrophobic ones. Observe also that in the case of HALO, more residues are active than for FENT ( Figure 4A,B). Figure 3 illustrates the contact analysis for the co-crystallized ligand PD144418. Similar to HALO, the most frequent contacts are found with the hydrophobic moieties (Ar, chair).
A comparison between the frequency of interactions of the FENT and HALO moieties indicates that the hydrophobic moieties of both molecules are definitely more active than the hydrophilic ones. However, a comparison of the frequency of contacts with the σ 1 receptor residues displays quite different patterns. For FENT, the majority of interactions occur with the hydrophobic residues.
For HALO, contacts with the hydrophobic residues are frequent but the hydrophilic contacts are important as well (see Figure 4B). It may be important that only these two most strongly binding ligands exhibit interactions with residues with labels exceeding 190-that is, Phe196, Leu199, and Phe200 (FENT) and Leu199, Tyr201, Thr202, and Leu203 (HALO) ( Figure 4B).
MORPH yields yet another interaction pattern. However, the hydrophobic moieties (Link, Ar(T) and chair) are responsible for the majority of molecular contacts (Figure 3c), while the hydrophilic residues Tyr120, Glu172, the hydrophobic Ile124 and Val162, and stacking forming Phe184 dominate in the frequency of contacts with the receptor. Therefore, the largest variety of different types of interactions are found in the case of MORPH, yet these interactions are focused on a much smaller number of residues ( Figure 4A,B).
Four ligands, Rand S-PTZ plus Rand S-PHZ, confirm the experimental observations about the importance of the ligand's stereochemistry. According to Sguazzine et al. [47], the stereochemistry of the fenpropimorph is not important for the binding to the σ 1 receptor, while the results found by Rossi et al. [29] suggest that in the case of the (R,S)-RC-33 ligands, there are preferences toward the S-conformation when the OH group is present in the chirality center. Hence, it is interesting how the PTZ and PHZ chirality affects the ligand-protein interactions.
In the case of R-PTZ, the most frequent contacts are with the hydrophobic moieties (Link, Me, Ar(T), chair, and CH 2 ), while the contacts with the hydrophilic moieties (NH + and OH(T)) are less frequent (Figure 3d). The highest peaks in Figure 4A are observed for the following residues: Tyr120, Ile124, Ala185, Tyr103, Leu182, Leu105, Ser117, and Glu172. Notice that this time the interactions with the hydrophobic (Leu, Ile, and Ala) and hydrophilic (Tyr, Glu, Ser, and Tyr) residues are, roughly, equally frequent. Interestingly, in agreement with the mutagenesis studies [4], frequent contacts of Tyr120 with R-PTZ can stabilize the ligand-receptor interaction and increase the potency of the R-form. A quite similar pattern is obtained for the S-PTZ (Figure 3f): the Link, CH 2 , Ar(T) and Me moieties form the most important hydrophobic interactions. However, two pairs of polar-hydrophobic contacts-Ile124 with Glu172 and Tyr120 with Ala185-seem to be promoting the leading contacts with high intensity suggesting the importance of the S-PTZ ionic interactions ( Figure 4B).
In contrast to R-PTZ, the R-PHZ has hydrophilic contacts with Trp89, Leu105, Phe184, Glu172, Thr181, Trp164, and Tyr120 residues ( Figure 4A), while the contacts with Leu105, Ile124, Phe184, and Ala185 residues have a hydrophobic character. However, the latter interactions are less frequent and play a minor role. In the case of S-PHZ ( Figure 4A), the hydrophobic contacts with Ile124, Ala185, Phe184, Phe107, and Leu186 have definitively higher frequencies. Nevertheless, less frequent hydrophilic constants with His154 and Glu172 have also been noticed ( Figure 4B). These findings suggest that Rand S-PHZ enantiomers display small but visible differences when binding to the σ 1 receptor. In conclusion, the chiral preferences of the σ 1 receptor seem to be clear yet they need further studies.   (A) (a-g) Normalized spectrum of the frequency of the contacts between receptor residues and the ligands; (B) On the left: Normalized spectrum of the frequency of the contacts between receptor residues and the ligands. On the right: Spectrum of the frequency of ligand contacts with receptor residues normalized according to the hydrophobic scale taken from Reference [46]. The positive bars correspond to hydrophobic contacts while negative bars correspond to hydrophilic contacts.

The Binding Energy Analysis
The mutagenesis and crystallographic experiments revealed that the two polar, anionic σ 1 residues Asp126 and Glu172 are obligatory for the ligand-receptor binding activity. In addition to those residues, several other hydrophobic residues (e.g., Ile128, Phe107, Lys105, Lys187, and Ala185) have been identified as playing a central role in the binding. Further insight into the origin of the ligand binding may be obtained from a decomposition of the total energy interactions (∆E tot , the enthalpic contribution) into the electrostatic (∆E elect ), the van der Waals term, as well as the interactions with polar (∆E polarsol ) and non-polar solvents (∆E nonpolarsol , within the SASA approximation). A comparison of the binding energy and its components for the series of the studied ligands within the MM/PBSA approach is presented in Table 1 and Figure 5.
First, it is noticeable that irrespective of the ligand's structure, the relative contribution of (a) the electrostatic term to the total attractive energy is ca. 75%; (b) the attractive van der Waals energy is ca. 22%, and (c) the only repulsive term, the interactions with the polar solvent, is also somewhat constant and is equal to ca. 60%. For all ligands, the attractive electrostatic term overcomes the repulsive term and is the main driving force for the ligand-receptor binding (Table 1, Figure 5a). In general, the attractive terms and even the electrostatic term alone overcome the repulsive one. Thus, it seems that the absolute attractive interaction contributions play a crucial role rather than the differences in the relative contributions of the attractive and repulsive forces. (A) Hence, the highest total binding energy is found for FENT, while the smallest is for S-PTZ. Moreover, even if the errors are taken into account and one considers only the low limits, the above conclusions remain true. (B) On the other hand, the values of total interaction energies of the most strongly interacting ligands-FENT, HALO, R-PHZ, and S-PHZ-are within the error bars of each other. The same holds true for the most weakly interacting ligands-S-PTZ, R-PTZ, and MORPH-for which the total interaction energies are also within the error bars of each other. (C) The high result for FENT is supported by both the highest electrostatic term and one of the highest van der Waals. FENT exhibits the largest attractive contribution and despite the fact that the repulsive term for it is also one of the two largest, the difference is still the biggest.
Notice that both the sum of the energy of the attractive interactions ∆E attr and the repulsive energy term ∆E polarsol linearly correlate with the number of contacts n cont (see Figure 5b). Thus, those two relatively independent computational measures of the ligand binding affinity are in close alignment. Let us now juxtapose the results of the calculations with the experimental characteristics, bearing in mind that those two refer to quite different conditions. We can remark that: (i) MORPH with the value of K i > 10000 nM [11] is endowed with the lowest affinity toward the σ 1 receptor among the series of studied compounds. This is consistent with both frequency of contacts and the most important terms of the interaction energy ( Figure 5b). (ii) Also, S-PHZ, which exhibits one of the highest K i in the studied series, 61.7 nM [23], is characterized by almost the smallest number of contacts and the smallest interaction energy terms in the presented calculations ( Figure 5b). (iii) On the other hand, the HALO ligand, which is the strongest σ 1 receptor antagonist considered here (K i = 0.65 nM [13]) exhibits the second highest number of contacts and the attractive energy term. (iv) Unfortunately, the computational predictions of the FENT binding energy were not successful.
For fentanyl, the experimental K i was estimated to be greater than 1000 nM [48,49], whereas in our calculations both the number of contacts and energy terms indicate FENT to interact very strongly ( Figure 5b). (v) The other studied ligands, R-PHZ, R-PTZ and S-PTZ, are located somewhere between the extreme interactions. However, they fit the predicted tendencies (Figure 5b).
At this moment we can state that there is in fact only one basic disagreement produced by the simulations with the FENT ligand. Although the experimental data can be thought to not be very exact, in the case of fentanyl the possibility of a mistake is very low. However, for now, we have no satisfactory explanation for the prediction of the very strong interaction between the FENT ligand and the σ 1 receptor. In fact, several factors can be responsible for such a disagreement. Still, most of these factors would refer to all groups of ligands. Nevertheless, out of these factors, the number of FENT contacts with the Leu105 residue is definitely the highest in the group of studied ligands (Figure 4Aa). The next highest number of such contacts is for HALO. Therefore, these contacts are probably overestimated for fentanyl.
Nevertheless, we think that the qualitative characteristics of the ligand interactions with the hydrophilic and hydrophobic residues of the receptor can differentiate the studied ligands in a good way. Indeed, among the set of studied ligands, both HALO and FENT share the most hydrophobic interactions with the receptor residues and they are essential for ligand binding as an antagonist. Our simulations confirm the important role of the interactions with the hydrophobic (Leu, Val and Leu) and hydrophilic (Trp, Thr, Ser, Glu and Tyr) residues, which are also identified in the experiment structure as crucial for the binding process [1]. In contrast, MORPH is the ligand without frequent non-polar interactions, which disqualifies it from playing the role of an antagonist. Instead, we found more frequent interactions with polar residues for PTZ and PHZ than antagonists such as HALO. Therefore, the important polar and non-polar contacts can be related to the agonistic and non-antagonistic properties of the ligands, respectively.

Conclusions
An interesting and intriguing aspect of the σ 1 receptor is that it can bind a wide spectrum of ligands of very different structural classes and diverse pharmacological applications, including, recently, a potential drug target for cancer therapy [7]. In contrast to the advances made for the family of G protein-coupled receptors (GPCRs) [50], the molecular mechanism governing the σ 1 receptor function remains obscure despite the profound insight obtained recently from the crystal structures [1]. Aiming to improve the understanding of the ligand binding of this receptor, we have performed, for the first time, MD simulations for seven ligand-receptor complexes (fentanyl, haloperidol, morphine, (R,S)-pentazocine, (R,S)-phenazocine) and additionally for PD144418 ligand, based on the published crystal structure of the σ 1 receptor (PDB ID: 5HK1). The combined approach (the MM/PBSA and the fingerprints analysis) unveiled the substantial role exerted by the polar and hydrophobic residue contacts in ligand binding so as to maintain the high affinity of the receptor. Ligand-protein interactions were categorized into two major types: hydrophilic and hydrophobic.
The results presented herein provide dynamic insight into the optimal binding of agonists and antagonists of the σ 1 receptor in atomic detail through computational modeling approaches. The altered interactions of agonists and antagonists with the σ 1 receptor have been known for a long time; however, the structural basis of the binding affinity was still elusive. The simulation results confirm the important role of particular protein residues-identified in the crystallographic σ 1 receptor-ligand data-which are crucial for the binding process. Analysis of the observed intermolecular contact frequencies together with the contact maps indicates that the hydrophilic interactions with Glu172, retained by all of the studied ligands, are a contributing factor to the affinity of the σ 1 receptor-ligand interaction. However, a high degree of contribution by hydrophobic interactions with Tyr120, Val162, Leu105, and Ile124 can clearly differentiate the studied agonist (MORPH) ligands from the antagonists (HALO). We anticipate that the results presented herein will pave the way for a deeper understanding of σ 1 receptor-ligand interactions, especially the fentanyl derivatives that have been a recent threat to public health in the form of illegal drugs.