Proposed Mode of Binding and Action of Positive Allosteric Modulators at Opioid Receptors.

Available crystal structures of opioid receptors provide a high-resolution picture of ligand binding at the primary ("orthosteric") site, that is, the site targeted by endogenous ligands. Recently, positive allosteric modulators of opioid receptors have also been discovered, but their modes of binding and action remain unknown. Here, we use a metadynamics-based strategy to efficiently sample the binding process of a recently discovered positive allosteric modulator of the δ-opioid receptor, BMS-986187, in the presence of the orthosteric agonist SNC-80, and with the receptor embedded in an explicit lipid-water environment. The dynamics of BMS-986187 were enhanced by biasing the potential acting on the ligand-receptor distance and ligand-receptor interaction contacts. Representative lowest-energy structures from the reconstructed free-energy landscape revealed two alternative ligand binding poses at an allosteric site delineated by transmembrane (TM) helices TM1, TM2, and TM7, with some participation of TM6. Mutations of amino acid residues at these proposed allosteric sites were found to either affect the binding of BMS-986187 or its ability to modulate the affinity and/or efficacy of SNC-80. Taken together, these combined experimental and computational studies provide the first atomic-level insight into the modulation of opioid receptor binding and signaling by allosteric modulators.

Opioid receptors belong to the superfamily of G protein-coupled receptors (GPCRs), the largest class of membrane proteins encoded by the human genome.[1] Like other GPCRs, they respond to extracellular stimuli such as small molecules, peptides, and ions, by undergoing ligand-specific conformational changes and consequently recruiting and activating accessory proteins such as G proteins and β-arrestins. Molecules targeting opioid receptors are efficacious therapeutic agents, especially against pain[2] or for controlling addiction,[3] and a growing body of evidence suggests that opioid receptor-targeted therapeutics may also have utility for the treatment of mood disorders.[4] However, the beneficial action of current opioid drugs is severely limited by the development of adverse effects, including respiratory depression, tolerance, dependence, constipation, and abuse liability.[5]

Most prescription opioid analgesics are small molecules (e.g., morphine, buprenorphine, codeine, fentanyl, etc.), and all bind to the primary “orthosteric” site on the receptor, which is targeted by endogenous peptides. The recent crystal structures of all four opioid receptor types, namely the δ, μ, κ, and nociceptin receptors,[6−10] provide high-resolution atomic insight into drug binding at opioid receptor orthosteric sites. This binding mostly involves amino acid residues in transmembrane (TM) helices TM3, TM6, and TM7, although interactions with TM2 have also been reported, depending on the specific ligand chemotype.[11,12]

In search of opioid therapeutics with reduced adverse effects, recent high-throughput screening campaigns have identified positive allosteric modulators (PAMs) of μ- and δ-opioid receptors,[13,14] that is, ligands that potentiate the orthosteric agonist-induced response of these receptors. This represents a significant discovery in view of the potential advantages of allosteric modulators compared to classical orthosteric ligands.[15] For instance, by targeting nonconserved, allosteric regions of the receptor, allosteric ligands may afford significant receptor subtype selectivity, resulting in limited off-target effects. This advantage of PAMs may be less important in the case of opioid receptors, where several selective opioid ligands are already available, and undesired responses often result from on-target effects of drugs binding the receptor in different tissues or brain regions. A potentially more important feature of opioid PAMs is that they would only act in the presence of the endogenous ligand, thus maintaining the temporal and spatial fidelity of opioid signaling in vivo. A direct consequence of this trait is potentially reduced receptor desensitization, resulting in the attenuation of the dependence and tolerance that are produced by classical orthosteric agonists. Additionally, opioid PAMs may be able to avoid opioid receptor-mediated adverse effects such as constipation and respiratory depression, again by virtue of acting only at cells/tissues where native opioid signaling is naturally occurring. Moreover, as the effect of allosteric modulators is limited by the degree of cooperativity with the orthosteric agonist, opioid PAMs are expected to exert fewer on-target overdosing risks, which represent a serious practical issue for current prescription painkillers. The reader is referred to a recent review[16] for additional information about potential advantages of opioid PAMs.

Using inferences from binding kinetic experiments, cannabidiol was previously proposed as a negative allosteric modulator of δ and μ receptors,[17] although it remains unclear whether this effect is mediated directly via the opioid receptors, or potentially via cooperative interactions between homo- or hetero-oligomeric receptors. In contrast, more recent studies have identified novel small molecules, e.g., BMS-986121 and BMS-986122, as first-in-class PAMs acting directly on the μ-opioid receptor.[18] Similarly, compounds BMS-986187 and BMS-986188 were recently shown to behave as potent PAMs at the δ-opioid receptor, modulating the affinity and/or efficacy of both peptidic (leu-enkephalin) and small-molecule agonists (SNC-80 and TAN-67) at several biological end points (receptor binding, G protein activation, β-arrestin recruitment, adenylyl cyclase inhibition, and extracellular signal-regulated kinases (ERK) activation).[13] Although these molecules are assumed to bind to sites that are topographically distinct from the orthosteric binding site, their actual binding pockets and modes are unknown. Through a combination of enhanced molecular dynamics simulations of binding of BMS-986187 to the δ-opioid receptor in the presence of an orthosteric ligand and experimental validation, the work reported here provides the first structural insights into the allosteric modulation of opioid receptors by small molecules.

