Abstract

In an attempt to guide experimental assessment of the identity of the minimal signaling unit for CBS, we conducted extensive coarse-grained (CG) molecular dynamics (MD) simulations of different combinations of the three major OR subtypes, i.e., u-OR, 6-OR, and K-OR, in an explicit lipid bilayer. Specifically, we ran multiple, independent MD simulations of each homomeric u-OR/u-OR, 6-OR/6-OR, and K-OR/K-OR complex, as well as two of the most studied heteromeric complexes, i.e., 6-OR/u-OR and 6-OR/K-OR, to derive the preferred supramolecular organization and dimer interfaces of ORs in a cell membrane model. These simulations yielded over 250 microseconds of accumulated data, which correspond to approximately 1 millisecond of effective simulated dynamics according to established scaling factors of the CG model we employed. Analysis of these data indicates similar preferred supramolecular organization and dimer interfaces of ORs across the different receptor subtypes, but also important differences in the kinetics of receptor association at specific dimer interfaces. We also investigated the kinetic properties of interfacial lipids, and explored their possible role in modulating the rate of receptor association and in promoting the formation of filiform aggregates, thus supporting a distinctive role of the membrane in OR oligomerization and, possibly, signaling.

Author Summary

Whether or not the different opioid receptor subtypes share a similar propensity to form di-/oligomers and specific receptor-receptor interactions is still unclear on the basis of the published data. This information, however, is necessary to predict stabilizing or de-stabilizing mutations and to design experiments to clarify the role of oligomerization in opioid receptor function. The inferences provided by the extensive molecular dynamics simulations reported herein constitute a first step in this direction.

Introduction

In particular, heteromerization between u-OR and 5-OR and between 5-OR and K-OR has been suggested to expand the repertoire of OR signaling by modulating ligand binding, receptor signaling, and/or trafficking properties [2]. This modulation has a direct translational relevance in view of the putative role that OR heteromers play in opioid-induced adverse effects (e.g., the development of analgesic tolerance [3,4]), and in the observed potentiation of morphine-induced analgesia by 5-OR-selective antagonists [5]. Thus, targeting OR oligomers directly may lead to novel drugs with potentially greater selectivity and reduced side effects compared to molecules targeting individual receptors.

To begin addressing these questions and eventually use these and other molecules as tools to elucidate the physiological relevance of OR oligomers, ligand-bound models of these receptor complexes are highly desirable. However, building such models is not feasible in the absence of reliable information about their interface of dimerization/oligomerization in a physiologically relevant environment. The recent X-ray crystal structures of the u-OR [12] and K-OR [13] have suggested specific receptor-receptor interactions involving transmembrane (TM) helices TM5 and TM6 or TM1, TM2, and helix 8 (H8). Although these interfaces are thermodynamically stable according to our recent free-energy calculations of u-OR and K-OR homo-dimers in an explicit lipid-water environment [14], the possibility cannot be ruled out that other receptor-receptor interactions are also feasible in the cell membrane, and perhaps more kinetically favorable than those inferred from crystallography. Here, we carried out extensive, unbiased coarse-grained (CG) molecular dynamics (MD) simulations of freely diffusing ORs in an explicit lipid-water environment to evaluate differences and similarities in the supramolecular organization and preferred dimeric interfaces of all three major receptor subtypes.

Results and Discussion

Specifically, to enhance the statistical significance of our description of the dimerization process, we ran five independent simulations for each of the homomeric (i.e., u-OR/u-OR, 5-OR/ 5-OR, and K-OR/K-OR) and the most studied heteromeric (i.e., 5-OR/u-OR and 5-OR/K-OR) complexes of ORs, starting from randomized initial orientations of CG models of sixteen individual receptors (8+8 in the case of heteromers) in an explicit CG palmitoyl-oleoyl-phos-phatidyl-choline (POPC)/ 10% cholesterol bilayer. Together, these simulations yielded over 250 us of accumulated data, which correspond to approximately 1 ms of effective simulated dynamics, according to the established (X4) scaling factor of the CG MARTINI model [16—18] we employed.

Preferred Supramolecular Organization of ORs

As an example, Fig. 1. shows the location of each of the 16 receptor molecules in one of the five simulations

Notably, linear arrays of receptors were observed in early atomic force microscopy images of the prototypic GPCR rhodopsin in native membranes [19], and were further supported by early CG simulations [20]. Below, we speculate that the formation of chains of receptors may be influenced by lipid dynamics (see the Dynamic Behavior of Lipid Molecule section for a rationale).

