Use of network model to explore dynamic and allosteric properties of three GPCR homodimers

Abstract

Recently, increasing experimental evidence has indicated that G-Protein Coupled Receptors (GPCRs) can form dimers, which are very possibly further potential functional units and new targets for drug development besides their monomeric units. However, knowledge about their structure and functional motion has been limited so far. Thus, we used an Elastic Network Model (ENM) and Protein Structure Network (PSN) to study three A GPCR homodimers (viz., CXCR4, κ-OR, β1AR) with two different interfaces based on their basic topologic structures. The low-frequency modes from ENM exhibit similarity to some extent, indicating similar functional motions shared by A GPCR dimers, such as asymmetric motion in the ECL2 and TM6 regions around the interface, which should contribute to the negative cooperation for ligand binding and asymmetric activation reported experimentally. The PSN results reveal that the dimerization can reduce the main informational flows from the extracellular to the intracellular domain and affect the contribution of TM regions to the allosteric paths. Some highly conserved residues were still observed to be hot residues in the meta-pathway, further confirming their conserved importance shared by A GPCR dimers; in particular for F6.44 and F6.48 residues and one non-conserved position X7.39. On the whole, dimerization plays a different role in influencing the dynamic motion of the protomer, dependent on the type of interface and contact area. Compared to the TM5–TM6 interface, TM1–TM2–H8 exhibits more a significant functional-role in influencing the dynamic behavior and allosteric paths.

---

Introduction

G-Protein coupled receptors (GPCRs) are the largest superfamily of membrane proteins, which share a similar architecture of seven alpha-helical transmembrane (TM) segments held together by tertiary contact. They are responsible for triggering many diverse cell responses and regulating their functions in the body, thus representing major targets for drug development.^1,2 GPCRs were originally regarded as a monomeric functional unit to couple to G-proteins and activate their signalling cascade.^3–5 But, recently, greatly increasing biochemical and biophysical evidence has indicated that GPCRs can form dimers or oligomers in living cells,^6–9 which may provide another main functional unit in signaling transduction and serve as novel targets for drugs. Thus, studies on GPCR dimers or oligomers have aroused considerable interest to date and have become another important subject in the GPCR field.

As revealed, GPCRs are allosteric proteins. The traditional mechanism of monomeric GPCRs is that the receptor is activated by ligands and undergoes a series of conformational changes to recruit the G-protein at a topographically distinct site and trigger signalling.¹⁰ In the context of GPCR dimers, similar allosteric modulations between the two protomers were observed in a few experiments.¹¹ The binding of ligands or proteins to one protomer could influence the binding to the other protomer through allosteric communications between the two protomers, leading to negative or positive cooperativity, which could modulate a GPCR’s functions such as ligand pharmacology and initiate different signalling pathways.^9,12 However, due to the few crystal structures of oligomers available and much fewer studies on the oligomers with respect to the monomeric GPCRs, there has been a lack of knowledge about the structure, dynamics and functional motions of GPCR oligomers, leading to many unanswered questions during experiments.¹³ For example, what are the functional interfaces? What is the functional mechanism of GPCR oligomers? Thus, in the case of the very limited oligomer structures available, it is highly desirable to introduce advanced computational methods to assist the experiments probing the above questions.

As known, molecular dynamics (MD) simulation is a powerful tool to provide microscopic information about structures and dynamics at the molecular level.¹⁴ Some all-atom molecular dynamics simulations were performed to study certain GPCR dimers in terms of their structures, dynamics and interaction with ligands within the nanosecond timescale.^15–19 However, the timescale of the functional motions of GPCRs are generally in the millisecond range, which is very difficult for all-atom MD due to the limited computing ability. Therefore, coarse-grained MDs were utilized to study GPCR oligomers at longer time scales, such as microseconds, which mainly focused on the self-assembly behaviors of the oligomers.^20–23 These MD studies provide valuable information for understanding the structures and dynamics of GPCR oligomers. The information, however, is still very limited, particularly in terms of the absence of dynamics information associated with important functional regions; for example, allosterism and the cooperative effect between the protomers. Therefore, it is still highly necessary to introduce some other effective computational techniques to obtain further more information to add to that from previous studies so that we can better understand the functional mechanism of GPCR dimers.

The elastic network model (ENM) is a harmonic vibration model based on structural contact topology and the statistical mechanical theory of polymer networks.²⁴ Although the method lacks atomic detail, it can provide insights into cooperative motion. As revealed, low-frequency modes in the cooperative motions predicted by ENM could well describe experimentally observed conformational change and are in general associated with functional motions.²⁵ Thus, ENMs have been broadly applied to unravel protein dynamics.^26–28

The protein structure network (PSN) method is a graph-based approach applied to protein structures.²⁹ The method could provide network features like nodes, hubs, and links, and gain insight into the holistic properties of protein dynamics, topological rearrangements and intra-molecular and inter-molecular communications, which are pivotal for proteins to execute their biological functions.^30,31 Consequently, it has been widely applied to study folding/unfolding, stability, conformational changes and allosterism in proteins.^32–35

In this study, we utilized ENM and PSN methods to study the three class A GPCR homodimers since the homodimers were considered to be a predominant species with respect to heterodimers.⁹ The available crystal structures of class A GPCR dimers are very limited, comprising only rhodopsin,³⁶ CXCR4,³⁷ the κ-opioid receptor³⁸ and β1AR,³⁹ in which there are two kinds of interfaces, viz., TM1–TM2–H8 for the dimers of rhodopsin, κ-OR and β1AR, TM5–TM6 for CXCR4 dimers. Herein, we select two homodimers with a TM1–TM2–H8 interface (viz., κ-OR and β1AR) and one homodimer comprised of CXCR4 with a TM5–TM6 interface to be representatives of the class A GPCR dimers to address some problems associated with functional mechanisms of the dimers; for example, how the dimerization influences the dynamic behaviors and the allosteric mechanism for the three GPCRs? How about correlation between the extent of impact and the two types of interface? The ENM method was utilized to predict large-amplitude dynamics while PSN was applied to explore the impact of dimerization on the allosteric structure communication. Common and different behaviors between the three GPCR homodimers were observed, dependent on the type of interface and contact area. Some observations provide further support for the experimental findings. More importantly, some novel observations are obtained, which provide valuable information for extending our understanding of the structure and functional mechanism of GPCR dimers.

Results and discussion

Collective motions characterized by the lowest frequency modes for the dimers

The Gaussian Network Model (GNM) is a kind of Elastic Network Model (ENM), which can provide information about N-dimensional properties such as mean-square fluctuations (MSFs) of residues, their cross-correlations, or movements along normal mode axes. Some studies already demonstrated that the collective motions of proteins derived from the lowest frequency modes are associated with functional conformation changes.^40,41 Thus, to gain insight into the functional significance of residue fluctuations, we used GNM to calculate the lowest-frequency modes for the three dimers.

Fig. 1 shows the first three non-trivial lowest-frequency GNM modes for the three dimers, in which the residue mobility given by eigenvectors of the most cooperative motion-mode are presented as curves of residue index i. In addition, Fig. 1 also displays ribbon diagrams of the three dimers color-coded according to the relative motions in the order of decreasing mobility: red (highest mobility), cyan, green, yellow, orange to blue (lowest mobility). It can be seen from Fig. 1 that the three dimers are to large extent similar for the variation trends of the first three GNM modes, implying that they should share to some extent similarity for some fundamental motions.