Results and Discussion

Predictions from multiple-walker well-tempered all-atom metadynamics simulations were used in combination with experimental structure–function analysis of binding and functional assays on both wild type (WT) and selected mutants of the δ-opioid receptor to identify the preferred modes of binding of the recently identified PAM BMS-986187[13] to the δ-opioid receptor in complex with the selective orthosteric agonist SNC-80 and embedded in an explicit lipid–water environment.

Free-energy Landscape for the Binding Process of BMS-986187 to the δ-Opioid Receptor

Automated docking algorithms using scoring functions of varying complexity and accuracy have successfully been applied to opioid receptor crystal structures over the past few years, allowing the identification of novel compounds that bind at the orthosteric binding site of these receptors (e.g., see refs (19), (20)). Our initial attempts of using automated docking algorithms to predict the binding mode of BMS-986187 at the δ-opioid receptor were hindered by the identification of multiple positions and orientations of the ligand in a relatively large binding pocket within the extracellular half of the receptor. The lack of crystal structures of opioid receptors bound to allosteric modulators and rough estimates of binding affinity from oversimplified scoring functions did not allow us to unambiguously favor a particular binding mode of the δ-PAM over another, thus motivating us to explore the predictive power of more sophisticated methods.

Recently, all-atom molecular dynamics (MD) simulations have successfully been applied to predict the spontaneous binding of allosteric modulators to GPCRs (e.g., see ref (21)). However, since ligand binding at a target site is a relatively rare event on microscopic time-scales, these simulations have required several million computing hours on special-purpose computational resources and remain of limited accessibility. We recently proposed the use of metadynamics[22] as a more efficient enhanced sampling method to study the long time scale process of ligand binding to GPCRs.[23] We first applied this method to successfully predict ligand binding to the δ-opioid receptor at a time when the crystal structure was not available yet[24] and have further validated it using crystal structures of various receptor types. Briefly, the method accelerates the conformational sampling of a system by adding to the potential energy a history-dependent term acting on a small number of collective variables (CVs) representing the slow degrees of freedom relevant to the process under study.

In the study reported here, we simulated the binding of the allosteric modulator, BMS-986187, to the δ-opioid receptor embedded in an explicit lipid–water environment and with SNC-80 bound to the orthosteric site. We carried out a total of 3.6 μs multiple-walker well-tempered metadynamics simulations, biasing the potential along the following two CVs: the BMS-986187–receptor distance and the number of BMS-986187–receptor interactions, herein labeled CV1 and CV2, respectively. To better characterize the different binding modes adopted by BMS-986187, the reweighting strategy described in ref (25) was used to calculate the free-energy as a function of the biased CVs (CV1 and CV2) and two additional degrees of freedom, namely the Z-component of the vector linking the ligand’s tricyclic structure with its ortho-substituted benzyl ring (CV3) and the XY projection of the BMS-986187-receptor distance (CV4; see the Methods section for more details).

A projection of the reconstructed four-dimensional free-energy landscape onto CV1 and CV2 is shown in Figure . This free-energy reveals a broad basin centered at a ∼2 nm distance between BMS-986187 and the receptor center (CV1), and encompassing binding states characterized by a different number of BMS-986187-receptor contacts (CV2), ranging approximately from 4 to 6. Free-energy convergence as a function of the different CVs was evaluated by plotting differences of its two-dimensional projections between the full simulation and its first three-quarters. The convergence of the reconstructed free-energy projected onto CV1 and CV2 is shown in Supporting Information Figure S1 as an example.

Figure 1

Reconstructed four-dimensional free-energy (in kcal/mol) landscape of the BMS-986187 binding to the δ-opioid receptor projected onto the ligand–receptor distance (CV1) and the number of contacts between the ligand and the receptor (CV2). Projections of the cluster medoids are indicated with circles whose areas are proportional to the cluster’s relative probability, with larger circles indicating the lowest free-energies.

Low free-energy states (within 2 × kBT of the minimum) in the four-dimensional phase-space defined by CV1–CV4 correspond to highly populated BMS-986187 binding poses. To characterize these poses, we calculated interaction fingerprints involving either the orthosteric or allosteric ligand and the δ-opioid receptor. We then clustered these poses based on fingerprint similarity and estimated the relative free-energy of the different states (see Methods section for details). This analysis revealed seven compact clusters characterized by different binding pockets and binding poses for the allosteric ligand. The positions of the cluster medoids are indicated by circles on the free-energy projection onto the CV1−CV2 plane shown in Figure . While clusters 2 and 3 partially overlap in the CV1–CV2 projection, their separation is evident when looking at the corresponding cluster medoids on the free-energy projection onto CV1 and CV3 (see Figure S2).

Predicted Lowest Free-Energy Binding Modes of BMS-986187 at the δ-Opioid Receptor

