Mechanistic insights into the allosteric modulation of opioid receptors by sodium ions.

The idea of sodium ions altering G-protein-coupled receptor (GPCR) ligand binding and signaling was first suggested for opioid receptors (ORs) in the 1970s and subsequently extended to other GPCRs. Recently published ultra-high-resolution crystal structures of GPCRs, including that of the δ-OR subtype, have started to shed light on the mechanism underlying sodium control in GPCR signaling by revealing details of the sodium binding site. Whether sodium accesses different receptor subtypes from the extra- or intracellular sides, following similar or different pathways, is still an open question. Earlier experiments in brain homogenates suggested a differential sodium regulation of ligand binding to the three major OR subtypes, in spite of their high degree of sequence similarity. Intrigued by this possibility, we explored the dynamic nature of sodium binding to δ-OR, μ-OR, and κ-OR by means of microsecond-scale, all-atom molecular dynamics (MD) simulations. Rapid sodium permeation was observed exclusively from the extracellular milieu, and following similar binding pathways in all three ligand-free OR systems, notwithstanding extra densities of sodium observed near nonconserved residues of κ-OR and δ-OR, but not in μ-OR. We speculate that these differences may be responsible for the differential increase in antagonist binding affinity of μ-OR by sodium resulting from specific ligand binding experiments in transfected cells. On the other hand, sodium reduced the level of binding of subtype-specific agonists to all OR subtypes. Additional biased and unbiased MD simulations were conducted using the δ-OR ultra-high-resolution crystal structure as a model system to provide a mechanistic explanation for this experimental observation.

Distinguished members of the G-protein-coupled receptor (GPCR) superfamily, opioid receptors (ORs) are the main targets for analgesics and play important roles in drug addiction. Because of pioneering studies using brain homogenates in the 1970s,[1−3] it has long been known that physiological concentrations of sodium decrease the level of binding of agonists, but not antagonists, to the μ-OR.[4] While similar allosteric effects were confirmed much later for several, albeit not all [e.g., the turkey β1-adrenergic receptor (B1AR)[5]], different family A GPCRs (see ref (6) for a recent review), the possibility that sodium differentially affects the binding of an agonist to the three major OR subtypes was also raised, with 65% agonist binding inhibition seen in μ-OR and δ-OR, but only 20% inhibition observed in κ-OR.[7] Moreover, a recent comparison of the effect of sodium, potassium, and lithium on δ-OR agonist binding suggested a differential modulation of δ-OR ligand binding parameters and G-protein coupling by monovalent ions, with sodium decreasing the level of δ-OR agonist binding more than the others.[8] Notably, treatment of membranes with reagents, particularly those attacking sulfhydryl groups, was shown to enhance the sodium effect.[9] On the other hand, divalent cations, and especially manganese ions (at 1 mM), were shown to almost restore full agonist binding in the presence of sodium at 100 mM in saturation studies,[10] while binding of antagonists remained unaffected.

At the molecular level, a possible explanation for the effect of ions on OR binding and signaling is that, like other molecules targeting allosteric sites, they affect the equilibrium between active and inactive states of the receptor, thus modulating the binding of native orthosteric ligands.[4,11] Mutagenesis studies in different GPCRs (e.g., see refs (12−18)) suggested a possible allosteric binding site for sodium, which involved a conserved aspartate in transmembrane helix 2 (TM2), namely D2.50 (the residue is labeled according to the Ballesteros–Weinstein generic numbering scheme,[19] which has been adopted throughout this work). Notably, pioneering molecular dynamics (MD) simulations of a model of the dopamine D2 receptor predicted a similar binding site for sodium ions diffusing freely from the extracellular side.[20] Similar conclusions were reached by a more recent MD study of the μ-OR, which also suggested the entry of sodium from the extracellular side, and binding to a site comprising residue D2.50.[21]

The first direct experimental evidence of binding of sodium to a GPCR came only very recently with the ultra-high-resolution crystallographic structure of the adenosine A2A receptor (A2AR; PDB entry 4EIY(22)), which was followed within just a few months by the very high (1.8–2.2 Å)-resolution crystallographic structures of the β1-adrenergic receptor (B1AR; PDB entry 4BVN(5)), protease-activated receptor 1 (PAR1; PDB entry 3VW7(23)), and the δ-OR (PDB entry 4N6H(24)). These structures revealed the precise location of the sodium ion (herein termed the “sodium crystallographic site” or “sodium allosteric site”), and its coordination by a salt bridge to D2.50, in addition to four polar interactions with receptor side chains and water molecules (S3.39, N3.35, and two water molecules in the case of δ-OR[24]).

To provide mechanistic details about the sodium control of GPCR binding and signaling, several MD simulations of the A2AR crystal structure with or without sodium at the allosteric site were recently conducted, and their results were interpreted in the context of radioligand binding, and thermostability experiments.[25] These studies further supported the idea that the binding of sodium and agonists is mutually exclusive by showing that the ion hampers possible activation-related conformational changes and rather gives preference to inactive conformations of the receptor. A similar conclusion was drawn from analysis of the agonist-bound crystallographic active-like state of the neurotensin receptor,[26] whose lack of sodium and collapsed sodium binding pocket implied that the high-affinity agonist binding and the presence of Na+ ions are mutually exclusive.

