Abbreviations

1. Introduction

G protein-coupled receptors (GPCRs) are the largest and most diverse class of membrane receptors in eukaryotes[1]; they bind many different types of ligands to initiate an array of intracellular signaling cascades. GPCRs primarily interact with membrane associated, heterotrimeric complexes called G proteins in order to transduce their extracellular signal into a cellular response. The three subunits, Gα, β, and γ, undergo conformational changes to interact with different protein binding partners along their signaling cycle in order to transmit the appropriate messages within the cell [1,2].

The most dynamic changes in structure and affinity can be seen in the Gα subunit, which mitigates each step of the complex’s signaling dynamics and function [3,4,5]. The affinity of the Gα subunit to each of its different binding partners is determined by the structural changes it undergoes within the signaling cycle [6]. Therefore one can think of the Gα subunit as the control center of this signal transducing machinery as it preferentially interacts with different proteins, complexes, and small molecules via conformational changes of its own structure to propagate the information to other signaling moieties within the cell (Figure 1A).

Figure 1. G protein coupled receptors (GPCRs) typically signal through interaction with membrane-associated heterotrimeric G proteins. These proteins become activated via the exchange of GDP for GTP, induced by the activated receptor (R*). Upon this nucleotide exchange, the heterotrimer dissociates into the monomer Gα(GTP) and Gβγ dimer which may then interact with downstream signal effector proteins (not shown for clarity) to propagate and amplify intracellular signaling. The cycle is complete when Gα hydrolyzes GTP to GDP + Pi which allows the trimer to reassemble into the basal, non-signaling state. A) Linear schematic of the G protein signaling cycle. B) ROSETTA-derived structural representations of the three Gα states examined herein; Gα_i1(GDP)β₁γ₁, R*-Gα_i1(empty)β₁γ₁, and Gα_i1(GTP).

DownLoad: Full-Size ImgPowerPoint

1.1. G protein signaling cycle

In its inactive state, the Gα subunit has a high affinity for the nucleotide GDP, possesses a closed helical domain, and interacts with the Gβγ subunits. Upon interaction with an activated GPCR, the Gα subunit undergoes conformational changes to accommodate binding the receptor (Figure 1B) [5,7,8]. This includes the rigid body rotation of its α5 helix up and into the receptor (Figure 2), as it moves along the hydrophobic β-sheets surrounding it to create new interactions sites within the GTPase domain and to the helical domain [9]. This rotation signals to the rest of the complex through an altered interaction network that the Gα subunit is bound to the receptor. In this receptor-bound conformation, the Gα subunit’s affinity for GDP is drastically reduced as its flexible helical domain opens to allow nucleotide escape [10]. Upon GDP release, the Gα subunit has a high affinity for GTP, though the nucleotides only differ in the addition of a single phosphate group. Once GTP is present in the binding pocket, the Gα subunit once again alters its conformation and affinity for both the activated receptor and the Gβγ subunits bound to it. Subsequent dissociation of the Gα subunit from this complex frees Gα, as well as Gβγ, to interact with downstream signaling effector proteins and regulator molecules in order to continue the signaling cascade [11]. In this GTP-bound, active conformation, the Gα subunit possesses different binding interfaces to interact with various effector moieties [3,12]. The intrinsic enzymatic ability of Gα hydrolyzes GTP back to GDP [2]. The rate of hydrolysis can be altered by interacting with accessory proteins which alter the enzyme’s catalytic efficiency [13]. Upon cleavage of the γ phosphate group, the Gα subunit structure returns to its basal state where its propensity to complex with Gβγ is once again higher than its affinity to interact with other signaling moieties; the reunion of the heterotrimer allows the signaling cycle to terminate or for the complex to begin additional rounds of signaling [1].

Figure 2. A) Representative Gα_i1(GDP) structure in the basal state (β₁γ₁removed for clarity) and rotated 180°. Basal nucleotide (GDP) is depicted in dark blue. B) The sequence is derived from rat Gα_i1. Secondary structure elements (red—helices, blue—sheets, and green—critical loop regions) are labeled as described in [3].

