Structural Determinants of Constitutive Activation of Gα Proteins

Dec 21, 2016 - Structural Determinants of Constitutive Activation of Gα Proteins: ... Dissecting intrinsic and ligand-induced structural communicatio...
0 downloads 0 Views 9MB Size
Download PDF
Article pubs.acs.org/JCTC

Structural Determinants of Constitutive Activation of Gα Proteins: Transducin as a Paradigm Angelo Felline, Simona Mariani, Francesco Raimondi,† Luca Bellucci,‡ and Francesca Fanelli* Department of Life Sciences, University of Modena and Reggio Emilia, via Campi 103, 41125 Modena, Italy S Supporting Information *

ABSTRACT: Heterotrimeric guanine nucleotide-binding proteins (Gα proteins) are intracellular nanomachines deputed to signal transduction. The switch-on process requires the release of bound GDP from a site at the interface between GTPase and helical domains. Nucleotide release is catalyzed by G protein Coupled Receptors (GPCRs). Here we investigate the functional dynamics of wild type (WT) and six constitutively active mutants (CAMs) of the Gα protein transducin (Gt) by combining atomistic molecular dynamics (MD) simulations with Maxwell-Demod discrete MD (MDdMD) simulations of the receptorcatalyzed transition between GDP-bound and nucleotide-free states. Compared to the WT, Gt CAMs increase the overall fluctuations of nucleotide and its binding site. This is accompanied by weakening of native links involving GDP, α1, the G boxes, β1-β3, and α5. Collectively, constitutive activation by the considered mutants seems to associate with weakening of the interfaces between α5 and the surrounding portions and the interface between GTPase (G) and helical (H) domains. These mutational effects associate with increases in the overall fluctuations of the G and H domains, which reflect on the collective motions of the protein. Gt CAMs, with prominence to G56P, T325A, and F332A, prioritize collective motions of the H domain overlapping with the collective motions associated with receptor-catalyzed nucleotide release. In spite of different local perturbations, the mechanisms of nucleotide exchange catalyzed by activating mutations and by receptor are expected to employ similar molecular switches in the nucleotide binding site and to share the detachment of the H domain from the G domain.

■

mechanism.6 The G domain holds a Rossmann fold, characterized by a 3-layer(αβα) sandwich architecture (Figure 1A and 1B). The nucleotide docks into a binding site contributed by the β1/α1, α1/β2 (αF/β2 in the Gα proteins), β3/α2, β5/α4, and β6/α5 loops (i.e., G boxes 1-5 (G1-G5), Figure 1A and 1B). Sequence conservation over the Ras superfamily resides just in the G boxes, also defined as ultraconserved regions (Figure 1A and 1B). G1, originally termed as the Walker A motif (motif GxxxxGKS/T), is also called P-loop (phosphate-binding loop). The latter, in fact, contacts the phosphates through main-chain NH groups and lysine side chain. G2 is also called switch I (swI or α1/β2 loop (αF/β2 loop in the Gα proteins)), and G3 is part of switch II (swII made by the β3/α2 loop plus the α2-helix).

INTRODUCTION

The Ras superfamily comprises many guanine nucleotidebinding proteins (G proteins) that are essential to intracellular signal transduction.1,2 These proteins act biologically as molecular switches cycling between ON and OFF states, thereby controlling a variety of processes.1 The switch-on process requires the release of the bound Guanosine DiPhosphate (GDP) and the subsequent binding of the Guanosine Tri-Phosphate (GTP), an intrinsically slow process catalyzed by Guanine nucleotide Exchange Factors (GEFs).1 In the inactive GDP-bound state, the G proteins of family α (Gα proteins) form membrane-associated αβγ heterotrimers, forming one of nature’s most important miniature (nano-) machines.3−5 In Ras GTPases nucleotide switch and hydrolysis are played by the conserved core, Ras-like domain or GTPase (G) domain, which represents the basic functional unit of G protein being deputed to GDP/GTP binding and the universal switch © XXXX American Chemical Society

Received: August 18, 2016 Published: December 21, 2016 A

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Primary sequence, crystal structure, and atomic fluctuations of Gt. A. The primary sequence of Gt is shown. Helices, strands, and loops are, respectively, violet, yellow, and white. The G boxes are delimited by green boxes. Orange asterisks indicate the positions of mutations considered in this study. Black numbers on the left side of the sequence refer to the sequential numbering, whereas black numbers above the sequences indicate the beginning of a secondary structure/G-box motif and serve to a position-dependent numbering valid for all members of the Ras-GTPase superfamily.7 Such numbering is characterized by the label of the secondary structure segment followed by the amino acid position in that segment (see the numbering in brackets here below). In those cases where the G-boxes overlap with the secondary structure segment, positions refer to the G-boxes. Loops other than G boxes maintain the sequential numbering. B. The cartoon representation of the Gt-GDP structure (PDB code: 1TAG) is shown in two different views (left and right panels). The Ras-like domain and the helical domain are colored according to secondary structure (i.e., helices, strands, and loops are, respectively, violet, yellow, and white). The GDP nucleotide is represented by sticks colored by atoms type. In the left panel, all the mutation sites are indicated by orange spheres centered on the Cα-atom. In the right panel, the mutated amino acids located in the linkers (G56 and G179(G2:6)) are indicated by an orange sphere centered on the Cα-atom, while mutated amino acids located in G5 (A322(G5:4)) and in α5 (T325(α5:1), V328A(α5:4), and F332(α5:8)) are represented by orange sticks. The residues involved in local perturbations of inter-residue interactions caused by the mutations are shown in sticks colored by atom type (namely, T44(α1:3), I45(α1:4), Q48(α1:7), and I52(α1:11), F185(β2:6), F187(β2:8), and F192(β3:3)). The β2/β3 hairpin in between swI and swII is also called interswitch. C. The Cα-Root Mean Square Fluctuation (RMSF) profiles from Molecular Dynamics (MD) trajectories of WT (cyan) and CAMs Gα (different colors) are shown. They refer to the 500000 frames constituting the 500 ns trajectory. The secondary structure elements are shown on the abscissa; G boxes are indicated in green on the abscissa. In all panels, the secondary structure elements are labeled following Noel’s nomenclature.8 B

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

effects overlap with the GEF action of rhodopsin. The questions were addressed by a number of analyses done on equilibrium MD simulations of WT and mutant Gt coupled with Maxwell-Demod discrete MD (MDdMD) simulations of the conformational transition between receptor-free GDPbound and receptor-bound nucleotide-free forms of the G protein.

