Structure-Function Relationships in the delta, kappa and mu Opioid Receptors: Modeling and Molecular Dynamics Simulations of the Receptor Structures

We have developed 3-D structural models of the opioid receptors (delta, kappa and mu), members of the G-protein coupled receptor (GPCR) superfamily using techniques of homology modeling ¹. Predictions of the structural elements are based on A) the consensus indicated by average properties (such as hydrophobicity, variability and conservation); B) the structural periodicities of these properties evaluated from a multiple sequence alignment of opioid, somatostatin, angiotensin and related receptors; and C) the resemblance of the tertiary structure of the opioid receptors to rhodopsin. The initial model structures of the three opioid receptor subtypes were refined with backbone constraints during energy minimization and molecular dynamics; these constraints were released in steps over 210 picoseconds. These relaxed structures were then simulated for an additional 1.8 nanoseconds at 300 °K.

The initial structures were generated using identical methodology, resulting in a backbone root mean square deviation (RMSD) between the minimized models of the various subtypes of < 2.0 Å. Several average structures were generated from the 1.8 nanosecond trajectory for each receptor subtype; these average structures have a backbone RMSD from the initial minimized structure of between 3 and 4 Å ; the range of inter-subtype backbone RMSD from comparisons among these average structures indicates the stability and similarity of these structures. During the simulations, structural elements that are common to the receptor subtypes and presumed to form the ligand binding site have diverged during the trajectories, resulting in different constructs that incorporate similar functional groups. These differences are likely to be responsible for differences in ligand-receptor complexes, highlighting the effect of the receptor micro-environment upon the binding site. In contrast, a common network of conserved interactions is observed in the intra-cellular portion of the simulated transmembrane helices in all three receptor subtype models, possibly indicating a core set of interactions related to a common process of signal transduction.

I. Introduction

Opiates are a large class of neurotransmitters and substances, mediating various behaviors and responses such as analgesia and euphoria; in clinical settings, these compounds are used for pain management, but often have undesirable side effects ^{2 3 4}. The diverse activities of various opiates led to the initial suggestion that there were divergent receptors for these compounds ⁵. These initial studies were initially extended and verified through the development of subtype-specific opiate compounds and their use in autoradiographical studies, thereby establishing clear differences in the pharmacological profiles of these receptor subtypes ⁶. The recent successful culmination of this research effort has been the cloning of the several opiate receptor subtypes and the discovery of new, closely related receptors and their endogenous ligands ^{7 8 9}.

The cloned opiate receptors are members of the heptahelical, G-protein coupled receptor family. Sequence alignments between individual opioid receptor subtypes confirm their relatedness, showing a high degree of sequence identity and homology (on the order of 61% and 75%, respectively). Moreover, the sequence identity in opiate receptor clones from different mammalian species generally exceeds 93%. The observed sequence differences between the opiate receptor subtypes are clustered in regions believed to fall within extra-cellular loops connecting the seven putative transmembrane helices and within the amino and carboxyl termini. Site-specific mutagenesis studies exploring conserved and divergent residues have demonstrated that both peptidyl and non-peptidyl opiate ligands form specific interactions in the transmembrane region, suggesting that a recognition domain common to all opiate ligands may lie within this region ^{10 11 12 13 14 15}. This concept of a recognition domain within the transmembrane region that participates in the binding of opiate ligands and which may cause a conformational change of the receptor leading to signal transduction is quite broadly related to the observed interaction of catecholamines within the transmembrane region ¹⁶. However, these same studies have also demonstrated that the known descriptions of subtype diversity can be accounted for by the transmembrane domains ^{10 11 12 13 14 15}.

We have constructed molecular models of the transmembrane region of the opiate receptors which is likely to contain key elements of the binding site for small non-peptidyl opiate ligands. The molecular models for each of the three major subtypes of opiate receptors, delta, kappa and mu, were based upon techniques of homology modeling, following the low-resolution three-dimensional structures of frog and bovine rhodopsin ^{17 18}. Molecular dynamics simulations of the final models showed them to be stable for very long simulation times, exceeding a total length of 2 nanoseconds. We have used these models to explore comparative features of the structures of the different opiate receptors and their ligand binding site. The closely related structural framework from which each model of the opiate receptor subtypes have been developed allows a direct structural comparison of both the dynamic and average properties of the models. This leads to the development of a description of the interactions both common to each subtype and unique to the subtype specific models, permitting a connection at a molecular level to the study of opiate receptor pharmacology.