This observation is in line with experimental estimates of dimer lifetime (a few seconds) obtained from recent single-molecule imaging studies of different, individually labeled, GPCRs in living cells [21]. Since dissociation is not observed in the simulations reported here, a direct calculation of free-energies is precluded. However, analysis of the executed twenty-five simulations instills confidence in that they might have captured the majority of fastest-forming dimer interfaces of ORs.

Preferred Homo-dimer and Hetero-dimer Interfaces of ORs

Since each protomer structure was kept fairly rigid by elastic network forces during the simulations (see Materials and Methods for details), no allosteric communication between protomers or interdependence between dimer interfaces were expected. Interfaces were defined based on the minimal number of residues on each receptor TM heliX within a certain distance cutoff from each other, and were clustered based on the similarity between inter-protomer contact maps (see details of the analysis under Materials and Methods).

The following interesting observations can be made on the basis of this analysis. First of all, not all possible combinations of TMs were found to be involved in dimer interfaces during the simulated timescale of 10 us. The only interfaces that formed in all studied homo and heteromeric systems are: TM1,2,H8/TM1,2,H8, TM1,2/

Other interfaces, such as TM4,5/TM4,5, TM4,5/TM5,6, TM5/TM5, and TM5/ TM1,2 did not form in at least one of the five studied u-OR/u-OR, 5-OR/ 5-OR, K-OR/K-OR, 5-OR/u-OR, and 5-OR/K-OR systems. In particular, the TM4,5/TM4,5 interface, which has been suggested to constitute a possible GPCR dimer interface by various experimental assays (reviewed in [22]), only formed during simulation of the 5-OR/5-OR system, and with a lower frequency with respect to other interfaces formed by 5-OR homomers. Notably, neither TM3 nor TM7 were ever found to be involved in a dimer interface, whether formed by the same or different receptor subtypes.

Specifically, the TM1,2 helices appear most frequently at the observed dimer interfaces of 5-OR and K-OR, followed by TM4,5 and TM5,6, with the latter heliX pair being more involved in a dimer interface in u-OR compared to 5-OR and K-OR. However, it is noted that the calculated confidence intervals for frequencies of the specific interfaces are quite broad and overlapping, and therefore the estimated differences between the three OR subtypes may not be as relevant as they appear to be.

Comparison with Putative Dimer Interfaces of GPCRS Inferred from Crystallography

To allow a quantitative comparison, we calculated the minimum Cor root mean square deviation (RMSD) distance between members of the cluster of dimeric complexes that formed during the simulations and each crystal structure listed in 82 Table. The cases where this distance resulted to be less than 10 A are reported in S3 Table for homo-dimers and in S4 Table for hetero-dimers.