The results of the enhanced simulations reported herein suggest that the allosteric modulator BMS-986187 adopts multiple metastable binding states during the binding process. The two most stable binding states of BMS-986187 (red and blue sticks in Figure ) occupy the same binding pocket within the δ-opioid receptor and have virtually indistinguishable free energy (i.e., there is just a 0.03 kcal/mol difference between them, which is smaller than the estimation error). Notably, there is no appreciable similarity between these predicted poses from simulations and the top-ranked docking conformations obtained using Glide XP version 6.9. As shown in Figure , representative structures of these two most stable BMS-986187 states share common interactions of the 2-methyl-benzyl group with the receptor but show different orientations of the fused tricyclic moiety within a binding pocket of the δ-opioid receptor surrounded by TM1, TM2, and TM7 helices, with some involvement of TM6. Figure and Table S1 show the details and probabilities of the BMS-986187-receptor and the SNC-80-receptor interactions for these representative, stable binding poses. In both states 1 and 2 (red and blue sticks, respectively, in Figure ), which are representative conformations of clusters 1 and 2 of the ligand–receptor interaction fingerprints, the allosteric ligand (see Figure a) forms direct polar, water-mediated polar, hydrophobic, and/or aromatic interactions with human δ-opioid receptor residues I52(1.35), Y56(1.39), L102(2.57), Y109(2.64), E112(2.67), V297(7.32), H301(7.36), I304(7.39), and Y308(7.43) (numbers in parentheses refer to two-digit numbers as per the Ballesteros–Weinstein generic numbering scheme[26]). The differences between these two poses (illustrated in separate panels in Figure S3) are mainly hydrophobic interactions between BMS-986187 and residues L46(1.29), L48(1.31), and L110(2.65) in conformations of cluster 1 but not cluster 2, and hydrophobic interactions between BMS-986187 and residues W284(6.58) and L300(7.35), and water-mediated interactions between BMS-986187 and residues Q105(2.60) and K108(2.63) in conformations of cluster 2, but not cluster 1. Note that the orientation of the Y109(2.64) side chain in state 1 is similar to that of the naltrindole-bound δ-opioid receptor crystal structure, but the Y109(2.64) side chain (shown in Figure a as red and blue transparent sticks for states 1 and 2, respectively) must rotate up and away from the center of the helical bundle for BMS-986187 to adopt the conformation of state 2.

Figure 2

Representative conformations of BMS-986187 states 1 and 2 (in red and blue, respectively). The left panel (a) illustrates a close-up of the binding region. In the right panel (b), below the 2D-diagrams of BMS-986187 and SNC-80, a top view of the δ-opioid receptor backbone in state 2 is shown in gray cartoon representation in addition to sections of the δ-opiod receptor and ligand densities in a plane normal to the membrane. Water molecules are shown as spheres. Receptor residues forming contacts with the allosteric ligand in state 2 are shown as sticks, colored by atom type, and labeled according to the Ballesteros–Weinstein numbering scheme. The conformers of residue Y109(2.64) in states 1 and 2 are highlighted in red and blue, respectively. The orthosteric agonist SNC-80 is shown in orange sticks. EL2 has been omitted to ease visualization.

Figure 3

Interactions formed in the top states 1 and 2 (in red and blue, respectively) by the allosteric ligand BMS-987187 and the orthosteric ligand SNC-80 in panels a and b, respectively. The opacity of the links is proportional to the probability of the interaction. Gray histograms represent the overall probability of forming an interaction with each residue.

Notably, during the whole simulation of BMS-986187 binding to the δ-opioid receptor, the orthosteric ligand SNC-80 maintains a stable orientation between helices TM3, TM5, TM6, and TM7 in the orthosteric binding pocket. A comparison of the interactions that the orthosteric ligand SNC-80 forms in either state 1 or state 2, is provided through the wheel plots of Figure b while a full list of the interactions that have a 20% probability for each pose is reported in Table S1. Specifically, the receptor residues involved in interaction with SNC-80 in both states are D128(3.32), Y129(3.33), M132(3.36), K214(5.39), V217(5.42), W274(6.48), I277(6.51), H278(6.52), F280(6.54), V281(6.55), W284(6.58), and L300(7.35), while the hydrophobic interactions between SNC-80 and I304(7.39) and Y308(7.43) are present in state 2 only.

Less stable binding modes of BMS-986187, with relative free energies higher than kBT, and thus significantly less populated at room temperature, included four additional binding states. These states, labeled states 3, 4, 5, and 7 in Figure S4, have free energies of 0.97, 0.98, 1.20, and 4.12 kcal/mol, respectively, relative to state 1. An additional state is observed in which the allosteric ligand is trapped in a location between the membrane and the receptor (state 6 in violet in Figure S4 has a relative free-energy of 3.62 kcal/mol relative to state 1). A comprehensive list of all the interactions established by both the allosteric and orthosteric ligands in each of these states is provided in Table S1.

Experimental Testing of the Predicted Binding Modes of BMS-986187 at the δ-Opioid Receptor