II. Methods

Construction of models: We have used a number of integrated techniques of homology modeling ¹ to construct molecular models of the delta, kappa and mu opioid receptors. A sequence alignment of opioid and related peptidyl G-protein coupled receptors (GPCRs) was constructed that included receptors for somatostatin, angiotensin, cholescystokinin, tachykinin, and interleukin-8. In the multiple alignment, the receptor sequences were adjusted to bring known highly-conserved residues within the transmembrane regions into alignment without allowing the introduction of insertions or deletions into the putative transmembrane helical regions. The multiple sequence alignment produced in the first step was used to calculate average property profiles of the sequence such as hydrophobicity, variability and conservation. These sequence-dependent property profiles were analyzed using Fourier transforms to generate an alpha-helical periodicity (AP) index. Consensus predictions from the AP index for each of the analyzed properties were used in the construction of the molecular models to rotate and orient the transmembrane alpha helices (Table 1).

Table 1 showing definitions of transmembrane helices

Molecular models of the transmembrane helix bundle for the opioid receptors without bound ligand were assembled on a rhodopsin-based template derived from Baldwin's analysis of conserved GPCR sequences using ideal helices, proline kink geometries and common, representative side chain rotamers ^{1 19 20 21 22}. Side chains were rotated only to introduce interactions, such as the established interaction between the conserved Asp 2.50 in Helix 2 and the conserved Asn 7.49 in Helix 7 or to relieve bad steric contacts between helices ²³.

In the final step of model construction, the structures were subjected to energy minimization and molecular dynamics (MD) using the CHARMm 23f3/24b2 package ²⁴. A distance-dependent dielectric constant of epsilon equal to 1/r was used throughout these simulations. The model for each receptor subtype (delta, mu, and kappa) was constructed using unified atom parameters including polar hydrogens. In the initial stage of the energy minimization, the backbone atoms of the models were held fixed during a minimization that used both steepest descents and adopted basis Newton-Raphson algorithms. The minimized models were imported into molecular dynamics simulations, using a Verlet algorithm, a step size of 1 femtosecond and updating nonbond interactions at 10 step intervals. The receptor models were gradually heated to 300 °K, using strong constraints on the backbone phi-psi angles. These constraints were gradually released over the first 210 picoseconds, producing "relaxed" structures for continued MD simulation. These trajectories were continued for another 1.8 nanoseconds after the constraints had been removed. Two-dimensional RMSD plots of the trajectory were used to monitor convergence and as a guide to extracting average, representative structures from homologous segments of the trajectory. Average structures were minimized for 500 steps using an adopted basis Newton-Raphson procedure.

Analysis of trajectories: In general, utilities present in the CHARMm 23f3/24b2 package were used for extraction of parameters. Hydrogen bonds between sidechains were determined by monitoring the distances between polar hydrogen donors and acceptors and the angles between donor antecedents, donor hydrogens and acceptors; if a distance was less then 3.0 Å and the acceptor angle was within the range 180° ± 30°, then a hydrogen bond was recorded. Stacking interactions were deduced using an algorithm developed for analyzing pi-pi interactions ²⁵. Electrostatic potentials were calculated using the Delphi software package; the parameters used included a salt concentration of 0.145 M, iterative focusing to develop high resolution grid mappings of the electrostatics and interior/exterior dielectrics of 2 and 80 respectively ²⁶.

III. Results

Stability of model structures: The opioid receptor models were quite stable during MD simulation. The structure of the models for the three receptor types clearly evolved and diverged during the course of the MD simulations; by the end of the simulation, the backbone of the simulated models developed an average RMSD greater than 3.5 Å from the initial models, even though the initial models were constructed using similar criteria and energy minimized to share backbone position deviations of less than 2 Å. This large structural deviation resulting from the development of alternative structures along each trajectory underscores a pertinent property of these MD simulations: the ability of the model structures, under the appropriate simulation conditions, to evolve into alternative structures with isopotential energy values; this property is largely absent from the energy minimization techniques applied to the construction of GPCR models and, in this case, would have been misleading about the structures available to each receptor model (see for example, ²⁷).