In particular, 5-OR and K-OR form TM1,2,H8/TM1,2,H8 interfaces in both homo and hetero-dimers that are very close (RMSDs below 4.3 A) to that seen in the crystal structure of K-OR (4DIH [13]). The closest crystal structure to the TM1,2,H8/TM1,2,H8 interface that forms during u-OR simulations is not the one inferred by the u-OR crystal structure (4DKL [12] ), but rather the one suggested by a [31 -adrenergic receptor (B lAR) crystal structure

Specifically, the interface herein termed TM1,2/TM5,6 deviated only 6.48 A from the asymmetric interface revealed by 30E8 (after overlapping the dimer from simulation with chains A and B of 3OE8) whereas the interface herein called TM4,5/TM5,6 deviated only 6.66 A from an interface inferred from SODU (after overlapping the dimer from simulation with chains A and B of 3ODU) and 6.62 A from an interface inferred from 30E8 (after overlapping the dimer from simulation with chains B and C of 3OE8). Fig. 2B and 2C show structural overlaps of the TM1,2/TM5,6 and TM4,5/TM5,6 interfaces of 5-OR dimers with 30E8. Larger RMSD values were obtained for the identified TM5/TM5 interface in both K-OR and u-OR simulations (see S3 and S4 Tables) after comparison with the available GPCR crystal structures of interacting parallel receptors that are listed in 82 Table. Specifically, in K-OR, this interface is 8.56 A apart from the putative TM5,6-TM5,6 dimer interface inferred from the u-OR crystal structure (after overlapping the dimer from simulation with chain A of 4DKL and its periodic image). Notably, u-OR simulations did not produce dimeric arrangements that were close enough to the crystallographic TM5,6-TM5,6 interface of u-OR, in spite of it being thermodynamically stable as we recently demonstrated through free-energy calculations [14]. The closest structure to the identified TM5/TM5 in u-OR simulations was the interface termed TM4,5/TM4,5 in the BlAR crystal structure corresponding to PDB code 4GPO (RMSD of 8.82 A, after overlapping the formed dimer with chain A of 4GPO and its periodic image). While the calculated RMSD of the identified TM5/TM5 interfaces of K-OR and u-OR homo and het-ero-dimers with respect to available crystal structures appear to be quite large, visual inspection of these overlaps (Figs. 3B and 4B for K-OR and u-OR, respectively) shows that the slight rotation of one protomer needed to match the two configurations would not dramatically change the nature of those interfaces. Thus, we speculate that the reason why the simulations reported herein are unable to reproduce the crystallographic TM5,6-TM5,6 interface of u-OR is that this interface may need longer times to form than its slight modifications (more on this below).

Insights into the Kinetics of OR Dimerization

The results are reported in Tables 3 and 4 for the simulated homo and hetero-dimers of ORs, respectively. The fastest forming homo-dimer interfaces are: TM1,2/TM4,5 and TM1,2/TM5,6 for 5-OR/5-OR, TM1,2,H8/TM1,2,H8 and TM1,2/TM4,5 for K-OR/K-OR, and TM5/TM5 for u-OR/u-OR. Notably, TM1,2,H8/TM1,2,H8, TM4,5/ TM1,2, and TM4,5/TM5,6 are the fastest interfaces for the 5-OR/K-OR hetero-dimer, Whereas the TM5/TM5 interface is the fastest forming for the 5-OR/u-OR hetero-dimer. This observation further stresses the higher propensity for K-OR to have the TM1,2 helices involved in a fast-forming dimer interface. In contrast, TM5 is more likely to be involved in a fast-forming dimer interface when one of the receptor partners is u-OR.

Dynamic Behavior of Lipid Molecules

Such a behavior, which is typical of glass-forming fluids and supercooled liquids [25] , and was also suggested for membrane proteins (e.g., see [26—28] ), has also been recently reported for water at the interface of globular proteins [29]. We speculate that this dynamical heterogeneity also controls membrane diffusion and viscosity near OR dimer interfaces, and plays an important role in modulating the rate of receptor association and the structure of the complex.

As seen in 81—85 Movies, which are provided as representative simulation trajectories of the 5-OR/5-OR, K-OR/ K-OR, u-OR/u-OR, 5-OR/K-OR, and 5-OR/u-OR systems, respectively, the density of slow lipids becomes higher (dark blue in 81—85 Movies) as protomers approach each other, suggesting that slow lipids may interfere with dimer formation at specific interfaces. To quantitatively investigate the role of lipid molecules in regulating OR dimerization, we calculated exchange and persistence time (tX and tp, respectively) distributions of the lipids at different positions relative to isolated receptors. These two quantities characterize the diffusion properties (D) of the lipids and the effective viscosity (1]) of the membrane around the protein, respectively (see details in the Methods section).

In the simulations reported here, the average exchange time of lipids increases up to (tX)%40—50 ns at specific regions near the OR surface, giving Dm2.5><10_7 cmz/s (or, effectively, Dm6.2><10_8 cm2/s). Notably, similar values of lipid diffusion constants have recently been reported in the literature [26,27] for comparable CG force fields, and a similar behavior was implied.