DownLoad: Full-Size ImgPowerPoint

1.2. Significance

Current progress in crystallization of GPCRs has greatly aided in our understanding of the G protein’s role within the ternary complex model. Recent work from the Kobilka laboratory has provided the first glimpse of an activated GPCR, the β2-Adrenergic receptor, in complex with a Gα_s heterotrimeric G protein [5]. However, the experimental structure does not provide information on the energetic interactions between amino acids critical for the signaling process. What are the energetic contributions of interactions, broken and newly made, that move the signal from the receptor to the nucleotide binding site? Such an analysis is complicated as the experimental structure presents a static image of interactions in a dynamic system. Crystal structures alone cannot show the conformational dynamics the Gα subunit must continue to undergo to propagate information to the rest of the complex. Further, the use of nanobodies, mutations, and various crystallization aids can alter physiologically relevant conformations of the protein to achieve the most energetically stable interactions for crystal formation.

To better understand the modulatory process the Gα subunit undergoes to propagate its signaling information, an energetic analysis of these conformational changes was performed. We introduce a new pairwise, residue-residue assessment of protein side chain and backbone interactions to describe tertiary topology. Using the available crystallographic structures of each conformation the Gα subunit progresses through during different signaling states, we have created interaction network “maps”. Specifically, we have chosen to investigate the heterotrimeric G protein α subunit in its basal, receptor-unbound Gα_i1(GDP)βγ state, the receptor-bound R*-Gα_i1(empty)βγ state, and the activated, monomeric Gα_i1(GTP) state using the protein software suite, ROSETTA (Figure 1B).

Understanding the mechanism of cellular signaling is a crucial step in understanding the biology of any living organism. This article analyzes changes in conformational and structural information by evaluating the predicted energy of interactions required to maintain function of the Gα subunit before, during, and after binding with the membrane-bound receptor. The ROSETTA protein modeling software allows interrogation of intra-protein and inter-protein interactions on the amino acid level. Using an established comparative modeling protocol [14,15] and binding interface analysis [16,17], we have created the first comprehensive framework for interrogation of pairwise amino acid interactions across each of the signaling states. This analysis has allowed us to create predictive communication maps between interacting side chain pairs throughout the Gα structure as the conformational shifts propagate.

2. Materials and Methods

2.1. Models

To create interaction networks within the different signaling states of the G protein α subunit, we have combined several methodologies. Using previously published comparative models of the GPCR-Gα_i1 heterotrimeric proteins [9], we have created an ensemble of structures for both the basal Gα_i1(GDP)βγ and the receptor-bound R*-Gα_i1(empty)βγ states. Likewise, we have utilized the available crystal structures of activated, monomeric Gα_i1 for a similar analysis (PDBIDs: 1GIA, 1GIL). Each structure of activated Gα_i1 was energy-minimized in the presence of its GTP-analogue. To ensure a robust sampling of the backbone and side chain conformational space consistent with low energy, 500 models were created based on a Gαi crystal structure (PDBID 1GIA). As more extensive sampling with 1000 poses was not shown to greatly increase model quality, generation of 500 models was used for all other structures. This is consistent with the findings of previous protocols [9]. Of these models, the ten lowest scoring models by ROSETTA score were shown to cover the spread of structural flexibility without allowing for larger structural deviations (Supplemental Figure 1). These ten models were employed for further analysis. For all analyses herein, each model possessed the appropriate nucleotide for the given signaling state during all calculations.

2.2. ΔΔG