Fig. 1 The first three lowest-frequency mode profiles for every protomer in the three dimers. Modes 1, 2, and 3 are shown in blue, red, and green, respectively. Black dots denote hinge sites. The ribbon diagrams of the three receptor dimers are colored by blue, cyan, green, yellow, orange and red, in accordance with their relative motions in the global mode profile from the lowest (blue) to the highest (red).

For the first mode (labelled as 1st), it is clear that the two protomers an exhibit antisymmetric trend with respect to the x-axis for the three dimers, indicating one anti-correlated motion of the two subunits surrounding the interfacial region. As clearly reflected by the blue region in the ribbon diagram in Fig. 1, the lowest mobility residues of 1st mode are located in the interfaces for the three dimers despite the different interfaces, probably because the motions of residues in the interfaces are restricted by interaction between the two subunits.

Different from the 1st mode, for the second (labelled as 2nd) mode, the two protomers present symmetric trends about the x-axis for the three dimers, suggesting identical motions for the two subunits. In addition, the profiles of the 2nd mode show that the eigenvector values fluctuate around the x-axis in the positive and negative directions along the y axis, implying correlated and anticorrelated movements between pairs of structural elements. Some crossover residues between the anticorrelated regions are found in the intersection site with the x axis (vide black dots in Fig. 1). The crossover residues exhibit zero mobility and were considered as hinge sites that are usually implicated in mechanical roles relevant to biological function.⁴³ The right ribbon diagrams in Fig. 1 further reveal that the hinge residues are located in the middle of the TM helices between the extracellular and cytoplasmic regions, thus dividing the protein into two halves and leading to their anti-correlated motions. Table 1 further lists the hinge residues for the three dimers. As shown in Table 1 and Fig. 1, these hinge residues distribute over each helix. Furthermore, most hinge residues locate or are close to the most conserved positions in each helix except for the residues in TM3; for example, F6.44 positions in CXCR4, κ-OR and β1AR were revealed to be conserved residues in TM6 and hubs in the signal transduction.⁴⁴ The hinge residues in TM2 are very close to the conserved position 2.50, which is involved in the stability of the transmembrane helix bundle.^45–47

Table 1 Global hinge sites^a (vide minima in global mode shapes of GNM in Fig. 1)

Receptor	CXCR4	β1AR	κ-OR
a The integers are the residue index, and the superscripts are the corresponding Ballesteros–Weinstein numbering.^42
TM1	51^1.45 to 54^1.48	54^1.45 to 57^1.48	72^1.45 to 78^1.51
TM2	85^2.51 to 89^2.55	86^2.49 to 89^2.52	102^2.47 to 106^2.50
TM3	114^3.30 to 117^3.33	124^3.35 to 127^3.38	143^3.37 to 146^3.40
TM4	167^4.56 to 169^4.58	166^4.50 to 170^4.54	181^4.48 to 185^4.52
TM5	202^5.41 to 204^5.43	214^5.45 to 217^5.48	233^5.45 to 236^5.48
TM6	247^6.43 to 249^6.45	299^6.44 to 304^6.49	283^6.44 to 286^6.47
TM7	294^7.45 to 296^7.47	334^7.45 to 338^7.48	323^7.44 to 326^7.48

---

Similar to the 1st mode but different from the 2nd mode, the third mode (labelled as 3rd) exhibits an antisymmetric trend with respect to the x-axis for the two protomers, also possessing an anti-correlation feather. The eigenvector variation trend of the 3rd mode is similar to that of the 2nd one, displaying antisymmetric movement about the x-axis and the same hinge residues. As shown in the ribbon diagram in Fig. 1, some of the lowest mobility residues are located in the interfacial region, similar to the 1st mode, while the other lowest mobility residues are distributed in the hinge site, similar to the 2nd mode. Thus, it can be assumed that the 3rd mode contains the features of the 1st and 2nd modes, in which the feature of the 2nd mode are more notable than the 1st one.

It was proposed by some experimental studies that the dimeric or oligomeric GPCRs are subject to some specific cooperative processes that regulate their biological function,^11,48,49 which should be associated with their asymmetric nature. In order to gain insight into the nature of low-frequency modes associated with the functional regions, we focused on the collective motions of some functional regions associated with the ligand binding and the activation of GPCR, as shown in Fig. 2. It can be seen for the CXCR4 dimer that the motion of residues associated with mode 8 presents a distinct difference between the two subunits, in particular for ECL2. ECL2 of one protomer exhibits conspicuous mobility and the other has very low mobility, showing anti-correlation between the two subunits. As revealed, ECL2 plays an important role in ligand binding, acting as a binding-gate.^50,51 The protomer with high flexibility should facilitate the entrance and binding of ligands while the other protomer with very low mobility should be disfavored.

Fig. 2 Some asymmetric modes associated with functional motions for the two subunits of the three dimers, derived from the first 20 lowest frequency modes. Ribbon diagram colored by blue, white and red, which denote the increasing mobility of the relative motions in the asymmetric mode profile from the lowest (blue) to the highest (red). The regions with marked different motions between the two protomers are highlighted by circles.

Some experimental studies⁵² also found that one protomer of β1AR dimer with a high ligand affinity conformation could promote the second protomer to have a low ligand affinity conformation, consistent with our observations. The asymmetric motions were considered to be associated with physiological significance. For example, the asymmetric motion involved in the ligand binding site would prevent the cells from excessive stimulation by the ligand through anticorrelation motion.⁹ In addition, Fig. 2 also clearly shows that TM6 of one protomer presents significant fluctuations, but the other protomer does not. As accepted, the most obvious feature of the activation state of class A GPCR is that TM6 in the intracellular region tilts significantly outward so that the receptor opens the G protein binding pocket.⁴⁴ The experimental study on leukotriene B4 receptor dimer also found that the important TM6 displays an asymmetric conformational change.⁵³ Similar to CXCR4, the two protomers of the β1AR and κ-OR dimers both show asymmetric feathers for the ligand binding domain and the G protein-binding domain, as shown in Fig. 2. The observations from the work provide support for the experimental assumption that GPCR dimers adjust their physiological signals in the cooperative processes through their asymmetric nature.^9,11,54 In fact, the cooperativity driven by asymmetric motions via chemical or physical effects like ligand binding and mechanics was considered to be a generalized manner of function for oligomers of many soluble and membrane proteins.^55,56 However, it should be noted that the asymmetric motions observed in the low-frequency modes do not preclude the existence of a symmetric motion in higher modes. In addition, as is known, the real motion is a linear combination of the individual modes and the motion of one mode may be partially cancelled by another.^24,40

Impact of different interfaces on the subunit dynamics

Structural models of GPCR dimers vary widely. Some crosslinking and mutagenesis studies have pointed to a number of potential interfaces.⁶ In these works, the three dimers present two kinds of interface, viz., TM1–TM2–H8 for β1AR and κ-OR, TM5–TM6 for CXCR4. As is known, ANM (Anisotropic Network Model) is another kind of ENM, which can give 3-D descriptions of 3N − 6 internal modes and in turn the directions of collective motions, so that it can permit us to evaluate directional preferences. Thus, we, herein, calculated ANM in order to gain insight into the role of dimerization in influencing the dynamics behavior of the protomer.