The presence of dynamical heterogeneity, with slow lipids close to the protein surface, corresponds to a breakdown of the Stokes-Einstein relation for lipid dynamics in this region. In other words, correlated lipid motion leads to an increased membrane effective viscosity, decoupled from the diffusional motion of the lipids, so that the inverse relationship between the diffusion coefficient and viscosity (nD o< constant), is no longer homogeneously valid. Analysis of the lipid dynamics around isolated ORs shows longer average persistence times <tp> (see Fig. 5, panels A, B, and C for 5-OR, K-OR, and u-OR, respectively) at specific locations of the protein surface, up to 100 ns. In general, persistence times increase more than exchange times (tX), so that the ratio <tp>/ (tX) is usually larger than 1 (see Fig. 5, panels D, E, and F for 5-OR, K-OR, and u-OR, respectively). Usually, lipid molecules in the region of helices TM1,2 display the shortest persistence times during simulation compared to the region of helices TM4,5 and TM5,6. This is interesting in view of the observed prevalence of interfaces that involve helices TM1,2 and the corresponding fastest on-rates, compared to the generally less frequent participation of TM4,5 and TM5,6 in dimeric interfaces of ORs. Based on this observation, it is tempting to speculate that this local effect on the membrane viscosity (i.e., <tp>/<tx>>1) and the presence of long-lived lipid microstates with long persistence times near the surface may delay the formation of specific interfaces, making them kinetically disfavored with respect to others. The Viscosity of the environment has a well-established, direct effect on the kinetics of biological processes as indicated by the expression of the rate constant in the widely applied Kra-mers’ framework [30]: where m is the effective mass associated with the order parameter used to describe the biological process, G;r is the free-energy of the transition state, to and (o;r are the curvatures of the free-energy profile at the bottom and top of the barrier, and n is the viscosity. According to this eXpression, the higher viscosity observed in regions of the environment with long persistence times results in slower rates, and therefore slower kinetics of the process. To quantitatively assess the relationship between regions of slow lipid dynamics around OR protomers and reduced kinetic rates, we calculated the position-dependent translational diffusion coefficients (DP) of u-OR protomers at the identified homo-dimeric interfaces using a Bayesian inference approach described in the literature [31]. According to these calculations, when the two protomers are within distances of a few nanometers. As reported in SS Table, average diffusion rates of ~5><10'7 cm2/s estimated when the protomers are far apart (d> 50 A) are reduced to ~l—2X10_7 cm2/s as the inter-protomer distance reaches values between 40 and 50 A. At shorter distances (d<40 A), the diffusion rates decreased to less than 10—7 cm2/ s for the slow-forming interfaces, but did not change much for the fast-forming interfaces (i.e, TM1,2/TM5,6 and TM5/TM5). This observation provides a possible explanation why the TM5,6/TM5,6 interface seen in the u-OR crystal structure did not form during the 10 us simulations, notwithstanding its thermodynamical stability [14]. As mentioned before, the two pro-tomers of the TM5/TM5 dimer identified by simulations need to rotate to obtain the crystallographic TM5,6/TM5,6 configuration. We estimated the free-energy associated with such a rotation using the results of steered MD simulations, and assuming no conformational changes occurring within the individual protomers based on the CG model we employed. These results, which are reported in S4 Fig., are consistent with a high free-energy barrier (~10 kcal/mol) between the identified TM5/TM5 and the crystallographic TM5,6/TM5,6 dimers, further supporting the hypothesis that longer time scales are needed for the latter to form. As also evident when viewing the 81—85 Movies, lipid regions of length scales of a few nano-meters between approaching receptors prior to dimer formation can also become locally trapped/restricted in motion (so-called “jammed” regions). While a complete kinetic analysis of the dimerization process that includes lipid dynamics cannot be achieved using the data reported herein, the simulations suggest a distinctive role for long-lived lipid microstructures as they appear to decrease the dimerization on-rate at specific interfaces, and kinetically select specific dimeric arrangements among all different possibilities.

Average values of cholesterol exchange times are reported in 82 Fig., panels AC, for 5-OR, K-OR, and u-OR, respectively, whereas average persistence times and persistence-to-exchange ratios of cholesterol molecules surrounding isolated 5-OR, K-OR, and u-OR are reported in Fig. 6A-C and 6D-F, respectively. Notably, regions with long cholesterol persistence times appear to be generally co-localized with regions with strong lipid molecule persistence, suggesting that both cholesterol-protein interactions and cholesterol-lipid interactions contribute to the kinetic selection of specific dimer interfaces.

For instance, a cholesterol binding pocket was identified in a groove characterized by highly conserved residues (so-called “consensus-motif’ residues) between the intracellular ends of helices TM2 and TM4 in two B2AR crystal structures, i.e., the carazolol-bound 2RH1 [32] and the timolol-bound 3D4S [33]. Cholesterol molecules were also observed in the ultrahigh resolution crystal structure of the A2A adenosine receptor corresponding to PDB code 4EIY [34]. While the “consensus motif’ residues identified in B2AR are conserved in A2A, no cholesterol was observed at the intracellular end of helices TM2 and TM4. In contrast, three cholesterol molecules were found at the extracellular sides of TM2,3, as well as TM5,6 and TM6,7. While no cholesterol molecules were resolved in the K-OR or 5-OR crystal structures, electron density was attributed to a cholesterol molecule in the u-OR crystal structure (4DKL), at the same location between TM6 and TM7 as seen in the A2A crystal structure 4EIY. Notably, the aforementioned “consensus motif ’ residues on TM2 (Y2.41, 82.45) are not conserved in members of the opioid receptor family. Moreover, the calculated high persistence time of cholesterol close to the extracellular ends of helices TM6 and TM7 from simulations of all three OR subtypes (see Fig. 6.) is consistent with the location of cholesterol molecules found in the ultra-high-resolution adenosine A2A crystal structure 4EIY and in the u-OR crystal structure 4DKL. Although the palmitoylation site C3.55 at the intracellular end of TM3 was proposed to constitute part of a cholesterol preferred binding site in the groove lining the intracellular region between TM4 and TM5 [35], this region does not show increased persistence times of cholesterol molecules in the simulations reported here.

Thus, for these specific calculations, we used the results of previous simulations [14] of u-OR crystallographic dimers, i.e., TM1,2,H8/ TM1,2,H8 and TM5,6/TM5,6, where the relative orientation and distance of the protomers in the dimer were maintained constant. As illustrated in Fig. 7, the protein surface adjacent to the TM1,2,H8/TM1,2,H8 and TM5,6/TM5,6 dimer interfaces is in contact with regions of the membrane with slower dynamics. These regions of long persistence lipids right at the dimer interface may provide a mechanistic explanation for the preferential arrangement of receptors into extended linear arrays rather than compact or branched aggregates. In summary, the simulations reported here suggest that both the formation of specific dimer interfaces and the overall topology of oligomeric aggregates depend on the kinetic kon rates. Through calculation of both persistence and exchange times of lipid and cholesterol molecules during simulation, we show the presence of ‘jammed’ lipid regions that exhibit long

Persistence times of lipid molecules around u-OR crystallographic interfaces using data from previously published simulations. Panels A and B show the persistence times of lipid molecules that are adjacent to the simulated TM1 ,2,H8/T M1 ,2,H8 and TM5,6/TM5,6 crystal dimers, respectively. Approximate locations of the centers of mass of the seven transmembrane helices (TM1-TM?) are indicated with solid and outlined red labels for each protomer, respectively. The lipid-exposed region of the protein surface adjacent to the interface is schematically highlighted with a red dashed line. persistence times and non-Poissonian dynamics, and we speculate that these regions play an important role in modulating the kinetics of GPCR di-/oligomerization.

Averaged measures of lipid diffusion have been reported (e.g., see [27]) as a function of the radial distance from the receptor, and have showed a generally slower lipid motion at small radial distances from the protein surface. Although informative, these measures fail to discriminate between the different behavior of lipids adjacent to different helices of the receptor, and are therefore of limited use. While reported measurements of the extent of local time-averaged lipid displacement (proportional to the local average velocity of lipid molecules) [26] have allowed to identify regions of the protein surface likely to be in contact with slower lipids, the average velocity is not a direct measure of the membrane viscosity, thus complicating the interpretation of the results. The present analysis provides rigorous information about the local behavior of the viscosity, thus allowing to analyze its relation to OR dimer formation at specific interfaces.

Our assessment of the rheological properties of the lipid bilayer also complements the analysis of membrane elasticity and mechanical properties Via macroscopic empirical models [36—38]. In particular, the results reported here provide direct evidence for the effect of the microscopic dynamics of lipids and cholesterol at dimeric interfaces on the mesoscopic length and time-scales of receptor interactions, further supporting an essential role of the lipid membrane in determining the identity of homomeric or heteromeric complexes of GPCRs.

Materials and Methods

System Preparation

Missing or unresolved residues (specifically, residues 263—270 in u-OR and residues 262 and 301—307 for K-OR) were built using ROSETTA [39]. The conformation for the intracellular loop 3 (IL3) was selected as the minimum energy structure reported by ROSETTA. For u-OR, the lowest energy conformation that also did not interfere (interlock) with the IL3 of the adjacent receptor at the TM5/ TM6 interface was selected. For 5-OR, the crystallographic IL3 was rebuilt using ROSETTA between residues 241 and 258. Notably, the root mean square deviation (RMSD) of this loop conformation from the resolved loop of the newest high-resolution crystal structure of 5-OR [40] is 0.46 A overall. Throughout this article we use the Ballesteros-Weinstein numbering scheme to facilitate comparison between the receptor subtype TM regions [41]. Accordingly, the first number of this scheme indicates the TM heliX in which the residue in question resides, and the second number indicates the position of the residue relative to the most conserved residue in the heliX, which is always numbered 50. The receptors were converted to a CG representation under the MARTINI force field (version 2.1) [16—18] and a modified elastic network was applied, as reported previously in the literature [42]. According to recent experimental findings [35], the receptors were also palmitoylated at the C355 position, with a palmitoylate chain consisting of 4 C1 beads, with a bond length of 0.47 nm, a force constant of 1250 k] mol-1 nm-2 and angles of 180°, with a force constant 25 k] mol_1.