A singularity of the Gα proteins is the presence of a third switch, swIII, that is the loop connecting β4 to α38 as well as the presence of an extra-Ras helical (H) domain forming an orthogonal bundle of a long central helix (αA) surrounded by five shorter helices (αB-αF, Figure 1A). The polypeptide chain leaves the G domain and enters the H domain through the α1/ αA loop or linker 1 (L1) and successively leaves the H domain and goes back into the G domain through the αF/β2 loop (swI or L2 (Figure 1A and 1B)). The interface between G and H domains constitutes the nucleotide binding cleft. Upon activation by extracellular signals, G protein-coupled receptors (GPCRs) catalyze the exchange of bound GDP for GTP, thus behaving as specialized GEFs for the Gα proteins. The latter recognize GPCRs by the C-term of α5 and the α/β loops, which are distal from the nucleotide binding site.9 Insights into the mechanism of GPCR-catalyzed GDP release arise from computational and biophysical experiments as well as from the crystallographic structure between agonist-bound β2adrenergic receptor (β2-AR) and nucleotide-free heterotrimeric Gs. These experiments indicate that, by acting on the orientation of α5 and on the conformation of the β6/α5 loop, GPCRs trigger the weakening of contacts between GDP and the G domain accompanied by a detachment of the H domain from the G domain, which ultimately leads to GDP release.9−16 The Gα protein transducin or Gαt (Gt) is the transducer in visual phototransduction, the process that generates a neuronal signal following light capture by visual pigments in photoreceptor cells (rods and cones). The visual pigment of rods is rhodopsin the cornerstone of family A GPCRs. Photoactivation of rhodopsin indeed catalyzes the exchange of GDP for GTP in Gt.17 A number of Gt mutants in the two linker regions (i.e., G56P and G179(G2:6)P; superscript brackets enclose the position-dependent amino acid numbering previously defined (see also the legend to Figure 1)7), in G5 (A322(G5:4)S), and in α5 (V328(α5:4)A, T325(α5:1)A and F332(α5:8)A), the fundamental recognition site for rhodopsin, have been shown to hold up to 165-fold increase in the intrinsic rate of nucleotide exchange.18−20 In detail, G59P and G179(G2:6)P, T325(α5:1)A and F332(α5:8)A, and A322(G5:4)S and V328(α5:4)A increased the exchange rate by >200%, 100−200%, and 50−100%, respectively. The localization of such constitutively activating mutations in G protein portions that likely exert a control on the intrinsic dynamics of the protein makes them useful tools to investigate the mechanism of nucleotide exchange at the atomic level. We recently used approaches like Principal Component Analysis (PCA) and Elastic Network Model-Normal Mode Analysis (ENM-NMA) to identify the important structural flexibilities that enable proteins in the Ras superfamily to switch between their active and inactive states.21,22 These analyses, carried out on the crystal structures of representative members, led to a hypothesis regarding the evolutionary adaptation of structural deformations by the individual members of the superfamily to fulfill their specialized function. Moreover, we employed atomistic Molecular Dynamics (MD) simulations also in combination with mathematical modeling of visual phototransduction to investigate nucleotide effects on the functional dynamics of the G protein7 and to infer the structural determinants of a spontaneous mutation of Gt associated with disease.23 Here we investigated if activating mutations exert common effects on the intrinsic dynamics of Gt and if some of those

■

RESULTS Constitutively Active Mutants Affect the Intrinsic Flexibility of Gt. The constitutively active mutants (CAMs) considered in this study, G56P, G179(G2:6)P, A322(G5:4)S, T325(α5:1)A, V328(α5:4)A, and F332(α5:8)A, are expected to affect the intrinsic flexibility of the protein. As for G56P and G179(G2:6)P, both location and character of the amino acid replacement prelude to conformational effects. Indeed, on one hand, the substitutions occur in either linkers that connect G and H domains; on the other hand, a residue intrinsically prone to enhance backbone’s degrees of freedom is replaced by a residue that instead confers rigidity to the main chain. The other four CAMs are associated with local perturbation in interresidue interactions (Figure 1B). In deep detail: a) serine substitution for A322(G5:4) establishes novel interactions with N265(G4:1); b) alanine substitution for T325(α5:1) weakens or abolishes the interactions found in the WT between the mutation site and both Q48(α1:7) and I52(α1:11); c) similarly, alanine substitution for V328(α5:4) weakens or disrupts the interactions found in the WT between the mutation site and both I45(α1:4) and I52(α1:11); finally, d) alanine substitution for F332(α5:8) disrupts the aromatic cluster formed by the mutation site and F185(β2:6), F187(β2:8), and F192(β3:2) (Figure 1B), and it weakens packing interactions with α1 amino acids such as M50(α1:8) and H54(α1:12). Remarkably, the three activating mutation sites in α5 belong to the so-called “N-terminal transmission module” transforming the structural information sent by receptor bound to the C-term of the helix (i.e., “Cterminal interface module”) into the switch of Gα protein activation.14 The separation of N-terminal and C-terminal portions of α5 in two modules bearing distinct functional roles emerged from a bioinformatics study combining the alignment of 973 Gα sequences with the analysis of noncovalent contacts between amino acid residues.14 Along the same line, the three activating mutation sites in α5 belong to the “activation cluster I” inferred from an alanine scanning mutagenesis study probing Gi1 activation at single-amino acid resolution (Supplementary Table 1 (Table S1).16 Collectively, the considered activating mutations in G5 or α5 perturb the G5 itself as well as packing interactions between α5 and a number of portions that precede and follow the two interdomain linkers. Mutational effects on Gt dynamics were investigated by comparing the MD trajectories. The analyses of the Cα-Root Main Square Fluctuation (Cα-RMSF) profiles highlight the high flexibility of linkers, switches, and interswitch, as an intrinsic feature of the G domain (Figure 1C). T325(α5:1)A holds an overflexible H domain, switches, and interswitch compared to the other mutants. Remarkably, an increased flexibility of the H domain, in particular, the αB/αC loop, distinguishes CAMs from the WT (Figure 1C). We also computed mutational effects on the overall fluctuation, Θ, of selected amino acid sets. The index Θ is a measure of the intrinsic flexibility of the whole protein (or of a given subset of residues) and is proportional to the extent of conformational space explored in a simulation.24 Θ is defined as C

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. Average mutation effect on the overall fluctuations of selected Gt portions. The average RelFluct% index concerning the six CAMs considered in this study is shown. Blue, green, yellow, orange, and red indicate, respectively, values ≥0%, ≥10%, ≥20%, ≥30%, and ≥40%. Here GD and HD stand for the G domain and the H domain, respectively.

increases in overall fluctuation with respect to the helices αDαF. The two switches, i.e. G2 and α2, are the portions of the G domain that undergo the widest increases in mean pairwise distance variance, both between each other or paired with portions of the G or H domains (Figure 2). Remarkably, alanine replacement for T325(α5:1) increases significantly (i.e., by 32%) the overall fluctuations of the Cαatoms in the H domain and the G domain compared to the WT. For the other mutants such increases are on average lower than 10%. In line with the analyses of protein flexibility, CAMs tend to share an increase in the distance between the δ-carbon atom of E112 in the αB/αC loop and the ε-nitrogen atom of K176(G2:3) compared to the WT. The highest and average differences in the E112-K176(G2:3) distance between WT and mutants are 11.3 and 5.7 Å, respectively (Figures S1 and S2). Collectively, Gt CAMs share significant increases in overall fluctuations of GDP and the Gt portions participating in the nucleotide binding site, with emphasis on the G boxes (including the two switches), and the helices with prominence to αD-αF, which may relate with the constitutive activity shared by these mutants. T325(α5:1) distinguishes from the WT and the other CAMs for a marked increase in the overall fluctuations of the G and H domains. Gt CAMs Weaken the Structure Network Involving Packing Interactions of α5, Nucleotide Binding Site, and Interdomain Interface. To infer similarities and differences likely related to function in spite of the heterogeneity of the local perturbations at the mutation site, we investigated longrange propagation of mutational effects by the Protein Structure Network (PSN) analysis, which is based on the application of the graph theory to protein structures.25 In a

root mean distance variance of each atom pair and is calculated by the following equation N

ΘAB =

M

∑i = 1 ∑ j = 1

F

∑k = 1(dijk − dij̅ )2 F

N×M

where A and B are two sets of residues, N and M are the total number of atoms in set A and set B, respectively, and F is the total number of trajectory frames. Furthermore, dkij is the distance between atom i from residue set A and atom j from residue set B in the kth frame, while d̅ij is the average distance between the same two atoms. Here we considered all Cα-atoms in each helix, strand, G box, and GDP, as well as the G and H domains to define the sets of amino acid residues for mean pairwise distance variance computation. The overall fluctuations of such residue sets in the six mutants were then compared with those in the WT thus computing a relative overall fluctuation index according to the following formula: RelFluct% = ((CAM/WT)*100)−100), where CAM and WT indicate the Θ index of sets A and B in the CAM and WT forms, respectively. The Gt portions undergoing increases of Θ in response to the six activating mutations include the G boxes taken singularly or as a whole, the single helices in the H domain, all strands and helices in the G domain, the G and H domains, and the nucleotide (Figure 2). On average, the highest RelFluct% increases involving GDP concern G2 (swI or L2), α2 (participating in swII), and αD-αF, with prominence to αF that is connected to the G domain via L2 (Figure 2). Remarkably, also the overall fluctuations of GDP and the H domain increase as an effect of mutation. In line with GDP, the single G boxes tend to give the most significant D

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 3. Native links weakened by the mutation. A. Native links undergoing an average reduction in frequency ≥10% in the six CAMs considered in this study are shown. Links are colored according to their frequency reduction, whereas nodes are colored according to the average frequency reduction of the links they are involved in. Colors light-cyan, cyan, pale-blue, and blue represent average reductions up to 25%, 50%, 75%, and 100%, respectively. The nucleotide is represented by a greater sphere centered on the N9 atom. B. Functionally relevant nodes involved in native links undergoing mutation-induced frequency reduction ≥30% are shown. Nodes are colored according to their belonging or location: ac1 cluster is magenta, sc is aquamarine, sc2 is violet, sc3 is light green, G boxes are dark green, and interdomain interfaces other than G boxes or sc are lemon green. The nucleotide is black. The diameter of the spheres is proportional to the number of mutants bearing the weakening of the considered link.