From these initial models we then probed for intra- and intermolecular interaction energies using the ROSETTA computer modeling software suite. Three signaling states of the G protein α subunit were addressed: Gα_i1(GDP)βγ, R*-Gα_i1(empty)βγ, and Gα_i1(GTP) (Figure 1B). For each state, the binding interface energy (ΔΔG) was calculated for various key inter-protein interfaces across the complex and within the GTPase and the helical domains of the Gα subunit. Regions for analysis were selected for their roles as protein-protein interfaces or for their apparent role in maintaining protein stability within each conformational state. Specifically, key secondary structure elements (Figure 2) were evaluated for their ability to contribute to overall protein stability by calculating the changes in free energy before and after removal from the structure. Note that all energies are given in ROSETTA Energy Units (REUs) and include predicted contributions of van der Waals interactions, desolvation effects, hydrogen bonds, and electrostatics. While the ROSETTA-predicted energy has been shown to correlate with the free energy in kcal/mol [16,17], it is important to highlight that inaccuracies in the structural models and simplifications in the ROSETTA energy function lead to deviations between predicted and experimentally observed energies. Furthermore, the internal energy of small molecules is assumed to be unaltered upon binding to the protein; the energy measurements herein reflect energy perturbations induced by the ligand when binding to the protein. All ΔΔG results are reported as the absolute value of REU scores for consistency with previously published data [9].

2.3. Pairwise interaction score analysis

Each of the three signaling states of the Gα_i1 subunit were then interrogated at the amino acid level utilizing ROSETTA’s pairwise score breakdown assessment. This feature calculates the interaction score for each possible amino acid pair. Note, that while this score is also measured in ROSETTA Energy Units (REUs) it is not a free energy in the thermodynamic sense. We therefore call these values consistently ‘interaction scores’. However, this analysis allows for intra-molecular probing of information flow across signaling states while creating a network of stabilizing amino acid interactions. A protocol capture for this application has been validated externally and is available for public use within the ROSETTA framework. Herein, we apply this method to the G protein α_i1 subunit to highlight the method’s effectiveness in predicting relevant structural nuances. Each of the signaling states of Gα_i1 were assessed by averaging the per-residue contribution of the top ten lowest scoring models. The appropriate nucleotides and subunits were present throughout all calculations.

2.4. Pairwise interaction score calculation

Pairwise interaction scores were calculated using the ROSETTA software suite. The per residue score breakdown was calculated on ten comparative models which were created as previously described [9].

/residue_energy_breakdown.linuxgccrelease -database /rosetta/main/database/ -in:files:s <list individual pdbs> -output:prefix <output file name> -restore_pre_talaris_2013_behavior

2.5. Protocols for pairwise interaction score analysis

Average per residue interaction pairs were calculated across ten models per signaling state in MATLAB using the following script:

file_1 = 'model_1.xlsm';

[~, ~, raw_1] = xlsread(file_1);

model_1 = zeros(5223, 3);

model_1(:, 1) = cell2mat(raw_1(1:end, 3));

model_1(:, 2) = cell2mat(raw_1(1:end, 5));

model_1(:, 3) = cell2mat(raw_1(1:end, 26));

new_matrix_1 = nan(354, 354);

for i = 1:size(model_1)

 new_matrix_1(model_1(i, 1), model_1(i, 2))= model_1(i, 3);

end

Continued for all models analyzed, then average scores across all models:

ave_matrix = nan(354, 354);

for ii = 1:354

 for jj = 1:354

  ave_matrix(ii, jj) = mean([new_matrix_1(ii, jj) new_matrix_2(ii, jj) new_matrix_3(ii, jj) etc.]);

 end

end

g = ave_matrix(~isnan(ave_matrix));

[i, j] = ind2sub(size(ave_matrix), find(~isnan(ave_matrix)));

fin = [i, j, g];

2.6. Protocol capture

For further breakdown of all computational methods utilized herein, please refer to the companion Protocol Capture. All in silico methods and calculations were graciously verified by an external reviewer, Dr. J. Koehler Leman, Chemical and Biomolecular Engineering, Johns Hopkins University, Baltimore, MD.

3. Results

3.1. Estimating free energy changes across protein-protein interfaces

Predicting free energy changes across protein-protein interfaces has been a staple in understanding the dynamics and kinetics of protein-protein interaction [16,18,19]. Used as a measure of binding efficiency, ΔΔG estimates are a useful means of probing the thermodynamic stability of a protein interface in the bound and unbound states [9,19,20]. For our purposes, we utilized this measure to assess the energetic contribution secondary structure elements possessed along intra-protein interfaces between the helical and GTPase domains as well as for inter-domain stability.

