|Home | About | Journals | Submit | Contact Us | Français|
Several G protein-coupled receptors (GPCRs), including opioid receptors δOR, μOR, and κOR, have been reported to form stable dimers or oligomers in lipid bilayers and cell membranes. This notion has been recently challenged by imaging data supporting a transient nature of GPCR association. Here we use umbrella sampling reconstructed free energies of δOR homodimers involving the fourth transmembrane helix to predict their association constant. The results of these simulations, combined with estimates of diffusion-limited association rates, suggest a short lifetime for δOR homodimers in the membrane, in agreement with recent trends.
Experimental studies using co-immunoprecipitation and bioluminescence resonance energy transfer suggest that the δ-opioid receptor (δOR)1 physically associates with itself, and with other members of the opioid receptor family (1−3). These methods, however, do not have the capability to determine the stability and mobility of these or other G protein-coupled receptor (GPCR) dimers or oligomers in the cell membrane.
A recent single-molecule study, using internal reflection fluorescence microscopy, has shown that it is possible to track the position of large numbers of individual molecules of a typical family A GPCR, the M1 muscarinic acetylcholine receptor, in living cells, over a period of several seconds (4). The results of this study point to a transient formation of M1 receptor dimers (lifetime of ~0.5 s at 23 °C) and to an estimate of ~30% total receptor molecules present in the cell as dimers. Transient association of GPCRs was also the main conclusion of recent fluorescence recovery after photobleaching (FRAP) studies of β1-adrenoceptors (5) and dopamine D2 receptors (6), wherein antibody-cross-linked receptors did not immobilize associated protomers, unless the association was established covalently by oxidative cross-linking.
The ability of GPCR dimers and oligomers to associate and dissociate rapidly suggests relatively small standard free energy differences between dimeric and monomeric GPCRs compared to those of protein complexes stabilized by multiple specific bonds. However, the nature of the interaction, whether transient or long-lasting, is unknown for the majority of GPCR dimers and oligomers, including opioid receptor complexes. Equally important is to understand the effect of the receptor sequence at the dimerization interface on the association−dissociation rate of the complex.
Here, we performed coarse-grained (CG) (7) umbrella sampling molecular dynamics (MD) simulations (8) of mouse δOR in an explicit palmitoyloleoylphosphatidylcholine (POPC) and 10% cholesterol/water environment to obtain estimates of the free energy difference between δOR monomers and homodimers involving the fourth transmembrane (TM) helix, from which to derive lifetime predictions. The focus on TM4 is justified by a number of published experimental studies on several GPCRs (9−20), including our early correlated mutation analyses of opioid receptor sequences (21,22), suggesting a direct primary association of lipid-exposed surfaces of these helices.
The TM region of mouse δOR was built by homology modeling with Modeler 9v3 (23), using the X-ray crystal structure of β2-adrenergic receptor (β2AR) at 2.4 Å resolution (Protein Data Bank entry 2RH1) as a structural template (24), and the β2AR-δOR sequence alignment deposited in the GPCR database (25), which is based on highly conserved functional residues in the TM segments. Specifically, the TM region of mouse δOR was defined by residues 461.29−761.59, 832.38−1112.66, 1183.22−1523.56, 1624.39−1854.62, 2105.35−2475.71, 2506.24−2876.61, and 2987.33−3187.53, with the superscript indicating the Ballesteros−Weinstein generic numbering scheme (26). δOR intracellular loops 1−3 (IL1−3, respectively) and extracellular loops 1−3 (EL1−3, respectively) were generated using the enhanced ab initio loop prediction approach implemented in the Rosetta 2.2 code (27). The protein N-terminus (residues 1−45) and C-terminus (residues 335−372) were not included in the receptor model. Initial configurations of δOR homodimers interacting across the TM4 interface were generated by manually positioning the protomers in a configuration compatible with symmetric contacts between residues at position 4.58.
The resulting δOR models were coarse grained (CG) according to the prescription of the Martini force field (7,28). An additional set of elastic potentials was included for beads lying within a distance cutoff (dCut) of each other in the receptor models, following the elastic network approach recently implemented by Periole and colleagues (29). Values for dCut and the elastic constants for helical and random coil regions were determined by comparison of the residue fluctuations of a monomeric δOR simulated for 50 ns using the Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA) force field in an explicit POPC/10% cholesterol bilayer with the same quantity obtained from a 50 ns simulation of the CG model. An extensive parameter search allowed us to fix the values as follows: dCut = 0.9 nm, kH = 1000 kJ mol−1 nm−2, and kL = 250 kJ mol−1 nm−2. In a slight deviation from the method of Periole and colleagues, the strength of the force constant of the elastic network was determined by the secondary structure of each of the receptor residues. If the residue was determined to have a defined secondary structure [by DSSP (30)], e.g., α-helix (including α- and 310-helix), as in the case of the TM regions of δOR, then a force constant of 1000 kJ mol−1 nm−2 was applied. For a sequence of more than two residues with undefined secondary structure (i.e., coil, bend, hydrogen bonded turn, or other undefined structure), a force constant of 250 kJ mol−1 nm−2 was applied. This allowed the secondary structure of the helices in the receptor to be maintained without compromising the flexibility of the loop regions.
A large membrane patch of 523 POPC lipids and 10% cholesterol [also described using the Martini prescription (7,28)] that could accommodate the dimeric configuration of δOR was equilibrated for 100 ns, prior to insertion of the dimer into the membrane following the protocol described in ref (31). Specifically, this protocol consists of subsequent compression and equilibration steps of the lipids following an initial expansion of the membrane in the x−y plane. The embedded receptor/membrane system was then hydrated, and counterions were added to neutralize the total charge. The final system (11.9 nm × 11.9 nm × 8.2 nm in size) consisting of 11052 beads was first equilibrated by progressively releasing constraints on the protein backbone beads and then by a final unconstrained run of 50 ns.
The reaction coordinates used to describe the relative position of interacting δOR protomers A and B during simulation are depicted in Figure Figure1.1. Specifically, the relative position of the two protomers was described by (i) the distance r between centers of mass CA and CB of the two TM regions; (ii) the rotational angle θA, defined by the projection onto the plane of the membrane of the center of mass of TM4A (defined by residues 1624.39−1854.62) CA of the protomer bearing TM4, and CB of the second protomer; and (iii) the equivalent rotational angle θB. To limit the exploration to the TM4 interface centered at position 4.58, the two rotational angles θA and θB were constrained to explore an ~25° interval around their centers by steep repulsive potentials. This interval was selected to allow the system to explore alternative homodimeric configurations of δOR exhibiting symmetric contacts between residues at position 4.58.
To take full advantage of the size of the simulation box, the x−y projection of center of mass CA of protomer A was kept fixed to its position, while the projection of center of mass CB of the second protomer was allowed to move along the line joining CA and CB, which is (approximately) parallel to the diagonal of the membrane x−y plane. These and the other constraints on the system were imposed using the Plumed plug-in (32) to the Gromacs 4.0.5 code (33).
During the simulations, the pressure was controlled with a semi-isotropic Berendsen barostat with a compressibility of 4.50 × 10−5 and a reference value of 1.0 atm. The temperature was controlled with the V-rescale thermostat to average 300 K, and a time step of 0.02 ps was used throughout the simulation.
For the umbrella sampling runs, 43 starting points at different distances (one point every 0.05 Å between 3.00 and 4.90 Å; to achieve uniform overlap of the probability distributions, six additional points were chosen at 3.125, 3.175, 3.225, 3.275, 3.325, and 3.775 Å) were extracted from unpublished CG well-tempered metadynamics simulations of δOR homodimers with TM4 at the interface, in which a bias potential had been applied to the distance between protomers while rotational angles were constrained to remain within an ~25° interval. These structures were equilibrated for 50 ns and simulated for additional 250 ns, harvesting the value of the distance each 2 ps. The resulting histograms showed continuous overlap, and the unbiased probability was reconstructed using the Weighted Histogram Analysis Method (WHAM) (34) and the WHAM code from the Grossfield Lab (University of Rochester, Rochester, NY) to calculate the free energy as a function of the distance around the transition states. To assess the accuracy of the calculated free energy, we performed commitment analysis (35) on the identified transition states at r values of 3.28 and 3.75. Specifically, we started 100 independent unbiased simulations from points at different distances, and verified that each transition state corresponded to equiprobable commitment to the adjacent basins.
Estimates of the dimerization constant for the symmetric TM4 interface of δOR homodimers in the lipid bilayer were obtained using an approach pioneered by Roux (36−38). In line with the formalism reported in the Supporting Information of ref (37), and using coordinates describing the interacting protomers, the dimerization constant was expressed as a function of the unconstrained free energy W, according to the following equation:
where r is the distance between the interacting protomers, Ω = (θA,θB) are the rotational angles describing their relative orientation (see Figure Figure1),1), (rM,ΩM) refer to a reference monomeric conformation, kB is the Boltzmann constant, T is the temperature, and the integration is restricted to dimeric states (D). When the relative orientation between protomers is constrained to a region Ω0, the free energy w(r) is given by
where the integral is limited to the constrained region Ω0 and C is an arbitrary constant. As in the work reported in ref (37), we define the w(r) to be zero in the bulk (i.e., at rM = 4.90 nm), such that the constant C in the equation above can be expressed as
where the integral is extended up to the maximum distance rD allowed for dimeric states.
For accurate calculation of a thermodynamically meaningful standard free energy (on the mole fraction scale) of δOR dimerization within the lipid bilayer, we applied the formalism described in detail in ref (39). Specifically, the binding free energy is expressed by the following equation:
where ΔGX° is the mole fraction standard state free energy change, R is the universal gas constant, T is the temperature in Kelvin, and KX is the association equilibrium constant on the mole fraction concentration scale. Thus, following the formalism described in ref (39), we expressed this equilibrium constant as a function of the mole fractions of dimeric (ND/NTot) and monomeric (NM/NTot) species in the membrane phase, according to the following equation:
Because the lipid concentration is much greater than the protein concentration, we can approximate the total number of molecules in the membrane phase of area A (NTot) with the number of lipids (NL) and calculate the KX according to the following equation:
is the dimerization constant.
We calculated the dissociation rate koff of δOR dimers involving the TM4 interface, according to the following equation:
where KD is the dimerization constant and kon is the association rate. The latter was approximated using its diffusion-limited value, following the Smoulchowski theory in two dimensions (40). Specifically, at long times, this rate is given by the following expression:
where Dc = DA + DB = 2DT is the sum of the diffusion constants of the two protomers A and B in the δOR dimer, R = 2R0 is the sum of the protein radii, γ is the Euler−Mascheroni constant, and t refers to typical experimental time scales explored to detect diffusion. If one starts from an initial value [D]0, for short times, the concentration of δOR dimers over time is given by the equation
Thus, the half-time of δOR dimers in the lipid bilayer can be calculated as
We recovered the atomic details of the structures representing metastable dimeric states of δOR by applying the protocol proposed in ref (41), which uses a simulated annealing algorithm. The resulting all-atom structures were then embedded into a pre-equilibrated all-atom membrane and solvated in a fashion similar to that described above for the CG systems. These systems were energy minimized and equilibrated by performing cycles of MD under successively relaxed position restraints (1000, 100, 10, and 0 kJ mol−1 nm−2) for a few picoseconds. After the embedding, the system was equilibrated with a 10 ns unbiased MD with the OPLS-AA force field, and contact maps were obtained by averaging the Cβ distance matrix over the equilibration run.
Estimates of the free energy difference between δOR monomers and homodimers involving TM4 interfaces, within an explicit POPC/10% cholesterol/water environment, were derived from CG (7) umbrella sampling MD simulations, according to the protocol described in Materials and Methods. Figure Figure22 shows the free energy w of the constrained dimeric system involving the TM4 interface, reconstructed as a function of r. Two different dimeric states of δOR (D1 and D2 in Figure Figure2)2) involving the TM4 interface were identified, which are separated from each other by a transition state at rTS1 = 3.28 nm, and from the monomeric state (r ≥ 4.90 nm) by a transition state at rTS2 = 3.75 nm. The insets in Figure Figure22 show the results of the commitment analysis (35) conducted on the two identified transition states to assess the accuracy of the calculated free energy. As these plots show, each transition state corresponded to equiprobable commitment to the two adjacent basins, leaning in favor of the results.
The two dimeric states, D1 and D2, exhibited a similar free energy, although D1 appears to be slightly more stable than D2 (by ~1 kcal/mol). These states are also very similar in structural terms. Figure Figure33 shows, as an example, a representative structure from the D1 basin, corresponding to a symmetric (average θA 20° and θB 23°) tight (r = 3.15 nm) δOR homodimer. As shown in this figure, TM4 from one protomer inserts into a groove on the opposite protomer formed by helix TM2, the C-terminal half of TM3, and TM4. Intracellular loop 2 (IL2) is also close to this interface. A list of the TM4 residues forming symmetric intermolecular contacts is presented in Table Table1.1. The structures in basin D2 (see a representative in Figure Figure4)4) are similar in the overall orientation of the protomers but correspond to a less compact (r = 3.40 nm), slightly asymmetrical interface, where θA 20° and θB 15°. Figure Figure44 also shows a comparison of the relative position of the TM4 helices of the two protomers in the two different conformations, D1 and D2, after alignment of one of the two protomers.
Using rD = rTS2 = 3.75 nm, the value of the free energy w(r) calculated with eq 4 of Materials and Methods from the curve reported in Figure Figure2,2, and the product of the allowed ranges for θA and θB in radians |Ω0| = 0.20, we obtain a KD of 1.02 μm2 for the δOR homodimer with TM4 at the interface. This value, added to eq 5 of Materials and Methods together with the surface density of the POPC/cholesterol patch used in the simulations (NL/A ~ 1.65 × 106 μm−2), allows us to calculate a standard state free energy change for the δOR homodimer with the TM4 interface [ΔGX° = −RT ln(1.68 × 106) = −8.54 kcal/mol].
Using a diffusion coefficient DT value of 0.08 μm2/s determined experimentally for μOR (42), we calculated the sum of the diffusion constants of the two protomers A and B in the δOR homodimers; i.e., Dc = DA + DB = 2DT ~ 0.16 μm2/s. Placing this value into eq 10 of Materials and Methods, together with the sum of the protein radii (R = 2R0 ~ 3 nm), the Euler−Mascheroni constant (γ 0.57), and using a mean value over the time interval of the experimental time scale explored to detect MOR diffusion [i.e., 2 ms ≤ t ≤ 30 s (42)], we obtained a mean association rate kon of 0.16 μm2/s. Using the calculated dimerization constants (KD) in combination with the experimental diffusion coefficient DT value of 0.08 μm2/s reported in ref (42), we calculated dissociation rates (koff) of 0.15 s−1 from eq 9 of Materials and Methods. Placing these values into eq 12, we calculated a half-time of 4.4 s.
The physical and chemical basis for opioid receptor interactions in the membrane is currently unknown. Here, we conducted CG umbrella sampling MD simulations to obtain theoretical estimates of the thermodynamics and kinetics of the dimerization of δOR involving the TM4 region, within an explicit POPC/10% cholesterol/water environment. These estimates were used to calculate theoretical dimerization constants that allowed us to obtain standard free energy of dimerization values on the mole fraction scale. According to Fleming (39), the latter corresponds to a thermodynamically meaningful free energy for the protein association within a membrane, since it can apply equivalently in all experimental systems. Notably, the calculated mole fraction standard state free energy of δOR dimerization (−8.54 kcal/mol) is very close to the corresponding experimental value of −7.0 kcal/mol for the association of glycophorin A transmembrane helices in C8E5 micelles at 25 °C (39). Using the total protein concentration of 1 μm−2 that resulted from the single-molecule study of M1 muscarinic receptors (42), the free energy difference between dimeric and monomeric δOR is much smaller (ΔG° 0.02 kcal/mol). If 1 μm−2 is a typical total protein concentration value, we could use it in combination with the calculated KD to calculate the mole fraction of dimers, fD. Thus, placing 2[D] = fDC° and [M] = (1 − fD)C° into the equation KD = [D]/[M]2 = fD/[2(1 − fD)2C°], we would obtain an fD of 0.50. Though not identical to the value recently proposed experimentally for M1 muscarinic receptors (30%) (4), this calculated value is remarkably close.
In combination with a diffusion coefficient determined experimentally for μOR (42), which is in line with values of 0.1 μm2/s determined experimentally for several other GPCRs, including rhodopsin (43), β-adrenoceptors (44−46), and M1 muscarinic receptor (4), our theoretical estimates of the dimerization constants of δOR homodimers involving TM4 helices allowed us to calculate short lifetimes of δOR dimers in the membrane (4.4 s), in line with the rapid association and dissociation of M1 muscarinic receptors assessed by single-molecule studies (4). Notably, the half-time of DOR would be even shorter (1.8 s) if we used the experimental time scale explored to detect diffusion of the M1 muscarinic receptors (i.e., 50 ms ≤ t ≤ 5 s) (4). Whether this suggested short lifetimes of δOR homodimers within the membrane have implications for the functional role of these receptor complexes and/or the specificity of their interactions remains to be established. Of particular interest is also the examination of the effect of different protein sequences and/or hydrophobic environments on the association rate of GPCRs, which are currently under investigation in our laboratory.
National Institutes of Health, United States
†This work was supported by National Institutes of Health Grants DA020032 and DA026434. The computations were supported in part by the National Science Foundation through TeraGrid advanced computing resources provided by TRAC MCB080077.