Although the lower-resolution crystal structures of μ-OR and κ-OR present chemically and conformationally conserved residues in their putative sodium allosteric sites compared to the ultra-high-resolution structure of δ-OR, the ion access from the bulk and binding to these sites may vary dramatically, because of differences in the number and location of negatively charged residues across the three receptor structures. Intrigued by the possibility that the mechanism of sodium control of OR binding and signaling may be different among the three receptors, we explored here the dynamic nature of sodium binding diffusing freely from the bulk to the interior of the ligand-free δ-OR, μ-OR, and κ-OR crystallographic conformations by means of unbiased, microsecond-scale MD simulations. The results of these simulations were coupled with those of radioligand binding assays in stably transfected cell lines to evaluate differences and similarities in the allosteric modulation of the three major OR subtypes by sodium ions. Insights from additional unbiased and biased MD simulations of antagonist-bound and agonist-bound, activated receptors using the δ-OR ultra-high-resolution crystal structure[24] as a model system further supported the mechanistic implications of both the ligand-free simulations and the experimental results.

Materials and Methods

Simulation System Setups and Force Field

The inactive crystal structures of mouse μ-OR, human κ-OR, and human δ-OR corresponding to PDB entries 4DKL,[27]4DJH,[28] and 4N6H,[24] respectively, were used as starting conformations for all-atom MD simulations in an explicit lipid–water environment. Crystallographic waters were kept if present, whereas other nonprotein atoms were deleted. Missing loops of μ-OR and κ-OR were modeled ab initio using ROSETTA,[29] whereas the missing C- and N-terminal segments were ignored. Missing side chains were inserted with the Swiss-Pdb Viewer.[30] To simulate the binding of free sodium to the ligand-free δ-OR crystal structure, the crystallographic sodium at D2.50 was removed at the start of the simulation. The δ-OR crystal structure was also simulated with the sodium bound at D2.50, with or without the antagonist naltrindole bound at the orthosteric site. Because no active OR crystal structure is available yet, all-atom, adiabatic biased molecular dynamics (ABMD) simulations[31] (see details in Molecular Dynamics Simulations) were employed to sample the conformational transition of the δ-OR TM bundle from its inactive crystal structure in the presence of a docked selective δ-OR agonist, i.e., (+)-4-{(αR)-α-[(2S,5R)-4-allyl-2,5-dimethyl-1-piperazinyl]-3-methoxybenzyl}-N,N-diethylbenzamide [SNC-80 (see the docking details in the following section)], instead of naltrindole, to the corresponding atomic coordinates of the TM region of the active β2-adrenergic receptor (B2AR) crystal structure in complex with the G-protein (PDB entry 3SN6(32)).

Each system was embedded into a pre-equilibrated patch of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and 10% cholesterol and solvated with 150 mM NaCl, corresponding to 0.00009 particle/Å3. All sodium ions were initially placed at a minimal distance of 3 Å from the protein, and additional chloride ions were added to neutralize the systems. Each complete system measured roughly 90 Å × 90 Å × 90 Å, consisting of the receptor, ∼20 cholesterol molecules, ∼200 POPC lipids, ∼11500 water molecules, ∼30 sodium ions, and ∼40 chloride ions, totaling close to ∼67000 atoms on average. The TIP3P water model and the CHARMM27 force field for the protein, lipids, and ions were used for the MD simulations. The CHARMM General Force Field[33] parameter set was used to simulate the ligands naltrindole and SNC-80.

Docking of SNC-80 to the Ultra-High-Resolution Crystal Structure of δ-OR

The system for docking of SNC-80 to the orthosteric binding pocket of the ultra-high-resolution crystal structure of δ-OR (PDB entry 4N6H(24)) instead of the crystallographic ligand naltrindole was prepared using Maestro within Schrödinger Suite 2014-1.[34] The SNC-80 molecule was built with the two-dimensional chemical sketcher tool and prepared with LigPrep at physiological pH 7.3, using Epik.[35] The ligand protonation state was assigned as suggested in the literature.[36] The receptor was prepared with the Protein Preparation Wizard tool to add hydrogens and assign bond orders. A grid box with dimensions of 26 Å × 26 Å × 30 Å with an inner box with dimensions of 10 Å × 10 Å × 14 Å centered on the crystallographic ligand naltrindole was used for flexible ligand docking calculations using Glide XP.[37] Ten poses of SNC-80 were selected for postdocking minimization, and the one exhibiting the best binding XP GScore (−5.860 kcal/mol), which incidentally overlapped nicely with the crystallographic ligand naltrindole, was chosen for further studies.

Molecular Dynamics Simulations