For our calculations, specific secondary structure elements (Figure 2) were evaluated for their ability to contribute to overall protein stability by calculating the changes in free energy before and after their removal from the subunit structure. For all calculations, the appropriate nucleotides were present. The top ten lowest scoring models for the Gα_i1(GDP)βγ, R*-Gα_i1(empty)βγ, and Gα_i1(GTP) states were each assessed, and their ROSETTA scores were averaged.

3.1.1. GDP vs. GTP-bound models

The Gαsubunit possesses similar energy in both the basal, Gα_i1(GDP)βγ, and activated, Gα_i1(GTP), states. This is expected as the two states differ only in the addition of a γ-phosphate ion. Though the Gβγ subunits were present for the basal calculations of the trimer, they do not significantly alter Gα’s energetics when evaluating regions such as the α1, α5, or αF helices (Figures 3-4,and Supplemental Tables 1-3). When evaluating these regions, the resulting energies highlight a consistency between these two states suggesting that any structural changes within these regions begin and end with similar energies of interaction.

Figure 3. Structural representation of predicted ΔΔG of the α1 helix across three states of Gα signaling—Gα_i1(GDP)β₁γ₁, R*-Gα_i1(empty)β₁γ₁, and Gα_i1(GTP). The α1 helix was defined as residues G45-E58 based on rat Gαi1 sequence. All calculations were averaged across the top ten scoring models for each state. Values reported here represent the absolute values of Rosetta Energy Units (REUs). REUs above 0.5 are considered significant. Residue contributions to the interface are color coded to indicate a greater contribution to stability. Lighter blue indicate a lower REU value relative to the darker shades.

DownLoad: Full-Size ImgPowerPoint

Figure 4. Structural representation of predicted ΔΔG of the α5 helix across three states of Gα signaling—Gα_i1(GDP)β₁γ₁, R*-Gα_i1(empty)β₁γ₁, and Gα_i1(GTP). The α5 helix was extended to include part of the TCAT motif and is defined as residues C325-F354 based on rat Gαi1 sequence. All calculations were averaged across the top ten scoring models for each state. Values reported here represent the absolute values of ROSETTA Energy Units (REUs). REUs above 0.5 are considered significant. Residue contributions to the interface are color coded to indicate a greater contribution to stability. Lighter blue indicate a lower REU value relative to the darker shades. *Note: Gα_i1(GTP) crystal structure only extends to residue I343 preventing analysis of the α5 helix beyond this residue in the activated, monomeric state.

DownLoad: Full-Size ImgPowerPoint

Noteworthy alterations in energies are seen around the nucleotide binding pocket and residues involved in stabilizing the Gβ interface between the basal and activated states. Examination of the P-loop and the variable Switch regions (I-III) (Figure 2) indicate more subtle ΔΔG changes across these regions (Supplemental Tables 4-7). In the basal, trimeric state, the Gβ subunit organizes the loop regions into a binding interface. In its absence, the activated monomeric models do not show significant changes as seen in ROSETTA energy scores overall, though specific amino acid positions are reported to modulate.

3.1.2. Receptor-induced conformational changes

In contrast to the basal and activated states, the R*-Gα_i1(empty)βγ models show a stark transition in the communication network across the secondary structure elements of the Gα subunit. During this phase of the signaling cycle, the Gα subunit undergoes a large conformational change which can be seen in the shifting of energetics around the α1 helix, the α5 helix (Figures 3-4,Supplemental Tables 1-2), and regions involved with nucleotide stability, namely the P-loop and the variable Switch (I-III) regions (Supplemental Table 4-7). It is during this stage of the signaling cycle that the receptor induces activation, the helical domain is opened, and the guanine nucleotide is allowed to exchange. The results from our models are consistent with experimental studies of these structural changes [9,10,21].

3.2. Predicting pairwise residue-residue contributions to protein stability

