Active-State Models of Ternary GPCR Complexes: Determinants of Selective Receptor-G-Protein Coupling

Based on the recently described crystal structure of the β2 adrenergic receptor - Gs-protein complex, we report the first molecular-dynamics simulations of ternary GPCR complexes designed to identify the selectivity determinants for receptor-G-protein binding. Long-term molecular dynamics simulations of agonist-bound β2AR-Gαs and D2R-Gαi complexes embedded in a hydrated bilayer environment and computational alanine-scanning mutagenesis identified distinct residues of the N-terminal region of intracellular loop 3 to be crucial for coupling selectivity. Within the G-protein, specific amino acids of the α5-helix, the C-terminus of the Gα-subunit and the regions around αN-β1 and α4-β6 were found to determine receptor recognition. Knowledge of these determinants of receptor-G-protein binding selectivity is essential for designing drugs that target specific receptor/G-protein combinations.


Introduction
G-protein-coupled receptors (GPCRs) are proteins that enable signal transduction through biological membranes. The more than 800 GPCRs (including receptors for olfaction and taste) constitute the largest family of membrane proteins in the human genome [1]. GPCRs show pronounced structural variety in their binding pocket and can thus be activated by diverse extracellular signals including photon-induced changes in ligand conformation, small molecules, peptides and proteins [2]. Agonist binding causes structural rearrangements in the intracellular part of the receptor [3][4][5][6][7][8][9] that enable binding of a heterotrimeric G-protein and thus formation of the ternary complex consisting of agonist, receptor and G-protein [10]. The ternary complex induces the transmission of signals that activate both distinct physiological processes involving sensory impressions such as vision, smell and taste and neurological, cardiovascular, endocrine and reproductive functions that make GPCRs (and G-proteins) important targets for drug design [11].
After the structural characterization of the b 2 -adrenergic receptor (b2AR) bound to an antagonist [12,13] and the first agonist-b2AR complexes [14,15], the crystal structure of the b2AR together with its signal-transducing G s -protein was determined by Brian Kobilka and his team [16]. This spectacular work offers important structural insights into the nucleotide-free ternary signaling complex that will be important for the rational, structurebased design of biochemical or computational studies to investigate ternary complexes. The G-protein as an intracellular binding partner has been shown to be a prerequisite for capturing the fullyactivated state of a GPCR in a crystal, since the recently determined structure of the b2AR bound to our agonist FAUC50 indicated a receptor conformation that was similar to the antagonist-bound form [14]. Only in the presence of a G-protein simulating nanobody [15] or the G-protein itself [16], could the rigid body movements described above be observed. Recently, NMR experiments investigating the dynamic behavior of b2AR emphasized the fundamental role of an intracellular binding partner in the stabilization of a fully-activated receptor conformation [17].
The crystal structure provides a physiological, atomistic template of a fully-activated G-protein-coupled receptor bound to and stabilizing a nucleotide-free G-protein. It represents a valuable template for homology modeling studies that explore high-affinity active-state binding sites of GPCR-G-protein complexes. Active-state homology models can be of great importance for identifying new agonist lead-structures, for example in docking campaigns [18]. Because many GPCRs can bind multiple Gprotein-subtypes, models of individual receptor-G-protein complexes are needed to design functionally selective drugs inducing the activation of a particular G-protein to a higher extend than coupling to alternative G-protein subtypes.
Herein, we describe the first active-state homology model of a G-protein-coupled receptor in complex with its preferred Gprotein based on the crystal structure of the b 2 -adrenergic receptor in complex with the G s -protein [16]. In order to identify the amino acids responsible for coupling selectivity between GPCRs and Gproteins, we examined the protein-protein interface of two different ternary complexes, the agonist-bound b2AR-Ga s crystal structure and, based on the b2AR-Ga s -structure, two homology models of the dopaminergic D 2 receptor (D2R), a drug target of particular interest for the treatment of neuropsychiatric disorders including Parkinson's disease and schizophrenia [19], in complex with dopamine and Ga i1 . We carried out one ms moleculardynamics simulations in a hydrated bilayer built of dioleoylphosphatidylcholine-lipids (DOPC) for each, and investigated the receptor G-protein interface by computational alanine scanning mutagenesis.