the number of links undergoing ≥10% increase in frequency (i.e., 134 and 46, respectively). This is in line with the active states being less stable than the inactive ones as previously shown for a GPCR.28 Thus, the analysis focused on link weakening. Remarkably, 52% of such links involve the clusters of amino acids found responsible for structural integrity and conformational changes associated with activation from alanine scanning mutagenesis on Gi1 (Tables S1 and S2, Figures 3 and S3).16 In particular, a) 25 weakened links involve nodes in the activation cluster I (ac1, Tables S1 and S3) that overlaps with the N-terminal transmission module inferred from the study by Flock et al.14 and that is made of amino acids in β1-β3, β5, α1, and α5 participating in the structural rearrangements that accompany receptor-catalyzed nucleotide exchange; b) 20 weakened links involve nodes in the stabilization cluster II (sc2, Tables S1 and S3) that is made of amino acids thought to stabilize the Rossmann fold of the G domain and located in αA, β4-β6, α3, αG, α4, and α4/β6 loop; c) 18 weakened links participate in the stabilization cluster III (sc3, Tables S1 and S2) made of amino acids in the H domain; finally, 7 weakened links involve nodes in the stabilization cluster at the interdomain interface (sc, Tables S1 and S2) and located in α1, αA, L1 (the loops that connects α1 and αA), and αF.16 The majority of weakened interportion links involving nodes in the activation cluster I connect α5, on one side, and either α1 (T325(α5:1)-Q48(α1:7), V338(α5:14)-I45(α1:4)), the interswitch (i.e., F332(α5:8)-F192(β3:3), or β4 (I339(α5:15)-C216(β4:1)), or β5 (V331(α5:7) -F263(β5:5) and V335 (α5:11) -V261(β5:3) ), or β6 (N327(α5:3)-H318(β6:4), F330(α5:6)-H318(β6:4), and I338(α5:14)Y316(β6:2)), G5 (i.e., V331(α5:7)-T320(G5:2)) on the other side. Other weakened links involving nodes in the activation cluster I include R28-L190(β3:1), L32(β1:4)-I218(β4:3), L34(β1:6)-K42(G1:7), L34(β1:6)-M194(β3:5), I45(α1:4)-A322(G5:4), M49(α1:8)-F192(β3:3), F187 (β2:8) -F192 (β3:3) , C216 (β4:1) -V261 (β5:3) , V261 (β5:3) -

Protein Structure Graph (PSG), each amino acid is represented as a node, and nodes are linked together based upon the strength of their noncovalent interaction. The method implemented in our Wordom software26 basically computes network features (e.g., nodes, hubs (i.e., hyperconnected nodes), links, etc.) and shortest communication pathways on MD trajectories. The employment of MD trajectories instead of a single structure serves to provide a dynamic description of the network as links break and form with atomic fluctuations. We have recently developed a strategy to infer a dynamic structure network even when dealing with a single structure rather than a trajectory. In this case, system’s dynamics is inferred from the coarse grained Elastic Network Model (ENM) paired with Normal Mode Analysis (NMA).27 The approach is herein defined as Protein Structure Networks-Elastic Network Model (PSN-ENM). Differently from the PSN-ENM method, the PSN applied to the MD trajectories is such that only links occurring in a given fraction of the trajectory frames, i.e. link frequency, enter in the PSG. Here, only links occurring in at least 30% of the trajectory frames were considered. The application of PSN to the trajectories of WT and CAMs was instrumental in computing differences in the PSG between mutants and WT. In particular we focused on those links occurring in ≥30% of the WT trajectory frames and undergoing frequency reduction ≥10% in the trajectories of the mutants. Pairwise comparisons were done between the PSG of WT and the PSGs of each mutant as well as between the PSG of WT and a consensus PSG obtained by averaging the frequency of each link over CAM trajectories. The latter comparison is the one we focus on here as it may provide additional structural/ dynamic signatures of mutational effects (Figure 3 and Table S2). The number of native links undergoing ≥10% reduction in frequency as an effect of mutation (i.e., in the consensus network compared to the WT network) is almost three times E

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 4. Overlap matrices of the first 3 eigenvectors obtained from PCA on the MD trajectories of WT and CAMs. On the right top corner, the porcupine representations of PC1 (A), PC2 (B), and PC3 (C) concerning the WT are shown. G and H domains are colored orange and green, respectively.

Y316(β6:2), R309-S259(β5:1), and H318(β6:4)-T320(G5:2) (Table S3). Furthermore, 31% of weakened links by mutation involve either the nucleotide, and/or the G boxes with emphasis on G5, and/or interswitch amino acids (i.e., shadowed in Table S3). Collectively, considering only intersecondary structure element links in the G domain and at the interface between G and H domains, in CAMs the structure network tends to weaken essentially at the interface between a) α1 and αF, G2, interswitch, G3, G5, and α5; b) G1 and β1, G2, αD/αE loop, and α3; c) G2 and G1, α1, αA, and αB/αC loop; d) G5 and α1, G4, β6, and β6/α5 loop, and α5; and e) α5 and α1, interswitch, β4, β5, β6, G5, and β6/α5 loop (Figure 3A and Table S2). Pairwise comparisons between WT and each CAM highlight differences in the effects of mutation (Figure S3). In this

respect, one of the less active mutants (A322(G5:4)S) shows a lower number of weakened links (97) compared to the others, i.e. 128 on average (Figure S3). The difference increases for links showing frequency reduction ≥30% (i.e., 35 for A322(G5:4)S and 72 on average for the other mutants, Figure S3). In the four most active mutants, the most significant link weakening concerns the nucleotide binding site (especially for G179P in L2), the interdomain interface (especially for G56P in L1 and T325(α5:1)A), and the interface between α5 and the interswitch (especially for F332(α5:5)A) (Figure S3). To gain mechanistic insight we analyzed native links undergoing mutation-induced reduction ≥30% and involving relevant amino acids, i.e. the components of the activation and stabilization clusters, the amino acids in the G boxes, and those lying at the interface between G and H domains not included in F

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. MDdMD simulation of the transition between GDP-bound and empty states of Gt catalyzed by the receptor. A. The input (A) and target (B) structures used for the MDdMD trajectory are shown. The G and H domains (plus αN) are orange and green, respectively. On the structures, the distances monitored during the trajectory are shown: a) distance 1 between T86(αA:28) and D234 (in the β4/α3 loop); b) distance 2 between S153(αE:7) and K329(α5:5); c) distance 3 between E167(αF:1) and E272(αG:6); d) distance 4 between D137(αD:8) and K329(α5:5); e) distance 5 between A134(αD:5) and E272(αG:6); and f) distance 6 between E112 (in the αC/αD loop) and K176(G2:3). Distances 1−5 were monitored between Cα-atoms, while distance 6 was monitored between the δ-carbon atom of E112 and the ε-nitrogen atom of K176(G2:3). C. The superimposed structures of the G domain plus αN concerning the GDP-bound (orange) and nucleotide-free states of Gt (gray) are shown. D. The porcupine representation of the first principal component from PCA on the first 3 ns of the MDdMD trajectory is shown. White arrows indicate maximum displacements of Cαatoms along the component. The G domain and the H domain are orange and green, respectively.