To interrogate the conformational changes that must occur at the amino acid network level between the signaling states, we devised a new application for the Rosetta modeling software’s per-residue assessment of predicted interactions (publically available); this application allowed us to evaluate individual amino acid contributions to stability and function. For each of the three signaling states, the top ten models were assessed for each amino acid pair contribution to stability. The average score across the ten models was then plotted for each state (Figures 5-7,Supplemental Figures 2-8).

Figure 5. Pairwise analysis of ROSETTA scores for individual amino acid interactions as a means to monitor information networks across signaling states. These interaction networks represent side chain and backbone atom contributions to the stability and functionality of the protein structure. These matrices compare across the protein signaling states to investigate predicted interaction (and therefore structural) changes between residue pairs. The x- and y-axes represent each residue position of the Gαsubunit compared across all other possible residue positions. Above the diagonal depicts the score difference (in REU) between the basal, heterotrimeric Gα_i1(GDP)β₁γ₁structure and the receptor-bound complex, R*Gα_i1(empty)β₁γ₁. The lower matrix below the diagonal depicts the score difference (in REU) between the basal, heterotrimeric Gα_i1(GDP)β₁γ₁structure and the activated Gα_i1(GTP) monomeric subunit after dissociation from Rhodopsin and β₁γ₁. Only the Gαsubunit’s amino acid contributions are shown for clarity. Residue-residue interactions in the range of −0.5 to 0.5 REU were removed to highlight more significant differences in contributions. Stabilizing residue interactions are depicted in red while a predicted loss of interaction scores are shown in blue. Note* The crystal structures used for the monomeric Gα_i1(GTP) models lack residues 1-34, and 343-354. These residues are therefore removed from the analysis.

DownLoad: Full-Size ImgPowerPoint

Figure 6. Pairwise analysis of individual amino acid interactions across three Gαsignaling states to investigate changes between residue pairs interacting with the α1 helix and Linker 1 region (K46-I56, & H57-S62, respectfully). A) The heterotrimeric Gα_i1(GDP)β₁γ₁structure in the basal state. B) The heterotrimeric Gα_i1(empty)β₁γ₁subunits interacting with the GPCR, Rhodopsin. C) The activated Gα_i1(GTP) monomeric subunit after dissociation from Rhodopsin and β₁γ₁. Only the Gαsubunit’s pairwise amino acid contributions are shown for A-C. REUs for individual pairs are color coded ranging from more stable predicted interaction scores (minimum −2.32) in blue to positive, repulsive scores terms (maximum 0.26) in red.

DownLoad: Full-Size ImgPowerPoint

Figure 7. Pairwise analysis of individual amino acid interactions across three Gαsignaling states to investigate changes between residue pairs interacting with the α5 helix (N331-I343). A) The heterotrimeric Gα_i1(GDP)β₁γ₁structure in the basal state. B) The heterotrimeric Gα_i1(empty)β₁γ₁subunits interacting with the GPCR, Rhodopsin. C) The activated Gα_i1(GTP) monomeric subunit after dissociation from Rhodopsin and β₁γ₁. Only the Gαsubunit’s amino acid contributions are shown for A-C. REUs for individual pairs are color coded ranging from more stable predicted interaction scores (minimum −2.30) in blue to positive, repulsive score terms (maximum 0.31) in red. *Note: Gα_i1(GTP) crystal structure lacks residues 344-354 preventing analysis of the full carboxyl terminus.

DownLoad: Full-Size ImgPowerPoint

To evaluate which interactions were made and broken between the different signaling states, we compared the basal, heterotrimeric Gα_i1(GDP)β₁γ₁ model scores to the receptor-bound, R* Gα_i1(empty)β₁γ₁ and to the monomeric, Gα_i1(GTP) active state (Figure 5). From this calculation we show the variability of the Switch regions, as interactions are lost, or are diminished during receptor binding (red, above the diagonal) and remade in the active state (blue, below the diagonal). Some of this variability may be due to the loss of β₁γ₁ binding upon activation.