We compared the ANM modes between the isolated monomer with the protomer, using the subsystem/environment coupling method described in the Methods section, in which one protomer is the subsystem and the other serves as the environment.^57,58 The overlaps between the six lowest-frequency ANM modes in the isolated monomer (vide y-axis) and those in the protomer (vide x-axis) are displayed in Fig. 3, in which the orange–red entries along the diagonal in each panel denote that the modes in the isolated GPCR are closely maintained in the dimeric receptor.

Fig. 3 Overlaps of the first six slowest ANM modes between the isolated monomer and the protomer in dimers (labelled as protomer^dimer) are shown in heat maps (top) for CXCR4 (A), β1AR (B) and κ-OR (C). Dynamic motions from several highly-overlapped modes are shown at the bottom. a, b and c denote the first three lowest-frequency modes of CXCR4, respectively, while d denotes the first lowest-frequency mode of β1AR. White and red arrows represent the directions of the collective motions of residues and the receptor, respectively.

As reflected by Fig. 3, there are high overlaps for the first four lowest-frequency ANM modes for CXCR4, indicating that dimerization does not play a significant role in the basic function of the protomer. Fig. 3 also shows the dynamic motions associated with the first three lowest-frequency modes. It is clear that the first three modes present a common antisymmetric rotation between the extracellular and intracellular regions, with minor differences in the extent and the direction. The anti-correlated motions were reported and confirmed by some combination studies of experiments and calculations on rhodopsin⁵⁹ and ghrelin receptors,⁶⁰ which were considered to play an important role in GPCR functions.

For β1AR, the extent of the overlap is significantly weakened and only the first mode is very similar between the isolated monomer and the protomer. The other low-frequency modes exhibit low overlap, displaying a more significant impact of the dimerization on the protomer dynamic motion for β1AR than CXCR4. Similar to the CXCR4, the first mode also presents to some extent anti-correlated motions between the two sides of the helix bundles for β1AR, as reflected by Fig. 3. The most significant difference is observed for κ-OR, in which case there is no overlap for the six low-frequency modes. κ-OR and β1AR have a similar interface (TM1–TM2–H8) but they present differences in the impact of the dimerization on the monomer dynamic motion. We calculated the inter subunit surface areas of CXCR4, β1AR and κ-OR, which are 754 Ų, 768 Ų and 944 Ų, respectively. The contact area of κ-OR dimer is the largest, which should be responsible for a larger perturbation of the intrinsic dynamics of the protomer upon dimerization than CXCR4 and β1AR, because the large interface area ensures closer association of the protomers. The contact areas of CXCR4 and β1AR are very close, but significant large differences between the monomer and the protomer are observed for β1AR with the TM1–TM2–H8 interface than CXCR4 with the interface of TM5–TM6. In fact, the TM5–TM6 interface of CXCR4 is only involved in the extracellular side, while the TM1–TM2–H8 interface of β1AR is not only involved in the extracellular but also the intercellular region of H8. The high overlap for the first four low-ANM modes in the CXCR4 indicates that TM5–TM6 dimerization interface plays a weaker role in the basic function of the monomer than the TM1–TM2–H8 interface, which should support the recent experimental assumption that the monomeric CXCR4 acts as a functional unit.³ The significant impact of the dimerization derived from the TM1–TM2–H8 interface suggests that the interface is more likely to be related to some specific function if the dimer can play a functional role. The latest quaternary structure study of the D₃ receptor showed that two dimers with the TM1–TM2–H8 interface could interact to form a ‘rhombic’ tetramer.⁶¹ In addition, a computational study of the self-assembly behavior of β1AR and β2AR also indicated that the interface involved in TM1/H8 is more stable than the other one, like TM4/3.⁶² Our work provides further support for this experimental assumption regarding the functional TM1–TM2–H8 interface.

Impact of the dimerization on the allosteric pathway from the extracellular domain to the intracellular one

GPCRs are allosteric proteins that transform extracellular signals into intracellular G proteins to promote nucleotide exchange. As is known, the allosteric effect mainly stems from the transformation of non-covalent interactions between residues. PSN can abstract the topological structure of a protein to a network, which is built on the non-covalent interactions between residues.³⁰ Since ENM analysis could obtain the correlation motions between the residues, we used the PSN–ENM method to calculate the structural allosteric pathways from the extracellular to the G-protein binding regions in order to explore the impact of the dimerization on their allosteric communication paths. The changes in the structure network derived from the PSN–ENM method could capture the reorganization of the structure communication between the dimer vs. the monomer, which may have implications in understanding the functional roles of the oligomeric states towards signaling.⁶³

Herein, the ligand-binding residues in the crystal structure were selected as the extracellular residues since the GPCRs are in general activated by ligands. In addition, residues in the G protein coupling, NPxxY and DRY motifs, which have the same B–W numbering, were chosen as the intracellular residues since residues with the same B–W numbering are in general considered to be involved in the analogous functions.⁶⁴ These selected residues are shown in Table S1 in the ESI.† We calculated the recurrence of nodes and links in the pathway pool and then constructed meta-pathways constituted by the recurrent nodes and links by means of statistical analysis, as shown in Fig. 4. Herein, we mainly focused on the main meta-pathways consisting of frequent nodes and links, with a recurrence of ≥0.2.^65,66 Table 2 lists the numbers of these main meta-pathways in the isolated monomer and the protomer in the corresponding dimers.

Fig. 4 Network representation of the meta-paths consisted of nodes and links with an rec value above 0.2 for the dimeric (A) and monomeric (B) types of the three GPCRs. Decimals near the residues are the Ballesteros–Weinstein numbering of these residues. The node sizes and the widths of the links are proportional to the recurrence (rec) values.

Table 2 The number of meta-paths with rec values above 0.2 and the number of nodes and links contained in these paths, as well as the similarity of the links between the isolated monomer and the protomer in the dimer, derived by PSN–ENM analysis

	System	Nodes	Links	Paths	Similarity^a
a Similarity = the number of common links shared by the isolated monomer and the protomer/the number of total links in the protomer.
CXCR4	Monomer	24	24	18	62.5%
	Protomer	21	21	3
β1AR	Monomer	28	29	49	51.7%
	Protomer	25	24	20
κ-OR	Monomer	21	18	20	44.4%
	Protomer	21	20	15

---