Both biased and unbiased MD simulations (see Table 1 for a summary) were conducted using the GROMACS 4.6 simulation package[38] in the NPT ensemble at 300 K and 1 bar, using a Nose-Hoover thermostat[39] and Parrinello–Rahman pressure coupling,[40] respectively. All bonds were constrained using the LINCS algorithm,[41] and a 10 Å cutoff was used for short-range nonbonded interactions. The receptors, lipids, and water and ions were defined as separate coupling groups to allow lateral diffusion. Prior to unbiased equilibration and production MD runs, energy minimizations were conducted using the steepest descent algorithm for 1000 steps on restrained heavy atoms, followed by 1000 additional steps on restrained Cα atoms. Equilibrations of receptors consisted of relaxations of 3 ns with decreasing positional restraints first on all the heavy atoms and crystal water and then on the Cα atoms only. Restraints were all removed during production runs of 1 μs and ten runs of 100 ns for each system, which used randomized atomic velocities according to the Maxwell distribution at 300 K.

Table 1

Summary of the MD Simulations Conducted on δ-OR, μ-OR, and κ-OR in the Presence of Physiological Concentrations of Sodium Chloride (150 mM)

system	simulation type	conformation	sodium at D2.50 in the starting structure	ligand	runs	simulation length (ns)	restraints
δ-OR	MD	inactive	no	none	1	1000	no
μ-OR	MD	inactive	no	none	1	1000	no
κ-OR	MD	inactive	no	none	1	1000	no
δ-OR	MD	inactive	no	none	10	100	no
μ-OR	MD	inactive	no	none	10	100	no
κ-OR	MD	inactive	no	none	10	100	no
δ-OR	MD	inactive	yes	none	1	500	no
δ-OR	MD	inactive	yes	naltrindole	1	500	no
δ-OR	ABMD	inactive to active	yes	SNC-80	1	115	no
δ-OR	MD	active	yes	SNC-80	3	100	no
δ-OR	MD	active	yes	SNC-80	1	100	on TM Cα
δ-OR	RAMD	active	yes	SNC-80	159	0.1	no

Adiabatic biased MD (ABMD) simulations[31] were performed using the PLUMED 1.3 plugin[42] within GROMACS to sample the conformational transition of the TM bundle of the inactive crystal structure of δ-OR, with the selective δ-OR agonist SNC-80 replacing naltrindole at the orthosteric binding pocket, to the corresponding atomic coordinates of the TM region of the active B2AR crystal structure in complex with the G-protein.[32] After the latter (residues 30–59, 66–96, 103–136, 147–171, 198–225, 266–298, and 306–328 in PDB entry 3SN6(32)) had been superposed onto the TM Cα atoms of the δ-OR inactive structure (residues 46–75, 82–112, 118–151, 162–186, 212–239, 254–286, and 298–320 in PDB entry 4N6H(24)), the TM Cα atoms of δ-OR with the sodium ion at the allosteric, crystallographic site were biased toward the coordinates of 3SN6 using the mean-square deviation (MSD) between these TM Cα atoms as a reaction coordinate and an increasing elastic constant (1000 kJ/Å2 for 70 ns, 2000 kJ/Å2 for 30 ns, and 10000 kJ/Å2 for 15 ns). The snapshot obtained at the end of the 115 ns ABMD simulations (MSD value of <0.7 Å) was used as a representative activated conformation of δ-OR and was further simulated with or without positional restraints on the TM Cα atoms.

Random accelerated molecular dynamics (RAMD[43]) implemented in the NAMD simulation package[44] was used to sample possible egress pathways of the sodium ion from the activated SNC-80-bound δ-OR conformation produced by ABMD, with random forces acting on the sodium at the allosteric site. GROMACS format topology/structure was converted to NAMD inputs using CHARMM.[45] Prior to these unrestrained RAMD simulations, the system was minimized for 1000 steps and then equilibrated for 20 ps with decreasing root-mean-square deviation (RMSD) restraints on the protein heavy atoms (10 kcal mol–1 Å–1 for 10 ps and then 1 kcal mol–1 Å–1 for 10 ps). The RAMD simulations used a 2 fs time step, and a NPT ensemble maintained at 300 K with Langevin dynamics. RAMD-specific parameters were chosen on the basis of the results of 192 RAMD trial simulations of 50 ps each, which used combinations of different accelerations (0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64, and 1.28 kcal mol–1 g–1), RAMD time steps (10, 20, 40, and 80 steps), and threshold distances (0.001, 0.002, 0.004, 0.008, 0.016, and 0.032 Å). In the end, 159 RAMD production simulations of 100 ps each were run with a 0.64 kcal mol–1 g–1 acceleration (no egress was observed at lower values of the acceleration constant), time steps of 10, combinations of 0.001, 0.002, and 0.004 Å threshold distances, and 53 random numbers. Simulations were terminated as soon as the ion moved more than 30 Å from the center of mass of the protein.

Analysis of Simulations

The distribution of sodium ions during simulations was examined using the Volmap tool of VMD[46] on each individual 1 μs simulation, or concatenated batch simulations of 100 ns for each system conducted to verify statistical significance. Snapshots at every 1 ns interval were extracted and fit onto starting structures for each system, using Cα atoms.