The predicted opened and closed conformations of the helical domain are also recognized when evaluating residue-residue interactions across states. The receptor binding induces structural rearrangements that ultimately lead to helical domain opening [5,9,10,21]; therefore, the upper matrix, above the diagonal, indicates the helical domain must break contacts for activation (red). The basal and active Gα_i1 subunits possess very similar secondary structure and tertiary fold. Therefore, fewer interactions are lost or diminished between the two states (below the diagonal). New interactions or more favorable interactions (blue) must be made to accommodate the GTP nucleotide and the lack of β₁γ₁ subunits.

However, this broad representation does not do justice to the nuanced alterations of residue-residue interaction. In addition the overall number of intraprotein Gα interactions is not expected to change across the different signaling states as all secondary structure elements and the global tertiary fold is maintained. Though there are technically fewer intraprotein contacts when the helical domain is opened during the receptor bound state, these differences are subtracted from interaction scores that are present in the GDP-bound trimer; the result is a change in magnitude from interaction to no interaction that is recorded in this matrix. Additionally changes in magnitude for the pairwise interactions in the range of −0.5 to 0.5 REU were removed to highlight more significant contributions.

3.2.1. The Switch Regions

Across the three signaling states, a pattern of interaction emerges. As seen with the ΔΔG calculations (Figures 3-4, Supplemental Tables 1-8), the basal and activated Gα subunits maintain similar amino acid interactions. However, the γ phosphate group present in the activated Gα monomer leads to shifts in the communication networks of the Switch I-III regions (Supplemental Figures 2-4). As implied by their name, these regions have been shown to alter their conformation in the presence of GTP instead of GDP in crystal structures [22,23].

The largest shift among these three elements is shown in the Switch II region (Supplemental Figure 3). This is expected as the Switch II region also interfaces with the Gβ subunit. Present in this analysis during the trimeric basal and receptor-bound states, the pairwise interaction map of Gβ with the Switch II region maintains interaction similarity and therefore structural similarity between these two states. Key differences can be seen around the Switch I and β2 elements with which the Switch II region interacts in the GDP-bound state, but not during the receptor bound state.

Alterations in the communication network along the Gβ subunit itself were not the primary focus of the current study; however, the Gβ/Gγ subunits were present for the analysis of the basal and receptor-bound states. Therefore, they are included as interaction partners along the corresponding interface residues. Interestingly, the Gβ subunit does show an altered conformation network between the receptor-unbound and bound states suggesting some flexibility between the two G protein subunits. This modulation of the Gα Switch II region does not seem to show similar intra-protein interaction flexibility within the Gα subunit itself, but rather it highlights relevant changes within the heterotrimeric complex which may contribute to the mechanisms of receptor induced activation.

3.2.2. Rearrangements from nucleotide exchange

With this detailed analysis, we show that more elements other than the Switch regions possess an altered communication network. Subtle changes in the α1 helix, α5 helix, αF helix and P-loop highlight structural alterations induced by the nucleotide (Figures 6-7, Supplemental Figures 5-6). Specifically, the GDP-bound versus GTP-bound Gα subunit show altered network intensities along the β6-α5 loop within the highly conserved TCAT motif (residues 326-329 in Figure 2B). Changes in the P-loop are also much more dynamic than we had originally predicted (Supplemental Figure 6). Interactions with the residues 147-150 of the αD-αE loop in the basal state are not recovered in the active state. Likewise Switch II and III, and β4 interact with variable degrees of binding intensity (as defined by ROSETTA Energy Units) with the P-loop suggesting more dynamic structural rearrangements in this region.

3.2.3. Receptor-induced network changes

As expected, the receptor-bound heterotrimer possesses an altered interaction network indicative of altered structure. These conformational changes are highlighted in interaction shifts along the α1 helix and the α5 helix (Figures 6-7) as these secondary structure elements move to make transient connections. Connections are also lost between Switch I and the Switch II/β3 interface during receptor binding, which are recovered upon Gα activation and dissociation.