The contribution of several δ-opioid receptor residues to the affinity of the allosteric ligand BMS-986187 and its ability to modulate the binding and/or the efficacy of the orthosteric agonist SNC-80 were probed experimentally via functional assays conducted at the WT and selected mutants of the human δ-opioid receptor. Specifically, the following mutants were generated and tested in an attempt to identify the molecular determinants responsible for BMS-986187 binding and/or modulatory action, as well as ways to discriminate between the two identified lowest-energy states of BMS-986187: Y56(1.39)A, Q105(2.60)A, K108(2.63)A, K108(2.63)N (as in μ-opioid receptor; a valine or an aspartic acid in κ-opioid receptor and nociceptin receptor, respectively), Y109(2.64)A, Y109(2.64)I (as in the nociceptin receptor), W284(6.58)K (as in μ-opioid receptor; a glutamic acid or a glutamine in κ-opioid receptor and nociceptin receptor, respectively), L300(7.35)W (as in μ-opioid receptor; a tyrosine in κ-opioid receptor), H301(7.36)R (as in nociceptin receptor), and H301(7.36)Y (as in κ-opioid receptor). Among them, W284(6.58)K and L300(7.35)W mutations had previously been shown to affect significantly the binding of SNC-80 to the δ-opioid receptor.[27,28]

First, cell surface expression levels of all δ-opioid receptor mutants were determined by anti-HA ELISA (see Figure S5); any significant variations relative to the expression of the WT were factored into the final parameter estimates of the signaling efficacy for both the orthosteric and allosteric ligands (see below).

Using ERK1/2 phosphorylation (pERK1/2) assays as a measure of δ-opioid receptor activity, concentration–response curves to the orthosteric agonist SNC-80 were constructed in the absence or presence of increasing concentrations of the PAM BMS-986187, at the WT and the aforementioned mutant δ-opioid receptors (Figure S6). The potency (EC50) and maximum response to SNC-80 (Emax) at the WT and mutant receptors are reported in Table S2. Analysis of the interaction between SNC-80 and BMS-986187 was performed using an operational model of allosterism and agonism.[29] Briefly, this model is obtained by combining a simple allosteric ternary complex model with the operational model of agonism. As such, it allows the description of the ability of both the orthosteric and the allosteric ligands to exhibit agonism (incorporating the intrinsic efficacy of each ligand, the total density of receptors, and the efficiency of the coupling of the ligand stimulus to the assay response), as well as capturing allosteric effects of BMS-986187 on both binding affinity and efficacy of SNC-80 (see Methods section for details).

This analysis yielded estimates of three key parameters: (i) the functional affinity of the allosteric modulator BMS-986187 for the unoccupied δ-opioid receptor (pKB, Table ), (ii) an overall cooperativity value, αβ, which incorporates modulatory effects on both orthosteric agonist affinity (α) and efficacy (β) (Table ), and (iii) the ability of each ligand to directly activate the receptor in its own right (τA for SNC-80 and τB for BMS-986187; Table S2 and Table S3). In addition, we performed whole cell radioligand competition binding experiments to allow us to discriminate, where possible, the modulatory effect of BMS-986187 solely on binding cooperativity at the WT and mutant receptors (i.e α). These experiments were performed using [3H]diprenorphine and increasing concentrations of SNC-80 in the absence or presence of increasing concentrations of the PAM BMS-986187 (Figure S7). The control SNC-80 competition curves were fitted to a one-site competitive binding model to derive affinity estimates (pKi) for the orthosteric agonist (Table S2), whereas the entire data set in the presence of an allosteric modulator was fitted to an allosteric ternary complex model[29] to derive estimates of the cooperativity between the PAM and the radioligand (α′; see Table S3) as well as the cooperativity between the PAM and SNC-80 (α). The radioligand interaction studies revealed that BMS-986187 displays negative cooperativity with [3H]diprenorphine (0 < α′ < 1, Table S3) and neutral binding cooperativity (α ∼ 1; see Table ) with SNC-80 at the WT receptor. This result extends from our initial functional characterization of BMS-986187 and reveals that this allosteric modulator is mainly an efficacy modulator of orthosteric agonists. A similar behavior of the PAM was observed at the Q105(2.60)A, Y109(2.64)A, and K108(2.63)A/N mutants, while it showed positive binding cooperativity with SNC-80 at Y109(2.64)I and H301(7.36)R/Y (Table , Figure S7). For the Y56(1.39)A, W284(6.58)K, and L300(7.35)W mutants, no modulation of SNC-80 or [3H]diprenorphine was observed, with the competition binding curves in the presence of PAM overlaying with the curve observed in the absence of PAM (Figure S7). Therefore, by combining our data from both functional and radioligand binding experiments, we were able to validate the predicted allosteric binding pocket at opioid receptors by relating mutations to changes in the affinity of BMS-986187 and its allosteric effects upon both orthosteric ligand affinity and efficacy. Figure represents a graphical illustration of the mutation-induced changes of the estimated allosteric parameters of BMS-986187.

Table 1

Estimates of BMS-986187 Binding Affinity (KB), and the Cooperativity Parameter αβ Which Incorporates Modulatory Effects on Orthosteric Agonist Affinity (α) and Efficacy (β), at the WT or Mutant δ-opioid Receptorsa