Results and Discussion
Active-state Homology Models of D2R-Ga i According to Kobilka et al. [16], the ''active state of a GPCR can be defined as that conformation that couples to and stabilizes a nucleotide-free Gprotein. '' We therefore used the crystal structure of the b2AR-Ga scomplex (PDB-ID: 3SN6) as a starting point for active-state homology models of D2R in complex with the nucleotide-free state of Ga i . We created alignments for the separated receptors and the G-proteins, combined them and subsequently started the modeling process using MODELLER 9v4. A more detailed description of the modeling process is provided in the Methods section. The models exhibited two different rotamer conformations of residue His393 6.55 in D2R with the side chain of histidine pointing either to the extracellular or to the intracellular part of the receptor (Figure 1). His393 6.55 has been shown to play a significant role in ligand binding and signaling bias at dopaminergic receptors [20][21][22] and that, in principal, both conformations are possible [23]. Therefore in the following studies, we decided to select two models of the D2R-Ga i -complex with both rotamer conformations of His393 6.55 , which are referred to in the following as D2 Up R-Ga i and D2 Down R-Ga i . The physiological agonist dopamine was docked manually into D2 Up R-Ga i and D2 Down R-Ga i in a way that the positively charged ammonium head group forms a salt bridge to Asp114 3.32 and that hydrogen bonds between the catechol moiety of dopamine and the side chains of Ser193 5.42 and Ser197 5.46 of D2R become feasible ( Figure 1). These serine residues, Ser193 5.42 and Ser197 5.46 , together with Ser194 5.43 , have been shown to be crucial for highaffinity catecholamine binding and for an effective receptor-Gprotein coupling [24,25].
Agonist binding of GPCRs leads to major structural changes within the receptors and the G-proteins that are consistent with the conformation of our active-state D2R-Ga i -complexes ( Figure  S1).

Molecular-dynamics Simulations
Three ternary complexes, b2AR-BI167107-Ga s and D2 Up R/ D2 Down R-dopamine-Ga i , were successfully embedded into a hydrated DOPC-bilayer. We cleared a space for the initial insertion of the protein structures into the bilayer by removing DOPC-molecules from the bulk of the membrane ( Figure S2a). A careful equilibration procedure was used to close the resulting gap between GPCRs and DOPC-residues ( Figure S2b, c) without water molecules flooding this gap. The resulting complexes were subsequently submitted to molecular-dynamics (MD) simulations for one ms each, with the interior of the DOPC-bilayer remaining free of water throughout the simulations ( Figure S3). The long simulation time of one ms for each complex was chosen to ensure the formation of sufficiently stable amino-acid contacts between the proteins in order to be able to elucidate amino acids that appear in the interface of GPCRs and G-proteins reliably.
All complexes remained very stable throughout the MD simulations showing low RMSD values for every member of the ternary complexes ( Figure S4). As the G-proteins were not stabilized by membrane lipids, they showed higher atomic fluctuations than the receptor moieties ( Figure S5). Substantial mobility was observed for the helical subunits of Ga s and Ga i , Ga s AH and Ga i AH, which have previously been shown to become highly flexible in their nucleotide-free state [26,27]. Comparing the atomic fluctuations of the two D2R-Ga i -complexes, we observed higher values for the D2 Up R-Ga i -simulation ( Figure S5), which were connected to a whole-body movement of Ga i starting at the lower part of the a5-helix, but leaving the majority of a5 and its C-terminus unaffected ( Figure S6a, b). The movement of Ga i originates in the enhanced flexibility of open ends in the N-terminal IL3, which is mainly associated with the absence of the bulk of IL3 ( Figure S5). This enhanced flexibility causes a loss of ionic interactions between residues from the Nterminal part of IL3 and residues from the area around a4-b6, which finally results in a displacement of Ga i around helix a4 in the D2 Up R-Ga i -simulation compared to the D2 Down R-Ga i -complex. As this conformation appeared to be stable for the remainder of the simulation and did not lead to the separation of D2R and Ga i ( Figure S4, S7), we continued investigating both D2R-Ga icomplexes. Additionally, our data give no indication for any displacements of GPCRs and G-proteins other than the one described for the D2 Up R-Ga i -complex.
The agonists BI167107 and dopamine in the b2AR-Ga scomplex and in the D2R-Ga i -complexes, respectively, are largely enclosed in their binding pockets. In the b2AR-Ga s -complex, BI167107 maintained its interactions with residues of TM2, TM3, TM5, TM6 and TM7, most of which were already present in the crystal structure (Figure 2a, b). In the case of the D2R-Ga icomplexes, dopamine showed a different orientation of its catechol moiety within the binding pockets. Whereas only the meta-hydroxy group of dopamine formed a hydrogen bond to Ser193 5.42 in the D2 Down R-Ga i -complex, both, the metaand para-hydroxy groups of dopamine were involved in the formation of hydrogen bonds to Ser193 5.42 and Ser197 5.46 of D2 Up R, respectively (Figure 2c, d).
This behavior may be associated with changes in the rotamer conformation of residue His393 6.55 throughout the D2R-Ga isimulations, where its side chain adopts three distinct dihedral angles, referred to as states 1, 2 and 3 ( Figure 3). In state 1, the side chain of His393 6.55 points towards the intracellular site of the receptor into the direction of TM7 (the initial conformation of the D2 Down R-Ga i -complex), where it is stabilized by an interaction to Figure 1. Initial conformation of dopamine in the D2R-Ga icomplexes. The backbone of D2R is shown as green ribbon, with important amino acids (indicated as green sticks) that stabilize the ligand dopamine in its initial conformation. Dopamine is represented as orange sticks and stabilized by ionic interactions to D114 3.32 and hydrogen bonds to S193 5.42 and S197 5.46 . The second conformation of residue H393 6.55 is shown as red sticks. doi:10.1371/journal.pone.0067244.g001 residue Tyr408 7.35 of upper TM7. State 2 shows the side chain of histidine pointing towards the extracellular part of the receptor (the initial conformation of the D2 Up R-Ga i -complex and the one observed in the crystal structure of D3R), where it regains spatial proximity to Tyr408 7.35 of TM7. The side chain is again oriented towards the intracellular site of the receptor in state 3, but now points in the direction of TM5, which enables a hydrogen bond to be formed to residue Ser193 5. 43 . We assume that the dihedral angle of His393 6.55 causes structural differences within the binding pocket of D2R, which lead to different conformations with respect to ligand binding. Structural connections between His 6.55 , Tyr 7.35 , TM5-serines and ligands that are able to discriminate between different downstream signaling pathways have been shown to be involved in biased signaling [23]. The agonist dopamine, which cannot cause functional selectivity, does not prevent the side chain of His393 6.55 from cycling between its possible rotamer conformations. Sterically more demanding ligands may lock His393 6.55 in one distinct rotamer conformation and thus trigger the activation of one distinct pathway. Therefore, further MDsimulations with selected ligands are necessary to elucidate the impact of His393 6.55 on functionally selective signaling.