Distance and dihedral measurements were obtained using either GROMACS or VMD tools.[46] Figures were rendered using Pymol.[47] Distance and angle measurements were plotted using the matplotlib package[48] in python.

RAMD trajectories were clustered on the basis of both the location of the ion at the end of each simulation and the putative egress pathway it follows. Clusters are labeled extracellular (EC), intracellular (IC), or membrane (M) depending on the location of the ion at the end of the simulation, with the specific egress pathway, i.e., the region through which they exit the bundle, specified in parentheses.

Membrane and Ligand Preparation

Chinese hamster ovary (CHO) cells stably transfected with mouse δ-OR (DOR-1), μ-OR (MOR-1), or κ-OR (KOR-1) were used for radioligand binding experiments. The membrane was prepared as follows: (1) rinse each plate with cold PBS buffer twice, (2) manually remove the cells adhered to the surface of the plates with cold PBS, (3) centrifuge at 3500 g for 5 min, (4) resuspend and homogenize the pellet for 15 s in tris buffer (50 mM Tris, 1 mM K+EDTA, and 100 mM NaCl) with 100 μM fresh PMSF (phenylmethanesulfonyl fluoride) at 4 °C, (5) incubate at 25 °C for 15 min, (6) centrifuge at 30000g for 45 min, (7) resuspend and homogenize the pellet for 15 s in 0.32 M sucrose at 4 °C, and (8) store in aliquots at −80 °C. Radiolabeled [3H]DAMGO [[3H]-[d-Ala2,NMe-Phe4,Gly-ol5]enkephalin, 46 Ci/mmol (Polypeptide Group, National Institute on Drug Abuse, Bethesda, MD)], [3H]U69,593 {(+)-(5α,7α,8β)-N-methyl-N-[7-(1-pyrrolidinyl)-1-oxaspiro[4.5]dec-8-yl]-benzeneacetamide, 43.6 Ci/mmol, PerkinElmer}, and [3H]DPDPE ([3H]-[d-Pen2,d-Pen5]rnkephalin, 46 Ci/mmol, Polypeptide Group, National Institute on Drug Abuse) were used as subtype-selective agonists for μ-OR, κ-OR, and δ-OR, respectively. Radiolabeled [3H]diprenorphine [17-(cyclopropylmethyl)-4,5-epoxy-18,19-dihydro-3-hydroxy-6-methoxy-α,α-dimethyl-6,14-ethenomorphinan-7-methanol, 57.4 Ci/mmol, PerkinElmer] was used as a nonselective antagonist for all receptor subtypes. Prior to ligand binding experiments, protein concentrations were determined by the Lowry method[49] as previously described, using bovine serum albumin (BSA) as the standard.

Radioligand Binding Experiments

All binding experiments were performed for 90 min in 25 mM Tris-HCl buffer, at pH 7.4 and 25 °C. The protein concentration was kept at approximately 150 μg/mL. For specific binding experiments, the specific binding was assessed by subtracting nonspecific binding from total binding. Levallorphan at 10 μM was used to measure nonspecific binding. Six sodium concentrations were examined: 5, 10, 25, 50, 100, and 150 nM. Each assay was repeated three times. For saturation experiments, six radioligand concentrations were examined for κ-OR and δ-OR: 0.125, 0.25, 0.5, 1, 2, and 4 nM. For μ-OR, 8 nM was also examined. Assays without sodium were repeated four times, and assays in the presence of 100 mM sodium were repeated three times. After incubation, the membranes were filtered through glass fiber filters and washed three times with ice-cold 5 mM Tris-HCl buffer on a Brandel cell harvester. Filters were then transferred into vials containing 3 mL of Liquiscent (National Diagnostics). The radioactivity in vials was determined using scintillation spectroscopy in a Packard TRI-CARB 2900TR counter. Kd and Bmax were calculated using nonlinear regression analysis (GraphPad Prism). Error bars were calculated as the standard error of the mean (SEM) of independent assays.

Results

The results of all the all-atom MD simulations we conducted on ligand-free crystallographic states of δ-OR, μ-OR, and κ-OR in the presence of physiological concentrations of sodium, as well as ligand-free or antagonist-bound, inactive δ-OR with sodium ion at the crystallographic site, and on an agonist-bound, activated model of δ-OR, are presented below in the various subsections. A summary of all simulations presented herein is provided in Table 1.

Preferred Receptor Locations Visited by Sodium