Molecular Dynamics Simulations

Temperature coupling at 300 K was achieved with the V-rescale algorithm, and pressure coupling at 1 bar was achieved with the Berendsen algorithm. Simulation were performed with GROMACS version 4.5.3 [43]. The systems were minimized and equilibrated with harmonic restraints of decreasing strength over 10 ns, on protein backbone beads. Water and counter ions were removed and sixteen of these membrane/protomer systems were randomly selected (repeat selections were permitted) and combined to create five 16-receptor setups for each of the studied homomeric (u-OR/u-OR, 5-OR/S-OR, and K-OR/K-OR) and heteromeric (5-OR/K-OR and u-OR/S-OR) systems, which were re-solvated, neutralized, and subsequently minimized. The heteromeric systems contained eight receptors of each subtype. Thus, twenty-five independently constructed 16-recep-tor systems, each with approximately 57,000 beads, were generated. Each of these systems was simulated for 10 us with a timestep of 20 fs, to give 50 us of pooled simulation time for each homo and heteromeric receptor system, and 250 microseconds in total, across all systems. Periodic boundary conditions were employed, and neighbor lists were updated every 10 steps. The Shift algorithm was used for electrostatic interactions. A single cutoff of 1.2 nm was used for Van der Waals interactions.

Interface Identification and Clustering