the preceding categories. Weakened links involving GDP were considered as well. The distribution of the nodes involved in such links onto the 3D structure clearly shows that activating mutations tend to affect nodes in the activation cluster I (at the interface between α5 and the surrounding portions like α1, the interswitch, and G5), in the G boxes, and in the stabilization cluster cluster at the interdomain interface (Figure 3B and Table S3). By considering that the majority of the G boxes, which make the nucleotide binding site, locate at the interface between G and H domains (especially G2 (swII) and G3 (swII)), it results that the relevant interactions weakened/lost as an effect of the activating mutation involve the interfaces between α5 and the surrounding portions and between G and H domains (Figure 3B and Table S3).

In summary, PSN analysis highlighted native link weakening in selected regions as possibly related to the constitutive activity of the six Gt mutants. The Essential Dynamics of Gt CAMs Overlaps in Part with the Principal Mode in Receptor-Catalyzed Transition from GDP-Bound to Nucleotide-Free States. We also inferred mutational effects on the essential dynamics of the G protein, focusing on the collective motions of the H domain with respect to the G domain. For all considered forms, the first three PCs cover ≥50% of the total variance defining the essential subspace. Therefore, PCA focused on them. For the WT, PC1 describes a concerted motion of the H domain around the main axis of αA, whereas PC2 and PC3 concerted detachments of the H domain from the G domain taking the C-terminal end of α1 as a hinge (Figure 4 and S4). For G59P, V328(α5:4)A, and G

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 6. Overlaps between MD and MDdMD trajectories. A. The overlaps between the first 3 eigenvectors from the MD trajectories of WT and mutants and the first eigenvector from the MDdMD trajectory are shown. The column “Sum” is the combination between the two eigenvectors that produce the highest overlap, i.e. PC2 and PC3, for the WT, PC1 and PC2 for G56P, PC2 and PC3 for G179(G2:6)P, PC1 and PC3 for A322(G5:4)S, PC1 and PC2 for T325(α5:1)A, PC1 and PC2 for T328(α5:4)A, and PC2 and PC3 for F332A(α5:8)A. B. For the WT and CAMs, the PCs either single or combined, which give the highest overlap with the PC1 of MDdMD, are shown in porcupine representation. In detail, combined eigenvectors are used for all forms but A322(G5:4)S and F332A(α5:8)A, for which no overlap improvement resulted from the pairwise combination of eigenvectors and the single PC1 and PC2 are shown, respectively. The G domain and the H domain are orange and green, respectively.

F332(α5:8)A, PC1 displays very high overlaps (i.e., ≥ 0.8) with the PC1 of WT (Figure 4 and S4). In contrast, the PC1 of the other three mutants, G179(G2:6)P, T325(α5:1)A, and A322(G5:4)S, overlaps between each other more than they do with the PC1 of WT (Figures 4 and S4). Pairwise comparisons of these three mutants display the highest overlaps along the diagonal (Figure 4). Highest overlaps along the diagonal concern also G56P and the WT, as well as V328(α5:4)A and F332(α5:8)A (Figure 4). To investigate if activating mutations would mimic some of the structural/dynamical effects exerted by the GEF on Gt, we employed MDdMD to simulate the transition of WT Gt from the receptor-free GDP-bound state to the receptor-bound nucleotide-free state. MDdMD is a method for estimating pathways for conformational transitions in macromolecules by means of discrete MD and biasing techniques based on a combination of essential dynamics and Maxwell-Demon sampling techniques.29 These computational experiments were based on the funded assumption that Gt and Gs, like the other Gα proteins, share the same mechanism of GDP release. The starting structure was the α-subunit extracted from the heterotrimeric GDP-bound state of the GtGi chimera, following mutation of the Gi portion into the corresponding sequence of bovine Gt. The last ten amino acids of the protein, missing in the crystal structure, were modeled according to the crystal structure of the C-terminal peptide of Gt in complex with opsin.30 The addition of the C-term was instrumental in using an input for MDdMD simulations with the same length as the target. In this case, the simulated transition from GDPbound to nucleotide-free forms of Gt did not include folding of the C-term of α5, which accompanies receptor binding.14,31 The choice of neglecting conformational changes in the Cterminus of α5 was also dictated by the absence of such portion in all G protein forms simulated by classical MD. The target structure representing nucleotide-free Gt was achieved by

comparative modeling (by means of Modeller32) by using the structure of Gs extracted from the crystallographic complex with the β2-AR (PDB code: 3SN69) as a template. The most surprising feature of the β2-AR-Gs complex is the dramatic displacement of the H domain with respect to the G domain. Indeed, the distance between the C-term of αA (i.e., the Cα-atom of N112(αA:28)) and swIII (i.e., the Cα-atom of T263) increases by 43 Å on going from the receptor-free GDPbound state to the receptor-bound nucleotide-free state of Gs.31 Analogously to Gs, the distance that increases the most for Gt following nucleotide release is that between T86(αA:28) and D234 (swIII) (Figure 5A,B and S5). Other structural differences between GDP-bound and the modeled nucleotide-free state of Gt concern α5, which is slightly misfolded at the N-terminus and is pulled by one turn as an effect of receptor binding, and the β6/α5 loop (including G5) whose conformational change is instrumental in GDP release. The pull of α5 is also evidenced by its position relative to the interswitch, whose turn interacts with F336(α5:12) in the inactive state while interacting with F332(α5:8) in the receptorcatalyzed active state (Figure 5C). Pulling of α5 upon receptor binding weakens the aromatic/hydrophobic residue-mediated packing of α5, α1, and β1-β3 while reinforcing packing of α5 against β5 and β6 (Figure 5C). Interestingly, the distance between E112 and K176(G2:3) that tends to increase in the majority of simulated CAMs of Gt increases as well, i.e. from 13.4 Å to about 30.0 Å, during the transition from GDP-bound and nucleotide-free states of Gt (Figure S5). The transition occurs in the first 3 ns of the MDdMD simulation, and it is essentially accounted for by the PC1 that describes a rotational outward motion of the H domain from the Ras-like domain, taking the C-terminal end of α1 as a hinge (Figures 5D and 6). H

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 7. Structure comparisons and global metapaths achieved by the PSN-ENM method. A, B. The superimposed crystal structures of WT Gt (1TAG) and G56P CAM (3V00, chain C) are shown. 1TAG is gray, whereas 3V00 is colored so as to distinguish the G domain (orange) and the H domain (green). C, D. The global metapaths computed on the crystal structures of WT Gt (C) and G56P CAM (D) are shown.

Remarkably, such essential motion computed over the first 3 ns of the MDdMD trajectory overlaps with either PC1 or PC2 from the 500 ns trajectories of the WT and CAMs (Figure 6A). In detail, for G56P, A322(G5:4)S, and T325(α5:1)A the highest overlap concerns PC1. Remarkably, the PC1 of T325(α5:1)A and that of MDdMD give an overlap equal to 0.7, which increases to 0.8 on combining PC1 and PC2 (Figure 6A,B). The same happens for G56P, even if with lower overlaps, i.e. 0.6 and 0.7, respectively. For A322(G5:4)S the highest overlap concerns only PC1, and no combination is able to improve it. The second highest overlap concerning a single PC (0.7) is reached by the PC2 of F332(α5:8). For the remaining forms, PC2 reaches the highest overlap, which however has to be combined with either PC1 or PC3 to surpass 0.6 (Figure 6A,B). Collectively, the PC1 of MDdMD gives significant overlaps (i.e., ≥ 0.6) with either PC1 or PC2 or a combination of the two for the majority of CAMs, suggesting that the majority of activating mutations considered in this study amplify a concerted outward motion of the H domain with respect to the G domain. Such motion overlaps with the mode describing the receptor-catalyzed transition associated with GDP release.