By revealing significant differences in the number and location of negatively charged residues, especially in the extracellular loop (EL) region, visual inspection of the crystal structures of δ-OR,[24] μ-OR,[27] and κ-OR[28] provided the rationale for studying whether sodium may follow different binding pathways in these receptors. Thus, we conducted all-atom, microsecond-scale MD simulations of ligand-free δ-OR, μ-OR, and κ-OR crystallographic structures (PDB entries 4N6H,[24]4DKL,[27] and 4DJH,[28] respectively) embedded in a hydrated POPC/10% cholesterol bilayer, and in the presence of physiological concentrations of sodium. Volumetric maps of sodium ions calculated using these microsecond trajectories and plotted onto representative structures from simulations of δ-OR, μ-OR, and κ-OR (panels A–C of Figure 1, respectively) show a common, highest sodium density at the D2.50 crystallographic sodium site. The representative conformations shown in Figure 1 correspond to the structures with the lowest average heavy atom RMSD from all other structures in the most populated cluster of conformations sampled during dynamics. Unlike μ-OR simulations, microsecond-scale simulations of the δ-OR and κ-OR revealed additional significant densities of the ion near the extracellular side. Specifically, these extra densities were seen at the nonconserved, negatively charged E6.58 (a W in δ-OR and a K in μ-OR) and EL2 D204 residues in κ-OR, as well as at the EL3 D293 and D288 residues in δ-OR. Notably, no significant sodium density was detected at the cytoplasmic side in any of the simulated ORs.

Figure 1

Vertical views of the representative structures from the microsecond simulations of (A) δ-OR (PDB entry 4N6H), (B) μ-OR (PDB entry 4DKL), and (C) κ-OR (PDB entry 4DJH) in cartoon representation with negatively charged (Asp and Glu) residues shown as sticks. Generic numbering for TM Asp/Glu residues and loop labels for other Asp/Glu residues are provided in parentheses. Sodium occupancy during the microsecond simulations at 0.12 and 0.24 particle/Å3 contour levels is shown in transparent and solid blue, respectively.

Sodium Binding Pathways

Sodium binding pathways from the bulk to the interior of the receptor were monitored during the microsecond simulations of δ-OR, μ-OR, and κ-OR and are illustrated in Figure 2A–C by connecting sodium positions at 1 ns intervals with lines colored by increasing simulation time using a blue–green–red color scale. It is indeed evident from this figure that at least one sodium ion accesses the receptor from the extracellular milieu, passes through the orthosteric binding pocket, establishing a direct interaction (distances of ≤2.5 Å) with residue D3.32 at the very beginning of the simulation (Figure 3A–C), and then accesses more or less rapidly [within the first ∼100 ns of simulation for δ-OR, ∼250 ns for μ-OR, and ∼20 ns for κ-OR (see Figure 3A–C)] the allosteric sodium site, forming a direct interaction with D2.50.

Figure 2

Sodium entry pathways for (A) δ-OR, (B) μ-OR, and (C) κ-OR microsecond simulations. The pathways are composed by connecting sodium positions at 1 ns intervals, colored by simulation time (blue to green to red). They are plotted on representative receptor simulation structures, i.e., the structures with the lowest average heavy atom RMSD from all other structures in the most populated cluster of conformations sampled during dynamics, which are shown in cartoon representation with negatively charged residues shown as sticks. Generic numbering for TM Asp/Glu residues and loop naming for non-TM Asp/Glu residues are included in parentheses.

Figure 3

Sodium coordination at the orthosteric ligand and allosteric sodium binding sites. Sodium minimum distance to side chain oxygens of D3.32 (at the orthosteric ligand site, gray) or D2.50 (at the allosteric sodium site, black) for (A) δ-OR, (B) μ-OR, and (C) κ-OR microsecond simulations. Representative conformational states are indicated as states 1 (extracted at 40, 145, and 6 ns for δ-OR, μ-OR, and κ-OR, respectively) and 2 (extracted at 768, 974, and 900 ns for δ-OR, μ-OR, and κ-OR, respectively), below the distance plots. Protein backbones from simulations are colored silver. Ion-interacting residues, bound sodium, and water from the simulation are shown as silver sticks, a blue sphere, and red spheres, respectively. In comparison, the sodium atom from the δ-OR crystal (PDB entry 4N6H) is shown as a cyan sphere, while protein and water molecules from the crystal are colored gray. Ion-interacting residues are labeled with generic numbering. Sodium coordination is indicated by black dashed lines.

To verify the statistical significance of our observation, i.e., that sodium ions access the interior of the receptor from the extracellular side and first stop at the D3.32 orthosteric ligand binding site before proceeding to the allosteric sodium site, we conducted ten 100 ns runs for each OR system and monitored the minimal distance of any sodium ion from both D3.32 and D2.50. As shown in Figure S1 of the Supporting Information, in general, ions access D3.32 before proceeding to D2.50. Notably, they do so quite rapidly in the κ-OR, compared to δ-OR and μ-OR. Moreover, unlike the δ-OR and μ-OR simulations, these relatively short κ-OR simulations show a few instances of simultaneous sodium occupancy of D3.32 and D2.50.

Sodium Coordination at the Orthosteric Ligand and Allosteric Sodium Binding Sites