Dimeric interfaces were clustered by a k-means algorithm using the distance induced by the Frobenius norm on the contact map matrices, and the clusters automatically labeled according to the TM regions involved in receptor-receptor interactions. Specifically, the dissimilarity of two interfaces k1 and k2 was defined as for hetero-dimers, While for homo-dimers the symmetrized form was used to account for the equivalence of the dimer With swapped protomers. The marginal contact map averages over the NC interfaces k in cluster C were calculated for one receptor as while for the interacting one the analogous expression with i and j swapped inside the sum was used. Helices for which at least 3 residues were involved in the interface, were included in the label of the cluster. We used a Bayesian inference framework to pool the information from the different trajectories and calculate estimates of the interface prevalence for each dimer and on-rates. The number of dimerization instances in a given trajectory i that yielded an interface belonging to cluster C can be described by a stochastic variable NC with a multinomial distribution where the probability of each cluster is pLC = wC/chc, wC being the overall weight of the cluster. We defined WC 2 exp(aC), used noninformative normal priors (with zero mean and large standard deviation) for the parameters ac, and employed Gibbs sampling to obtain the posterior distributions of pC, for which we report the average as well as (2.5%,97.5%) confidence intervals. While best estimates of the probabilities pC are well approximated by the pooled fractions, the larger confidence intervals obtained by Bayesian inference reflect the variations in the different trajectories. Estimates of the on-rate for dimerization k0n(C) at each interface C were obtained assuming that associations are independent of each other, and that the number of interface-specific di-merization events 11 follows a Poisson distribution: Where t is the time and c1- the concentration of monomers in trajectory 1'. Again, we defined the parameters 190- so that the Poisson intensity is (t kon ci) 2 exp(bCi), where 190- have non-infor-mative normal priors, and sampled the posterior distribution of kon using Gibbs sampling. For all estimates, 104 samples from the posterior distributions were obtained after a 5><103 burn-in phase using Markov-chain Monte Carlo techniques [44].

In order to ignore the structural differences in the monomeric structures, and capture only the degree of similarity of the OR dimer interfaces from simulation with those inferred by crystal structures, we aligned the individual CG ORs to the receptors in each crystal dimer. For the evaluation of the identified hetero-dimeric interfaces, both pair combinations (e.g., u-OR/S-OR and 5-OR/u-OR) were considered for the superposition onto crystal structures, and the RMSD was defined as the minimum between the two individual RMSD values.