Thus, rhodopsin and selected CAMs of Gt, with emphasis on G56P, T325(α5:1)A, and F332A(α5:8)A, prompt for the detachment of the G and H domains accompanying nucleotide release. The Inferences from Computational Experiments Match with X-ray Crystallographic Determinations. The X-ray crystal structure has been recently released for the GDPbound state of the G56P CAM (PDB code 3V00).20 The structure diverges from the WT (PDB code: 1TAG) in the three switches and in the α/β loops, in particular the α4/β6 loop (Figure 8A). Similarly to the simulated structures of Gt CAMs achieved in this study, the crystal structure of the G56P mutant does not show a significant alteration in the relative juxtaposition of the G domain relative to the H domain, presumably because the G56P mutant has retained bound GDP, which helps maintaining a “closed” conformation rather than an “open” state.20 This suggests that, similarly to our simulated CAMs, the X-ray structure of GDP-bound G56P might mimic an intermediate state along the receptor-mediated nucleotide exchange pathway that would form prior to the actual release of GDP.20 Remarkably, and analogously to the simulated structures of Gt CAMs, no pull of α5 is observed in the crystal I

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 8. Comparisons of the collective motions of the G56P mutant from MD simulations of mutated 1TAG structure and the ENM-NMA on the crystal structure of the mutant (PDB code: 3V00, chain C). The overlap matrix between the first 3 eigenvectors obtained by PCA on G56P-Gt MD trajectory (rows) with the first 3 normal modes obtained from ENM-NMA on the crystal structure of the G56P mutant (columns) are shown. On the right, the porcupine representation of the PC1 from the G56P-Gt MD trajectory and the first normal mode from the ENM-NMA analysis of 3V00 are shown. Both structures are colored according to the secondary structure. PCA was done on the Cα-atoms 27-340 (i.e., both the G and H domains).

(Figure 7C). This suggests that, in the inactive state, the α5 of Gt communicates with αA in the H domain through α1 and L1. Such communication involving hydrophobic amino acids may confer stiffness to the region that holds the hinge point for the rotation of the H domain away from the G domain. Interestingly, this communication does not characterize the G56P CAM (Figure 7D), in line with the increase of the overall fluctuations of the two domains and with the results of PCA, highlighting collective motions of the H domain as well as of linkers and switches, as intrinsic features of Gα amplified by activating mutations.

structure of the G56P mutant compared to the WT form (Figure 7A,B), which is suggestive of divergent local perturbations as productive for GDP release catalyzed by receptor and activating mutations. The average Cα-RMSD between the crystal structure of G56P and the 500000 frames constituting the trajectory achieved following transformation of WT Gt in the G56P mutant is 2.19 ± 0.18 Å over the whole structure. Consistent with the results of simulations of the G56P mutant, the L1 in the crystal structure is characterized by an increase in distance between E112 and K176(G2:3) compared to the WT. Indeed such distance is 10.9 Å in the structure of the WT (PDB code: 1TAG) and 16.8 in the structure of the mutant (PDB code 3V00, chain C), thus corroborating the hypothesis that this is a structural feature of the CAM states. We also compared the principal components inferred from the MD trajectory of the mutant with the normal modes achieved by applying ENM-NMA to the crystal structure. Interestingly, the PC1 computed on the MD trajectory overlaps with the first normal mode from ENM-NMA (overlap 0.8; Figure 8), thus strengthening the concerted motions of switches, linkers, and helical domain as a feature of the Gα structure. We then compared the structural communication features of WT and G56P mutant Gt by applying the mixed PSN-ENM approach to the 1TAG and 3V00 crystal structures, respectively. The approach combines ENM-NMA with the analysis of Protein Structure Networks (PSN) to compute communication pathways on single structural models from Xray crystallography, NMR, or structure predictions.27,33,34 For each of the two proteins, a coarse view of the global communication propensity was inferred through global metapaths, i.e. coarse paths made of the most recurrent nodes and links in the path pool. In this study, metapaths are made of nodes present in ≥50% of the path pool (i.e., “frequent nodes”) and of links satisfying both conditions of being present in ≥50% of the paths and of connecting “frequent nodes”. The global metapath of WT Gt is made of nodes at the interface between α1 and α5, in L1, and along the main axis of αA

■

DISCUSSION The details of the early events in receptor-induced nucleotide exchange in heterotrimeric G proteins and the role of the βγ dimer in such process still await elucidation. Different mechanisms of receptor-catalyzed GDP release have been proposed (reviewed in ref 31). Significant insights in that respect were provided by biophysical experiments on Gi1,11,16,35−37 which identified a possible allosteric pathway propagated along swI to αF, that, like α5 and the β6/α5 loop, forms part of a putative GDP exit route.36 Collectively, the results suggested that the receptor uses contacts with the extreme C-term of Gα to communicate structural changes through α5 presumably to the β6/α5 loop and induce the release of GDP. The studies provided also evidence that GDP release is associated with significant displacement of the H domain with respect to the G domain,11 which found a high resolution picture in the crystal structure of the ternary complex between agonist-activated β2-AR and nucleotide-depleted heterotrimeric Gs, representing a late state in the pathway of receptor-catalyzed GDP release. In such a complex, the 6 Å displacement toward the receptor and rotation of α5 and the shift of the β6/α5 loop outward from the nucleotide binding cleft, as well as changes in conformation of β1 and the P loop, generated the hypothesis that α5 and β1 mediate receptorinduced signal propagation from the G domain to the H domain.9 The dramatic displacement of the H domain with respect to the G domain in the β2-AR-Gs complex is J

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

increase the overall fluctuations of GDP as well. In fact, the mean pairwise distance variance of the nucleotide on one side and G2 (swI or L2), α2 (participating in swII), αD-αF, with prominence to αF, and the whole H domain increases significantly as an effect of mutation. The same is observed for the G domain with respect to the H domain, whose mutation-induced increase in overall fluctuation is maximal for the T325(α5:1). These effects on protein flexibility relate with weakening of links in the native protein structure network. Remarkably, 31% of weakened links by mutation involve either the nucleotide, and/or the G boxes with emphasis on G5, and/ or interswitch amino acids. Moreover, 52% of weakened links involve the clusters of amino acids found responsible for structural integrity and conformational changes associated with activation in an alanine scanning mutagenesis study on Gi1.16 Collectively, considering only intersecondary structure element links in the G domain and at the interface between G and H domains, in CAMs the structure network tends to weaken essentially at the interface between a) α1 and αF, G2, interswitch, G3, G5, and α5; b) G1 and β1, G2, αD/αE loop, and α3; c) G2 and G1, α1, αA, and αB/αC loop; d) G5 and α1, G4, β6, and β6/α5 loop, and α5; and e) α5 and α1, interswitch, β4, β5, β6, G5, and β6/α5 loop. With the analysis of weakened native links involving functionally relevant nodes, i.e. the nucleotide, the members of the activation and stabilization clusters, the nodes lying on the G boxes and at the interface between G and H domains, we could infer that the relevant interactions weakened/lost as an effect of activating mutation involve the interfaces between α5 and the surrounding portions and between G and H domains. Starting from the funded assumption that Gtα and Gsα share the same mechanism of receptor-catalyzed GDP release, we modeled the receptor-induced nucleotide-free state of Gt by using the crystal structure of nucleotide-free Gs as a template. MDdMD simulations of the transition between GDP-bound and GDP-depleted Gt allowed us to infer hypotheses on the collective motions associated with the receptor-catalyzed uncoupling of H and G domains. Interestingly, the major mode over the 3 ns of the MDdMD simulation, which consists in a displacement of the H domain employing α1 as a hinge, gives overlaps higher than 0.6 with either PC1 or PC2 or a combination of the two for the majority of CAMs (G56P, T325(α5:1)A, V328(α5:4)A, and F332(α5:5)A). For A322(G5:4)S the highest overlap is modest (0.5) but involves PC1. In contrast, for the WT and G179(G2:6)P, an overlap with the PC1 from MDdMD is obtained only by the combination of less relevant eigenvectors, PC2 and PC3. T325(α5:1)A, a mutation characterized by release/weakening of interactions at the interface between the beginning of α5 and α1, as well as at the interface between G and H domains shows the highest overlap. These results are in line with mutation-induced increases in interdomain mean pairwise distance variance, T325(α5:1)A, showing the most marked effect. These results suggest that the uncoupling of H and G domains is an intrinsic feature of the Gα structure and receptors or certain activating mutations prioritize such motion. To induce these effects, the receptor interacts with and folds the C-term of α5, pulls α5, and induces long-distance misfolding of the first turn of the helix, thus allosterically destabilizing the conserved packing interactions between the Nterminal half of α5 and both α1 and the interswitch. In contrast, G56P and G179(G2:6)P affect the flexibility of L1 and L2, respectively, A322(G5:4)S perturbs the interactions between G4