The Receptor-G-protein Interface
Our ms MD-simulations were carried out in order to identify stable amino-acid contact sites between the receptors and their Gproteins that are maintained for long periods. Early experimental work, which focused on elucidating the interface between rhodopsin and its G-protein transducin using synthetic peptides that correspond to different regions of rhodopsin and transducin, identified the intracellular loops 2 and 3, the junction between TM7 and helix 8 of rhodopsin [28] and the area around a4-b6 and the C-terminal helix of transducin's Ga subunit, Ga t [29], as important contact sites between the two binding partners. These contact areas were further strengthened by a disulfide cross-linking study using the muscarinic M3 receptor and Ga q [30]. A first structural glimpse of the amino acids involved in binding GPCRs to G-proteins was provided by crystallizing light-activated opsin together with a synthetic peptide (GaCT, residues ILENLKDCGLF) derived from the C-terminus of Ga t [31]. By mutating the residues in GaCT into the corresponding amino acids of Ga s , we were able to delineate its interactions with b2AR [9]. Now, with the crystal structure of an entire ternary b2AR-Ga s complex at hand, we have an excellent framework for investigating active-state models of structurally unknown ternary GPCRcomplexes via computational methods. The trajectories of the MD simulations were therefore screened for amino-acid contacts between the receptors and the appropriate G-proteins. The receptor-G-protein interfaces are shown in Figure 4 as individual alignments for the receptors and for the G-proteins. Amino acids are highlighted in the alignment when at least one atom of an amino acid approaches at least one atom of another amino acid closer than 3.5 Å and when this interaction is found in more than 50% of the simulation. Detailed connection tables are provided in the (Table S1, S2, S3).
The receptor-G-protein interface of these fully-activated, nucleotide-free ternary complexes is comprised of homologous regions within the b2AR-Ga s -complex and the D2 Up/Down R-Ga icomplexes. The amino-acid contacts within the two D2R-Ga i simulations were found to be highly congruent, despite the differences concerning the displacement of Ga i discussed above ( Figure S6). GPCR contacts include the area around IL2, the Nand C-terminal parts of IL3 and the junction of TM7 and helix 8. The latter area only emerged as a contact region during the MD simulations and is not visible in the crystal structure of the b2AR-G s -complex. This observation underlines the importance of dynamic techniques such as MD simulation, which are not limited to a static snapshot of the protein. The G-protein contact regions consist of the aNb1-loop, the area around b2-b3, the area around a4-b6 (with different distributions of the contact residues for the b2AR-Ga s -and the D2 Up/Down R-Ga i -complexes) and the Cterminal a5-helix together with its C-terminus.
Additional information about the receptor-G-protein interfaces is given by highlighting the individual residues that appear in these interfaces with different colors that show the number of individual contacts from one residue to others: the darker the color (from yellow over green to blue) the more neighbors an amino acid has and the more important it is likely to be for receptor-G-protein coupling. Thus, the C-terminal domain of Ga, where high densities of tightly packed amino acids occur, can be assigned an outstanding role for complex stabilization and coupling selectivity arising from the G-protein. This is because the C-terminal a5helix together with its extreme C-terminus is incorporated in the cavity formed by the outward movement of TM6 during receptor activation, which enables pronounced interactions with all of the contact regions of the GPCRs depicted. On the side of the receptors, we observed pronounced interactions for residues belonging to the areas around IL2 and the junction of the distal part of TM5 connected to the N-terminal part of IL3.