On the backside of the Gα subunit, the P-loop, which has also been implicated in nucleotide stability and possible mechanisms of release [9,24,25] shows a drastic structural rearrangement and transition during receptor binding (Supplemental Figure 6). As observed above, the P-loop possesses a structural alteration resulting in a loss of interaction with the αD-αE loop that is not present during receptor binding nor is it recovered post-dissociation in the monomeric, active state.

The linker 1 region connecting the helical domain to the GTPase domain via the α1 to αA helices also possesses a shift in conformation (Figure 6,Supplemental Figure 7). This element was hypothesized to be an important mechanistic feature to allow domain opening for nucleotide escape [10,24]. However some movement is expected as it does not possess any secondary structure elements.

3.2.4. The Helical Domain as a rigid body

Interesting secondary structural elements within the helical domain, such as the αA helix, do not drastically alter their interaction networks across the three signaling states. This suggests that these elements move together while maintaining a similar tertiary fold (Supplemental Table 8,Supplemental Figure 8). These results are in agreement with DEER, EPR, NMR and crystallographic data [5,10,25], which suggests the helical domain moves as rigid body away from the nucleotide binding pocket [5,10,21].

4. Discussion and Conclusions

The heterotrimeric G protein undergoes dynamic changes in its structure and its binding affinity throughout the stages of the signaling cycle. We utilized structural models of these conformational states to analyze the energetic contributions that stabilize intra- and inter-molecular interactions that define these states, specifically within the Gα_i1 subunit. This new analysis application predicts key amino acids to be nodes within the information network that propagate the signal across the complex upon interaction with the receptor.

Utilizing the ROSETTA software suite, we computed energy values for residue interactions along different binding interfaces. This benchmarked computational technique has been shown to provide useful insight in the following studies [16,17,26]. Likewise, ROSETTA was used to compute pairwise interactions between individual amino acids within the Gα subunit of the heterotrimeric G protein. This technique allowed us to compare the predicted thermodynamically stabilizing interactions between the basal, receptor-bound and activated conformations of the Gα subunit. Through this analysis we were able to detect intra-protein differences in amino acid interaction networks important for propagating conformational changes.

In a previous analysis [9] the Gα_i1 subunit was evaluated in the basal, GDP-bound trimeric state and in the receptor-bound state through the use of ΔΔG analysis. Our current study expands on this progress by also including the activated monomeric state for comparison of energy contributions made by key secondary structure elements to evaluate critical regions for G protein activation. In addition we have modeled all three signaling states to evaluate changes in residue pair contributions during signaling.

4.1. GDP- vs GTP-bound models

From this analysis, we have highlighted the similarity of the Gα_i1 GDP- versus GTP-bound structures. By excising specific structural elements, a broad map of protein stability can be painted. Regions important for interfacing with other proteins, such as the α5 helix and the Switch II domain show the most altered energy changes between these two states (Figure 4,Supplemental Table 2,6). This is expected as the binding partners contribute to the relative energy of the system and inhibit interface flexibility. Regions not involved in protein-protein interactions or large structural rearrangements, such as the αA helix (Supplemental Table 8), remain more energetically stable and consistent across the different models in the GDP and GTP-bound Gα subunit. This result is in agreement with other structural studies that suggest the helical domain moves as a rigid body throughout G protein activation [4,10,21,27].

4.2. Receptor-induced activation

The ΔΔG calculations serve to highlight the role of key secondary structure elements as well as specific non-structured linker regions in G protein activation via receptor coupling, R*-Gα_i1(empty)βγ. During this structural transition state in which Gα must undergo a dynamic conformational change, the ΔΔG analysis shows a shift in interaction partners for the α1 helix, α5 helix, and P-loop (Figures 3-4,Supplemental Tables 1-2,4). This conformation must therefore propagate from the receptor to the helical domain of the Gα subunit in order to disrupt binding of GDP. Each of these elements has been implicated in the mechanism of nucleotide escape and G protein activation [22,24,25,28]. From this analysis alone, however, no direct conclusions could be made on the order or dynamics of conformational propagation across the subunit.