The RMSD for the comparisons of backbone atoms between average structures derived from 1.8 nanoseconds of trajectory (Table 2) is generally less than 2.0 Å, suggesting that the structures generated within a trajectory have stabilized (Figure 1).

Figure 1: Two dimensional root mean square deviation (RMSD) plot of kappa opioid receptor model. The trajectory shown consists of structures sampled every 5 picoseconds over a length of 1.98 nanoseconds. The upper left triangle is a comparison of the structures RMSD when considering all atoms; the lower right triangle is a comparison of the structures RMSD when considering only backbone heavy atoms. The color scale for interpreting the plot is shown at right.

The RMSD between the comparisons of the averaged structures of the subtype models lies between 2.7 and 4.2 Å, illustrating a second property of these MD simulations: the intrinsic ability of sequence-specific variation in the structural models to overcome the biases imposed by homology requirements. This property of the simulations has allowed the development of divergent model structures from a common initial topology, dictated solely by the intrinsic energy requirements of these structures. Although these opioid receptor models possess divergent structural qualities, it is of especial interest that these models retain faithfully the known, incorporated properties of G-protein coupled receptor structure.

General features of model structures: The model structures consist of seven alpha helices, extracted from the cloned sequences of the delta, mu and kappa opioid receptors (Figure 2). The structures of the opioid receptors have been modeled using a counter-clockwise orientation of the seven helices (when viewed extracellularly) since this is presumed to lead to a more probable arrangement of the model ²⁰.

Figure 2: Backbones of delta, kappa and mu opioid receptor models. The structures shown are average, energy minimized structures from the final 400 picoseconds of the trajectory. The structures in the lower half have been rotated by 90° around the horizontal axis in the plane of the figure from the structures in the upper half. Structures have been superimposed to minimize the RMSD. The average isotropic fluctuation per residue is represented in the ribbon structures as the width of the ribbon; the width is linearly related to the fluctuations with a slope of approximately 1.0. The transmembrane helices are identified only for the delta opioid receptor; the other structures follow in a similar fashion. In the top half of the figure, the view of the models is presented from the extracellular side; in the bottom half of the figure, the extracellular side of the models is oriented towards the top of the figure.

The interhelix crossing angles between the interacting helices of the opioid receptor models fall well within the range observed in crystallographic studies (Table 2) ²⁸. Although the observed interhelical angles vary among models, none of the values falls significantly beyond the comparable range in the other receptor subtype models. Within the simulation of a single model, the widest variation in the crossing angles is seen with helix 1; this helices lie on the 'edge' of the receptor topology, sharing few interactions with the other helices and are consequently much more flexible ²⁰. Large values for the crossing angle, ranging from ~20 to ~30°, are observed for the angles that helices 1, 2, 3, 5 make with helices 4, 6 and 7, reflecting the assignment of the helices to the low resolution structures of rhodopsin and the consequent assumption of tilt in the construction of the model structures ^{1 20}.

	 Hlx 1 - 2     Hlx 1 - 7     Hlx 2 - 3      Hlx 2 -7    Hlx 3 - 4
delta	14.2 ± 4.0     8.7 ± 4.5     6.7 ± 2.3    17.0 ± 2.6    9.7 ± 1.4   
kappa	 8.3 ± 4.1    14.0 ± 4.4     9.0 ± 3.2    20.2 ± 2.1    7.9 ± 2.2   
mu	14.7 ± 4.2    24.3 ± 3.6    24.1 ± 2.4    28.2 ± 2.6    8.8 ± 4.1   

	 Hlx 3 - 7     Hlx 4 - 5     Hlx 5 - 6     Hlx 6 - 7
delta	12.5 ± 2.5     5.7 ± 1.8    16.5 ± 3.3    18.0 ± 4.7
kappa	25.3 ± 1.7    19.3 ± 2.9    10.4 ± 2.9     5.4 ± 3.0
mu	30.3 ± 1.5    16.4 ± 5.4    25.0 ± 3.6    18.4 ± 1.8