mutant	pKB	log αβ	log α	log β
WT	6.24 ± 0.20	1.08 ± 0.09	= 0	1.08 ± 0.09
Y56(1.39)A	5.64 ± 0.25	1.28 ± 0.31	= 0 (NC)	1.28 ± 0.31
Q105(2.60)A	5.09 ± 0.30*	2.08 ± 0.67	= 0	2.08 ± 0.67
K108(2.63)A	5.34 ± 0.40	1.77 ± 0.41	= 0	1.77 ± 0.41
K108(2.63)N	5.70 ± 0.34	0.55 ± 0.32*	= 0	0.55 ± 0.32*
Y109(2.64)A	4.93 ± 0.77	2.22 ± 0.72	= 0	2.22 ± 0.72
Y109(2.64)I	4.82 ± 0.53**	2.46 ± 0.36	0.54 ± 0.21*	1.92 ± 0.42**
W284(6.58)K	5.49 ± 0.19*	1.36 ± 0.24	= 0 (NC)	1.36 ± 0.24
L300(7.35)W	4.74 ± 0.11***	2.14 ± 0.12**	= 0 (NC)	2.14 ± 0.12**
H301(7.36)R	5.27 ± 0.31	1.65 ± 0.34	0.90 ± 0.29**	0.75 ± 0.45
H301(7.36)Y	4.95 ± 0.30***	1.96 ± 0.36	0.24 ± 0.27	1.72 ± 0.45

Parameters are obtained by operational model fitting of data from [3H]diprenorphine binding and pERK phosphorylation assays. Equal sign before a value indicates that it was constrained during the fitting of the operational model; “NC” indicates that the value could not be obtained from the binding studies. Statistical significance levels from Dunnet’s test p values are indicated as stars (respectively p < 0.05, 0.01, 0.001).

Figure 4

Effects of δ-opioid receptor mutations on BMS-986187 binding affinity (KB), and on the orthosteric agonist affinity (α) and efficacy (β). Bars represent the difference in pKB (top panel), binding cooperativity (log α, middle panel), and efficacy cooperativity (log β, bottom panel) of BMS-986167 relative to WT as derived from functional and binding interaction experiments with SNC-80 (Table ). Data represent the mean ± SE of three experiments performed in duplicate. Statistical significance levels from Dunnet’s test p values are indicated as stars (respectively p < 0.05, 0.01, 0.001).

At the WT receptor, BMS-986187 exhibited an affinity of ∼0.6 μM and positive allosteric modulation with a functional cooperativity factor of β ∼ 12 (Table ). Five of the 10 tested mutants, specifically Q105(2.60)A, Y109(2.64)I, W284(6.58)K, L300(7.35)W, and H301(7.36)Y, display a significant reduction of the estimated BMS-986187 binding affinity (KB) with respect to the WT receptor (see Table , Figure ), suggesting the involvement of the corresponding wild-type δ-opioid receptor residues in a direct or water-mediated interaction with the allosteric ligand. Two of these five residues (W284(6.58) and L300(7.35)) are involved in interactions with SNC-80 in the two states (see Table S1), and the pronounced effect of their mutations on affinity and efficacy of SNC-80 (see Table S2) is in line with the hypothesis of a direct effect on the binding of the orthosteric ligand as well as previously published experimental data.[27,28]

Mutation of the Y109(2.64) residue to isoleucine did not only affect BMS-986187 binding affinity but it also resulted in a significant increase in both its binding and functional cooperativity (α and β values in Table , respectively). Interestingly, the L300(7.35)W mutation also displayed a significant increase in the functional cooperativity of BMS-986187 (βL300W ∼ 140 vs βWT ∼ 12), despite negligible effects in the binding interaction experiments (αL300W ∼ 1). Finally, two mutations that had no effect on BMS-986187 binding were found to affect the binding cooperativity (H301(7.36)R) or functional cooperativity (K108(2.63)N) between the modulator and the orthosteric agonist.

Collectively, our data support BMS-986187 state 2 (blue sticks in Figure ) as the most substantiated state, as it is the only one where all the δ opioid-receptor residues that are inferred to contribute to BMS-986187 binding affinity (e.g., Q105(2.60), Y109(2.64), W284(6.58), L300(7.35), and H301(7.36)) have a high probability to be involved in interaction with the allosteric ligand (see Table S1). In particular, residues Q105(2.60), W284(6.58), and L300(7.35) are found to form interactions with BMS-986187 in state 2 only. In spite of these differences, our computational predictions and experimental validation are currently insufficient to discriminate unambiguously between the ligand’s states 1 and 2, and the current assumption is that they would be equally populated at room temperature. Nonetheless, both identified states and experimental data support the first identification of an opioid receptor allosteric pocket in proximity, but more extracellular, to the orthosteric binding site. Notably, this binding pocket, mostly defined by TM1, TM2, and TM7, does not coincide with that seen in the crystal structure of the peptide-bound δ-opioid receptor,[30] where the DIPP-NH2 tetra-peptide extends toward the receptor’s TM3 and extracellular loop (EL) 2. Although this peptide binding region does not overlap with states 1 or 2, it does so with the less stable state 3 identified in our simulations (see Figure S8).

Conclusion