4.3. Residue-Residue changes within the network

To better address this, a more detailed analysis of the structural differences was performed. The informational network mapping through the per-residue pairwise analysis highlighted subtle changes in G protein side chains induced by the γ phosphate group of the nucleotide. These altered interaction scores are indicative of altered structures which may prove to be important for interaction with downstream signaling and regulator moieties. However, we do not predict that all changes seen between these two states contribute to effector selectivity and interaction, as some of the altered network must be involved in maintaining the stability of the new structure without contributing to function.

Our pairwise analysis provides insight into possible routes of this information flow from the receptor to the nucleotide binding pocket. Through examination of the α1 helix, the α5 helix, and the P-loop, extreme displacement of the interaction pairs predicts the importance of these structural elements in allowing nucleotide exchange and G protein activation (Figures 6-7,Supplemental Figure 6). To further test and validate these predictions, additional experiments must be performed to further elucidate the mechanism of G protein activation.

From these analyses, we have created full, downloadable interaction matrices of our results to provide further understanding of G protein structure and modulation (Supplemental Material). By including pairwise score information across several signaling states, we hope this data will prompt new and unique questions on G protein activation and its signaling mechanics through investigation of these interactive communication maps. The values represent averaged relative interaction scores within these protein complexes as derived from comparative modeling based on published crystal structures. Future studies will be required to investigate the true predictive power of these results in vitro.

4.4. Method development

The use of ΔΔG calculations in evaluating protein-protein interfaces has long been an important application within ROSETTA [9,16,17,26]. We utilized this analysis not only for evaluating changes along known protein-protein interfaces, but also along key secondary structure elements thought to be important for propagating conformational changes across the protein subunit or necessary for stability. By mapping the ΔΔG of critical structures across multiple models, we were able to compare relative energy contributions as described by the ROSETTA score function for multiple structural snapshots.

One of the primary purposes for the creation of these energy calculations was to apply and validate a new method of interaction analysis available in the ROSETTA modeling software suite. Here we introduce a new methodology for evaluating the communication networks underlying three dimensional protein topology. By evaluating the residue-residue contributions to protein structure, we have created a technique to map interaction partners necessary for structural stability and conformation transmission. The ROSETTA score term for each contributing residue pair provides a roadmap for amino acid interactions necessary for both structure and function. This pairwise analysis also highlights key nodes of information flow when calculated across multiple protein structural states. The protocol utilized herein has been externally validated and made available for academic and public use with the ROSETTA software suite.

5. Downloadable Communication Maps

From these analyses, we have created downloadable interaction matrices available as supplementary material. They combine secondary structure stability with individual ROSETTA scores of interactions on a residue-residue level. This novel perspective has allowed us to begin to probe regionally specific interactions required for GPCR-G protein interaction, residues required to propagate intra-domain conformational changes, and stabilize the basal, receptor-bound, and activated Gα states. The download also possesses general features about the regional selection such as secondary structure elements, relative evolutionary conservation, amino acid composition etc. as specific to the Gα_i1 subunit sequence.

Acknowledgements

The authors would like to thank Dr. Rocco Moretti for creating the residue_energy_breakdown application within ROSETTA as well as Dr. Ali Kaya for reviewing the manuscript and for helpful discussion.

Funding for this work was provided by the NIH (R01 EY006062, R01 GM080403, R01 GM099842, R01 DK097376, R01 HL122010, R01 GM073151, U19 AI117905) and NSF (CHE 1305874).

Conflicts of Interests

All authors declare no conflicts of interest in this manuscript.

Author Contributions

AL performed the analyses, wrote the manuscript, created models for the signaling states, and wrote some of the scripts needed for analysis. JLK edited the manuscript, wrote some of the scripts needed for analysis, and validated the Protocol Capture. KK wrote some of the scripts needed for analysis. NA created models for the signaling states and wrote some of the scripts needed for analysis. HH edited the manuscript and provided insight and discussion. JM planned the research, edited the manuscript, and provided guidance in protocol development.