The illustrations in Figure 3A–C show the coordination that sodium acquires at the orthosteric ligand or allosteric sodium binding sites during the δ-OR, μ-OR, and κ-OR simulations. Specifically, these figures depict representative conformational states, herein termed 1 and 2, extracted at 40 ns, 145 ns, and 6 ns or 768 ns, 974 ns, and 900 ns to show sodium coordination at the orthosteric ligand or allosteric sodium binding sites of δ-OR, μ-OR, and κ-OR, respectively. For comparison, these structures are overlapped onto the ultra-high-resolution crystal structure of δ-OR corresponding to PDB entry 4N6H. As shown in the illustrations of state 1 for all three simulated receptors, no receptor residue other than D3.32 coordinates sodium ions at the orthosteric ligand binding pocket. Indeed, the first coordination shell is made of a bidentate interaction of D3.32 and interactions with four different water molecules. As the simulations evolve, the sodium at D3.32 moves toward the allosteric sodium binding site where it engages in coordination with D2.50. During the course of the microsecond simulations, the ion is expected to acquire the coordination seen in the crystal structure of PDB entry 4N6H, which involves residues D2.50, S3.39, and N3.35 and two water molecules. Indeed, this occurs for all the OR systems, as shown by the state 2 illustrations in Figure 3A–C, as well as the distance plots in Figure S2 of the Supporting Information. However, while stable interactions between the ion and residues D2.50, S3.39, and N3.35 are established by the end of the μ-OR and κ-OR simulations, the interactions with S3.39 and N3.35 keep fluctuating during the microsecond simulations of δ-OR, with the exact crystal-like coordination of the ion at the allosteric sodium site last seen at around 780 ns (Figure 3A and Figure S2A of the Supporting Information). Figure S3 of the Supporting Information shows the residues in the first and second coordination shells of the sodium ion at the allosteric sodium site for each opioid receptor subtype. While D2.50 and S3.39 are still coordinating the ion at the end of the δ-OR microsecond simulation, the N3.35 side chain rotates away [χ2 changes from ∼50° to −25° (see Figure S2A of the Supporting Information)] to interact with L2.52 and A2.53. This conformational behavior of N3.35 is also noticed in the μ-OR simulation (see the χ2 plot in Figure S2B of the Supporting Information) but is not observed in the microsecond κ-OR simulations (see the χ2 plot in Figure S2C of the Supporting Information). Notably, the δ-OR A2.53 residue with which N3.35 engages in interaction is the same as in μ-OR but is substituted with V2.53 in κ-OR.

Sodium Coordination in Antagonist-Bound δ-OR Simulations

To test the hypothesis that the sodium coordination seen in the ultra-high-resolution crystal structure of δ-OR is stably maintained only in the presence of the ligand bound at the orthosteric pocket, we conducted additional δ-OR simulations with sodium bound at D2.50, with or without the antagonist naltrindole. The simulations, 500 ns each, show significant differences in the dynamic behavior of the N3.35 and S3.39 side chains, and consequently the coordination of the sodium ion, compared to that observed in the ligand-free δ-OR simulations. In the ligand-free δ-OR simulations with sodium bound at D2.50 as seen in the ultra-high-resolution crystal structure, sodium coordination by N3.35 is lost almost immediately (Figure S4A of the Supporting Information) as the residue χ2 dihedral changes from 50° to −25°. In contrast, in the presence of the antagonist naltrindole, sodium coordination remains the same as that seen in the ultra-high-resolution crystal structure for the entire duration of the simulation (Figure S4B of the Supporting Information). Notably, the H-bonding network that was found in the ultra-high-resolution crystal structure of δ-OR to link naltrindole to sodium is preserved in these simulations (Figure S5 of the Supporting Information). Part of this network is the charge interaction between the naltrindole amine and D3.32, which is also involved in stabilizing N3.35 through water-mediated hydrogen bonding.

As expected, in the simulations of δ-OR with sodium at the allosteric D2.50 site with or without naltrindole at the orthosteric site, the sodium second-coordination shell residue W6.48, which is also called the “toggle switch” in GPCR activation[3,50] because of its side chain dihedral change hypothesized earlier to occur upon activation, remains close to the conformation seen in the inactive crystal structure. The values of the W6.48 χ2 angle monitored during the ligand-free or antagonist-bound simulations, as well as the ligand-free microsecond simulations of δ-OR, μ-OR, and κ-OR, are reported in Figure S6 of the Supporting Information. As shown in this figure, the W6.48 χ2 angle is rather stable throughout all simulations, in spite of some fluctuations recorded in the κ-OR system.

Whereabouts of Sodium in a Fully Activated Receptor Conformation