Computational Alanine-scanning Mutagenesis
To elucidate the importance of each amino acid that appears in the interface between receptors and G-proteins, we carried out computational alanine-scanning mutagenesis of the b2AR-Ga sand the D2 Up R/D2 Down R-Ga i -interfaces. This approach has been shown to be a valuable tool for estimating the contribution of individual amino acids to the stabilization of protein-protein interactions [32] and to be able to reproduce experimental investigations qualitatively [33]. We therefore used the MM-GBSA-method (Molecular Mechanics-Generalized Born Surface Area) [34], implemented in MMPBSA.py [35], to calculate the relative binding free energy changes (DDG) between alaninemutant complexes and the corresponding wild-type complexes in order to identify so-called hot-spot residues within the GPCR-Gprotein interfaces that contribute to both coupling affinity and selectivity. (atoms: C-CA-CB-CG) is depicted as green and red lines for the D2 Down R-Ga i -and the D2 Up R-Ga i -simulations, respectively. The right column shows representative snapshots taken from the D2R-Ga i -simulations and visualizes the interactions of residue His393 6.55 with amino acids S193 5.43 and Y408 7.35 depending on its dihedral angle (orange: state 1; purple: state 2; dark-cyan: state 3). Helices 5, 6 and 7 are shown as ribbons, whereas the amino acids are represented as sticks. Additionally, state 2 shows the conformation of residue His 6.55 taken from the crystal structure of the dopaminergic D 3 receptor, as grey sticks. doi:10.1371/journal.pone.0067244.g003 In a first step, we omitted water and membrane molecules and calculated the binding free energies (DG) of the b2AR-Ga s -and the D2 Up R/D2 Down R-Ga i -interfaces using the GBSA-method within MMPBSA.py in order to prove that the complex is energetically favorable and that the energy values remain generally consistent over the time scales investigated. Conformationally stable time periods within the three ternary complex trajectories were identified based on RMS deviations ( Figure S4) and used to generate the required trajectories for the receptor-and the Gprotein-parts with intervals of 500 ps between snapshots. Our calculations showed consistently negative DG-values for the systems on the time scales investigated, which indicates energetically favorable interactions between receptors and G-proteins ( Figure S8). We subsequently performed computational alaninescanning mutagenesis for the amino acid residues within the receptor-G-protein interfaces of b2AR-Ga s -and the D2 Up R/ D2 Down R-Ga i -complexes that are highlighted in Figure 4, except for alanine-, glycine-and the C-terminal residues L380 and F354 from Ga s and Ga i , respectively. In cases where only one amino acid of the D2 Up R/D2 Down R-Ga i -complexes constitutes a contact residue, we nevertheless performed alanine scanning on both amino acids. Important results of the alanine scan are shown in Figure 5, the complete results are provided in the (Table S4). In general, a positive value for the binding free energy change (DDG) is associated with an amino acid that contributes to stabilizing the ternary complex, and vice versa.
For the b2AR-Ga s -system, we found that residues R131, I135, F139, Q229, K232, I233, E237, K270 and R333 from b2AR and H41, Y344, D367, I369, Q370, R371, H373, L374, R375, Y477, E378 and L379 from Ga s stabilize the receptor-G-protein interface. Among these residues, F139 from IL2, Q229 and E237 from TM5-IL3 and R333 from helix 8 have been found to Figure 4. Alignment of the amino-acid contacts between receptors and G-proteins. Individual alignments for the receptors and the Gproteins are shown. A colored background indicates that the residue forms contacts to other amino acids (yellow: 1 or 2 contacts; green: 3 or 4 contacts; blue: at least 5 contacts). Red letters indicate residues involved in ionic interactions, whereas dotted underlines indicate contacts present in the crystal structure of b2AR-Ga s . doi:10.1371/journal.pone.0067244.g004 Figure 5. Summary of selectivity determining amino acids within the b2AR-Ga s -and the D2R-Ga i -complexes and representative values of the alanine scanning mutagenesis. The grey columns in the middle refer to the regions within GPCRs and G-proteins, to which the mentioned amino acids belong. Amino acids in italic letters have not been mutated in the computational alanine scanning (n.d.). Blue, green and red bars show the binding free energy differences of the alanine scanning mutagenesis for the b2AR-Ga s complex and the D2 Down R-Ga i and the D2 Up R-Ga i -complexes, respectively. The orange and red rectangles besides the amino acids correspond to hydrophobic or polar interactions to other residues, respectively. doi:10.1371/journal.pone.0067244.g005 be of major importance for b2AR, whereas D367 and R375 from a5 and E378 from the C-terminus of the Ga-subunit are important for Ga s . F139 interacts tightly with a hydrophobic pocket comprised of residues H41, V203, F205, F362, C365, R366 and I369 from Ga s (Figure 6f, Table S1) and thus stabilizes the interface of IL2, aN-b1, b2-b3 and a5. It has been shown that mutating F139 to alanine in b2AR prevents activation of adenylyl cyclase by Ga s and that, in general, a bulky amino acid is necessary in this position for effective receptor-G-protein coupling [36]. Residue Q229 from the N-terminal IL3 forms the center of an extended hydrogen-bond network to residues D367, Q370 and R371 of the a5-helix of Ga s and K232 from TM5 of b2AR (Figure 6e, Table S1). Residue E237 from IL3 and R333 from H8 of b2AR are involved in salt bridges to residues R375 from the a5helix and E378 from the C-terminus of Ga s , respectively (Figure 6c, Table S1).
For the two D2R-Ga i -simulation systems, amino acids important for receptor-G-protein-binding were found to be, in general, qualitatively comparable between D2 Down R-Ga i and D2 Up R-Ga i . The main difference is caused by the movement of Ga i within the D2 Up R-Ga i -complex discussed above, which weakens the interactions between residues from the extreme N-terminal IL3 and the area around a4-b6. Taken together, residues R132, V136, M140, Y142, R145, R150, R219, R222, K226 Down , R227 Down , K367 and K370 from D2R and residues E25, E28, E308 Down , D315, E318, D341, I344, L348, D350 and L353 from Ga i were revealed to be important for the stability of the complexes. The most interesting observation within the interface of D2R and Ga i is the density of positively charged amino acids from the receptor and of negatively charged amino acids from the G-protein, which mainly form salt bridges to each other. Salt bridges involve residues from IL2/TM4 (R145, R150) and TM5/IL3 (R219, R222, K226 Down , R227 Down ) of D2R, which are connected to residues from aN-b1 (E25, E28) and a4-b6 (E308 Down , D315 Down , E318 Down )/a5 (D341) of Ga i , respectively (Figure 6b, g). The importance of basic amino acids of D2R, which interact with negatively charged residues from Ga i , is emphasized by the observation that the alanine scanning mutagenesis for basic amino acids from Ga i (R24, R32, K192, K345, K349) finds a destabilizing effect on the receptor-G-protein-interface. Our results indicate that basic residues from TM6 (K367, K370), despite not forming contacts to acidic amino acids from Ga i (Table S2, S3), participate in stabilizing the receptor-G-protein interface. This can be attributed to interactions with C-terminal residues of Ga i , especially F354, where a cation-p-interaction can be formed. As MMPBSA.py does not allow alanine scans for terminal residues, it was not possible to perform an alanine scan for this C-terminal residue, but as the corresponding amino acid to F354 is a leucine in Ga s and the amount of direct interactions to surrounding amino acids suggest a great importance for this residue, cation-p-interactions seem to constitute an additional determinant of coupling selectivity. Comparable to residue F139 from b2AR, M140 of D2R is stabilized by a hydrophobic pocket comprised of different amino acids within the two D2R-Ga i -simulations (K192, L194, F336 and T340 in D2 Down R-Ga i and R32, V34, L193 and I343 in D2 Up R-Ga i , Figure 6g, h). These differences are likely to be caused by the movement of Ga i within the D2 Up R-Ga i -simulation ( Figure S6). A significant difference between the D2 Down R-Ga i and D2 Up R-Ga icomplexes lies in the conformation of residue R132 from TM3 (Figure 6a, d). Whereas the side chain of R132 points ''downwards'' in the direction of the C-terminal a5-helix of Ga i in the D2 Up R-Ga i -complex, its side chain reaches out directly towards the junction of TM7/H8 in the D2 Down R-Ga i -complex. R132 forms a salt bridge to residue D350 from the C-terminus of Ga i in D2 Up R-Ga i . In contrast, R132 and D350 do not show direct D2 Down R-Ga i -interactions. Thus, the conformation of R132 is stabilized by residue F429 from H8 of D2R and D350 of Ga i forms a hydrogen bond to residue N430 of D2R.