Divergent features of model structures: The divergent nature of the models is driven by the sequence variability at several levels. The individual helices of the model structures have divergent backbone topologies, driven by nonlocal sequence variation. For example, the backbone comparisons of the average structure of helix 3 from the last 400 picoseconds of the simulations indicates an RMSD on the order of 2 Å (Table 4). However, helix 3 has minimal sequence variability in all three receptors (Table 1), indicating that the structural variation in other helices has necessarily driven the structural divergence of this helix in the model. This divergence may be important since identical residues in helix 3 have been demonstrated to display different phenotypes when mutated in the different subtypes of receptors ¹². Helices 1, 4 and 6, which share the least sequence homology in the opioid receptor family, all show strong backbone diversity. In the experimental studies, subtype diversity is correlated with both mutants and chimeras of these helices, suggesting a correlation with the structures from these simulations ^{10 12 14 15}. Lower RMSDs are seen with helices 2, 5 and 7, suggesting that the sequence variation in the receptors generates only subtle local changes.

Helix 1  kappa   mu	Helix 2  kappa   mu	Helix 3  kappa   mu
delta    1.5 Å  2.0 Å	delta    1.0 Å  0.6 Å	delta    1.5 Å  2.6 Å
kappa           1.3 Å	kappa           1.0 Å	kappa           1.5 Å

Helix 4  kappa   mu	Helix 5  kappa   mu	Helix 6  kappa   mu
delta    1.4 Å  2.2 Å	delta    1.0 Å  1.8 Å	delta    2.2 Å  2.3 Å
kappa           2.6 Å	kappa           1.5 Å	kappa           2.4 Å

Helix 7  kappa   mu
delta    1.4 Å  1.6 Å
kappa           1.1 Å

A finer detailing of the structural divergence of the receptor models can be examined through plots of the backbone RMSD on the basis of individual residues (Figure 3). The lower RMSD in helices 2 and 7 is observed to be a common backbone topology stretching over several helical turns in the central portion of the receptor, while the backbone topologies at the termini of these helices is quite different. Helix 5, while retaining significant backbone homology among the different models (Table 4), is seen in the kappa and mu models to have undergone apparent rigid body motions around a central pivot, relative to the remaining receptor helices; in the delta model, this helix appears to have been both translated and rotated in the vicinity of the conserved Pro 5.50(225).

Figure 3: RMSD of superimposed structures of delta, kappa and mu opioid receptor models. The structures shown are average, energy minimized structures from the final 400 picoseconds of the trajectory. The structures in the right half have been rotated by 90° around the vertical axis in plane of the figure from the structures in the upper half. Structures have been superimposed to minimize the RMSD. The average RMSD for the backbone atoms in each residue is represented in the ribbon structures as the color of the ribbon; the blending of the colors is linearly related to the distance between the three models. The color scale used is shown at lower left. If the backbone atoms of a receptor model are more than 2 Å apart from the homologous atoms in both of the other two models, the backbone is colored with a primary color, as shown in the vertices of the color scale. If the backbone atoms of one model approach only one other model within 2 Å , the primary colors representing each model are blended as demonstrated by the centers of the sides of the color scale. If all three models lie within 2 Å of each other, all three primary colors are blended and the backbone is colored white.

Significant backbone similarity between the three models is observed in helix 3, stretching approximately from residue 3.29 to 3.39; this indicates that the divergence of the backbone RMSD observed in Table 4 is concentrated at the extracellular and intracellular termini of this helix. Helix 4 of delta and mu is seen to share a common backbone topology through the central portion of this helix; the structural variation in this helix is quite visible at the amino and carboxyl termini. Helix 1 in delta and mu is seen to follow similar backbone paths in a central turn; kappa follows a different path than either of the other two models, although it is structurally related to both (Table 4). Helix 6 simultaneously shows both structural integrity and diversity among the receptor subtypes; significant structural homologies are observed over half of the helix between delta and kappa (Figure 3), although the backbone of this helices are more than 2 Å distant between the model structures. In contrast, the backbone of helix 6 in mu is both structurally diverse and follows a quite distinct path at all points.