Lipid Velocity Distribution and Residence Times

We characterized the dynamical properties of the lipid (POPC) and cholesterol molecules by calculating their exchange and persistence time distributions (tX and tp, respectively) at different positions relative to the protein. Specifically, we applied the equations reported in [25] , and calculated local averages of lipid exchange and persistence time distributions, (tX), and

5—6). Formally, this approach consists of viewing the lipid dynamics as a continuous time random walk (CTRW) model [45], defined by a coarse-graining length d, i.e. as a sequence of displacements of length d occurring at different times. The reported results for d = 10 A do not change significantly for other choices of d. Increased persistence times tp reflect long-lived microscopic arrangements of lipid molecules, thus representing a measure of structural correlation of the lipid motion or the membrane local viscosity noc(tp). The exchange time is the time elapsing between subsequent displacements of length d. Thus, the first exchange time tXJ is the minimum time for which ||x(tX,1)-x(tp) | | >d, where x is the position of the lipid, and the following values of the exchange time ti are defined similarly by the condition | |x(tx,i) -x(tx,(i_1)) I I >d.

The latter characterizes the dynamics of the lipid molecules by a diffusion constant Do< 1/ (tX). The distributions of persistence and exchange times, pp and pX, respectively, are related by: so that in normal bulk conditions, when pX is Poissonian, so is pp, and (tX)~(tP).

The effective viscosity 1] increases and is de-coupled from the diffusional motion of the lipids. We calculated the distributions pp and pX in bins of IAX 1A using the dynamics of the phosphate group in the first third of each trajectory, and averaging over upper and lower leaflet lipids, over all homomeric trajectories for a given receptor subtype, and over all protomers in the simulation, after alignment to a common reference protein molecule. Although the heteromeric simulations yielded comparable results, these are not reported here because of their reduced statistical significance.

Position-Dependent Diffusion Coefficients

Continuous trajectories of pairs of interacting proteins were extracted from the original trajectories containing 16 copies of the receptor. To avoid problems with the interpretation of the results, frames in which a given interacting pair of protomers was in contact with a third protomer were discarded. This was achieved by checking that the distance between protomers in a given pair was smaller than distances from any other nearby protomer (see S3 Fig). The time dependence of the two angles defining the relative orientation of the protomers (4)1, (p2) and their distance (d) was calculated for all the frames of the resulting trajectories of protomer pairs. Trajectories relative to the formation of specific interfaces were selected by restricting the analysis to regions I = {ocl $431361 0623432362}, and the distances {3 nm 3 d g 9 nm} binned to obtain a set of unbiased binned trajectories X = {dj(r)}. The diffusion coefficient D was calculated by sampling the posterior distribution of the rate matriX from the posterior distribution p(K|X) = p(X|K) p(K), assuming a uniform prior p(K), and the likelihood Where ND» is the number of observed transitions between distance bin j and bin j’. From the distribution of the rate matrix K, the distribution of the diffusion coefficients was obtained from Average values of these quantities and their (25%, 95%) credible intervals are reported in SS Table.

Free-Energy Barrier Estimation from Steered MD

The free-energy barrier separating the u-OR TM5/TM5 interface from simulations and the crystallographic TM5,6/TM5,6 interface was estimated using the Iarzynski equality [46] and multiple runs in which the angles (p1 2 (p2 were steered from 1.2 degrees (corresponding to the TM5/TM5 interface) to 2.1 degrees (corresponding to the TM5,6/TM5,6 interface), with a velocity of 0.05 rad/ns and an elastic force of k = 800 (kcal mol_1)/rad2. The Iarzynski expression was applied separately to 5 sets of 5 independent runs to obtain the potential of mean force (PMF) and error estimates.

Supporting Information

Approximate locations of the center of mass of the seven transmem-brane (TM1-TM7) helices are indicated with red labels.

Approximate locations of the center of mass of the seven

Only frames in which no nearby protomers occupy the gray shaded region are included in the analysis. The relative position of the protomers is defined by the COM-COM distance d, and the angles of the COM distance with respect to the direction of the

PMF of the transition from the observed TM5/TM5 interface (4)1 = ¢2~1.2) and the crystallographic TM5,6/TM5,6 interface (4)1 = ¢2~2.4) of p-OR homo-dimers. Solid line is the average PMF obtained from the application of the Iarzynski equation to 5 independent sets of 5 simulations, and errors are the corresponding standard deviations. (TIF)