significantly greater than the maximal displacement estimated on the Gi-rhodopsin complex11 raising the question whether it may represent a physiological condition or an artifact of crystallization. In this respect, it has been argued that the position of the H domain in the crystallographic β2-AR-Gs complex may reflect only one of an ensemble of conformations that it can adopt under physiological conditions, which has been stabilized by crystal packing interactions. Collectively, biophysical and computational experiments together with X-ray crystal structure determinations provide funded evidence that the receptor-catalyzed nucleotideexchange on heterotrimeric G proteins requires a significant displacement of the H domain away from the G domain.9−11 The triggers of such displacement would include motions of α5, which represents a fundamental recognition point for the receptor, and changes in the conformation of the P-loop and the β6/α5 loop, which favors GDP exit.9,10,13 Uncoupling of the H and G domains would be a consequence of nucleotide release or at least a coincident event.9 This evidence was recently strengthened by microsecond MD simulations on receptor-bound nucleotide-free heterotrimeric Gs and on receptor-free GDP-bound/GDP-unbound heterotrimeric Gt, which showed that the H and G domains separate spontaneously even in the absence of the receptor and that domain separation is necessary but not sufficient for rapid nucleotide release.12 Rather, receptors catalyze nucleotide release by favoring an internal structural rearrangement of the Ras domain involving α5 and the β6/α5 loop that weakens its nucleotide affinity. The inferences above agree with the outcome of a computational study that inferred a universal allosteric mechanism for Gα activation by GPCRs based on the analysis of functional state-dependent consensus contact networks in 80 crystallographic Gα protein structures coupled with the information on amino acid conservation from sequence alignment.14 The analysis, which emphasized the role of α5, suggested that GPCRs interact and activate Gα proteins through a highly conserved mechanism, in which the interruption of the contacts between α1 and α5 is a key step of GDP release. It was inferred that while α1 is the molecular switch for GDP release, α5 is the distal trigger that is pulled upon receptor binding.14 Along this line, an alanine scanning mutagenesis study on Gi16 found an activation cluster I consisting of residues in α1 and α5 packed against residues in β1-β3 in the nucleotide-bound states. In the receptor-bound state, the interactions between α5-α1 and β1-β3 are weakened and are compensated for by a new set of interactions between α5 and strands β4-β6. Overall, the study suggested that the most consequential event in activation of Gi1 is the destabilization of α1 caused by a rearrangement on activation cluster I. This leads to a perturbation and weakening of the interdomain interface, dissociation of the helical domain from the GTPase domain in a rigid-body movement and release of the GDP.16 In the present study, comparative MD simulations on WT Gt and on six CAMs located in the two interdomain linkers as well as in G5 and α5 served to investigate the structural determinants of constitutive activation of Gt by mutation. Possible fingerprints of mutational effects arose from the overall fluctuation and PSN analyses. Indeed, compared to WT, Gt CAMs share the increase in mean pairwise distance variance of selected portions, participating in the nucleotide binding site, with emphasis on the G boxes, including the two switches and the helices with prominence to αD-αF. Activating mutations K

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation and G5, V328(α5:4)A and T325(α5:1)A perturb the interactions between α5 and α1, and F332(α5:8)A perturbs the interactions of α5 with β1, α1, and interswitch. Thus, switching packing interactions are perturbed directly by activating mutations and allosterically by the GPCR.

■

constrain all bond lengths, allowing for an integration time step of 2 fs through the leapfrog algorithm. The v-rescale thermostat43 was employed to keep the system at a constant temperature of 300 K, by using a coupling constant (τt) of 0.1 ps. The pressure of the system was kept fixed at 1 atm, using the Berendsen weak coupling algorithm44 with a coupling constant (τp) of 1 ps. The pre-equilibrated systems were then subjected to 500 ns of unrestrained isothermal−isobaric (T = 300 K, P = 1 atm) MD simulations (the RMSD time series are shown in Figure S6). MD Analyses of the Intrinsic Flexibility. MD trajectories were subjected to a variety of analyses aimed at inferring the time series of a number of structural descriptors (e.g., interatomic distances, Cα-RMSD), the intrinsic flexibility (e.g., Cα-RMSF, overall fluctuations,24 and Essential Dynamics (ED) or PCA45) of the systems. All these analyses were performed by means of the Wordom software.26 As for PCA, the covariance matrices were built on the Cα-atoms of the isolated MD trajectories, by fitting the Cα-atoms of the G domain. ENM-NMA on the crystal structure of the G56P mutant was carried out by means of Wordom as well.26 MDdMD Simulations. MDdMD simulations were carried out by using a Web-based application accessible at http://mmb. irbbarcelona.org/MDdMD. The starting structure was the αsubunit extracted from the heterotrimeric GDP-bound state of the GtGi chimera, following mutation of the Gi portion into the corresponding sequence of bovine Gt. The last ten amino acids of the protein, missing in the crystal structure, were modeled according to the crystal structure of the C-terminal peptide of Gt in complex with opsin.30 The target structure representing nucleotide-free Gt was achieved by comparative modeling (by means of Modeller32) by using the structure of Gsα extracted from the crystallographic complex with β2-AR (PDB code 3SN6)9 as a template. PSN Analysis on the MD Trajectories of WT and CAMs. Building of the PSG is carried out by means of the PSN module implemented in the Wordom software.26 PSN analysis is a product of graph theory applied to protein structures.46 A graph is defined by a set of vertices (nodes) and connections (edges) between them. In a PSG, each amino acid residue is represented as a node, and these nodes are connected by edges based on the strength of noncovalent interactions between residues.47 The strength of interaction between residues i and j (Iij) is evaluated as a percentage given by the following equation

CONCLUSIONS

The present study allowed us to infer a number of structural determinants of the constitutive activity held by six mutants of Gt. In spite of the differences in the local perturbations induced by the mutation, the six activating mutations of Gt increase the overall fluctuations of the nucleotide and its binding site. This is accompanied by the weakening of native links involving GDP, as well as the interfaces between α5 and the surrounding portions and the interface between G and H domains, which involve the G boxes. More than 50% of the links weakened in the consensus network of the six mutants involve amino acids found to contribute to the structural stability and the activation process of Gi1.16 These effects of mutations associate with increases in the overall fluctuations of the G and H domains, which reflect on the collective motions of the proteins. Indeed, compared to the WT, activating mutations of Gt, with prominence to G56P, T325(α5:1)A, and F332(α5:8)A, prioritize also collective motions of the H domain overlapping with the collective motions associated with GPCR-catalyzed GDP release. In spite of different local perturbations, the mechanisms of nucleotide exchange catalyzed by activating mutations and by receptor are expected to employ similar molecular switches in the nucleotide binding site and to share the detachment of the H domain from the G domain.

■