The allosteric sodium pocket has been observed to collapse in a fully activated conformation of the receptor.[6] We sampled the conformational transition from the δ-OR inactive crystal structure, in the presence of sodium at the allosteric site and the δ-OR-selective agonist SNC-80 at the orthosteric ligand binding site, to an active-like form using the G-protein-bound active crystal structure of B2AR as a target structure with ABMD (see Materials and Methods for details). At the end of the ABMD simulation, the MSD of the TM Cα atoms reached a value of <0.7 Å2 with respect to the target structure, yet the sodium ion did not leave the allosteric pocket but rather kept its coordination with residues D2.50, N3.35, and S3.39. Notably, a bidentate coordination by D2.50 was observed at the end of ABMD simulation, resulting in the six coordination of the ion (Figure S7A of the Supporting Information) compared to the five coordination seen in the δ-OR ultra-high-resolution crystal structure (Figure S5 of the Supporting Information). The hydrogen bond network around sodium is otherwise similar to that observed in the δ-OR crystal structure and ligand-free simulations (Figure S5 of the Supporting Information). The structure obtained at the end of the ABMD simulation was then subjected to 100 ns standard MD simulations with the TM Cα atoms restrained to verify that the side chains would not rearrange when no biasing forces are applied. Indeed, no significant changes were recorded after 100 ns, as shown in Figure S7B of the Supporting Information. Three additional unrestrained standard MD simulations of 100 ns each were conducted starting from the structure obtained at the end of the ABMD simulation to verify possible changes in the size of the allosteric sodium pocket assessed by monitoring the Cα distance between residues I3.40 and N7.45. This distance is 12.0 Å in the ultra-high-resolution inactive crystal structure of δ-OR but was reduced to 10.2 Å by the end of the ABMD simulation. During the three independent MD simulations without any restraint, the I3.40–N7.45 distance gradually returned to around 12 Å in two of the three runs (Figure S8 of the Supporting Information).

To explore the possible egress of the ion from the allosteric sodium site at much longer time scales, we used random accelerated molecular dynamics simulations (RAMD) (see Materials and Methods for details). From clustering of the RAMD trajectories (see Table S1 of the Supporting Information for details), three major egress paths emerged (shown in Figure 4). This figure shows that sodium could exit the receptor TM bundle through either the extracellular side or the intracellular side. Sodium interacted with N3.35 and D3.32 side chains along one of the most favorable pathways through the extracellular side (EC [TM3, EL2], 16% of total egress trajectories). When the sodium ion engaged in interaction with D3.32, the interaction of the ligand amine with D3.32 momentarily broke to promptly re-form soon after the ion left the site. The alternative most favorable pathway through the extracellular side (EC [TM2, EL2], 14% of total egress trajectories) differs from EC [TM3, EL2] in that the ion is coordinated by N3.35 and S7.46 after it leaves crystal coordination and never interacts with D3.32. On the other hand, the major sodium egress pathway through the intracellular side (IC [TM2, TM7], 15% of total egress trajectories) is characterized by an interaction of sodium with D7.57 at the cytoplasmic end of TM7. Prior to engaging in this interaction, the ion is mostly coordinated by backbone oxygen atoms.

Figure 4

Representative main ion egress pathways derived from RAMD simulations. Pathways are shown as colored tubes, while charged δ-OR side chains that interact with the ion along the pathways are shown as sticks. The bound agonist SNC-80 is shown in ball-and-stick representation.

Radioligand Binding Experiments Confirm a Similar Effect of Sodium on Ligand Binding at All Major ORs

Using stable cell lines expressing KOR-1, DOR-1, and MOR-1, we examined the effects of sodium ions on agonist binding using [3H]U69,593, [3H]DPDPE, and [3H]DAMGO, respectively, and [3H]diprenorphine to measure antagonist binding (Figure S9 of the Supporting Information). Binding of the radioligand to the cloned receptors was greatly influenced by sodium ions. Binding of the agonist to all three receptors was sensitive to sodium ions, with a loss of 50% of specific binding seen at approximately 25 mM NaCl. Of the three, binding of [3H]U69,593 to the κ-OR was slightly more sensitive. Sodium ions had an opposite effect on the binding of the antagonist [3H]diprenorphine. There was a modest increase in the level of binding to the κ-OR and the δ-OR with increasing sodium concentrations. However, the level of binding to the μ-OR was increased over 40%.

Saturation studies with [3H]DAMGO were best fit with a two-site model (Table S2 of the Supporting Information). When the experiments were performed in the presence of sodium ions at 100 mM, the high-affinity Bmax was almost 95% lower with little effect on the Bmax of the lower-affinity site, similar to results previously reported in brain membranes.[1] These observations are consistent with the ability of sodium ions to shift the conformation of the receptor from a high-affinity agonist state to a low-affinity agonist state. Analysis of the saturation results for the κ-OR and δ-OR also showed decreases in binding levels. Unlike those of the μ-OR, these data were best fit with a single site. Our inability to detect a lower-affinity binding site may reflect the technical difficulties involved upon examination of binding at the high radioligand concentrations needed to define the low-affinity binding site.

Discussion

Allosteric effects of sodium on GPCR agonist binding and activation have been known for years but were recently revived by new insights provided by ultra-high-resolution (1.8–2.2 Å) crystallographic structures of the A2AR (PDB entry 4EIY(22)), the B1AR (PDB entry 4BVN(5)), the PAR1 (PDB entry 3VW7(23)), and the δ-OR (PDB entry 4N6H(24)), which have revealed important details of binding of sodium to GPCRs (recently reviewed in ref (6)). Not only have these structures shown that the ion is coordinated by chemically and conformationally conserved residues such as D2.50 and S3.39 in the middle of the 7TM bundle, but they have also revealed an extended network of conserved water molecules forming hydrogen bonds with a number of highly conserved receptor residues across family A GPCRs, including L2.46, A2.49, W6.48, N7.45, and N7.49.[22]