We have presented here results from computational predictions and experimental validation that provide the first atomic-level insight into the modulation of opioid receptor binding and/or signaling by allosteric modulators. Specific binding sites and conformational states were identified for BMS-986187, a recently discovered positive allosteric modulator of the δ-opioid receptor, based on inferences from free-energy calculations. This information led to the identification of specific molecular determinants that are responsible for binding of the allosteric modulator to the receptor and/or its ability to modulate the binding affinity and/or efficacy of the orthosteric agonist SCN-80. Specifically, mutations of residues Q105(2.60), K108(2.63), Y109(2.64), W284(6.58), L300(7.35), and H301(7.36) are shown to have an effect on either the allosteric modulator binding affinity, or cooperativity, or both. While Q105(2.60)A, W284(6.58)K, and H301(7.36)Y are shown to only affect the binding affinity of the allosteric ligand, Y109(2.64)I and L300(7.35)W appear to affect both the allosteric binding affinity and its cooperativity. In contrast, mutants such as H301(7.36)R and K108(2.63)N are found to only affect allosteric cooperativity. The results presented here represent the initial step in the identification of the binding sites of allosteric modulators at the opioid receptor family. Future work will assess whether increased PAM affinity or cooperativity can be achieved through alternative mutations within these sites. For instance, the replacement of Y109(2.64) with a positively charged amino acid may increase PAM affinity and/or cooperativity based on the results of our simulations. Together, these data will inform future rational design and drug discovery efforts toward novel allosteric modulators of the δ-opioid receptor.

Methods

System Setup and Equilibration

The ultrahigh resolution crystal structure of the δ-opioid receptor (PDB ID: 4N6H(6)) without the flexible N-terminal amino acid residues (G36SPGA40) was used as a starting conformation. With the exception of crystallographic water molecules, all other non-protein atoms, as well as the cytochrome b562 RIL insert, were removed. The δ-selective agonist SNC-80 (4-[(R)-[(2S,5R)-4-allyl-2,5-dimethylpiperazin-1-yl](3-methoxyphenyl)methyl]-N,N-diethylbenzamide) was docked into the ligand-free receptor using the Schrödinger Suite 2014-1 and the strategy we previously reported in the literature.[31] Briefly, the ligand was prepared using LigPrep version 2.9[32] at physiological pH while the Protein Preparation Wizard tool was used to add hydrogen atoms and assign bond orders. Flexible ligand docking of SNC-80 to the δ-receptor was performed using extra-precision Glide XP version 6.2.[33]

The complex formed by SNC-80 and the δ-opioid receptor was embedded in a pre-equilibrated patch of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and 10% cholesterol, further solvated with water and NaCl at 150 mM. Additional chloride ions were added to neutralize the system. The entire simulation box measured 83 × 83 × 110 Å3 and consisted of the δ-receptor, ∼20 cholesterol molecules, ∼200 POPC molecules, ∼15 000 water molecules, ∼35 sodium ions, and ∼50 chloride ions, totaling ∼77 600 atoms. This system was simulated using the TIP3P water model and the Charmm36 force field for the protein, lipids, and ions. Ligand parameters were obtained from the CHARMM General Force Field, validated according to the published guidelines,[34] and are available upon request.

All MD simulations reported herein (restrained, unrestrained, and metadynamics simulations) were performed using the Gromacs 4.6 simulation package[35] in the NPT ensemble at 300 K and 1 bar, with a Nose–Hoover thermostat,[36] Parrinello–Rahman pressure coupling,[37] and a 2 fs time step. All bonds were constrained using the LINCS algorithm,[38] and short-range, nonbonded interactions were cut off at a 10 Å distance.

The SNC-80-bound δ-receptor in the membrane mimetic environment was equilibrated for 6 ns with decreasing positional restraints first set on all heavy atoms and then on Cα atoms only, followed by a 50 ns unrestrained equilibration run. The pose of SNC-80 remained stable during the equilibration run as judged by a heavy atom root-mean-square deviation (RMSD) of 0.9 Å between the initial docking pose and the equilibrated one. An initial pose of BMS-986187 (3,3,6,6-tetramethyl-9-(4-((2-methylbenzyl)oxy)phenyl)-3,4,5,6,7,9-hexahydro-1H-xanthene-1,8(2H)-dione) inside the δ-receptor was obtained by docking the ligand into the equilibrated SNC-80-bound δ-receptor using the same procedure used to dock SNC-80. This system was equilibrated for 6 ns with decreasing positional restraints set first on all heavy atoms and then on Cα atoms only. To obtain additional, different starting structures of BMS-986187 for the multiple-walker metadynamics procedure, both inside the receptor and in the bulk solvent, a 40 ns metadynamics run was carried out using PLUMED version 1.3[39] and the GROMACS version 4.6 package as well as one CV describing the distance between the centers of mass of the BMS-986187 heavy atoms and the transmembrane helical bundle of the δ-receptor. A deposition rate of 10 ps, a bias factor of 10, a Gaussian hill height of 0.8 kJ/mol, and a Gaussian width of 0.125 Å were used in these metadynamics simulations. Lower and upper limits of 8 and 37 Å, respectively, on the aforementioned CV were imposed using steep, harmonic restraints with an elastic constant of 4000 kJ/nm2. It must be noted that starting positions of the different walkers could have been provided by other methods (e.g., automated docking algorithms). The use of a short metadynamics run to generate initial conformations, has, however, some technical advantages. First and foremost, it provides a validation of the ability of the chosen CVs to effectively enhance the sampling of different ligand conformations. Moreover, it generates solvated and equilibrated starting systems, whereas independently docked ligands would have to undergo additional preparation steps before any MD-based simulation could be performed.

Multiple Walker Metadynamics