Selectivity Determinants
Selectivity of a GPCR for a distinct G-protein (or vice versa) arises from structural differences at the interacting epitopes. Figure 5 provides a direct comparison between residues of the b2AR-Ga s -and the D2R-Ga i -complexes that participate in stabilizing the receptor-G-protein interfaces while showing sequence differences at the same time. Highlighted amino acids of b2AR and D2R are suggested to be crucial for coupling to Ga sand Ga i -proteins, respectively, as they exhibit a high degree of sequence conservation within the subfamily of aminergic Ga s and Ga i coupled receptors, which is depicted in Figure 7.
Significant amino acids that control selective receptor-G-protein coupling are located mainly at the intracellular end of TM5 and the N-terminal region of IL3, which comprise a coupling domain for the C-terminal part of Ga and the a4/b6 domain ( Figure 5). Interactions of b2AR with the C-terminus of Ga s are supported by residues from TM3-IL2, TM6 and TM7-H8. Among these residues, I135, A226, Q229, I233, E237, T274 and I278 represent a strongly conserved feature of aminergic GPCRs that couple preferentially to Ga s (Figure 7). The equivalent of A226 in TM5 is represented by an alanine residue for every Ga s -coupling amine receptor, but differs within the Ga i -coupled subfamily. The Cterminal parts of Ga differ significantly (Figure 4, 6). Residues Y377, E378 and L380 as well as D350, C351, G352 and F354 in Ga s and Ga i , respectively, are differently stabilized within their GPCR-pockets and lead to a different orientation of their Ctermini ( Figure 8). Together with residues from the lower parts of the a5-helix (Q370, R371, H373 and R375 in Ga s and I344 and N347 in Ga i ), which interact with amino acids from the Nterminal part of IL3 of the receptors, they constitute, in general, the main determinant of coupling selectivity of G-proteins. The importance of these regions is supported by mutational studies [37][38][39][40]. In agreement with functional experiments with artificial model proteins indicating the importance of the N-terminal part of IL3 for D2R coupling [41,42], the selectivity-determining areas of the D2R-Ga i -complexes were found to be located in the intracellular TM5/N-terminal IL3-region of D2R and the Cterminal part of Ga i . Selective coupling is supported by the junction of TM3 and IL2, the C-terminal TM6 and the junction of TM7 and helix 8 ( Figure 5) when the major amino acids of GPCRs that couple mainly to Ga i were shown to be a valine residue (V136 in D2R) in TM3 and two residues from TM7/H8, F429 and N430 in D2R (Figure 7). The most striking difference between the b2AR-Ga s -complex and the D2R-Ga i -complexes was identified for the interaction of the GPCRs' intracellular loop 2 and the domains aN/b1 and a4/ b6 of the G-proteins. Thus, the intracellular loop 2 of b2AR presents a phenylalanine (F139) interacting with a hydrophobic pocket formed by residues from aN-b1, b2-b3 and a5 of Ga s (Figure 6f). Especially the aromatic amino acids H41 and F205 from Ga s are suggested to enable a highly efficient stabilization of the aromatic residue F139 or, to a lesser extent, other bulky, hydrophobic residues in the equivalent position of IL2 in other GPCRs ( Figure 5, 7). In contrast to the hydrophobic interaction of b2AR and Ga s , D2R and Ga i form ionic interactions between basic amino acids of D2R and negatively charged amino acids of Ga i . Ionic contacts involve arginine residues of IL2/TM3 (R145, R150) and TM5/IL3 (R219, K226, R227) of D2R and glutamate residues of aN-b1 (E25, E28) and a4-b6 (E308, E318) of Ga i . The importance of these basic amino acids in D2R is proposed to be a general determinant of coupling selectivity towards G i , as the structures of G i preferring aminergic GPCRs exhibit homologous residues in the corresponding positions (Figure 7).