METHODS MD Simulations. The structural models of the constitutively active mutants of Gt, G56P, G179(G2:6)P, A322(G5:4)S, V328(α5:4)A, T325(α5:1)A, and F332(α5:8)A were achieved by amino acid replacements in the crystal structure of the GDP state (PDB code: 1TAG). All the simulated systems hold the Mg2+ ion together with the coordinating water molecules. Due to missing amino acids at the two terminations, the N-terminus was acetylated and the C-terminus was N-methyl-amidated. MD simulations were carried out by using the GROMACS4 simulation package38 with the AMBER03 all atoms force field.39,40 The TIP3P water model was employed to describe the solvent. AMBER parameters to describe the GDP and GTP molecules were taken from the literature.41 Periodic Boundary Conditions (PBC) were applied by using an octahedric box as a unit cell, imposing a minimum distance of 12 Å between the solute and the box boundaries. A 0.1 M concentration of NaCl was employed while using a slightly higher number of Na+ ions (i.e., Na+ = 47, Cl− = 38) to neutralize the system. MD simulation setup is the same as the one recently employed to simulate a number of Ras GTPases.7 All the input crystallographic structures were subjected to energy minimization keeping restricted the positions of the main chain atoms, the nucleotide, the Mg2+ cation, and the coordinating water molecules. The systems were then equilibrated at 300 K for 4 ns of backbone restricted MD simulations. The Particle Mesh Ewald (PME) method was employed to compute the electrostatic interactions. Short range repulsive and attractive interactions were computed by using a Lennard-Jones potential with a cutoff of 10 Å. The LINCS algorithm42 was used to

Iij =

nij NN i j

× 100

where Iij is the percentage interaction between residues i and j; nij is the number of atom−atom pairs between the side chains of residues i and j within a distance cutoff (4.5 Å); Ni and Nj are normalization factors for residue types i and j, which account for the differences in size of the amino acid side chains and their propensity to make the maximum number of contacts with other amino acids in protein structures. The normalization factors for the 20 amino acids were derived from the work by Kannan and Vishveshwara.48 The normalization value for GDP, 220.19, was derived from a database of 55 G proteins. The normalization factors of Mg2+ were 14.65 in the GDP-bound form (based on 41 structures). Finally, the normalization factor for water, 27, was computed on the crystal structures of rhodopsin and four Gα proteins (PDB codes: 1GZM, 1CIP, 1CUL, 1TAG, and 1TND). L

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

■

Thus, Iij is calculated for all nodes, excluding i ± n, where n is a given neighbor cutoff, 2. An interaction strength cutoff Imin is then chosen and any residue pair ij for which Iij ≥ Imin is considered to be interacting and hence connected in the PSG. Node interconnectivity is finally used to highlight clusterforming nodes, where a cluster is a set of connected amino acids in a graph. The node clustering procedure is such that nodes are iteratively assigned to a cluster if they could establish a link with at least one node in such cluster. A node not linkable to existing clusters initiates a novel cluster and so on until the node list is exhausted. Cluster size, defined as the number of nodes, varies as a function of the Imin, and the size of the largest cluster is used to calculate the Icritic value. The latter is defined as the Imin at which the size of the largest cluster is half the size of the largest cluster at Imin = 0.0%. At Imin = Icritic weak node interactions are discarded, emphasizing the effects of stronger interactions on PSN properties. We set the Imin equal to the Icritic approximated to the second decimal place. Therefore, it is possible to obtain different PSGs for the same protein structure depending on the selected Imin, and, consequently, Imin can be varied to obtain graphs with strong or weak interactions forming the edges between the residues. Only edges occurring in a given fraction of the trajectory frames, i.e. link frequency, enter in the PSG. A link-frequency cutoff of 30% was employed here. Analysis of the Communication Pathways by Means of the PSN-ENM Method. The structural communication features of WT and G56P mutant Gt were investigated by applying the mixed PSN-ENM approach to the 1TAG and 3V00 crystal structures, respectively. The approach combines ENM-NMA with PSN analysis to compute communication pathways on single structural models from X-ray crystallography, NMR, or structure predictions.27,33,34 The structural communication between two distal sites is inferred from the computation of the shortest paths (i.e., chains of pairwise linked nodes) between pairs of nonlinked nodes. The first step consists in performing the PSN analysis on a single high resolution structure, which means building the Protein Structure Graph (PSG). The latter is searched for the shortest paths between pairs of nodes. The algorithm first defines all possible communication paths between selected node pairs, and then it filters the results according to crosscorrelation of atomic motions, as derived from ENM-NMA. Filtering consists in retaining only the shortest path(s) that contains at least one residue correlated (e.g., with a crosscorrelation ≥0.8 in this study) with either one of the two extremities (i.e., the first and last amino acids in the path). Metapaths made of the most recurrent nodes and links in the path pool (i.e., global metapaths) are worth computing to infer a coarse/global picture of the structural communication in the considered system.27 In this study, metapaths were made of nodes present in ≥50% of the considered path pool (i.e., “frequent nodes”) and of links satisfying both conditions of being present in ≥50% of the paths and of connecting “frequent nodes”.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Francesca Fanelli: 0000-0002-7620-6895 Present Addresses †

University of Heidelberg, Im Neuenheimer Feld 267, 69120 Heidelberg, Germany. ‡ NEST, Istituto Nanoscienze-CNR, Piazza San Silvestro 12, 56127 Pisa, Italy. Author Contributions

A.F. and S.M. contributed equally. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding

This work was supported by Telethon-Italy grants [GGP11210 and GGP13227A] to F.F. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The employment of the PyMOL 0.99rc6 software for the realization of all drawings is acknowledged.

■

REFERENCES