The most striking differences in the way the ion binds were observed in the PAR1[23] and δ-OR[24] ultra-high-resolution crystal structures, because of the involvement of normally nonconserved residues (e.g., D7.49 in PAR1 and N3.35 in δ-OR) in the first coordination shell of the ion. Notably, the nitrogen atom of the N3.35 side chain linked the sodium allosteric site to the orthosteric ligand binding site in the δ-OR crystal structure through a water-mediated hydrogen bond with the main chain carbonyl atom of D3.32, which is the residue forming a salt bridge with the ligand nitrogen group. Mutating N3.35 to alanine or valine abrogated or reduced the allosteric effect of sodium on ligand binding, in addition to abolishing the G-protein signaling while producing high levels of constitutive activity at the noncanonical β-arrestin-mediated signaling pathway.[24] Similarly, mutating D2.50, N7.45, and N7.49 to alanine transformed classical δ-OR antagonists into potent β-arrestin-biased agonists,[24] highlighting a role for sodium-coordinating residues as efficacy switches between G-protein- and β-arrestin-mediated signaling pathways within δ-OR.

Similar to the results of simulations reported in the literature,[20,21] the simulations presented here, in which sodium ions diffuse freely from the bulk to the interior of the δ-OR, μ-OR, or κ-OR bundles, show that sodium always enters the receptor from the extracellular milieu and reaches the allosteric site seen crystallographically quite rapidly (nanosecond time scale), following similar pathways. Specifically, before coordinating D2.50, S3.39, N7.35, and two water molecules at its allosteric crystallographic site, sodium stops at the orthosteric site where it is found to coordinate D3.32. However, the ion density at D3.32 is not found to be as pronounced as at nonconserved E6.58 and EL2 D204 residues in κ-OR, or EL3 D293 and D288 residues in δ-OR. These molecular determinants appear to be responsible for the different tendency of a second ion to enter the receptor from the extracellular side. In the absence of these molecular determinants, sodium appears to take longer to access the μ-OR orthosteric site, and we speculate that the possible reduced level of competition between sodium and the ligand for this site may be responsible for the differential increased antagonist binding affinity of μ-OR observed by the specific binding experiments conducted in transfected cells.

On the other hand, specific ligand binding experiments using subtype-selective agonists in cells transfected with the individual ORs showed that, unlike antagonist binding, binding of agonists to either receptor is consistently affected by sodium. Decreases in agonist affinity and in the maximal number of binding sites were recorded for all three receptor subtypes, although a slightly larger decrease in affinity was measured for κ-OR than for μ-OR and δ-OR. We sought a mechanistic explanation for the decrease in agonist affinity and the maximal number of binding sites by adding to the analysis of the ligand-free simulations that of additional biased or unbiased MD simulations of antagonist-bound and agonist-bound, activated receptors conducted using the δ-OR ultra-high-resolution crystal structure as a model system. A few differences in the dynamics of the three OR systems were noted in our simulations, most notably the higher flexibility of the S3.39 and N3.35 side chains of δ-OR compared to that of the side chains of μ-OR and κ-OR during the simulated time scales. The flexibility of these two residues, as well as that of the so-called rotamer toggle switch (W6.48), which is a sodium second-coordination shell residue, is impaired in the presence of an antagonist, implying that an intact network linking the allosteric sodium binding site to the orthosteric ligand binding site is what stabilizes the first and second coordination shells of the ion and perhaps prevents the large-scale movement of TM6 away from TM3 leading to activation.

Finally, in light of recent suggestions that the collapse of the allosteric sodium site accompanying activation (reviewed in ref (6)) may be triggering the departure of the ion from the molecule, we investigated the whereabouts of the ion following the transition from an inactive to active conformation of the receptor. In contrast to previous suggestions[6,25] that the ion leaves the allosteric binding pocket upon activation and does so through the cytoplasmic side, we find that undergoing a transition from an inactive to an active conformation in the presence of sodium does not produce the spontaneous release of the ion from the bundle, at least not on the simulated time scales. Reasoning that the ion may take much longer to leave the receptor, we forced its departure from the bundle with random accelerated MD and studied its possible egress pathways. Notably, if the ion ever had to leave the allosteric binding site in the middle of the bundle, it would do so through either the extracellular side or the intracellular side. We do not know what egress pathway is preferable on the basis of the results of these simulations, but the latter suggests residues along the pathways (e.g., N3.35, D3.32, S7.46, and D7.57) that may play a role in the egress of the ion from the bundle, and hence be worthy of experimental testing. Notably, this group does not include residue D3.49 of the so-called DRY motif, which had been hypothesized to participate in the egress of the ion from the bundle.[6]