Conclusions
To evaluate receptor binding and activation of unexplored GPCR subtypes and to understand the variety of functionally relevant conformations better, the recent crystal structure of the ternary b2AR-Ga s -complex must be complemented by dynamic techniques such as molecular-dynamics simulations, NMR or mass spectroscopy. Because many GPCRs are able to bind more than one G-protein-subtype, models of individual receptor-G-protein complexes will facilitate the rational design of functionally selective drugs inducing the activation of a particular G-protein to a higher extend than coupling to alternative G-protein subtypes. Activation of multiple G-protein dependent and independent pathways and the existence of functionally biased ligands have been demon-strated for the pharmacologically relevant D2R [43][44][45]. The different coupling characteristics of the Ga-subunits Ga i and Ga o towards D2R are associated with subtle sequence differences within their GPCR binding interfaces that involve ionic interactions in Ga i (E28, D315 and D350) missing in Ga o ( Figure S9).
Examination of functionally biased ligands in previous studies attributed a major significance to His393 6.55 in TM6 [21,23], whose distinct rotamer conformations were herein shown to stabilize different conformations of the ligand-binding pocket of D2R ( Figure 3). Thus, His393 6.55 can act as a switch that connects the behavior of ligands to distinct conformational ensembles on the intracellular side of the receptor. Further MD-simulations with selected ligands and/or G-proteins are therefore necessary to elucidate the impact of His393 6.55 on functionally selective signaling on a molecular level.
We exploited the crystal structure of the ternary b2AR-Ga scomplex to establish an active-state model of the pharmacologically highly relevant dopaminergic D 2 receptor in complex with  Figure 7. Alignment of contacts areas to G-proteins of aminergic GPCRs. Amino acids of the receptors supposed to determine selective coupling between b2AR-Ga s and D2R-Ga i are highlighted in dark-blue and dark-green, respectively. A brighter color, light-blue or light-green, is attributed to amino acids, which show an identical sequence compared to b2AR and D2R, or, in the case of arginine and lysine residues, a similar sequence, whereas a grey color points to sequence differences. Amino acids, which appear in the interface of b2AR-Ga s and D2R-Ga i , but are not supposed to determine selective coupling, are colored in yellow. doi:10.1371/journal.pone.0067244.g007 the G-protein subunit Ga i and the endogenous ligand dopamine. Different computational methods including molecular-dynamics simulations and computational alanine-scanning mutagenesis were used to identify distinct hot-spot residues that determine receptor-G-protein selectivity (Figure 4, 6). Additionally, we transferred our results to closely related aminergic GPCRs and found highly conserved amino acids of receptor subtypes preferentially coupling to Ga s or Ga i (Figure 7). As structural information for most GPCR-G-protein complexes is still missing, the computational approach described here is of general importance for investigating protein-protein interfaces of ternary complexes and understanding the determinants of functionally selective signaling.
Our computational approach provides firm predictions with respect to amino acids determining selectivity between GPCRs and G-proteins that can now be confirmed experimentally. The impact of water molecules and possible entropic contributions to selective receptor-G-protein coupling were neglected. In the near future, increasing computational power may give the modeling community the opportunity to visualize the activation of a GPCR and its binding to a G-protein in ''real time'' and to perform such investigations on a higher level of accuracy. A detailed knowledge of the distinct conformational steps involved in receptor activation upon ligand binding and receptor-G-protein coupling will be a prerequisite on the way to fully reveal the secrets of GPCRsignaling.