Metadynamics[22] enhances the sampling of selected CVs by applying to the system’s dynamics a history-dependent, periodically updated bias that discourages the system from revisiting conformations that have already been sampled. In its well-tempered version,[40] the bias contributions are scaled so that the CV can overcome free-energy barriers of controlled height, allowing assessment of the convergence more directly. To improve the sampling efficiency, multiple independent trajectories (i.e., “walkers”) that contribute to the same bias can be simulated in parallel.[41] The multiple-walker well-tempered metadynamics protocol employed in this study involved 15 walkers and resulted in a total simulation time of ∼3.6 μs. Initial height of the Gaussian bias contributions of 0.8 kJ/mol, a deposition rate of 5 ps, and a bias factor of 15 were used for these simulations. The sampling was biased along two CVs: CV1 represented the distance between the centers of mass of the BMS-986187 heavy atoms and the TM helical bundle of the δ-receptor, and CV2 corresponded to the total number of polar and hydrophobic contacts formed at any time between BMS-986187 and the δ receptor. More specifically, this was defined aswhere r is the distance between the center-of-mass (COM) of groups i and j, and r0 was set to 5.0 Å. The receptor polar group comprises all the side-chain heavy atoms of charged and polar residues in the binding funnel; the receptor hydrophobic group comprises all the side-chain heavy atoms of hydrophobic residues in the binding funnel; ligand polar groups comprise each of the ligand’s four oxygen atoms (see Figure S9); the ligand hydrophobic groups comprise heavy atoms in the methyl substituents on the fused tricyclic ring (defined as the COM of atoms C20–C21 for one group, and COM of C22–C23 for the other, see Figure S9), as well each of the of two six-membered rings (defined as the COM of atoms C2 and C5 for one ring, and COM of C29 and C31 for the other, see Figure S9). Gaussian widths of 0.125 Å and 0.15 were selected for CV1 and CV2, respectively, based on inspection of the initial dynamics of the system. Lower and upper limits of 8 and 47 Å, respectively, for the values of CV1 were enforced using steep harmonic potentials with an elastic constant of 4000 kJ/nm2. To avoid the ligand sampling of the membrane region, additional harmonic restraint (elastic constant of 4000 kJ/nm2) was enforced on the XY component of CV1, with an upper limit of 32 Å.

Clustering Based on Interaction Fingerprints

Because of the applied bias, the simulation does not sample the Boltzmann ensemble and the trajectory must therefore be reweighted before further analysis. To this end, the unbiasing technique described in ref (25) was used to calculate the free-energy of microstates defined by the values of four CVs optimally describing the binding pose of the ligand in the allosteric pocket. Specifically, in addition to the two CVs described above (CV1 and CV2), CV3 was defined as the Z component of the vector linking the center of mass of the BMS-986187 tricyclic moiety (carbons C7–C19 in Figure S9) with the center of mass of the ortho-substituted benzyl ring (carbons C25–30 in Figure S9), and CV4 was defined as the XY component of CV1. The two additional CVs were introduced to discriminate poses with similar values of CV1 and CV2, but different orientations and positions in the helix bundle.

Microstates were defined by dividing the range of each CV into 75 bins. Microstates with free-energy below 5 kJ/mol (∼2 kBT) were considered for further analysis, and the binding pose in each microstate s was described by the probability δ(s) of any specific contact i being formed in the microstate itself. Any interaction between the ligands and the receptor, or between SNC-80 and BMS-986187 ligands, was considered and classified either as direct contact (further specified as hydrophobic interaction; H-bond with the ligand as a donor, H-bond with the ligand as acceptor; aromatic π-cation, edge-to-face or face-to-face), or a water-mediated interaction (i.e., one water molecule simultaneously forming H-bonds with one ligand and the protein, or with both ligands, simultaneously). Direct interactions with the receptor were further classified based on whether the backbone or the side chain of a specific residue was involved in the interaction. Microstates were clustered into macrostates using the dissimilarity between fingerprints defined asand a Ward agglomerative algorithm. The total number of clusters, which needs to be specified by the user, was determined so that the average Tanimoto coefficient (that ranges from 0 to 1, with 1 being the highest similarity) between clusters is less than 0.5. The free-energy of each cluster α was calculated aswhere the integration is extended to all the microstates s comprising cluster α, kBT is the thermal energy, Z is the partition function, and t is the simulation time.

Convergence of the free-energy estimates was assessed by plotting the difference of the free-energy between 75% and 100% of the simulation and checking that the difference was less than 1 kJ/mol in absolute value in the region of the minima.

After extracting medoid structures for each cluster, the probability of any specific contact being formed in each macrostate was obtained as a weighted average of the microstate probabilities, that is

Contacts formed with probability higher than 20% are reported in Table S1.

Cell Culture and Receptor Mutagenesis

Mutation of the 3xHA-hDOR sequence was achieved using the QuikChange Site-Directed Mutagenesis kit (Agilent Technologies, La Jolla, CA) following the manufacturer’s instructions. All mutations were confirmed by DNA sequencing (AGRF, Melbourne, Australia). Mutant 3xHA-hDOR pEF5/FRT/V5-DEST constructs were transfected into FlpIn CHO cells (Life Technologies, Melbourne, Australia) and selected using hygromycin B (Roche, Sydney, Australia) for stable expression. Cells were maintained in Dulbecco’s modified Eagle medium (DMEM, Life Technologies) with 5% v/v FBS and 700 μg/mL hygromycin B in a humidified incubator with 5% CO2 at 37 °C.