Networks of interactions in opioid receptor models: Two forms of networks are observed in the simulations of the opioid receptor models. In one, aromatic residues are observed to cluster while maintaining geometries resembling favorable stacking interactions in known protein structures (Table 5). While all of these aromatic residues are conserved among members of the opioid receptor family, others are conserved in the larger family of GPCRs. Residues conserved within the opioid receptor family that are present in this cluster are Tyr 3.33, Phe 5.43 and His 6.52. The aromatic residues that are highly conserved in most GPCRs and are present in this cluster include Trp 4.50, Phe 5.47, Phe 6.44 and Trp 6.48 (Figure 4). Five of these aromatic residues have been examined by site-directed mutagenesis in the delta opioid receptor ¹³. The most severe effects were observed for the mutation of Tyr 3.33(129) and Trp 4.50(173), separately, to alanine. In the model, these residues flank the important residue Asp 3.32(128), a probable key component of the binding site (see Electrostatics), and are positioned in the upper half of the transmembrane helices in the putative binding site, suggesting a ready explanation for the effect of these mutations. The presence of both family specific and highly conserved aromatic residues in a conserved network of interactions suggests that many of these could be elements of a conserved signal transduction pathway.

		Aromatic-aromatic Interactions observed

			Tyr 3.33 - Trp 4.50
			Tyr 3.33 - Phe 5.43
			Trp 4.50 - Phe 5.43
			Trp 4.50 - Phe 5.47
			Trp 4.50 - Trp 6.48
			Phe 5.43 - Phe 5.47
			Phe 5.47 - His 6.52
			Phe 5.47 - Trp 6.48
			Trp 6.48 - Phe 6.44
			Trp 6.48 - His 6.52

Table 5: Interactions between aromatic residues observed in opioid receptor simulations. Interactions were defined as the common intersection in the three simulations of all aromatic-aromatic interactions that were present for at least a stretch of 100 picoseconds, truncating from the list interactions that neither formed beyond one additional stacking partner nor with the main aromatic cluster. Stacking interactions were deduced by applying an algorithm developed for examining pi-pi interactions in phenylalanine to all aromatic residues in the simulations ²⁵. 3960 structures were analyzed, sampled at a spacing of 0.5 picoseconds.

Figure 4: Structure of a network of aromatic residues in the delta opioid receptor model. The structure shown is the average, energy minimized structure from the final 400 picoseconds of the delta model opioid receptor trajectory. Aromatic residues associated with a network of aromatic-aromatic interactions in the simulated models are shown. Helices 1, 2 and 7 are not shown.

A second network observed in the model simulations is formed of hydrogen bonds between residues conserved across the opioid receptor subtype. This network includes several residues that are very highly conserved in all known GPCRs; these are Asn 1.50, Asp 2.50, Trp 6.48, Asn 7.45, Ser 7.46, Asn 7.49 and Tyr 7.53. Several residues, conserved only in the opioid receptors, are also hydrogen bonded within this network; these are Asn 3.35, Ser 3.39 and Thr 3.42 (Figure 5). The interactions among Asp 2.50, Asn 7.49 and Asn 1.50 were designed into the initial model based on deductions from experimental studies ^{1 23} and persisted throughout the simulations. Interestingly, these initial modeled interactions were joined by the remaining hydrogen bond elements of the network that developed during the simulation procedure. The roles of Asp 2.50 and Trp 6.48 have been tested through mutagenesis in the mu and delta receptors ^{7 10 14}, but the remaining residue interactions have not been examined. It is noteworthy that several residues in this network, most notably Asn 7.45 and Asn 3.35, appear to play a role in connecting the hydrogen bond network in the receptor models to the binding site, thereby linking regions known to be closely involved with activation and the sodium response of the opioid receptors to structural regions known to be important for binding ^{9 10}.

Figure 5: Structure of a network of hydrogen bonded residues in the mu opioid receptor model. The structure shown is the average, energy minimized structure from the final 400 picoseconds of the delta model opioid receptor trajectory. Residues associated with a persistent network of hydrogen bonded interactions in the simulated models are shown. Hydrogen bonds are represented as pink balls at a spacing of 0.5 Å connecting polar moieties.