Homology Modeling
We used the crystal structure of the b 2 -adrenergic receptor (b2AR) together with a heterotrimeric G-protein [16] (PDB-ID: 3SN6) as a starting point for our calculations. The coordinates of the b2AR and the GaRAS-part of the Ga s -protein were used as a template to create a homology model of the dopaminergic D 2 receptor (D2R) in complex with a Ga i1 -protein. We omitted the bc-subunit because it has been shown that the (acylated) a-subunit is sufficient to interact with a G-protein coupled receptor [46]. Three amino acids in the extracellular loop 2 (EL2) of b2AR that are not resolved in the crystal structure were taken from a nanobody-stabilized active-state structure of the b2AR [15] (PDB-ID: 3P0G), the 16 residues missing in the area around a4 of Ga s RAS were modeled manually according to the structure of the GTPcS-bound Ga s -protein [47] (PDB-ID: 1AZT). The aminoacid sequences for GPCRs and G-proteins were retrieved from the SWISS-PROT database [48]. b2AR and D2R sequences (together with 16 additional sequences of family A GPCRs) as well as Ga s and Ga i1 sequences (together with 4 additional Ga protein sequences) were aligned using ClustalX [49] (Gonnet series matrix with a gap open penalty of 10 and a gap extension penalty of 0.2). The initial sequence alignment was manually refined where necessary by means of BioEdit [50] in order to achieve a perfect alignment of the highly conserved amino acids. Absent parts of the b2AR-Ga s -complex structure (i.e. intracellular loop 3 of b2AR and Ga s AH of Ga s ) were omitted in the alignment. It has been shown experimentally that removing the bulk of IL3 within the b2AR does not prevent the receptor from coupling to its G-protein [51]. In addition, constructs of the muscarinic receptors M2 and M3, in which the central region of IL3 (more than 100 amino acids) was omitted, were still able to bind their G-proteins selectively and with near wild type efficacy [30,38]. Therefore, we assume that the truncated D2R used in our investigations is still able to bind to the G i -protein selectively, especially as the important N-and C-terminal portions of IL3 are present.
Based on the final alignment and the b2AR-Ga s RAS-complex structure as a template, we created 50 models of the D2R-Ga i RAS-complex using MODELLER 9v4 [52]. We observed two different rotamer conformations of residue His393 6.55 in the D2R models with the side chain of His393 6.55 pointing to the extracellular and intracellular part of the receptor, respectively. We selected two models of the D2R-Ga i RAS-complex (referred to as D2 Up R and D2 Down R) for further investigation. The models showed the canonical disulfide bond between residue Cys107 3.25 of transmembrane helix 3 (TM3) and residue Cys182 of extracellular loop 2 (EL2). A second disulfide bond between residues Cys399 and Cys401 of EL3 was attributed to the models because of the spatial proximity of the cysteine residues involved and the observation that the highly homologous dopaminergic D 3 receptor exhibits a second disulfide bond in an equivalent position [53].

Structure Refinement and Modification
The two D2R-Ga i RAS-complexes were submitted to energy minimization in order to remove bad van der Waals contacts of the amino-acid side chains. The SANDER classic module of the AMBER10 program package was used [54]. We applied 500 steps of steepest descent minimization, followed by 4,500 steps of conjugate gradient minimization. The minimization steps were carried out in a water box with periodic boundary conditions and a nonbonded cutoff of 10.0 Å . The all-atom force field ff99SB [55] was used.
In order to avoid unnecessarily high flexibility during the simulation process caused by open ends in the Ga part of the complexes, we completed the structure of Ga by modeling the missing helical part of Ga i (Ga i AH) manually according to the crystal structure of a GDP-bound heterotrimeric Ga i1 b 1 c 2 protein [56] (PDB-ID: 1GP2) and submitted both complexes to energy minimization (see procedure described above). Dopamine was manually docked into D2 Up R-Ga i and D2 Down R-Ga i RAS to obtain agonist-bound ternary GPCR-G-protein systems. The two nucleotide-free ternary D2R-complexes were minimized with SANDER according to the procedure described above using the general AMBER force field (GAFF) [57] for the dopamine atoms and ff99SB for protein residues. Parameters for dopamine were assigned using antechamber [54] and charges were calculated using Gaussian 09 [58] at the HF/6-31(d,p) level and the RESP procedure according to the literature [59]. A formal charge of +1 was defined for dopamine.
The structural information for the majority of the missing Ga s AH in the b2AR-Ga s RAS-complex was taken from the crystal structure of the GTPcS-bound Ga s -protein (PDB-ID: 1AZT). A small loop of Ga s that connects the Ga s RAS-and the Ga s AHsubunits, the a1/aA-loop, still not resolved, was modeled manually according to the crystal structure of 1GP2 (residues I55 to K70). Non-conserved residues between Ga s and Ga i were mutated by means of PyMOL [60]. The final structure, comprised of the agonist BI167107, b2AR and the nucleotide-free Ga s , was submitted to energy minimization using the procedure described above for the D2R-Ga i RAS-systems. Parameters and charges for the ligand BI167107 were used as described above and a formal charge of +1 was attributed to BI167107.