Summary of all the molecular dynamics simulations of coarse-grained representations of opioid receptors reported in this manuscript. (DOCX)

Putative dimer interfaces inferred from recent GPCR crystal structures. (DOCX)

Minimum RMSD distances (below 10 A) between homo-dimeric complexes formed during simulation and selected crystal structures. For each cluster of dimeric complexes that formed during simulation, the configuration corresponding to the lowest RMSD (highlighted in bold) is depicted in Figs. 2, 3, and 4 of the manuscript for 5-OR/5-OR, K-OR/ K-OR, and u-OR/u-OR, respectively.

Minimum RMSD distances (below 10 A) between hetero-dimeric complexes formed during simulation and selected crystal structures. (DOCX)

Position-dependent diffusion rates calculated for the identified homo-dimeric configurations of u-OR from simulations. Protein translational diffusion coefficients (in 10—7 cmZ/s) for different interfaces and different distance ranges.

Example of a simulation run for the homomeric 6-OR/6-OR system. The instantaneous velocity of lipid molecules was calculated after interpolating the motion of the corresponding head groups over 20 ns. The density of lipids with velocity within the first quartile was calculated using a Gaussian kernel (with bandwidth 5 A) and plotted in colors (blue, high density; white low density). The location of protein molecules is indicated approximately by plotting one single isodensity contour of the protein beads obtained with the same Gaussian kernel. The location of the center of mass of helices TMl, 4 and 6 is indicated with white, gray, 82 Movie. Example of a simulation run for the homomeric K-OR/K-OR system. The instantaneous velocity of lipid molecules was calculated after interpolating the motion of the corresponding head groups over 20 ns. The density of lipids with velocity within the first quartile was calculated using a Gaussian kernel (with bandwidth 5 A) and plotted in colors (blue, high density; White low density). The location of protein molecules is indicated approximately by plotting one single isodensity contour of the protein beads obtained with the same Gaussian kernel. The location of the center of mass of helices TM1, 4 and 6 is indicated with White, gray, and black dots, respectively.

Example of a simulation run for the homomeric u-OR/u-OR system. The instantaneous velocity of lipid molecules was calculated after interpolating the motion of the corresponding head groups over 20 ns. The density of lipids with velocity within the first quartile was calculated using a Gaussian kernel (with bandwidth 5 A) and plotted in colors (blue, high density; white low density). The location of protein molecules is indicated approximately by plotting one single isodensity contour of the protein beads obtained with the same Gaussian kernel. The location of the center of mass of helices TMl, 4 and 6 is indicated with white, gray,

Example of a simulation run for the heteromeric 6-OR/K-OR system. The instantaneous velocity of lipid molecules was calculated after interpolating the motion of the corresponding head groups over 20 ns. The density of lipids with velocity within the first quartile was calculated using a Gaussian kernel (with bandwidth 5 A) and plotted in colors (blue, high density; white low density). The location of protein molecules (red for 5-OR and black for K-OR) is indicated approximately by plotting one single isodensity contour of the protein beads obtained with the same Gaussian kernel. The location of the center of mass of helices TM1, 4 and 6 is indicated with white, gray, and black dots, respectively. (MP4)

Example of a simulation run for the heteromeric 6-OR/u-OR system. The instantaneous velocity of lipid molecules was calculated after interpolating the motion of the corresponding head groups over 20 ns. The density of lipids with velocity within the first quartile was calculated using a Gaussian kernel (with bandwidth 5 A) and plotted in colors (blue, high density; white low density). The location of protein molecules (red for 5-OR and black for [.1-OR) is indicated approximately by plotting one single isodensity contour of the protein beads obtained with the same Gaussian kernel. The location of the center of mass of helices TM1, 4 and 6 is indicated with white, gray, and black dots, respectively.

Acknowledgments

Computations relied on the Extreme Science and Engineering Discovery Environment (XSEDE) under MCB080109N, and on the computational resources and staff expertise provided by the Scientific Computing Facility at the Icahn School of Medicine at Mount Sinai.

Author Contributions

Performed the experiments: DP MBB IM]. Analyzed the data: DP IM] MF. Contributed reagents/materials/analysis tools: DP. Wrote the paper: DP IM] MF.