Electrostatic fields in opioid receptor models: Because a protonated amine group is known to be a conserved functional group of opioid ligands, the conserved residue Asp 3.32 in the opioid receptors has long been considered by analogy with catecholaminergic receptors to be the complementary binding site in the opioid receptors ^{16 29}. While mutations of Asp 3.32(147) to alanine and asparagine in the mu receptor have severely impaired ligand binding affinity, the complementary mutations of Asp 3.32(128) in the delta receptor have yielded equivocal information, suggesting alternate roles for this residue ^{10 12}. The presumed interaction site for the protonated amine functional group, Asp 3.32, is not the only negatively charged residue within the transmembrane domain of the opioid receptors. Other residues include Asp 2.50, Asp 3.49 and, in the kappa receptor, Glu 6.58(297). This latter residue has been observed to be a binding determinant in the kappa receptor for ligands such as norbinaltorphimine, suggesting that this glutamate may represent the interaction site for the second amine group of this bivalent ligand ¹¹. Analysis of the electrostatic potential fields in the opioid receptor models can help clarify the roles of these residues and their possible function in the receptors.

The electrical potential perceived by ligands arising from electrostatic charges is known to be heavily modulated by both the interior dielectric and the exterior shape of the associated protein ²⁶. The electrostatic potential generated by each of the receptor subtypes should possess both common and divergent properties which are modulated by the receptor tertiary structure. Important determinants in these properties are likely to be conserved charged residues that are buried deeply within the receptor (e.g. Asp 2.50), residues exposed near the surface that have variable conformations, (e.g. Asp 3.32) or residues representing sequence variants of a particular receptor (e.g. Glu 6.58).

The electrostatic properties of each receptor subtype model were evaluated through calculation of the electrostatic potential using the Delphi software package ²⁶. The results for various two-dimensional slices through the receptor models are shown in Figure 6. Several common features are evident for each receptor. In all three subtypes, the central portion of the receptor model that incorporates Asp 3.32 has a strong negative potential. The negative potential both deepens and widens as the section descends through the plane of the carboxylate atoms of Asp 3.32 towards the intracellular half of the receptor and weakens as the section moves towards the extracellular surface. The negative potential has a 'bubble' which extends through helix 3 and 4. The prime determinant of the electrostatic potential in the portion of the transmembrane domain towards the electrostatic side is Asp 3.32 in the delta and mu models, with the addition of Glu 6.58 in the kappa model. The low dielectric interior of the model structure acts as a 'guide' to the potential, projecting the electric field of the proximal charged residue, Asp 3.32, through helix 3 and 4. A similar effect can be observed in the calculated potential for the kappa model with Glu 6.58 and helix 6. The strengthening of the negative potential with the descent into the receptor model suggests that Asp 2.50 and Asp 3.49 are generating this strong potential. This is consistent with a path of a sodium cation penetrating into the receptor towards its probable interaction with Asp 2.50 ¹⁰.

Figure 6: Electrostatic potentials in the delta, kappa and mu opioid receptor models. The structures shown are the average, energy minimized structure from the final 400 picoseconds of the trajectories. Electrostatic contours were drawn at ±3.0 and ±0.5 kT. The backbone of the transmembrane helices are shown as stick figures. Negative potentials are drawn in red, positive potentials are drawn in blue.

The model of the kappa receptor possesses the strongest negative electrostatic potential which, at all levels, has a more widespread projection than do the negative potentials in the other subtype models, most likely due to Glu 6.58(297) which is present only in this receptor. The conformation of Asp 3.32 in the opioid receptor model simulations can also be seen as a determinant of the electrostatic potential, since the delta and mu receptor have comparable numbers of charged anionic residues, and different conformations of this aspartate. The preferred trans conformations of this residue in the mu and kappa models appears to shift the projection of the electrostatic potential towards an area including helices 1 and 7 at the plane of the putative binding site and to reduce the projections through helices 3 and 4. An interesting observation is that the potential in the region near helices 1 and 7 in the plane of Asp 3.32 in the delta model appears to possess a strong positive character, unlike the potentials in this region in the kappa and mu models. At the three levels analyzed here, the delta model appears to have a stronger negative potential than the mu model, suggesting an additional explanation for the unusual effects observed with the mutation of Asp 3.32 to alanine in the delta opioid receptor ¹².

IV. Discussion