Preparation of the Simulation Systems
Parameter topology and coordinate files for the minimized complexes (BI167107-b2AR-Ga s , dopamine-D2 Up R-Ga i and dopamine-D2 Down R-Ga i ) were build up using the tleap module of AMBER10 and subsequently converted into GROMACS input files [61,62].
Each complex was inserted into a dioleoylphosphatidylcholine (DOPC) membrane according to a procedure applied successfully earlier [9].
A pre-equilibrated system bearing a hydrated membrane with 72 DOPC lipids [63] was used. This system had to be enlarged in the x, y and z dimensions in order to surround the ternary complexes fully using a method described earlier [9]. The resulting membrane contained 460 DOPC lipids. According to the density profiles of the membrane, the distribution of all components was confirmed to be as expected without water invading the lipophilic parts of the membrane ( Figure S3).
The charges of the simulation systems were neutralized by adding 3 sodium and 8 chloride atoms to the b2AR and the D2R complexes, respectively. In total, the BI167107-b2AR-Ga s system consisted of 223,264 atoms (659 amino acids, 49,661 water molecules), the dopamine-D2 Up R-Ga i system of 227,641 atoms (624 amino acids, 51,333 water molecules) and the dopamine-D2 Down R-Ga i system of 224,760 atoms (624 amino acids, 50,188 water molecules).

Membrane Simulations
For all simulations, GAFF was used for the ligands and the DOPC molecules and the force field ff99SB for the protein residues. The SPC/E water model [64] was applied.
After insertion into the prepared membrane, the simulation systems were submitted to energy minimization, equilibration (100 ns) and production molecular-dynamics simulation runs (1 ms) at 310 K using the GROMACS simulation package [61,62] as described earlier [9]. Initial gaps between GPCRs and DOPC-lipids were shown close perfectly during the equilibration ( Figure S2).
Throughout the productive simulations a force of 1.0 kcal mol 21 Å 22 was applied to the N-terminal part of the G-protein's aN-helix. In vivo, the aN-helix is anchored to the membrane via acylation with fatty acids and further stabilized by the bc-subunit when the G-protein is nucleotide-free or bound to GDP [46,65]. The aim of the applied force is to avoid spurious conformations caused by the high flexibility of the aN-helix in the absence of both the bc-subunit and the stabilizing acylations because the amino acids that could potentially be acylated are not resolved in the crystal structure of the ternary complex (PDB-ID: 3SN6).

Data Analysis
The analysis of the trajectories was performed with the PTRAJ module of AMBER10. Calculation of the binding free energies and computational alanine scanning mutagenesis was accomplished using the script MMPBSA.py [35]. As our simulations systems are very large, water molecules had to be deleted from the trajectories before analyzing the data in order to reduce the computational demand of the calculations. Therefore, we cannot preclude the existence of further interactions between GPCRs and G-proteins mediated by water molecules. At least for the interactions revealed by our contact analysis, the interacting amino acids are close enough to each other to form stable interactions, even without water molecules.
Figures were prepared using PyMOL [60].  Figure S8 Free energies of binding for the ternary complexes. The free energies of binding for the b2AR-Ga s system (A), for the D2 Down R-Ga i system (B) and for the D2 Up R-Ga i system (C) are shown. Here, the free energy of binding consists of a molecular mechanics energy term (internal energy of bonds, angles and dihedrals), the polar contribution and the nonpolar contribution of the solvation free energy (polar contribution calculated using the Generalized Born equation and the nonpolar contribution using the molecular solvent-accessible surface area). The curves exhibit a best fit line with a positive gradient for (A) and (B) (0.012 and 0.021 for the b2AR-Ga s -and the D2 Down R-Ga i -system, respectively), and a negative gradient for curve (C) (20.021 for the D2 Up R-Ga i -system). As these gradients are very small, we expect that the values will converge to zero for longer simulation times. (TIFF) Figure S9 Alignment of contact areas of chosen Gasubunits. Amino acids within the Ga s and Ga i sequences forming stable contacts to receptor residues are highlighted with a blue and green background, respectively (according to Figure 4). Red backgrounds point to sequence differences between Ga i and Ga o subunits. Red letters indicate residues involved in ionic interactions.

Supporting Information
(TIFF) Table S1 Amino-acid contacts within the b2AR-Ga ssimulation. The occurrence for each amino-acid contact throughout the MD simulation is shown in the grey columns. (DOC)