Whole Cell Radioligand Binding

Cells were seeded at 50 000 per well in a 96-well Isoplate (PerkinElmer, Melbourne, Australia), allowed to adhere for 6 h, and then serum starved overnight. Plates were washed once with ice-cold assay buffer (146 mM NaCl, 10 mM d-glucose, 5 mM KCl, 1 mM MgSO4, 2 mM CaCl2, 1.5 mM NaHCO3, 10 mM HEPES, pH 7.4). Cells were incubated with increasing concentrations of SNC-80 (in the absence or presence of increasing concentrations of allosteric modulator) for 4 h at 4 °C in the presence of 0.3 nM [3H]-diprenorphine (PerkinElmer, specific activity 36.1 Ci/mmol). Nonspecific binding was determined by the coaddition of 100 μM naloxone. After washing in cold saline, cells were solubilized in Optiphase scintillant, and radioactivity was measured in a MicroBeta counter (PerkinElmer).

ELISA

Cells were seeded at 125 000 per well in a 48-well cell culture plate and allowed to adhere overnight. Plates were washed with Tris-buffered saline (TBS; 50 mM Tris-HCl, 150 mM NaCl, pH 7.5), fixed with 4% w/v paraformaldehyde for 30 min at RT, and incubated with blocking buffer (1% w/v skim milk powder in 100 mM NaHCO3, pH 8.6) for 4 h at RT. Surface Human influenza hemagglutinin (HA)-tagged receptors were detected using the HA-7 mouse anti-HA antibody (1:1000, Sigma-Aldrich, Sydney, Australia), followed by HRP-conjugated goat antimouse IgG (1:2000, Sigma-Aldrich). After washing with TBS, the peroxidase substrate SIGMAFAST OPD was added, and the reaction was terminated by the addition of 1 M HCl. The colored reaction product was detected at 490 nm in a multilabel plate reader (EnVision, PerkinElmer). The absorbance values for stably expressing cells were normalized to those of untransfected cells.

ERK1/2 Phosphorylation

Cells were seeded into transparent 96-well plates at 50 000 per well, allowed to adhere for 6 h, and then serum starved overnight. Previous studies have shown that maximal stimulation of DOR by SNC-80 is achieved after 5 min.[13] Thus, cells were stimulated with ligands for 5 min at 37 °C in 5% CO2. For interaction studies with BMS-986187, increasing concentrations of SNC-80 and allosteric ligand were added simultaneously. The reaction was terminated by removal of media and ligands, and the samples were processed using the AlphaScreen SureFire p-ERK1/2 kit (PerkinElmer) as per the manufacturer’s instructions. The fluorescence signal was measured using a Fusion-α plate reader (PerkinElmer). Data were normalized to the maximal response elicited by 10% v/v FBS at the same time point.

Data Analysis

All data were analyzed using Prism 6.0f (GraphPad Software, San Diego, CA). Competition binding curves between [3H]-diprenorphine and an unlabeled ligand were fitted to a one-site binding model.[42] Binding interaction studies with allosteric ligands were fitted to the following allosteric ternary complex model, eq :[43]where Y is percentage (vehicle control) binding; Bmax is the total number of receptors; [A], [B], and [I] are the concentrations of radioligand, allosteric modulator, and orthosteric ligand, respectively; and K, K, and K are the equilibrium dissociation constants of the radioligand, allosteric modulator, and orthosteric ligand, respectively. α′ and α are the binding cooperativities between the allosteric modulator and radioligand and the allosteric ligand and orthosteric ligand, respectively. For the WT, Y56(1.39)A, Q105(2.60)A, K108(2.63)A/N, Y109(2.64)A, W284(6.58)K, and L300(7.35)W, no modulatory effect upon SNC-80 affinity was observed; consequently, log α was constrained to 0.

SNC-80 concentration response curves were analyzed using a three-parameter logistic equation as previously described.[44]where “Bottom” represents the E response value in the absence of SNC-80, “Top” represents the maximal stimulation in the presence of SNC-80, [A] is the molar concentration of the ligand, and EC50 represents the molar concentration of ligand required to generate a 50% response between minimal and maximal receptor activation.

Concentration–response curves for the interaction between the allosteric and orthosteric ligand in the pERK1/2 functional assay were globally fitted to the operational model of allosterism and agonism, eq :[29]where Emax is the maximum possible cellular response, [A] and [B] are the concentrations of orthosteric and allosteric ligands, respectively, K and K are the equilibrium dissociation constant of the orthosteric and allosteric ligands, respectively, τA and τB are operational measures of orthosteric and allosteric ligand efficacy, respectively, α is the binding cooperativity parameter between the orthosteric and allosteric ligand, and β denotes the functional cooperativity between the orthosteric and allosteric ligand. n is a transducer slope factor linking occupancy to response, for which no significant variation from unity was observed.

All affinity, potency, and cooperativity values were estimated as logarithms,[45] and statistical comparisons between values were by one-way analysis of variance followed by Dunnett’s multiple comparison post hoc test to determine significant differences between mutant receptors and the WT 3xHA-hDOR. A value of p < 0.05 was considered statistically significant.