Edited by Michael L. Klein, Temple University, Philadelphia, PA, and approved March 28, 2012 (received for review December 14, 2011)
Protein allosteric pathways are investigated in the imidazole glycerol phosphate synthase heterodimer in an effort to elucidate how the effector (PRFAR, N′-[(5′-phosphoribulosyl)formimino]-5-aminoimidazole-4-carboxamide ribonucleotide) activates glutaminase catalysis at a distance of 25 Å from the glutamine-binding site. We apply solution NMR techniques and community analysis of dynamical networks, based on mutual information of correlated protein motions in the active and inactive enzymes. We find evidence that the allosteric pathways in the PRFAR bound enzyme involve conserved residues that correlate motion of the PRFAR binding loop to motion at the protein-protein interface, and ultimately at the glutaminase active site. The imidazole glycerol phosphate synthase bienzyme is an important branch point for the histidine and nucleotide biosynthetic pathways and represents a potential therapeutic target against microbes. The proposed allosteric mechanism and the underlying allosteric pathways provide fundamental insights for the design of new allosteric drugs and/or alternative herbicides.
Allostery is a fundamental property that allows for the regulation of function and dynamic adaptability of enzymes and proteins. Allosteric enzymes contain at least two distant binding sites, including the active site responsible for catalytic activity, which binds the substrate, and the allosteric site, which binds the effector and initiates the allosteric signal propagation to the active site. In V-type systems, substrate binding is not affected by the presence of the effector but if the effector is not bound, the allosteric protein is usually catalytically inactive (or poorly active), indicating that the effector binding is coupled to the kinetic and/or thermodynamic parameters of the biochemical reaction in the active site. Allosteric information transfer can range from large, enthalpically driven conformational changes to purely entropically driven motions or a combination of both enthalpic and entropic effects, but in each case the kinetic parameters of the catalyzed reaction at the substrate binding site are altered. At the heart of allosterism there is intramolecular thermodynamic coupling over long distances (> 10 Å), between the active and allosteric sites. An important challenge for fundamental studies is the elucidation of the allosteric pathways that connect the two ligand-binding sites.
In this work, we combine community network analysis based on molecular dynamics (MD) simulations and NMR studies of protein motion based on relaxation dispersion techniques and chemical shift titrations experiments to provide an atomistic description of allostery in the V-type allosteric enzyme imidazole glycerol phosphate synthase (IGPS) from the thermophile Thermotoga maritima (Fig. 1). IGPS is a tightly associated heterodimeric enzyme in which each monomer enzyme catalyzes a different reaction (1⇓–3). The 23 kDa HisH enzyme is a member of the glutamine amidotransferase family that catalyzes the hydrolysis of glutamine (Gln) to produce ammonia. The ammonia travels to the other monomer, HisF, passing through the (β/α)8 barrel and reacting with the effector PRFAR; i.e., N′-[(5′-phosphoribulosyl)formimino]-5-aminoimidazole-4-carboxamide ribonucleotide, to form imidazole glycerol phosphate (IGP) and 5-aminoimidazole-4-carboxamide ribotide (1, 4, 5). These two products enter the histidine and purine biosynthetic pathways, respectively. IGPS is not present in mammals, but it is involved in essential biosynthetic pathways of microorganisms. In particular, many plant pathogens, and importantly opportunistic human pathogens such as Cryptococcus, Candida, and Ajellomyces, have an IGPS that is highly homologous to the T. maritima enzyme. The allosteric IGPS complex is, thus, a potential target for antifungal, antibiotic, and herbicide development (2, 6, 7).
Structural model of IGPS heterodimer from T. maritima, including protein subunits HisH and HisF, glutamine substrate (green) and allosteric effector, PRFAR (red), and flexible loop1 in the cyclase domain.
IGPS is a V-type allosteric enzyme because the Gln binding affinity is unaffected by the presence of PRFAR. However, hydrolysis of Gln is accelerated 5,000-fold upon PRFAR binding (8). Here, we explore the effect of PRFAR on correlated motions of amino acid residues and we use this information to evaluate how the effector binding affects the protein network of the apo enzyme. Our methodology is inspired by studies of signaling pathways in tRNA:protein complexes (9). We combine solution NMR techniques to probe the motion of amino acid residues in the ms time scale and chemical shift titrations data with MD simulations for the analysis of mutual information based on Shannon’s entropy to quantify the extent of correlated motion of amino acid residues (10). The resulting analysis shows evidence that PRFAR binding alters the community structure of the inactive IGPS protein network, along with the collective pathways linking the effector and active sites, providing valuable insight on the mechanism of information transfer in this V-type allosteric enzyme.
The paper is organized as follows. We first analyze motions in the apo and the effector bound complexes, then compare correlated motions among amino acid residues by measuring the mutual information. We use this measurement to build dynamical protein networks and compare community structures of apo and PRFAR bound IGPS complexes. From the comparison of community structures and the analysis of information paths, IGPS allosteric pathways are obtained. We then propose an allosteric mechanism for the activations of the glutaminase reaction induced by PRFAR binding. Conclusions are summarized in the final section.
MD simulations show differences in motion of the flexible loop1 (HisF domain) when comparing apo and PRFAR bound IGPS complexes (Fig. 2). We observed that loop1 in HisF is more ordered and that the two-stranded β-sheet is further extended to the beginning of fα1 upon PRFAR binding, stabilized by hydrogen bonding between the highly conserved amino acid residues fV17 and fG30.
(A) Root mean square fluctuations (rmsf) of Cα atoms for apo (blue) and PRFAR bound (red) IGPS complexes, calculated with respect to the average structures of the four independent 100 ns MD simulations. (B) HisF average structures of apo (blue) and PRFAR bound (red) enzymes, showing differences induced by PRFAR binding in loop1 and fα1, fα2, and fα3 helices. PRFAR dynamics in the cyclase binding site (colored sticks from red to green for 0–100 ns) are shown. (C) and (D) Hydrogen-bonds in (HisF) loop1 in apo (C, in blue) and PRFAR bound (D, in red) enzymes at 50 ns, showing the extension of the small two-stranded β-sheet via the formation of the fV17-fG30 hydrogen-bond at the loop1-fα1 junction. (E) and (F) Breathing motion in apo (blue) and PRFAR bound (red) IGPS complexes. Instantaneous configuration at t = 10 ns (E), showing the larger angle θ between Cα atoms (F) of fF120, hW123, hG52 (red and blue balls) in PRFAR bound (red) with respect to the apo complex (blue). Thermal fluctuations of θ for representative trajectories of apo and PRFAR-bound complexes (F). Note that θ fluctuates back and forth between open and closed configurations 4 times/100 ns for PRFAR bound and 1/100ns for apo.
PRFAR binding also alters the breathing motion at the interdomain interface. Fig. 2 shows the low-frequency interdomain motions observed along representative 100 ns trajectories of both apo and PRFAR-bound IGPS complexes. This breathing motion allows for fluctuations between two states (open and closed IGPS heterodimer) that are accessible at thermal equilibrium in both the apo and PRFAR complexes. However, opening-closing events are observed more frequently when PRFAR is bound to the cyclase domain. The red trace in Fig. 2F shows that the angle θ [i.e., the angle between Cα atoms of fF120, hW123, hG52] opens and closes 4 times/100 ns when PRFAR is bound. Note the maximum at approximately 15 ns followed by a minimum at approximately 30 ns, then a maximum at approximately 40 ns followed by a minimum at approximately 50 ns, a maximum at approximately 60 ns followed by a minimum at approximately 70 ns, and then a maximum at approximately 85 ns followed by a minimum at 100 ns. In contrast, for apo-IGPS (blue trace) θ opens and closes only once/100 ns. Similar results are shown for other simulations in SI Text. These motions are observed in several independent 100 ns simulations performed, suggesting that PRFAR binding alters the protein-protein interactions at the HisF-HisH interface, consistent with the hinge-like displacement observed in the comparative analysis of crystal structures for the IGPS apo complex from T. maritima (5) and the PRFAR bound IGPS from yeast (11) (SI Text).
These observations corroborate the hypothesis that the motion of loop1 is coupled to the interdomain “breathing” motion, as previously indicated by steered MD simulations of the effector undocking (12). While a network of conserved residues, including hW123-fR249 (hinge) and hK181-fD98 (salt-bridge anchor), has been proposed to play a crucial role (12), the allosteric pathway has yet to be established. In the remainder of this section we address the allosteric mechanism in terms of the analysis of correlations in the motion of amino acid residues, with emphasis on the effect of PRFAR binding on those correlations. We find evidence of the allosteric pathway from the analysis of dynamical networks obtained by weighting the graphs with the mutual information that results from correlations of nuclear displacements in apo and PRFAR bound complexes.
The analysis of correlations of fluctuations in the positions of Cα-atoms is based on the generalized correlation coefficient rMI 
computed in terms of the mutual information (10), I[xi,xj] = H[xi] + H[xj] - H[xi,xj]. Here, H[xi], H[xj] and H[xi,xj] are the marginal and joint (Shannon) entropies (13) for atomic vector displacements (xi,xj), computed as ensemble averages over multiple 50 ns windows along eight independent 100 ns MD simulations for both apo and PRFAR-bound IGPS complexes (see Materials and Methods).
Fig. 3 shows our results for rMI of apo-IGPS and differences with respect to the corresponding values for the PRFAR bound complex. Large values (rMI > 0.8) are obtained for highly coupled motions in the apo enzyme, including residues with common secondary structure elements (near off-diagonal elements in Fig. 3A). Amino acid residues at the HisF-HisH interface also show relatively high correlations in the apo complex (rMI around 0.5–0.6), including residues in the hα1 helix and the Ω-loop (HisH loop1 connecting hβ1 and hα1) at one side of the complex (Fig. 3A, red box), and interface residues in the fβ4, fβ5 and hβ7 sheets and in the fα4 helix at the opposite side of IGPS (Fig. 3A, blue box). In contrast, amino acid residues in the range 201–253 (belonging to fα7, fβ8, fα8′-fα8 and the C-terminal loop of HisF) are poorly correlated with other residues in the HisH and HisF domains (Fig. 3A, dashed lines).
Comparison of generalized correlation coefficients rMI for apo (A) and PRFAR-minus-apo (B) IGPS complexes. Note low correlations for HisF amino acid residues 200–250 (A, bands within dotted lines), high correlation (rMI > 0.50) between HisF residues and HisH residues in hα1) and Ω-loop (A, red box) and interface residues in the fβ4, fβ5 and hβ7 sheets and in the fα4 helix (A, blue box) in apo IGPS. Decreased correlation for residues 100–220 in HisF (B, green features within dotted bands) and for residues in the 49-PGVG sequence (orange box) and around the active site (blue boxes) in HisH, increased correlation within residues 50–100 (B, magenta ellipse), between residues 50–100 and HisF loop1 (B, orange ellipse), between hβ1-hβ3 sheets, hα1 helix, Ω-loop, and HisF loop1 (B, blue ellipse), and between residues 222–230 (PRFAR gP-binding site) and residues in the sideR of HisF (B, black ellipse) upon PRFAR binding.
Comparison of correlated motions between apo and PRFAR bound enzymes reveals a reduction of correlations, induced by PRFAR binding, in agreement with the large, positive entropy change for PRFAR binding observed in isothermal titration calorimetry (ITC) experiments (14). The correlation matrix difference between PRFAR and apo complexes, shown in Fig. 3B, is indeed prevalently characterized by negative values (Fig. 3B, green spots). In particular, the largest reduction of correlations is observed for residues in the fα4′-fα7 helices and in the fβ5-fβ7 β-sheets on one side of the HisF domain (residues 100–220) (Fig. 3B, dashed lines). Note that the features highlighted in Fig. 3B result from thermal averaging, including all the possible combinations of PRFAR-minus-apo differences available from multiple MD simulations (SI Text). The green color code indicates features with rMI(PRFAR)‐rMI(apo)= -0.1, including a decrease of correlation for residues 100‐220 (Fig. 3B, green bands). These changes are 5‐10 times larger than the rmsds of the average matrix elements (SI Text). Magenta features, indicating an increase of correlation upon PRFAR binding, are also much larger than the rmsds.
The depletion of correlations indicates a certain degree of isolation of the decorrelated HisF amino acid residues from the rest of the HisF domain and from most of the residues in the HisH enzyme. The reduction of correlations in the PRFAR bound complex is also indicative of a weakening of interface interactions in this region upon effector binding. Changes in correlations at the interface are associated with the previously discussed changes in breathing motions, observed in MD simulations of apo and PRFAR bound complexes (Fig. 2). To understand the molecular origin of these changes in motion, we characterized the detailed PRFAR-HisF interactions and how these interactions affect the structure of correlations.
The pose of PRFAR in the cyclase active sites extends across the bottom of the (β/α)8 barrel, lying on a long cleft that separates the underside of HisF in two regions (Fig. 4).
PRFAR pose in the cyclase domain (HisF) and binding sites of the PRFAR phosphate groups. (A) Phosphate groups at the PRFAR ribose (rP) and glycerol (gP) sides are located between fβ3-fα3 and fβ4-fα4′ turns (green) and fβ7-fα7 and fβ8-fα8′ turns (green), respectively. (B) and (C) H-bonds (white dashed lines) at the interstitial site between fβ3-fα3 and fβ4-fα4′ turns (green), showing direct interactions between rP and fT104 (fβ4-fα4′ turn) and water-mediated interactions between rP and fG83 (fβ3-fα3 turn) in the PRFAR bound complex (B), water oxygen atoms shown as red balls. In the apo complex (C) only water molecules (blue balls) are present in the interstitial zone forming an hydrogen-bonding network. (D–F) H-bonds at the interstitial site between fβ7-fα7 and fβ8-fα8′ turns (green), showing hydrogen bonds between gP and fA224-S225 (fβ8-fα8′ turn) and water-mediated interactions between gP and fG203 (fβ7-fα7 turn) in the PRFAR bound complex (D). The extended hydrogen-bonding network in the apo complex at the fβ7-fα7 and fβ8-fα8′ site is shown in (E). (F) H-bonds formed by the hydroxyl groups of the PRFAR glycerol side (gOH) with amino acid residues in the fβ7-fα7 turn, including fS201, fG202 and fG203.
Phosphate groups at the glycerol side (gP) and at the ribose side (rP) of PRFAR are bound to the cyclase domain, occupying interstitial sites between the fβ7-fα7 and fβ8-fα8′ turns and between the fβ3-fα3 and fβ4-fα4′ turns, respectively. Hydrogen bonds between the rP group and both the side chain and the backbone of the highly conserved fT104 (in the fβ4-fα4′ turn) are the strongest interactions between PRFAR and HisF, as also suggested by undocking simulations of the effector (12). PRFAR-HisF interactions via the rP-fT104 hydrogen bonds induces desolvation in the interstitial region of the fβ3-fα3 and fβ4-fα4′ turns. Water-mediated hydrogen bond interactions present in the apo complex at this site are replaced by PRFAR-mediated interactions (Fig. 4). On the glycerol side of PRFAR, hydrogen bonds of the gP group with the backbone of fA224 and hydroxyl group of fS225 side chain are the first contacts to be broken during undocking of the effector (12), indicating a formation of weaker interactions with respect to interactions in the ribose side. Our MD simulations also indicate that hydroxyl groups of the PRFAR glycerol moiety create a network of hydrogen bonds with highly conserved fS201 and fG202 residues at the end of fβ7 (Fig. 4), contributing to the consistent number of protein-effector interactions in the interstitial region between the fβ7-fα7 and fβ8-fα8′ turns. These effector-protein contacts substitute for protein-solvent interactions that modulate interresidue communication in this region of the apo complex. The observed desolvation upon PRFAR binding is consistent with an entropically driven binding of the effector, as observed in the ITC experiments (14). Desolvation around fβ3-fα3, fβ4-fα4′, fβ7-fα7, and fβ8-fα8′ turns (residues 100–220) and formation of selective effector-protein interactions at residues fT104, fS201-fG202, and fA224-fS225, upon PRFAR binding, induce the observed decrease of correlations within amino acid residues in the range 100–220.
Notably, PRFAR binding at the bottom of the HisF barrel induces a reduction of correlations in the HisH domain, which is more than 20 Å away from the effector binding site. In particular, amino acid residues of the IGPS catalytic triad (hC84, hH178, hE180) (1, 4, 5) and their closest neighborhoods show lower rMI values in the PRFAR bound complex than in the apo complex (Fig. 2B, blue boxes). Motions of amino acid residues in the hP49-hG50-hV51-hG52 sequence (49-PGVG, also called the oxyanion strand) are also affected by PRFAR binding, with a significant reduction of correlations induced by effector binding (Fig. 3B, orange box). This result is consistent with the previously proposed allosteric mechanism invoking conformational changes of the 49-PGVG sequence induced by PRFAR binding (2, 14), and more generally, with the notion that effector binding triggers specific motions in the (distant) active site that favor catalytically active conformations.
In contrast to the decrease in correlations observed for some groups of amino acid residues upon PRFAR binding, increased correlations (Fig. 3B, magenta features) are observed for residues on the opposite side of the HisF protein, sideR (Fig. 1, Right). This side of the heterodimer includes the flexible loop1 at the bottom of the (β/α)8 barrel, the fα2-fα3 helices and the fβ2-fβ4 sheets in HisF, the hα1 helix and the Ω-loop at the HisF-HisH interface, and the hβ1-hβ3 sheets in HisH. Interestingly, an increment of correlated motions is observed between residues belonging to fα2-fα3 helices and fβ2-fβ4 sheets (Fig. 3B, magenta ellipse) and residues in the flexible loop1, including residues fN22 and fF23 (Fig. 3B, orange ellipse). Moreover, an increase of correlations due to effector binding is observed between HisF and HisH residues at the sideR (Fig. 3B, blue ellipse). In particular, motions of residues in the hβ1-hβ3 sheets in HisH are correlated with motions of the HisF loop1 through interface residues in the hα1 helix and the Ω-loop. Finally, larger rMI values between residues in the PRFAR gP-binding site (at the fβ8-fα8′ turn, residues 222–230) and residues in the sideR of HisF are obtained in the presence of PRFAR (Fig. 3B, black ellipse), indicating the presence of other correlated motions involving the sideR of the (β/α)8 barrel.
In summary, the mutual information analysis indicates that PRFAR binding induces changes in correlations of motions that are heterogeneously distributed within the heterodimer complex. On one side of the IGPS complex, sideR, the effector binding induces an increase of correlations for residues within the HisF subunit and between HisF and HisH residues. Because amino acid residues in sideR such as the molecular hinge and the salt bridge anchor and provide structural stability to the IGPS heterodimer, the increase of correlations in this region upon PRFAR binding indicates that effector binding would not destabilize the heterodimeric form of IGPS during Gln hydrolysis. On the other hand, a significant decrease of correlations is observed in the opposite side of the IGPS complex (hereafter named sideL) indicating that the entropically driven PRFAR binding promotes the separation of the HisF-HisH domains in the sideL. More importantly, PRFAR binding induces a significant reduction of correlations for motion of residues belonging to the Gln active site and adjacent 49-PGVG β-strand, more than 25 Å away from the effector-binding site.
The correlation analysis does show evidence of significant changes in the protein motion induced by effector binding. However, a more complete analysis of the structure of correlations is required to elucidate the allosteric pathways and residues that are crucial for signal propagation from the allosteric to the active sites. For this purpose, we use a network theory method that provides understanding of the underlying allosteric mechanism in the IGPS heterodimer.
Dynamical weighted networks for both apo and PRFAR bound complexes were built from MD calculations of generalized correlation coefficients (rMI) (see Materials and Methods). The weighted graphs are obtained using the average values of rMI calculated for all the MD simulations, as reported in Fig. 2, representing the average protein networks. Analysis of the time-dependent evolution of the weighted networks along 100 ns indicates that the time-averaged networks adequately describe the dynamical networks, capturing the protein conformational changes induced by effector binding during the 50 ns dynamics of the IGPS complexes. While rearrangement of residues in the communities and changes in the total number of communities are observed among different community structures, due to limitations of the modularity methodology (SI Text), many common features are present in different community networks and consequently in the average community structures. As previously reported for the analysis of allosteric pathways in tRNA:protein complexes (9), this methodology can provide valuable insight not readily available from the measures of correlation alone. Here, we sensibly improve the evaluation of the correlated motions by using the generalized correlation coefficient instead of the Pearson coefficient, overcoming the underestimation of correlation associated with the conventional covariance analysis (10). We focus on the generalized correlation coefficients based on mutual information because these are well defined measures of information content that can be exchanged between two amino acid residues. Moreover, the mutual information defined in terms of the Shannon entropy (10) allows for a direct estimation of the entropy redistribution within the protein complex, induced by correlations established upon effector binding.
Fig. 5 shows the optimum structure of communities obtained for apo-IGPS, including h1 (dark blue), h2 (green), h3 (cyan), and h4 (light blue) in HisH; f1 (orange), f2 (black), f3 (magenta), and f4 (gray) in HisF, and the network of the PRFAR bound complex, including h1′ (dark blue), h2′ (green), and h3′ (cyan) in HisH, and f1′ (orange), f2′ (black), f3′ (magenta), and f4′ (red) in HisF. A coarse-grained picture of the networks of communities of apo and PRFAR bound complexes is also shown in Fig. 5 (Bottom). The width of the bonds connecting communities is proportional to the sum of edge betweennesses of the edges connecting them (i.e., the intercommunities total betweenness, ITB). The optimum apo community structure has more communities than that of the PRFAR-bound complex (eight vs. seven communities, respectively). PRFAR binding induces changes in both the HisH and HisF domains with a reassignment of a few residues from h2 and h3 into the newly formed communities h2′ and h3′ and the merging of communities h1 and h4 into the new community h1′. In the HisF domain, communities f1 and f4 merge into f1′ and a new community f4′ is formed by some residue coming off from f3.
Color coded optimal community network of the apo (Left) and PRFAR bound (Right) IGPS complexes. Color Scheme: h1, h1′ (dark blue), h2, h2′ (green), h3, h3′ (cyan), h4 (light blue), f1, f1′ (orange), f2, f2′ (black), f3, f3′ (purple), f4 (gray), and f4′ (red), with predominant contributions from HisH (h) and HisF (f) subunits. The apo communities reported in brackets in the PRFAR network indicate those communities that contribute the most to the newly formed PRFAR communities. Widths of connectivity lines are proportional to the ITB with connections smaller than 10% of the maximum ITB value not shown, for clarity.
The analysis of optimal community structures typically reveals communities of nearby residues with common secondary structure elements although they might be distant in sequence, as shown in Fig. 5 (Top) and Table 1, where the distribution of secondary structure elements in the optimal community structures for apo and PRFAR bound complexes is reported.
Distribution of secondary structure elements in optimal community structures of apo and PRFAR bound complexes
In the apo complex, f3 is an important community located at sideR and is strongly connected to the other four communities in the IGPS heterodimer. Such a community includes residues at the interface belonging to different secondary structures of both HisF and HisH domains, including interface residues of the fα1, fα2, fα3, and hα4 helices and the fβ1, fβ2, fβ3, and hβ7. Interestingly, this community contains the conserved residues previously identified as the most stable contacts between the HisF and HisH domains (12), including the molecular hinge hW123-fR249 and the salt bridge anchor hK181-fD98 directly adjacent to the glutaminase active site. The grouping of these residues in a special “interface” community (f3) confirms the high level of correlation between elements in the IGPS sideR, providing structural stability and flexibility of the IGPS heterodimer.
The community f3 also contains the residues from the flexible loop1 and other residues at bottom of the HisF barrel. Upon effector binding, f3 is mostly conserved and forms the new community f3′. Interestingly, while the community is mostly preserved, we observed the removal of loop1, fα1 and residue hW123 and contemporary inclusion of some residues from hα1 in f3′. This redistribution of residues in f3' is indicative of a change in correlations at the molecular hinge (hW123-fR249) upon effector binding, consistent with the changes in breathing motion, and increased communication between the PRFAR binding site and the HisH active site. In the HisF domain, we observed the clustering of communities f1 and f4 forming community f1′ (residue 200–250) upon PRFAR binding, with partial reassignment of a few residues of f2 into f3′ at the sideL. The formation of community f1' and the rearrangement of residues in community f2' are due to the general lowering of correlation among residues of sideL (residues 100–200) in HisF upon effector binding (Fig. 3B) and to the contemporary increase of local correlations at the gP binding site in the region of the fβ7-fα7 and fβ8-fα8′ turns (Fig. 3B).
Our community analysis indicates that the glutaminase domain (HisH) of the IGPS apo complex is divided into four communities. Upon PRFAR binding, the HisH protein is instead divided into three communities, due to clustering of communities h1 and h4 into the new community h1′, which now includes residues hH178, hE180 from the catalytic triad. Community h3′ contains the 49-PGVG sequence and hC84 from the catalytic triad. Notably, effector binding in the HisF domain alters the community network of the HisH protein, affecting the correlations of residues in the catalytic active site, at more than 25 Å from the PRFAR binding site. Therefore, it is important to analyze the communication pathways between the two protein domains and the implications on the allosteric mechanism of IGPS induced by effector binding.
Fig. 1 shows the active and allosteric sites located approximately 25 Å apart in the HisH and HisF domains, respectively. Communication through these long distances is thus essential for the allosteric activity. The community network analysis shows that a large amount of the information flux must pass through the HisF-HisH interface. In particular, pathways connecting the HisF community f3 with HisH communities h2 and h4, and community f1 with community h1, involve most of the information flow in the IGPS apo complex (Fig. 5), while the f3′-h2′ channel is the most correlated interface route in the PRFAR bound complex.
In the apo complex, the information flow between HisF-f3 and HisH-h2 passes through two conserved interactions (12), including the hN15-fD98 and the hN15-hK181 hydrogen bonds that are part of the network of interdomain interactions adjacent to the glutaminase active site. These results are consistent with hN13-alanine mutagenesis studies in yeast (a site structurally equivalent to hN15 in T. maritima), showing the key role of this residue in providing structural rigidity and stability to the architecture of the IGPS heterodimer (15). Our community analysis consistently indicates that hN15 plays an important role in the communication network of the apo complex, providing stable contacts between the hα1 helix and the fα2-fβ3 and fβ3-fα2 turns at the sideR of IGPS. Interestingly, hN15 is connected to the rest of the h2 community via hN12 in the Ω-loop.
PRFAR binding induces an increase in the communication flow between communities f3′ and h2′ with respect to the flow between the apo communities f3 and h2. This increment in information flow is due to the enhanced correlations in the sideR of the PRFAR-bound complex, with increased communication between the helix hα1 and the Ω-loop in HisH and the HisF interface residues belonging to the fα2-fα3 helices. Interestingly, PRFAR binding induces significant changes in the connections within communities h2′ and h3′ in the glutaminase domain, with the hN12-hP10 connection in the Ω-loop having the highest betweenness in the new h2′ community, and connections between the Ω-loop and the Gln binding site (including the conserved 49-PGVG sequence) being loosened with respect to that in the apo community h2. These observations reflect a decrease in correlation, upon PRFAR binding, for the 49-PGVG sequence and the HisH active site (Fig. 2, blue and orange boxes) and the concomitant increase of correlations for residues at the sideR of HisF, including the fα2-fα3 helices (Fig. 2, magenta circle).
To analyze the interactions that generate the changes in the community network, we must focus on specific hydrophobic, ionic, and hydrogen-bonding interactions associated with the amino acid residues identified by the community network analysis. Such study indicates that the changes in the communication pathways induced by PRFAR binding are triggered by alteration of hydrophobic and ionic interactions at the sideR of the IGPS complex. Fig. 6 shows the time-dependent distances between hydrophobic residues at the bottom of the HisF (β/α)8 barrel. There are clear differences between the apo and PRFAR bound complexes. The residues primarily affected are those in the fβ2 strand that line the ammonia tunnel. The residues closest to the effector binding site, including fL50, fI52, and fF23, become more compact upon PRFAR binding, forming a hydrophobic cluster. The fL50-fI52-fF23 cluster formation is primarily driven by interactions with the π-system of the imidazolecarboxamide group of PRFAR (Fig. 6). The distance between fV48 and fL50 residues elongates considerably in the PRFAR bound complex while both fL50-fF23 and fL50-fI52 shorten upon PRFAR binding.
Changes in hydrophobic interactions induced by PRFAR in the HisF β∕α8 barrel. (A–D) show specific interaction distances as a function of simulation time for apo (blue) and PRFAR (red), for the hydrophobic residues shown in the bottom box. All highlighted residues but fF23 (in loop1) belong to the fβ2 strand, making a direct link between loop 1 and the ammonia tunnel. The representative MD snapshot (t = 40 ns) shown in the bottom box includes only HisF for simplicity. Interresidues distances (a, b, c, and d) are indicated in analogy to time-dependent distances reported in (A–D).
Our NMR titration experiments (Fig. 7) show prominent changes, induced by PRFAR binding, in the 1H and 13C chemical shifts for the methyl groups of key residues in the fβ2 strand. Fig. 7 shows that PRFAR induces a large change in both the chemical shifts and the line widths for the resonances of fL50 and fV48, indicating a significant change in the local chemical environment. This result is consistent with the conformational changes due to formation of the fL50-fI52-fF23 cluster. The disruption of the fV48-fL50 interaction on the ns timescale may represent the initial loosening of the protein core that precedes the wholesale enhancement of ms motions observed by solution NMR. Moreover, the alteration of hydrophobic interactions in the fβ2 strand, upon PRFAR binding, can also affect the opening of the ammonia tunnel gate favoring the NH3 transfer from the glutaminase subunit to the bottom of the (β/α)8 barrel. More relevantly, formation of the fL50-fI52-fF23 cluster in the fβ2 region also involves extension of the β-hairpin in loop1 and rearrangement of the adjacent secondary structures elements, such as the fα1 and the fα2 helices. Some of these rearrangements involve ionic interactions at the HisF/HisH interface, consistent with previous studies (12, 16).
The 1H-13C HSQC spectral overlays of the titration of PRFAR into IGPS (0.89 mM). Spectra were collected for 0 mM (red), 0.56 mM (yellow), 0.74 mM (green), 1.81 mM (blue), and 3.47 mM (purple). Large chemical shift perturbations and significant line broadening are observed for V48 and L50 of HisF.
Our MD simulations show that upon binding of PRFAR there is an alteration of the network of salt-bridges belonging to the hα1, fα2, and fα3 helices. As shown in Fig. 8, these α-helices are located in the sideR of the IGPS complex, with charged side chains placed at the protein surface and fα2 packing against the hydrophobic fβ2 region. As with hydrophobic interactions, these Coulombic contacts span over the PRFAR and Gln binding sites.
(Top) Changes in the network of salt-bridges at sideR of IGPS, induced by PRFAR binding. (Left) The fα2, fα3, and hα1 helices are highlighted in gray, with positively and negatively charged residues in green and yellow sticks, respectively. (Right) Schematic representations of ionic interactions in apo (blue) and PRFAR bound (red) complexes with the most stable ion-pair interactions highlighted with solid black lines. (Bottom) Thermal fluctuations of the Cζ-Cδ distance between the charged side chains of fR59 and fE91, respectively, along a representative 50 ns MD simulation of both apo (blue) and PRFAR bound (red) IGPS complexes.
MD simulations indicate that the network of salt bridges formed by residues fR59, fE67, fE71, fE91, fR95, and hR18 is significantly altered by PRFAR binding (Fig. 8). In particular, the fR59-fE91 salt bridge located at the bottom of HisF breaks in less than 20 ns in the apo complex while remaining stable for the entire simulation (50 ns) when PRFAR is bound (Fig. 8). Due to the weak fR59-fE91 interaction in the apo complex, fα2 and fα3 helices can be disrupted by thermal fluctuations, indicating a certain degree of flexibility of this region, adjacent to loop1. In the presence of PRFAR, instead, fα2 and fα3 helices are less flexible and their motion is more tightly coupled, as indicated by the increased correlations of these residues in the PRFAR bound complex (Fig. 3). Considering the proximity of loop1 to the fβ2 and that of fβ2 to fα2, the formation of the hydrophobic cluster fL50-fI52-fF23 and the extension of the small two-stranded β-sheet in loop1 (fV17-fG30 H-bond formation, SI Text) should favor a conformation of HisF with fα2 near fα3, promoting the stabilization of the fR59-fE91 salt bridge.
The alteration of the fR59-fE91 interaction, due to PRFAR binding, propagates through the adjacent fE91-fR95 and fR95-fE67 ion-pairs and in turn through the hR18-fE67 and the hR18-fE71 salt bridges at the HisH-HisF interface. Breaking of the fR59-fE91 salt bridge in the apo complex stabilizes interactions between charged side chains with alternating sign, composed by fE91-fR95-fE67-hR18-fE71 (Fig. 8), which remains stable for > 35 ns along the MD simulation. As shown in Fig. 8, this chain of ionic interactions is altered upon PRFAR binding, with a stable fR59-fE91 salt bridge limiting the interaction between fE91 and fR95 residues and stabilizing specific ion pairs such as fR95-fE67 and hR18-fE71. This rearrangement of ionic interactions upon PRFAR binding induces the increased communication flow between sideR communities f3′ and h2′ by coupling more tightly the motions of α-helices hα1, fα2, and fα3. Overall, the conformational rearrangement of helices hα1, fα2, fα3, and Ω-loop, upon effector binding, do not damage the structural stability of the IGPS heterodimer as indicated by preservation of clustering of interface residues in the f3 and f3′ communities. However, the reorganization of ionic interactions at the HisF-HisH interface is coupled to the breathing motion of the IGPS complex, with the PRFAR bound complex assuming a conformation that facilitates motions at the heterodimer interface. The breathing motion present in the PRFAR bound complex exposes more frequently the interface region to the solvent water molecules, which reduces the correlations in motions at the bottom of the HisH domain, including correlations between motions of the Ω-loop and the oxyanion strand at sideR, and motions between hβ7 and hβ8 β-sheets at the sideL of HisH. As we describe below, these results have significant implications on the allosteric mechanism of IGPS.
A recurring feature among enzymes with the catalytic triad motif Cys-His-Glu, such as guanine monphosphate synthetase (GMPS) (17), carbamoyl phosphate synthetase (CPS) (18), and anthranilate synthase (AS) (19), is the so-called “oxyanion strand,” a β-strand adjacent to the Gln binding site (the nucleophilic cysteine, hC84, of the catalytic triad). The importance of the oxyanion strand is due to its crucial role in forming an “oxyanion hole” that stabilizes a tetrahedral intermediate during the Gln hydrolysis, assisting the enzymatic activity. The oxyanion strand in the crystal structures of GMPS, CPS, and AS (17⇓–19) has been resolved in its “active conformation.” The orientation of the backbone of amino acid residues hG50 and hV51 (T. maritima IGPS numbering) leaves the NH group of hV51 pointing into the catalytic site, thus becoming part of the oxyanion hole. In IGPS from T. maritima, the oxyanion strand is constituted by the highly conserved sequence hP49-hG50-hV51-hG52 (49-PGVG) in the HisH subunit. The crystal structures of IGPS from T. maritima (5) indicates that the PGVG β-strand has an improper (inactive) conformation, with the carbonyl group of hG50 pointing into the Gln binding site (Fig. 9). Consequently, the NH group of hV51 is not able to form the hydrogen bond with the negatively charged oxygen of the tetrahedral intermediate and thus is not in optimal position for taking part in the oxyanion hole. This “inactive conformation” of the PGVG backbone has also been observed in the IGPS from Saccharomyces cerevisiae (2). These observations led to the hypothesis that a conformational change of the IGPS oxyanion strand, associated with a 180° rotation of the 49-PGVG backbone, is induced by effector binding (Fig. 9). This hypothesis is recently supported by our 1H, 15N NMR analysis of the 49-PGVG sequence, revealing conformational exchange motions of hG50 on the μs-ms timescale, upon titration with PRFAR (14).
(A) and (B) MD representative snapshots (at t = 50 ns) of apo (A) and PRFAR bound (B) IGPS complexes, showing H-bonds between hP10 (Ω-loop, in green) and hV51 (oxyanion strand, in blue). (C) and (D) Partial reorientation of the oxyanion strand triggered by PRFAR binding. Thermal fluctuations of the 49-PGVG backbone are shown (sticks colored from red to blue for 0–50 ns) highlighting the NH amide group of hV51 (larger sticks). Apo (C) and PRFAR bound (D) average structures (colored cartoons) and Gln binding site (hC84, colored sticks) are also shown. (E) Schematic representation of the oxyanion hole formation during the Gln hydrolysis reaction.
In the inactive conformation of the 49-PGVG backbone, the orientation of the NH group of hV51 is stabilized by a hydrogen bond with the carbonyl oxygen of hP10 in the Ω-loop. Disruption of this interaction is, therefore, a prerequisite for allowing the conformational change associated with the flip of the oxyanion strand backbone. All MD simulations indicate that this hydrogen bond is stable during the entire 100 ns trajectories (SI Text) in the apo IGPS complex. On the other hand, this important hydrogen bond breaks within 50 ns upon PRFAR binding, as indicated by all four independent MD simulations of the PRFAR-bound IGPS complex (SI Text). The hV51-hP10 hydrogen bond is responsible for connecting the Ω-loop with the HisH active site, including the 49-PGVG strand and the catalytic triad. The community network analysis reveals that breaking of the hV51-hP10 hydrogen bond is coupled with changes in the interface communications at the IGPS sideR and between the Ω-loop and the glutaminase active site. Fig. 9 A and B show the difference in the hV51-hP10 interaction as observed in a representative MD simulation of the apo and PRFAR bound IGPS complexes. The hV51-hP10 hydrogen bond begins to weaken after 12 ns and it completely breaks after 27 ns. Interestingly, we observed that the 49-PGVG strand partially rotates in the ns time scale upon breaking of the hV51-hP10 H-bond (Fig. 9 C and D). The displacement of the 49-PGVG sequence from the Ω loop confers more rotational degrees of freedom to the oxyanion strand, enabling it to rotate with respect to glutaminase binding site. The complete rotation of the oxyanion strand requires overcoming a large activation barrier > 20 kcal mol-1 on the ms time scale, as suggested by NMR solution experiments (14).
In summary, PRFAR binding modifies the protein network of IGPS, affecting the two sides of the heterodimer complex in two different ways. At sideL, entropically driven PRFAR binding decreases the degree of correlation in the motion of residues, inducing the splitting of the large HisF community f2 and favoring the transfer of entropy that increases the overall breathing motion of the IGPS complex. At sideR, PRFAR binding affects hydrophobic interactions at the bottom of the (β/α)8 barrel, inducing conformational changes that propagate through ionic interactions of α-helices fα2, fα3, and hα1 and, finally, affecting hydrogen-bonding interactions between the Ω-loop and the oxyanion strand adjacent to the HisH active site.
The allosteric pathway of the V-type IGPS heterodimer has been elucidated by combining NMR, MD simulations, and the correlation analysis of protein motions based network theory methods. PRFAR binding changes the structure of interactions that are essential for the allosteric mechanism. Dynamical rearrangements involve the flexible loop1 and the hydrophobic region (fβ2) at the HisF allosteric site, ionic interactions in the interdomain region next to the loop1 (sideR), including residues in the fα2 and the fα3 helices in HisF, and in the helix hα1 in HisH, as well as hydrogen bonding between the Ω-loop and the conserved 49-PGVG sequence (oxyanion strand) adjacent to the glutaminase active site. The entropically driven effector binding enhances the overall HisF-HisH interdomain motion, with the sideR portion of the protein complex remaining more rigidly closed than the sideL part. The resulting protein motions partially promote rotation of the conserved oxyanion strand responsible for the inactive-to-active allosteric transition. Single amino acid residues, such as fL50, fR59, fE91, and hP10, play an important role in the allosteric mechanism and should be good candidates for unique mutagenesis studies. Therefore, the proposed allosteric pathway provides useful information for the design of new molecules that might hamper glutaminase activity in plant or human pathogens by interrupting the allosteric signal propagation of IGPS.
The computational structural models for apo and PRFAR bound IGPS complexes are based on the crystal structure of the bienzyme complex from T. maritima at 2.4 Å resolution (Protein Data Bank code 1GPW) (1, 4, 5). The HisH-HisF apo-complex having the C-terminal loop of the HisF unit resolved in a closed (active) conformation (i.e., chains C and D) was chosen as starting structure. All crystal waters belonging to this protein complex were kept. The active-site mutation (N11D) present in the original crystal structure was mutated back to its wild-type form. The positions of the phosphate groups in the T. maritima IGPS crystal structure and in the PRFAR bound complex from yeast (PDB code 1OX5) (11) match nicely after superposition of the two X-ray structures. Therefore, the initial pose of the PRFAR effector in the HisF cavity was chosen according to its position in the yeast structure (SI Text). MD simulations of the apo and PRFAR bound structures of IGPS are based on the AMBER-ff99SB (20) force field for the protein and Generalized Amber Force Field (21) for the PRFAR ligand (see SI Text), as implemented in the NAMD2 software package (22). Four independent MD simulations have been performed for each complex (apo and PRFAR-bound) for a total simulation time of 0.8 μs. Further details of preequilibration procedure and MD production runs are described in SI Text.
The cross-correlation analysis is an established method used for quantifying correlated motion in proteins from MD simulations (23, 24). According to this approach, correlations between residues motions are obtained from MD trajectories by calculating the Pearson correlation coefficient 
[1]where xi,j are the fluctuation vectors of atoms i and j, defined as the deviation from the mean; i.e., 
, of the atomic position vector v, with 
denoting the ensemble average. Estimate of correlations from the Pearson coefficient is not independent from the relative orientation of the atomic fluctuations and is limited to the detection of linear correlations only, with relevant artifacts in the quantification of motion correlations of amino acid residues (10). Lange, et al. (10) recently proposed a method to measure generalized correlations in biological systems, which overcomes the problems related to the use of the Person coefficient, as coefficient of nondetermination. The joint probability distribution p(x) of (N) random variables xi, with i = 1, ..,N, equals the product of their marginal distribution, pi(xi), if and only if the components xi are independent; i.e., uncorrelated. Therefore, the ensemble-averaged deviation from the uncorrelated distribution is given by the mutual information (MI), defined as 
[2]MI vanishes in the case of independent random variables; i.e., it vanishes in the bivariate case of the fluctuation vectors xi and xj for fully uncorrelated motions of residues i and j. The above defined MI yields values in the range [0,..,∞) leading to a less intuitive interpretation of the correlations with respect to the Pearson coefficient, which values vary from -1 to 1. For an easier interpretation of the MI, it is thus convenient to define a generalized correlation coefficient, rMI, as the Pearson-like coefficient derived from the mutual information. Considering colinear Gaussian distributions of unit variance (d = 3), the generalized correlation coefficient rMI is defined as 
[3]The coefficient rMI is zero for fully uncorrelated variables and assumes values up to 1 for fully correlated variables.
MI is closely related to the Shannon entropy (13), H[x], which is defined as the expected value of the information content of a discrete random variable x, having probability mass function p(x), 
[4]The mutual information in the case of bivariate random variables can be calculated as 
[5]where the marginal, H[xi] and H[xj], and the joint, H[xi,xj], entropies could be estimated with different methods (25).
In this work, the information content H[x] is calculated using the atomic positions fluctuations derived from the MD simulations by means of the k-nearest neighbor distances algorithm (25) as implemented in GROMACS (10, 26).
The generalized correlation coefficient, rMI, for each pair of amino acid residues of the IGPS complexes are, then, obtained according to Eq. 5 considering the motion of Cα atoms, averaging for each independent 100 ns MD simulation over six 50 ns windows, and then averaging over four independent MD simulations for both apo and PRFAR-bound complexes (for a total simulation time of 0.8 μs).
The IGPS network is defined as the set of amino acid residues; i.e., the nodes (or vertices), connected by edges whose length is weighted using the generalized correlation coefficient rMI. For a fair comparison between the apo and the PRFAR bound IGPS networks, PRFAR is not considered to be a node, with apo and PRFAR protein networks having the same number of nodes. Two nodes are considered connected if any heavy atoms from the two residues are within 5.0 Å of each other (distance cutoff) for at least 75% of the frames (percentage cutoff) analyzed (SI Text).
The edge lengths in the network are obtained using the generalized correlation coefficients associated to the mutual information of each pair of residues. We first obtained the matrix of direct communication ‘distances’ 
between all pairs of residues (ij). The matrix of shortest distances 
was determined by using the Floyd-Warshall algorithm (27), considering direct distances 
as well as up to N possible intermediate residues mediating indirect communication pathways (where N is the number of residues in the system). The matrix of shortest pathways 
provides the edge-betweenness matrix, with elements b(ij) defined as the number of shortest paths that include pair (ij) as one of its communication segments. Therefore, the edge betweenness matrix b(ij) is a measure of ‘traffic’ passing through edges within the network and can be used for partitioning the communication network into communities. We use the Girvan-Newmann algorithm (27, 28), maximizing the modularity measure (Q) (29) to optimize the quality of the community structure (SI Text).
HisF and HisH (with a C-terminal His6 tag) codon optimized (Genscript) T. maritima genes for Escherichia coli were individual cloned into a pET43.1b vector between the NdeI-HindIII and NdeI-Xhol sites, respectively. 2H, 15N, and 13C methyl groups of Ile, Leu, and Val residues (ILV) HisF-2H HisH isotope enriched samples were prepared by transforming the HisF and HisH plasmids into separate E. coli BL21(DE3) cells, which were grown in 100% D2O M9 minimal media with 
as the nitrogen source and the isoleucine, leucine, and valine precursors 2-ketobutyric acid-4-13C (Isotec) and 2-keto-3-(methyl-D3)-butyric acid-4-13C (Isotec) for HisF. The cells were grown to an OD600 of 0.9–1.0 at 37 °C and were induced with a final concentration of 1 mM IPTG. To achieve 
isoleucine and leucine as well as 
valine, 60 mg 2-ketobutyric acid-4-13C and 100 mg 2-keto-3-(methyl-D3)-butyric acid-4-13C per liter of M9 were added 1 h prior to induction (30). Carbon-3 bound hydrogens for 2-ketobutyric acid-4-13C and 2-keto-3-(methyl-D3)-butyric acid-4-13C were exchanged with deuterium by dissolving the compounds in 30 mL and 20 mL of D2O, respectively, with one KOH pellet. The high pH mixture was then incubated overnight at 37 °C to facilitate exchange. Before addition to the culture, the pH was lowered to 7.0. The cells were harvested after 8 h of induction. HisF and HisH cells were resuspended in 5 mL 100 mM potassium phosphate, 100 mM Tris-HCl (pH 9.0), and 1.0 M NaCl per gram of cell. The resuspended cells were colysed by sonication with a 1.5 times growth culture volume excess of deuterated HisH. Cell waste was removed though centrifugation and the supernatant was mixed with Ni/NTA resin. The Ni/NTA resin was washed with 100 mM potassium phosphate (pH 7.5), 10 mM Tris-Base (pH 9.5), and 10 mM imidazole. The HisF-HisH complex was eluted in 10 mM N-cyclohexyl-3-aminopropanesulfonic acid (CAPS), 10 mM Tris-Base (pH 9.5) and 250 mM imidazole and dialyzed against 4.0 L of 10 mM Hepes (pH 7.3), 10 mM KCl, 1.0 mM EDTA and 1.0 mM DTT. PRFAR was dissolved in 10 mM Hepes (pH 7.3), 10 mM KCl, 1.0 mM EDTA, 1.0 mM DTT and the pH was adjusted to 7.3. 1H-13C spectra were collected for 0, 0.19, 0.56, 0.74, 1.81, and 3.47 mM PRFAR with a starting IGPS concentration of 0.89 mM.
We thank J. Lipchock for acquisition of the NMR titration data, and O. Lange and V. Limongelli for very useful discussions. V.S.B. acknowledges supercomputer time from the National Energy Research Scientific Computing Center (NERSC) and support from the National Institutes of Health Grant 1R01-GM-084267-01 for the developments of methods and models implemented in this study. G.M. acknowledges support from National Institutes of Health training Grant T32GM008283.
Author contributions: I.R., J.P.L., and V.S.B. designed research; I.R., M.M.S., N.-S.L., and G.A.M. performed research; I.R., M.M.S., N.-S.L., G.A.M., J.P.L., and V.S.B. analyzed data; and I.R., J.P.L., and V.S.B. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
See Author Summary on page 8366 (volume 109, number 22).
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1120536109/-/DCSupplemental.