In the following discussions, we are not only concerned with the specific residues but also their positions in the meta-paths, in order to provide instruction for the other class A GPCR dimers. Venkatakrishnan analyzed seventeen class A GPCR crystal structures and identified a conserved network involved in 36 topologically equivalent positions, which are maintained across these class A GPCRs.⁶⁴ Although differences in some specific residues located in these topological positions exist between different GPCRs, mutation studies indicated that most of the 36 positions identified could represent structurally important positions in the GPCR receptors. These observations clearly showed the conserved nature of these important positions for class A GPCRs. Thus, we can surmise that key positions identified by PSN from the three GPCR dimers should be important for the other A GPCR dimers if they are also located in the 36 conserved positions. Fig. S1 in the ESI† shows the sequence alignments of the three GPCRs and reports Ballesteros–Weinstein numbering, TM annotation, identity, conservation and the positions in the conserved interaction network.⁶⁴ The sequence alignments are processed by GPCRdb⁶⁷ (http://www.gpcrdb.org/). Fig. 4 displays the meta-paths consisting of nodes and links with an rec value above 0.2 for the dimeric and monomeric types of the three GPCRs. In addition, Fig. 5 further shows the meta-paths embedded in the structural architectures of the three GPCRs, in which the specific paths to the dimer/monomer and the common paths shared by the monomer and the protomer are highlighted by different colors.

Fig. 5 Pathways consisting of links with rec > 0.2 for the monomers and the dimers. The nodes and links colored orange and purple are specific to the monomer and the protomer, respectively, while those shared by the monomer and dimer are colored green. The size of nodes and width of links are proportional to their recurrence.

For CXCR4, the PSN result shows that the difference is not significant between the monomer and the protomer. The protomer in the dimer (labelled as protomer^dimer) shares 62.5 percent of the links of the isolated monomer (vide Table 2). There are 18 pathways composed of links with rec > 0.2 in the monomeric CXCR4. These pathways are mainly from the ligand binding site involved in TM2 and TM7 and the cytoplasmic side of TM3, TM5 and TM6 associated with the G-protein binding site. For the protomer in the dimer, there are only three pathways, indicating that the dimerization weakens the main structural communications. The decrease mainly stems from the weakness in the links between N119^3.35, H294^7.45 and C295^7.46, as reflected by Fig. 4A. The three paths of the protomer in the dimer are same as the three pathways in the monomer, which start at E288^7.39 and end at T241^6.37, I222^5.61 and R134^3.50 located in the cytoplasmic side of TM6, TM5 and TM3, respectively. The nodes and links in the three pathways have a high rec, as revealed by Fig. 6, indicating that they are essential channels for the allosteric pathways from the ligand-binding site to the G-protein binding one. In addition, the three pathways shared by the monomer and the protomer contain the same segment, comprised of Y255^6.51–W252^6.48–F248^6.44–L127^3.43–I245^6.41–Y219^5.58. This segment locates in the middle of the helix bundle and consists of many highly conserved residues like W252^6.48, F248^6.44, L127^3.43, I245^6.41 and Y219^5.58, which were reported to play an important role in GPCR activation.⁴⁴ For example, the importance of highly conserved F6.44 and W6.48 were confirmed by many experiments.^44,64 The latest mutation experiment of CXCR1, which is highly homologous with CXCR4, indicated that F6.44 is spatially located in a “hot spot” and is essential for CXCR1 activation.⁶⁸ The W6.48 was also proposed to be the water-gating micro-switch residue⁶⁹ and participate in controlling the relative position between TM3 and TM6 in the activation process through forming inter helical interactions.⁷⁰ Compared to the inactive state of rhodopsin, Y5.58 in the active state was observed to have a significant conformational change that repositions it close to the ionic lock region in direct hydrophilic interaction with R3.50.⁷¹ In addition, it was reported that Y5.58 could form a hydrogen bonding network with the Y7.53 of the NPxxY motif through a water molecule in the activation process of rhodopsin. Residues L3.43, X6.40/X6.41, and F6.44 were revealed to participate in the “hydrophobic hindering mechanism” (HHM), which could hinder the formation of the water channel, facilitating the formation of the activated state for β2AR.⁴⁴ On the whole, although the dimerization of CXCR4 reduces the main communication flow, the three essential pathways consisting of these important residues above are still retained.

Fig. 6 Rec values of some key residues identified by PSN for the isolated monomer and the protomer in the dimer for the three GPCRs. The height of the histogram represents the rec.

Compared to CXCR4, the impact of dimerization on the allosteric pathway is more significant for β1AR. The protomer of β1AR dimer shares 51.7 percent of the links of the isolated monomer. The monomer contains 49 main paths, which are from the ligand binding sites involved in TM3, TM5, TM6 and TM7 to the G-protein binding sites involved in TM6 and TM7, as observed from Fig. 4. The protomer contains 20 paths, which are mainly from the ligand binding sites of TM3, TM6 and TM7 to the G-protein binding sites associated with TM2, TM3, TM6 and TM7. The observation also shows that the dimerization of β1AR can restrict its main information flow due to the weakness in some links. For example, there is a branch from F299^6.44 to the cytoplasmic sides of TM6 and TM7 in the monomer, acting as one main link. However, the dimerization weakens the link, leading to its disappearance from the main meta-pathway presented in Fig. 5. The protomer in the β1AR dimer has eight pathways commonly shared by the monomer. The eight common paths are from W303^6.48, N329^7.39, Y333^7.43, D121^3.32 to P340^7.50, Y343^7.53, respectively. The remaining 12 paths in the protomer are different from the monomer, in which the differences mainly distribute over the start and end residues and are less involved in the middle region of the helix bundle. All pathways in the protomer and six paths in the monomer contain the same segment of N329^7.39–W303^6.48–F299^6.44–S128^3.39–D87^2.50–N59^1.50, which have a high rec, as revealed by Fig. 5. These pathways also include the highly conserved W6.48 and F6.44, which were considered to be the active switch. The D2.50 can be bound to sodium ion,⁴⁶ which plays a role in participating the regulation of signal transduction of Rh-GPCRs.⁷² The N1.50 is a totally conserved residue⁴⁴ and mutation study on NK1R (A GPCR) showed that N1.50 is essential for three signaling pathways involved in Gq, Gs, and β-arrestin mobilization.⁷³

For the κ-OR, the most significant difference is observed. The protomer only share 44.4% percent of the links of the monomer. The monomer contains 20 main paths, which come from the ligand binding sits of TM6, TM7 to the H8 and cytoplasmic side of TM3, TM6 and TM7. The protomer has 15 paths, which come from the ligand binding sites of TM6 and TM7 to cytoplasmic side of TM5 and TM7, as shown in Fig. 4. The monomer and the protomer share about six paths, which are from W287^6.48, I290^6.51 and I316^7.39 to N326^7.49 and Y330^7.53, respectively. The residues in these common pathways mostly have relatively low rec values, indicating that these common pathways are not important. The six common pathways share one segment consisting of W287^6.48–N322^7.45–N326^7.49–Y330^7.53, which also contains the water-gating W6.48, N7.49 and Y7.53 of NPxxY motif. The relatively conserved N322^7.45 was found to be sensitive in Gα_i2-overexpressed cells in a mutation experiment.⁷⁴ However, the conserved N322^7.45, N326^7.49 and Y330^7.53 in the common segment also present relatively low rec values, as reflected by Fig. 5. The protomer has 9 different paths from its monomer, and their differences mainly locate in the middle link region and the cytoplasmic side, as shown in Fig. 4 and 5. For example, the protomer has one branch from I146^3.40 to the cytoplasmic side of TM5, while in the monomer, the branch is from I146^3.40 to the cytoplasmic side of TM3 and TM6. The residues that are involved in the different paths hold higher rec values than the common path, in particular for the residues S153^3.47 and I242^5.54, as observed from Fig. 6. S120^3.47 in CXCR1 was revealed by mutation experiment to be involved in G-protein binding and the activation of the receptor⁶⁸ while previous MD study⁷⁵ indicates that I290^5.54 plays an important role in breaking the key ionic lock involved in CB1 receptor. Combined with our observation, it can be assumed that the two positions are also important for the dimer of κ-OR. The high rec values of some important residues in the different paths between the monomer and the protomer indicate that the dimerization plays significant role in influencing the communication paths in κ-OR.

As observed above, meta-paths identified by PSN do indeed consist of some important residues or positions and TMs revealed in the signaling. Table S2 in the ESI† lists the positions located in these meta-paths, which simultaneously belong to the 36 positions of the conserved interaction network revealed by Venkatakrishnan.⁶⁴ As can be seen from Table S2,† seven positions identified by the meta-paths belong to the 36 of the monomer and the protomer of CXCR4, in which only one position is different between the monomer and the protomer. For β1AR, 12 positions of the meta-paths belong to the conserved network, one of which is different between the protomer and the monomer. Similarly, the largest difference between the protomer and the monomer is observed for κ-OR. There are 10 and 15 positions belonging to the conserved network for its monomer and protomer, respectively. Only eight positions are same between the monomer and the protomer for κ-OR. This observation indicates that the positions located in the reported conserved interaction network also contribute to the meta-path, although the number of the positions and the type of the residues are different, dependent on the type of the receptors. However, the W6.48, F6.44 and X7.39 positions always exist in the meta-paths and present high rec values, regardless of the monomer or the dimer or the different GPCRs. As reported, the highly-conserved W6.48 and F6.44 serve as the active switch in the activation process.^44,64 Whereas for X7.39, although it is the one non-conserved position, some experiments^76,77 and computational models^78,79 indicated that the Glu288^7.39 of CXCR4 participates in polar interactions with the novel HIV-1 inhibitor and DV1 as well as a salt bridge with the charged N-terminal end of the ligand CXCL12. In addition, it was observed from the crystal structure that Asn^7.39 of β1AR forms H-bonds with the hydroxyl oxygen of the inverse agonist.⁸⁰ Studies on the T355N^7.39 mutant of the 5HT1D receptor demonstrated the role of this Asn residue in recognizing the hydroxyl group of these antagonists.⁸¹ Substitutions of Asn312^7.39 could modify the ligand binding specificity of the β2AR, resulting in a reduction in binding affinity for typical antagonists in the N312T mutant of the 2-adrenergic receptor.⁸² For the three GPCRs in this work, there are very different residues in the non-conserved position X7.39 (E288, N329, I316 for CXCR4, β1AR, κ-OR, respectively) while it was observed that ligands of the three receptors are very different.^83,84 Taken together, we may deduce that the X7.39 position should play, to some extent, a role in ligand specificity for the class A GPCRs. Thus, it is also reasonable to surmise that the meta-pathways very possibly catch the most crucial information flow in the network and these key positions identified should contribute to the allosteric paths in other A GPCR dimers.

On the other hand, in order to gain clear insight to the contribution from TMs, we counted the number of links between TMs involved in the meta-pathways above and the statistical results are illustrated by Fig. 7. For the CXCR4, the distribution of the links in the isolated monomer was similar to those in the protomer, except for the link between TM3 and TM7, further displaying the minor impact of the dimerization on the monomer of CXCR4. For the β1AR, the main distribution regions of the links are also similar between the monomer and the protomer, excepting the formation of H8–TM6 and H8–TM1 as well as the disappearance of TM3–TM2 upon dimerization. Significantly different from CXCR4, however, the marked variations in the number of links between the TMs occur upon dimerization for β1AR, like a significantly weakened link of TM6–TM7 and the enhanced link of TM3–TM2. The results show that the interface of TM1–TM2–H8 indeed influences the links associated with the three regions. Similar to the observations above, the most significant difference between the monomer and the protomer is observed for κ-OR in either the distribution region of the links or the number of the links. The links between TM2 and the other TMs existing in the monomer are almost completely broken upon dimerization, in particular for TM2–TM4. Two new links are formed. One is a link between TM3–TM5, which should be strong since the number of the links is up to four. The other is a single link of TM5–TM6. Different from β1AR with the same TM1–TM2–H8 interface, however, there is no link observed between the TM1, TM2 and H8 for κ-OR. In other words, the formation of the TM1–TM2–H8 interface does not increase the interaction between TM1–TM2–H8 like β1AR. In contrast, the link TM2–H8 is broken upon dimerization.

Fig. 7 A comparison of the connection among TMs in the meta-paths between the isolated monomer and the protomer in the dimer for the three GPCRs. The digit in the line denotes the number of the links and the width of the line is proportional to the number of links.

A comparison of the three types of GPCRs indicates that TM3 is a communication hub among TMs, either for the dimer or the monomer. In fact, it was already revealed that mutations in many positions in TM3 would cause inactivation or constitutive activation of the receptors,⁶⁴ in line with our observations. In addition, although the dimerization enhances some interactions between some helices, more contacts are weakened, which may be attributed to the reason that the dimerization restricts the flexibility of the protomer.

Conclusions

Since it is a new important topic for GPCR oligomers in the GPCR fields and a very limited number of GPCR dimer crystal structures are available so far, there has been an severe lack of knowledge of the structures, dynamics and functional motions of GPCR dimers. In this work, we used the elastic network model (ENM) and the protein structure network (PSN) to study the effect of the dimerization on the dynamics behavior and the allosteric communication path associated with functional mechanism for the three class A GPCR homodimers (viz., CXCR4, β1AR and κ-OR) with the two different interfaces (viz., TM1–TM2–H8 for β1AR and κ-OR, and the TM5–TM6 for CXCR4).

For the three A GPCR dimers similarities, to some extent, were observed for some low-frequency modes derived from NMA analysis, implying similar functional motions for these A GPCR dimers, such as asymmetric motions between the two protomers associated with certain important regions like ECL2 involved in the ligand binding and TM6 involved in G-protein binding. This observation provides support for the ligand binding and asymmetric activation of GPCR dimers through negative cooperation. However, some differences are observed between the monomer and the protomer, dependent on the interface and the contact area. The important anticorrelation motion between the intracellular and extracellular regions in the monomer is retained for the CXCR4 to large extent and β1AR to some extent, but not for κ-OR owing to the limitation of the largest contact area.

The PSN–ENM analysis indicates that some important and conserved positions reported for the monomeric GPCRs still retain their important role in the dimeric unit irrespective of different GPCR types and interfaces, in particular for W6.48 and F6.44 as well as the one non-conserved position X7.39. The dimerization significantly weakens the information flow of the protomer for the three dimers. The weakened extent also follows the order of κ-OR > β1AR > CXCR4. For CXCR4, the dimerization plays a minor role in influencing the communication pathway and TMs participating in the allosteric pathway, in which TM3, TM6 and TM7 still retain important contributions in the dimer and the minor differences induced by the dimerization mainly locate in the cytoplasmic side. For β1AR, however, a more significant impact on the allosteric path is observed, in which the contributions from TM6 and TM7 are significantly weakened upon dimerization, and the dimerization-induced difference mainly distributes in the start and end residues, with a minor fraction in the middle link region. The most significant impact of the dimerization is observed for κ-OR, in which the contribution from TM2 to the main communication paths is significantly weakened, while the contributions from TM3 and TM5 are enhanced. The difference upon dimerization mainly involves in the middle link region and the cytoplasmic side for κ-OR.

On the whole, the effect of dimerization on the dynamics behavior of GPCRs is dependent on the type of interface and the contact area. Significant differences are observed for κ-OR and β1AR with the interface of TM1–TM2–H8, in particular for κ-OR with a greater contact area. The smallest difference is observed for CXCR4 with a TM5–TM6 interface despite its contact area being close to that of β1AR, implying a more functional impact from the TM1–TM2–H8 interface with respect to that of TM5–TM6.

Methods

Structural data

Crystal structures of the three homodimers were obtained from Protein Data Bank (PDBID 3ODU for CXCR4,³⁷ PDBID 4DJH for κ-OR³⁸ and PDBID 4GPO for β1AR³⁹). The monomers are chosen from their corresponding dimers. To prevent the receptor suffering from the “tip effect”, we truncated the unstructured C-terminus in the CXCR4 chemokine receptor at residue 306.⁸⁵

Elastic network model (ENM)

The Elastic Network Model (ENM) has two variants, viz., the Gaussian Network Model (GNM)^86,87 and the Anisotropic Network Model (ANM).⁸⁸ Nodes in the two network models are given by the positions of C_a atoms in the residues and the edges are gained by the distances between the residues i and j. The GNM normal modes are acquired by pseudo-inversion and eigenvalue decomposition of a Kirchhoff matrix Γ (N × N) given by the potentials V_GNM in eqn (1)where R⁰_ij and R_ij are original and instantaneous distance vectors between residues i and j based on their C_α-atom positions, h(x) denotes the Heaviside function, which is equal to 1 if x is positive, and zero otherwise, and R^GNM_cut is the cutoff distance for inter-residue edges, adopted here to be 10 Å.⁵⁹ The establishment of an elastic network model is based on a harmonic approximation of perturbation in the system near the equilibrium position. It can be assumed that the potential V in the initial coordinate is zero and there exists a minimum value. Accordingly, the potential energy can be expanded as a power series, in which the first and second terms are zero; then, a first order approximation is made to gain Γ. The eigenvalues and eigenvectors of Γ are representative of the frequencies of the individual modes and define the shapes of the modes (see ref. 24 for details). The GNM modes are obtained from the transformation of the Kirchhoff matrix, which provide N-dimensional information about the displacement of residues along each mode axis,⁸⁹ but cannot provide the motion direction.

In principle, ANM analysis is also one NMA applied to an ENM, the potential of which is defined as V_ANM in terms of eqn (2).

Herein, R^ANM_cut is adopted to be 15 Å.⁵⁹ The equation is very similar to that of GNM and the major difference is that the distances (scalars), |R_ij| and|R⁰_ij|, replace the vectors R_ij and R⁰_ij in GNM. The 3N × 3N Hessian matrix H can be obtained by deducing V_ANM as it is in GNM. The ANM modes can be deduced from Hessian matrix, which could provide 3-D descriptions of the 3N − 6 internal modes, thus permitting us to evaluate directional preferences.²⁵ Herein, the GNM and ANM calculations were performed on oGNM⁹⁰ and ANM2.0 web servers.⁹¹

PSN–ENM path analysis

In order to perform PSN–ENM analysis, we first need to build a PSG (Protein Structure Graph). In the graph, each amino acid in the protein structure is considered as a node based on its C_α position, and an edge is constructed based on the strength of interaction (I_ij) between residues i and j, which is calculated in terms of eqn (3).⁹²

In eqn (3), n_ij is the number of atom–atom pairs between the side chains of residues i and j within a distance cutoff (4.5 Å); and N_i and N_j are normalization factors for residue types i and j, which were taken from the work of Kannan and Vishveshwara⁹³ on 20 different amino acids in a non-redundant set of protein structures. I_ij is calculated for all node pairs. If I_ij of any residue pair is more than a given interaction strength cutoff I_min, this residue pair will be considered to be interacting and is connected in the PSG. After that, the shortest paths between pairs of nodes in the PSG are searched using Dijkstra’s algorithm.⁹⁴ Then, correlation matrix derived from ENM was adopted to filter these shortest paths. Paths in which at least one node holds a correlation motion with either one of the two extremities based on the 0.6 cutoff value of the correlation coefficient were retained to comprise a pool of the paths. We selected the frequent nodes and links that recur in 20% of the considered path pool to build the global meta-paths. More information about the approach can be found in ref. 63. The PSN–ENM path approach was carried out using WebPSN, a freely available web server built by Seeber and coworkers.⁹⁵

NMA of a subsystem coupled to a dynamic environment

In order to study the impact of the elastic couplings between the protomers in the GPCR dimers modeled as a C_α-only elastic network, the GPCR dimer can be divided into two components (viz., subsystem and environment). One of the two protomers is considered as the subsystem (S) and the other one serves as the environment (E). The Hessian of the whole system is partitioned into four submatrices,⁹⁶ as shown in eqn (4).where H_SS is the Hessian submatrix for the subsystem, H_EE is that of the environment and H_SE (or H_ES) refers to the coupling between the subsystem and the environment. An effective Hessian for the subsystem (H^eff_SS) can be obtained in terms of eqn (5).

H^eff_SS = H_SS − H_SEH_EE⁻¹H_ES (5)

The NMA derived from H^eff_SS could effectively describe the collective dynamics of the subsystem in the presence of elastic coupling to the environment. The approach is implemented by Prody.⁹⁷

Subspace overlap

The overlap between the dynamic conformational spaces can quantify similarities of ANM modes from the two different systems, which can be described by two sets of m and n modes in the two different systems, u_k and v_l, respectively. For all m × n pairs of modes, the overlap between every pair of the modes is calculated by the correlation cosine, |cos(u_k·v_l)|,^57,58 in which k and l denote kth and lth eigenvectors, respectively.

Acknowledgements

This project is supported by the National Science Foundation of China (Grant No. 21273154, 21573151) and the Sichuan Province Science and Technology Support Program (Grant No. 2015GZ0193).

Notes and references

1. D. M. Rosenbaum, S. G. Rasmussen and B. K. Kobilka, Nature, 2009, 459, 356–363.
2. C. Zhang, T. Hou, Z. Feng and Y. Li, Curr. Comput.-Aided Drug Des., 2013, 9, 60–75.
3. I. Kufareva, B. S. Stephens, L. G. Holden, L. Qin, C. Zhao, T. Kawamura, R. Abagyan and T. M. Handel, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, E5363–E5372.
4. J. F. White, J. Grodnitzky, J. M. Louis, L. B. Trinh, J. Shiloach, J. Gutierrez, J. K. Northup and R. Grisshammer, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 12199–12204.
5. M. R. Whorton, M. P. Bokoch, S. G. Rasmussen, B. Huang, R. N. Zare, B. Kobilka and R. K. Sunahara, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 7682–7687.
6. G. Khelashvili, K. Dorff, J. Shan, M. Camacho-Artacho, L. Skrabanek, B. Vroling, M. Bouvier, L. A. Devi, S. R. George, J. A. Javitch, M. J. Lohse, G. Milligan, R. R. Neubig, K. Palczewski, M. Parmentier, J. P. Pin, G. Vriend, F. Campagne and M. Filizola, Bioinformatics, 2010, 26, 1804–1805.
7. K. Palczewski, Trends Biochem. Sci., 2010, 35, 595–600.
8. S. R. George, B. F. O’Dowd and S. P. Lee, Nat. Rev. Drug Discovery, 2002, 1, 808–820.
9. S. Ferre, V. Casado, L. A. Devi, M. Filizola, R. Jockers, M. J. Lohse, G. Milligan, J. P. Pin and X. Guitart, Pharmacol. Rev., 2014, 66, 413–434.
10. V. C. Vsevolod Katritch and R. C. Stevens, Annu. Rev. Pharmacol. Toxicol., 2013, 53, 531–556.
11. X. Rovira, J. P. Pin and J. Giraldo, Trends Pharmacol. Sci., 2010, 31, 15–21.
12. L. F. Agnati, K. Fuxe and S. Ferré, Trends Biochem. Sci., 2005, 30, 188–193.
13. N. S. Schonenbach, S. Hussain and M. A. O’Malley, Wiley Interdiscip. Rev.: Nanomed. Nanobiotechnol., 2015, 7, 408–427.
14. F. Fanelli and P. G. De Benedetti, Chem. Rev., 2011, 111, PR438–PR535.
15. Q. Bai and X. Yao, Sci. Rep., 2016, 6, 21763.
16. M. Neri, S. Vanni, I. Tavernelli and U. Rothlisberger, Biochemistry, 2010, 49, 4827–4832.
17. D. Rodriguez and H. Gutierrez-de-Teran, Proteins, 2012, 80, 1919–1928.
18. A. Bruno, A. E. Guadix and G. Costantino, J. Chem. Inf. Model., 2009, 49, 1602–1616.
19. A. Bruno, C. Beato and G. Costantino, Future Med. Chem., 2011, 3, 665–681.
20. S. Mondal, J. M. Johnston, H. Wang, G. Khelashvili, M. Filizola and H. Weinstein, Sci. Rep., 2013, 3, 233–239.
21. X. Periole, A. M. Knepp, T. P. Sakmar, S. J. Marrink and T. Huber, J. Am. Chem. Soc., 2012, 134, 10959–10965.
22. D. Provasi, M. B. Boz, J. M. Johnston and M. Filizola, PLoS Comput. Biol., 2015, 11(3), e1004148.
23. X. Prasanna, A. Chattopadhyay and D. Sengupta, Biophys. J., 2014, 106, 1290–1300.
24. I. Bahar, T. R. Lezon, A. Bakan and I. H. Shrivastava, Chem. Rev., 2010, 110, 1463–1497.
25. E. Eyal, A. Dutta and I. Bahar, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 426–439.
26. U. Bastolla, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 488–503.
27. A. Uyar, N. Kantarci-Carsibasi, T. Haliloglu and P. Doruker, Biophys. J., 2014, 106, 2656–2666.
28. E. Fuglebakk, S. P. Tiwari and N. Reuter, Biochim. Biophys. Acta, 2015, 1850, 911–922.
29. K. V. Brinda and S. Vishveshwara, Biophys. J., 2005, 89, 4159–4170.
30. S. Vishveshwara, A. Ghosh and P. Hansia, Curr. Protein Pept. Sci., 2009, 10, 146–160.
31. B. M. S. Ghosha and S. Vishveshwaraa, Curr. Protein Pept. Sci., 2015, 16, 193.
32. V. P. Souza, C. M. Ikegami, G. M. Arantes and S. R. Marana, FEBS J., 2016, 6, 1124–1138.
33. F. Fanelli, A. Felline, F. Raimondi and M. Seeber, Biochem. Soc. Trans., 2016, 44, 613–618.
34. M. S. Vijayabaskar and S. Vishveshwara, Biophys. J., 2010, 99, 3704–3715.
35. S. Mariani, D. Dell’Orco, A. Felline, F. Raimondi and F. Fanelli, PLoS Comput. Biol., 2013, 9, e1003207.
36. J. H. Park, P. Scheerer, K. P. Hofmann, H. W. Choe and O. P. Ernst, Nature, 2008, 454, 183–187.
37. B. Wu, E. Y. Chien, C. D. Mol, G. Fenalti, W. Liu, V. Katritch, R. Abagyan, A. Brooun, P. Wells, F. C. Bi, D. J. Hamel, P. Kuhn, T. M. Handel, V. Cherezov and R. C. Stevens, Science, 2010, 330, 1066–1071.
38. H. Wu, D. Wacker, M. Mileni, V. Katritch, G. W. Han, E. Vardy, W. Liu, A. A. Thompson, X. P. Huang, F. I. Carroll, S. W. Mascarella, R. B. Westkaemper, P. D. Mosier, B. L. Roth, V. Cherezov and R. C. Stevens, Nature, 2012, 485, 327–332.
39. J. Huang, S. Chen, J. J. Zhang and X. Y. Huang, Nat. Struct. Mol. Biol., 2013, 20, 419–425.
40. P. Petrone and V. S. Pande, Biophys. J., 2006, 90, 1583–1593.
41. S. Mahajan and Y. H. Sanejouand, Arch. Biochem. Biophys., 2015, 567, 59–65.
42. J. A. Ballesteros and H. Weinstein, Methods Neurosci., 1995, 25, 366–428.
43. L. W. Yang and I. Bahar, Structure, 2005, 13, 893–904.
44. B. G. Tehan, A. Bortolato, F. E. Blaney, M. P. Weir and J. S. Mason, Pharmacol. Ther., 2014, 143, 51–60.
45. J. L. Millergallacher, R. Nehmé, T. Warne, P. C. Edwards, G. F. Schertler, A. G. Leslie and C. G. Tate, PLoS One, 2014, 9, e92727.
46. W. Liu, E. Chun, A. A. Thompson, P. Chubukov, F. Xu, V. Katritch, G. W. Han, C. B. Roth, L. H. Heitman and A. P. Ijzerman, Science, 2012, 337, 232–236.
47. J. H. Yun, M. Kim, K. Kim, D. Lee, Y. Jung, D. Oh, Y. J. Ko, A. E. Cho, H. S. Cho and W. Lee, Biochim. Biophys. Acta, Biomembr., 2015, 1848, 1294–1302.
48. Y. Han, I. S. Moreira, E. Urizar, H. Weinstein and J. A. Javitch, Nat. Chem. Biol., 2009, 5, 688–695.
49. J. R. Lane, P. Donthamsetti, J. Shonberg, C. J. Draper-Joyce, S. Dentry, M. Michino, L. Shi, L. López, P. J. Scammells, B. Capuano, P. M. Sexton, J. A. Javitch and A. Christopoulos, Nat. Chem. Biol., 2014, 10, 745–752.
50. A. González, T. Perezacle, L. Pardo and X. Deupi, PLoS One, 2011, 6, e23815.
51. R. O. Dror, A. C. Pan, D. H. Arlow, D. W. Borhani, P. Maragakis, Y. Shan, H. Xu and D. E. Shaw, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 13118–13123.
52. K. Gherbi, L. T. May, J. G. Baker, S. J. Briddon and S. J. Hill, FASEB J., 2015, 29, 2859–2871.
53. M. Damian, A. Martin, D. Mesnier, J. P. Pin and J. L. Banères, EMBO J., 2006, 25, 5693–5702.
54. N. J. Smith and G. Milligan, Pharmacol. Rev., 2010, 62, 701–725.
55. D. S. Goodsell and A. J. Olson, Annu. Rev. Biophys. Biomol. Struct., 2000, 29, 105–153.
56. G. Maksay and O. Tőke, Prog. Biophys. Mol. Biol., 2014, 114, 153–169.
57. E. Marcos, R. Crehuet and I. Bahar, PLoS Comput. Biol., 2010, 6, e1000738.
58. E. Marcos, R. Crehuet and I. Bahar, PLoS Comput. Biol., 2011, 7, e1002201.
59. B. Isin, A. J. Rader, H. K. Dhiman, J. Klein-Seetharaman and I. Bahar, Proteins, 2006, 65, 970–983.
60. M. Damian, S. Mary, M. Maingot, C. M’Kadmi, D. Gagne, J. P. Leyris, S. Denoyelle, G. Gaibelet, L. Gavara, M. Garcia de Souza Costa, D. Perahia, E. Trinquet, B. Mouillac, S. Galandrin, C. Gales, J. A. Fehrentz, N. Floquet, J. Martinez, J. Marie and J. L. Baneres, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 1601–1606.
61. S. Marsango, G. Caltabiano, C. Pou, M. J. Varela Liste and G. Milligan, J. Biol. Chem., 2015, 290, 15146–15162.
62. J. M. Johnston, H. Wang, D. Provasi and M. Filizola, PLoS Comput. Biol., 2012, 8, e1002649.
63. F. Raimondi, A. Felline, M. Seeber, S. Mariani and F. Fanelli, J. Chem. Theory Comput., 2013, 9, 2504–2518.
64. A. J. Venkatakrishnan, X. Deupi, G. Lebon, C. G. Tate, G. F. Schertler and M. M. Babu, Nature, 2013, 494, 185–194.
65. F. Raimondi, A. Felline and F. Fanelli, J. Chem. Inf. Model., 2015, 55, 1878–1893.
66. F. Fanelli, A. Felline and F. Raimondi, Methods Cell Biol., 2013, 117, 43–61.
67. V. Isberg, S. Mordalski, C. Munk, K. Rataj, K. Harpsoe, A. S. Hauser, B. Vroling, A. J. Bojarski, G. Vriend and D. E. Gloriam, Nucleic Acids Res., 2016, 44, D356–D364.
68. X. Han, Y. Feng, X. Chen, C. Gerard and W. A. Boisvert, FEBS Open Bio, 2015, 5, 182–190.
69. S. Yuan, Z. Hu, S. Filipek and H. Vogel, Angew. Chem., Int. Ed. Engl., 2014, 54, 556–559.
70. T. W. Schwartz, T. M. Frimurer, B. Holst, M. M. Rosenkilde and C. E. Elling, Annu. Rev. Pharmacol., 2006, 46, 481–519.
71. X. Deupi, P. Edwards, A. Singhal, B. Nickle, D. Oprian, G. Schertler and J. Standfuss, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 119–124.
72. M. S. Parker, Y. Y. Wong and S. L. Parker, Amino Acids, 2008, 35, 1–15.
73. L. Valentinhansen, T. M. Frimurer, J. Mokrosinski, N. D. Holliday and T. W. Schwartz, J. Biol. Chem., 2015, 290, 24495–24508.
74. F. Yan, P. D. Mosier, R. B. Westkaemper and B. L. Roth, Biochemistry, 2008, 47, 1567–1578.
75. J. Y. Shim and L. Padgett, in G-Protein Coupled Receptors: Structure, ed. P. M. Conn, 2013, vol. 520, pp. 337–355.
76. W. T. Choi, S. M. Tian, C. Z. Dong, S. Kumar, D. X. Liu, N. Madani, J. An, J. G. Sodroski and Z. W. Huang, J. Virol., 2005, 79, 15398–15404.
77. S. M. Tian, W. T. Choi, D. X. Liu, J. Pesavento, Y. L. Wang, J. An, J. G. Sodroski and Z. W. Huang, J. Virol., 2005, 79, 12667–12673.
78. P. Tamamis and C. A. Floudas, Biophys. J., 2013, 105, 1502–1514.
79. P. Tamamis and C. A. Floudas, J. Chem. Inf. Model., 2014, 54, 1174–1188.
80. S. Vanni, M. Neri, I. Tavernelli and U. Rothlisberger, PLoS Comput. Biol., 2011, 7(1), e1001053.
81. R. A. Glennon, M. Dukat, R. B. Westkaemper, A. M. Ismaiel, D. G. Izzarelli and E. M. Parker, Mol. Pharmacol., 1996, 49, 198–206.
82. S. Suryanarayana and B. K. Kobilka, Mol. Pharmacol., 1993, 44, 111–114.
83. T. Warne, R. Moukhametzianov, J. G. Baker, R. Nehmé, P. C. Edwards, A. G. Leslie, G. F. Schertler and C. G. Tate, Nature, 2011, 469, 241–244.
84. M. Filizola and L. A. Devi, Nature, 2012, 485, 314–317.
85. D. Kolan, G. Fonar and A. O. Samson, Proteins, 2014, 82, 579–586.
86. I. Bahar, A. R. Atilgan and B. Erman, Folding Des., 1997, 2, 173–181.
87. T. Haliloglu, I. Bahar and B. Erman, Phys. Rev. Lett., 1997, 79, 3090–3093.
88. A. R. Atilgan, S. R. Durell, R. L. Jernigan, M. C. Demirel, O. Keskin and I. Bahar, Biophys. J., 2001, 80, 505–515.
89. D. Riccardi, C. Qiang and G. N. Phillips, Biophys. J., 2009, 96, 464–475.
90. L. W. Yang, A. J. Rader, X. Liu, C. J. Jursa, S. C. Chen, H. A. Karimi and I. Bahar, Nucleic Acids Res., 2006, 34, 24–31.
91. E. Eyal, G. Lum and I. Bahar, Bioinformatics, 2015, 31, 1487–1489.
92. M. Bhattacharyya and S. Vishveshwara, Biochemistry, 2011, 50, 6225–6236.
93. N. Kannan and S. Vishveshwara, J. Mol. Biol., 1999, 292, 441–464.
94. E. W. Dijkstra, Numer. Math., 1959, 1, 269–271.
95. M. Seeber, A. Felline, F. Raimondi, S. Mariani and F. Fanelli, Bioinformatics, 2015, 31, 779–781.
96. W. Zheng and B. R. Brooks, Biophys. J., 2005, 89, 167–178.
97. A. Bakan, L. M. Meireles and I. Bahar, Bioinformatics, 2011, 27, 1575–1577.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra18243g

---