(1) Vetter, I. R.; Wittinghofer, A. The guanine nucleotide-binding switch in three dimensions. Science 2001, 294, 1299−1304. (2) Colicelli, J. Human RAS superfamily proteins and related GTPases. Sci. Signaling 2004, 2004, RE13. (3) Gilman, A. G. G proteins: transducers of receptor-generated signals. Annu. Rev. Biochem. 1987, 56, 615−649. (4) Clapham, D. E. The G-protein nanomachine. Nature 1996, 379, 297−299. (5) Oldham, W. M.; Hamm, H. E. Structural basis of function in heterotrimeric G proteins. Q. Rev. Biophys. 2006, 39, 117−166. (6) Wittinghofer, A.; Vetter, I. R. Structure-function relationships of the G domain, a canonical switch motif. Annu. Rev. Biochem. 2011, 80, 943−971. (7) Raimondi, F.; Portella, G.; Orozco, M.; Fanelli, F. Nucleotide binding switches the information flow in ras GTPases. PLoS Comput. Biol. 2011, 7, e1001098. (8) Noel, J. P.; Hamm, H. E.; Sigler, P. B. The 2.2 A crystal structure of transducin-alpha complexed with GTP gamma S. Nature 1993, 366, 654−663. (9) Rasmussen, S. G.; Devree, B. T.; Zou, Y.; Kruse, A. C.; Chung, K. Y.; Kobilka, T. S.; Thian, F. S.; Chae, P. S.; Pardon, E.; Calinski, D.; Mathiesen, J. M.; Shah, S. T.; Lyons, J. A.; Caffrey, M.; Gellman, S. H.; Steyaert, J.; Skiniotis, G.; Weis, W. I.; Sunahara, R. K.; Kobilka, B. K. Crystal structure of the beta(2) adrenergic receptor-Gs protein complex. Nature 2011, 477, 549−555. (10) Raimondi, F.; Seeber, M.; De Benedetti, P. G.; Fanelli, F. Mechanisms of inter- and intramolecular communication in GPCRs and G proteins. J. Am. Chem. Soc. 2008, 130, 4310−4325. (11) Van Eps, N.; Preininger, A. M.; Alexander, N.; Kaya, A. I.; Meier, S.; Meiler, J.; Hamm, H. E.; Hubbell, W. L. Interaction of a G protein with an activated receptor opens the interdomain interface in the alpha subunit. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 9420−9424. (12) Dror, R. O.; Mildorf, T. J.; Hilger, D.; Manglik, A.; Borhani, D. W.; Arlow, D. H.; Philippsen, A.; Villanueva, N.; Yang, Z.; Lerch, M. T.; Hubbell, W. L.; Kobilka, B. K.; Sunahara, R. K.; Shaw, D. E. SIGNAL TRANSDUCTION. Structural basis for nucleotide exchange in heterotrimeric G proteins. Science 2015, 348, 1361−1365. (13) Scheerer, P.; Heck, M.; Goede, A.; Park, J. H.; Choe, H. W.; Ernst, O. P.; Hofmann, K. P.; Hildebrand, P. W. Structural and kinetic

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00813. Tables S1−S3 and Figures S1−S4 (PDF) M

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation modeling of an activating helix switch in the rhodopsin-transducin interface. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 10660−10665. (14) Flock, T.; Ravarani, C. N.; Sun, D.; Venkatakrishnan, A. J.; Kayikci, M.; Tate, C. G.; Veprintsev, D. B.; Babu, M. M. Universal allosteric mechanism for Galpha activation by GPCRs. Nature 2015, 524, 173−179. (15) Goricanec, D.; Stehle, R.; Egloff, P.; Grigoriu, S.; Pluckthun, A.; Wagner, G.; Hagn, F. Conformational dynamics of a G-protein alpha subunit is tightly regulated by nucleotide binding. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, E3629−3638. (16) Sun, D.; Flock, T.; Deupi, X.; Maeda, S.; Matkovic, M.; Mendieta, S.; Mayer, D.; Dawson, R. J.; Schertler, G. F.; Babu, M. M.; Veprintsev, D. B. Probing Galphai1 protein activation at single-amino acid resolution. Nat. Struct. Mol. Biol. 2015, 22, 686−694. (17) Hofmann, K. P.; Scheerer, P.; Hildebrand, P. W.; Choe, H. W.; Park, J. H.; Heck, M.; Ernst, O. P. A G protein-coupled receptor at work: the rhodopsin model. Trends Biochem. Sci. 2009, 34, 540−552. (18) Marin, E. P.; Krishna, A. G.; Sakmar, T. P. Rapid activation of transducin by mutations distant from the nucleotide-binding site: evidence for a mechanistic model of receptor-catalyzed nucleotide exchange by G proteins. J. Biol. Chem. 2001, 276, 27400−27405. (19) Majumdar, S.; Ramachandran, S.; Cerione, R. A. Perturbing the linker regions of the alpha-subunit of transducin: a new class of constitutively active GTP-binding proteins. J. Biol. Chem. 2004, 279, 40137−40145. (20) Singh, G.; Ramachandran, S.; Cerione, R. A. A constitutively active Galpha subunit provides insights into the mechanism of G protein activation. Biochemistry 2012, 51, 3232−3240. (21) Raimondi, F.; Orozco, M.; Fanelli, F. Deciphering the deformation modes associated with function retention and specialization in members of the Ras superfamily. Structure 2010, 18, 402− 414. (22) Fanelli, F.; Raimondi, F. Nucleotide binding affects intrinsic dynamics and structural communication in Ras GTPases. Curr. Pharm. Des. 2013, 19, 4214−4225. (23) Mariani, S.; Dell’Orco, D.; Felline, A.; Raimondi, F.; Fanelli, F. Network and atomistic simulations unveil the structural determinants of mutations linked to retinal diseases. PLoS Comput. Biol. 2013, 9, e1003207. (24) Munz, M.; Hein, J.; Biggin, P. C. The role of flexibility and conformational selection in the binding promiscuity of PDZ. domains. PLoS Comput. Biol. 2012, 8, e1002749. (25) Brinda, K. V.; Vishveshwara, S. A network representation of protein structures: implications for protein stability. Biophys. J. 2005, 89, 4159−4170. (26) Seeber, M.; Felline, A.; Raimondi, F.; Muff, S.; Friedman, R.; Rao, F.; Caflisch, A.; Fanelli, F. Wordom: A user-friendly program for the analysis of molecular structures, trajectories, and free energy surfaces. J. Comput. Chem. 2011, 32, 1183−1194. (27) Raimondi, F.; Felline, A.; Seeber, M.; Mariani, S.; Fanelli, F. A Mixed Protein Structure Network and Elastic Network Model Approach to Predict the Structural Communication in Biomolecular Systems: The PDZ2 Domain from Tyrosine Phosphatase 1E As a Case Study. J. Chem. Theory Comput. 2013, 9, 2504−2518. (28) Gether, U.; Ballesteros, J. A.; Seifert, R.; Sanders-Bush, E.; Weinstein, H.; Kobilka, B. K. Structural instability of a constitutively active G protein-coupled receptor. Agonist-independent activation due to conformational flexibility. J. Biol. Chem. 1997, 272, 2587−2590. (29) Sfriso, P.; Emperador, A.; Orellana, L.; Hospital, A.; Gelpi, J. L.; Orozco, M. Finding Conformational Transition Pathways from Discrete Molecular Dynamics Simulations. J. Chem. Theory Comput. 2012, 8, 4707−4718. (30) Scheerer, P.; Park, J. H.; Hildebrand, P. W.; Kim, Y. J.; Krauss, N.; Choe, H. W.; Hofmann, K. P.; Ernst, O. P. Crystal structure of opsin in its G-protein-interacting conformation. Nature 2008, 455, 497−502. (31) Fanelli, F.; De Benedetti, P. G. Update 1 of: Computational Modeling Approaches to Structure-Function Analysis of G ProteinCoupled Receptors. Chem. Rev. 2011, 111, PR438−535.

(32) Sali, A.; Blundell, T. L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993, 234, 779−815. (33) Seeber, M.; Felline, A.; Raimondi, F.; Mariani, S.; Fanelli, F. WebPSN: a web server for high-throughput investigation of structural communication in biomacromolecules. Bioinformatics 2015, 31, 779− 781. (34) Raimondi, F.; Felline, A.; Fanelli, F. Catching Functional Modes and Structural Communication in Dbl Family Rho Guanine Nucleotide Exchange Factors. J. Chem. Inf. Model. 2015, 55, 1878− 1893. (35) Oldham, W. M.; Van Eps, N.; Preininger, A. M.; Hubbell, W. L.; Hamm, H. E. Mechanism of the receptor-catalyzed activation of heterotrimeric G proteins. Nat. Struct. Mol. Biol. 2006, 13, 772−777. (36) Oldham, W. M.; Van Eps, N.; Preininger, A. M.; Hubbell, W. L.; Hamm, H. E. Mapping allosteric connections from the receptor to the nucleotide-binding pocket of heterotrimeric G proteins. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 7927−7932. (37) Van Eps, N.; Oldham, W. M.; Hamm, H. E.; Hubbell, W. L. Structural and dynamical changes in an alpha-subunit of a heterotrimeric G protein along the activation pathway. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 16194−16199. (38) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (39) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668−1688. (40) Sorin, E. J.; Pande, V. S. Exploring the helix-coil transition via all-atom equilibrium ensemble simulations. Biophys. J. 2005, 88, 2472− 2493. (41) Meagher, K. L.; Redman, L. T.; Carlson, H. A. Development of polyphosphate parameters for use with the AMBER force field. J. Comput. Chem. 2003, 24, 1016−1025. (42) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (43) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (44) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (45) Amadei, A.; Linssen, A. B.; Berendsen, H. J. Essential dynamics of proteins. Proteins: Struct., Funct., Genet. 1993, 17, 412−425. (46) Vishveshwara, S.; Brinda, K. V.; Kannan, N. Protein structure: insights from graph theory. J. Theor. Comput. Chem. 2002, 01, 187− 211. (47) Vishveshwara, S.; Ghosh, A.; Hansia, P. Intra and intermolecular communications through protein structure network. Curr. Protein Pept. Sci. 2009, 10, 146−160. (48) Kannan, N.; Vishveshwara, S. Identification of side-chain clusters in protein structures by a graph spectral method. J. Mol. Biol. 1999, 292, 441−464.

N

DOI: 10.1021/acs.jctc.6b00813 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX