Conformational Changes of Glutamine 5′-Phosphoribosylpyrophosphate Amidotransferase for Two Substrates Analogue Binding: Insight from Conventional Molecular Dynamics and Accelerated Molecular Dynamics Simulations

Introduction

GPATase catalyzes the first step of de novo purine nucleotide synthesis and takes the sum of two ‘half reactions’ at different active sites in separate domains (i.e., the Gln and the PRTase domains) (Krahn et al., 1997; Wolan et al., 2004; Wang et al., 2009; Fan et al., 2019; Cao et al., 2020). It is a key enzyme in the total synthetic pathway, and it transfers nitrogen from glutamine to phosphoribosylpyrophosphate (PRPP) (Chittur et al., 2001; Zhu et al., 2017; Semmelmann et al., 2019). This enzyme consists of both a glutaminase domain and a synthase domain. The glutamine domain in its N-terminal is an Ntn hydrolase in conserved enzyme groups (Brannigan et al., 1995; Teplyakov et al., 1999). Its C-terminal domain is a phosphoribosyltransferase (PRTase) (Musick and Nyhan, 1981) (Figure 1).

FIGURE 1

FIGURE 1. The overview of GPATase structure and the ligands. (A) The stereochemical structure of DON and cPRPP; (B) GPATase subunit structure and topological structure of the protein in the glutamine domain and PRTase domain made by Pro-origami.

GPATase is a type I PRTase (Grubmeyer et al., 2012). Several crystal structures of type I PRTases (Krahn et al., 1997; Singh and Christendat, 2006; González-Segura et al., 2007; Fan et al., 2009; Xiang et al., 2009; Lund et al., 2010; Kang et al., 2012) have been determined. Ammonia is transferred by the ammonia channel, and this process depends on the active sites of the two domains. However, two questions must be answered. How do the two domains mutually effect each other? What causes conformational changes in GPATase? Recently, enhanced sampling molecular simulations have been used to successfully probe the effect of multisite PTMs on the mobility of several regions for HMGB1-PtDNA complexes (Lyu et al., 2018). Consequently, cMD simulations and aMD simulations were performed to elucidate the changes in protein structures that occur when substrate analogues bind to GPATase. The systems and all molecular dynamic simulations performed in this work are summarized in Table 1.

TABLE 1

TABLE 1. Summary of molecular dynamics simulations performed in this work.

cMD simulations as an efficient tool are increasingly used to probe the structural biology of biomolecules (Wang et al., 2018; Zhu et al., 2020). However, cMD simulations have difficulty capturing numerous conformational changes because of several limitations. Accelerated molecular dynamics (aMD) (Hamelberg et al., 2007; Kappel et al., 2015) applies a bias potential so that researchers can not only monitor conformational changes more efficiently than cMD but also use fewer computing resources (Tsai et al., 1999; Gasper et al., 2012; Spindel, 2012). More importantly, this technology has achieved success in probing proteins conformational changes induced by ligand bindings (Wang et al., 2016; Kokh and Amaral, 2018; Peng et al., 2018; Wang and Miao, 2019). For example, Grant et al. employed aMD simulations to probe the GTP and GDP modified structure changes of Ras proteins, their findings suggest that aMD can capture nucleotide dependent conformational switching but cMD simulations cannot, showing the advantages of aMD simulations. In this study, we first performed cMD to observe the effect of two substrate analogues (i.e., DON and cPRPP) on the structural stability of GPATase. Despite the impact of different ligand binding modes on the conformational change of two active sites being studied by cMD simulations, it is hard to capture how the conformational transition process induced by two domains mutually effect each other. Therefore, to realize our aims, performing aMD simulations could provide detailed pictures of the conformational changes in the Glutamine loop and PRTase flexible loop of the enzymes. Moreover, it can describe the mechanism of the synergistic actions of these two domains.

Materials and Methods

System Setup

The coordinates of E. coli GPATase (PDB code 1ECC) (Krahn et al., 1997) complexed with cPRPP and DON were downloaded from the Protein Data Bank (Figure 1). As the initiating structure, cPRPP was used as a nonhydrolyzable analogue of PRPP substrate. In this crystalline structure, the side chain of Cys1 is covalently modified by DON, which is a glutamine affinity analog (Krahn et al., 1997; Parry et al., 1997). Four models were constructed: 1) GPATase without ligands, 2) only DON remained in the GPATase (Mn²⁺ was coordinated by O@water in Model 1 and Model 2), 3) only cPRPP remained in the GPATase (Mn²⁺ was ligated by O1, O2, O3, O1β@cPRPP, and O@water), and 4) GPATase with two ligands (this model is a combination of Model 2 and Model 3). MCPB. py (Li and Merz, 2016) (version 3.0, released in AmberTools 17) was utilized to build force fields based on quantum calculations and simulate metal ions containing coordination compounds via the bonded model approach (Cornell et al., 1995). The Mn²⁺ parameters used in this program were from Li et al. (Li et al., 2013). All models were prepared using the Leap module from the AMBER16 suite and were characterized by employing the ff14SB force field (Maier et al., 2015). All histidine residues were protonated at the ε position. All missing hydrogen atoms were added to cPRPP using the Leap program from the AMBER16 package (Götz et al., 2012). The geometries of cPRPP were refined with Gaussian16 at the B3LYP/6-31G* level (Slanina et al., 2004). The optimized geometries were used to calculate charges derived from electrostatic potential derived (ESP) following the RESP methodology (Cornell et al., 1993) implemented in the antechamber module in Amber16. The neutral charge of the systems was maintained by adding 9, 11, 12, and 15 Na⁺ ions respectively to Model 1, Model 2, Model 3, and Model 4. Each model was solvated in a cubic box with a buffer of no less than 12 Å and used TIP3P water molecules to generate water model parameters (Paschek et al., 2011).

Molecular Dynamics Simulations

Four cMD simulations were performed using the pmemd. cuda module of AMBER16 (Götz et al., 2012). Electrostatic interactions were calculated via the smooth particle mesh Ewald method, and the cut-off distance of non-bonded real space interactions was 10 Å (Darden et al., 1993). The bonds containing hydrogen atoms were restrained, using the SHAKE algorithm, and the integration time step for four simulations was 2 ps (Ryckaert et al., 1977). Inappropriate geometry and spatial conflicts were avoided by minimizing the energy before the cMD simulations via the steepest descent algorithm and conjugate gradient minimizer method (Jie and Liu, 2003). The cMD production for all four systems were ran under constant temperature (310 K) and 1 atm constant pressure conditions (Berendsen et al., 1998). Subsequently, 200 ns cMD simulations were performed for each system with periodic boundary conditions. Saving atomic coordinates every 100 ps for further analysis.

Accelerated Molecular Dynamics Simulations

aMD is a means of enhancing sampling technology by bridging the potential energy surface. It can overcome the drawbacks of cMD, which has difficulty monitoring domain structure variation and long-distance correlated motions that occur over very long timescales (Hamelberg et al., 2004; Hamelberg et al., 2007; Markwick and Mccammon, 2011). In this study, four dual-boost aMD simulations of GPATase were initiated from the final structure of the corresponding cMD simulations. Before running dual-boost aMD, the following aMD parameters were calculated: 1) EthreshP ($E_{tot}$), which denotes the average total potential energy threshold; 2) alphaP ($Alph{a}_{tot}$), which is the inverse strength boost factor for the total potential energy; 3) EthreshD ($E_{dih}$), which represents the average dihedral energy threshold; and 4) alphaD ($Alph{a}_{dih}$), which is the inverse strength boost factor for the dihedral energy. These parameters were calculated using the following formula (Pierce et al., 2012):

$E_{tot}=E_{tot avg}+N_{atoms}\times 0.16(1)$

$Alph{a}_{tot}=N_{atoms}\times 0.16(2)$

$E_{dih}=E_{dih avg}+N_{resid}\times 4(3)$

$Alph{a}_{dih}=\frac{1}{5}\times \left(4\times N_{resid}\right)(4)$

where $E_{tot avg}$ and $E_{dih avg}$ are obtained from cMD simulations. $N_{atoms}$ is the number of the whole system atoms, and $N_{resid}$ is the total number of residues. For each system, 400 ns aMD simulations were performed and aMD parameters settings are based on cMD simulation results.

Subnetwork Analysis of Protein-Ligand Complexes

Protein structure networks offer invaluable insights into the structural properties, functions, and stability of proteins (Atilgan et al., 2004; Brinda and Vishveshwara, 2005). In this study, RINerator software was utilized to obtain the RINs of GPATase, GPATase–DON, GPATase–cPRPP, and GPATase–DON–cPRPP. RINalyzer was used to perform the RIN analysis (Doncheva et al., 2011), The results were then visualized using Cytoscape (Shannon et al., 2003). Three-dimensional structures were visualized using the PyMOL software (Delano, 2002) and VMD software (Humphrey et al., 1996). Cluster analysis was performed using the K-means clustering algorithm available in CPPTRAJ (Shao et al., 2007). Detailed clustering results for the top six clusters were summarized in Supplementary Tables S1–S4. The centroid structure of the most populated cluster that had the lowest cumulative distance to every other point was selected as the representative structure. The representative structures obtained from cluster analysis of the four complexes were applied to conduct protein-ligand network analysis.

Cross-Correlation Analysis and Principal Component Analysis

Dynamic cross correlation analysis (DCCA) was conducted using the Bio3D R package version 2.3.0 (Hünenberger et al., 1995; Grant et al., 2006; Skjærven et al., 2014). The scales of the coordinate axis represent the atomic number in the covariance matrix map, whereas the gradation of color displays the extent of atomic motion correlation. The red regions denote the positive correlation between atoms reflecting motion along the same direction, whereas the blue regions show that the corresponding two atomic motion modes are negatively correlated. PCA was performed using Bio3D (version 2.3.0) and visualized by VMD (version 1.9.3) (Grant et al., 2006; Skjærven et al., 2014). The energy of protein conformations for the first two principal components (PC1 and PC2) were characterized in Free Energy Landscape (Hünenberger et al., 1995). The trajectories to PC1 and PC2 of motion were projected by employing the program converting dot distribution to probability distribution (ddtpd) (Nicolaï et al., 2013; Iida et al., 2016).

Results and Discussion

Structural Stability Analysis during cMD Simulations

The structural stability of the four complexes was analyzed. The root-mean-square deviation (RMSD) of the protein Cα atom backbone was calculated to illustrate conformational changes. The radius of gyration (Rg) can reflect the plasticity potential of protein structures. The free energy landscapes for both the Rg and RMSD of the four systems are shown in Figures 2A–D. The DON–cPRPP–GPATase system exhibited the main energy basin located at the RMSD value of ∼2.2 Å and the R_g value of ∼23.4 Å, both of which were smaller than those of the three other systems. The corresponding probability of RMSD and R_g of the four systems is displayed in Figures 2E,F. The four systems achieved convergence during the cMD simulations (Supplementary Figure S1A). The RMSDs of the GPATase without ligand and the GPATase with two substrate analogue bonds reached equilibrium at 2.6 Å after 130 ns. By comparison, the RMSDs of the GPATase with DON reached equilibrium at 3.2 Å after 175 ns. The Rg value of the GPATase-DON-cPRPP system stabilized at about 23.8 Å, whereas the R_g value of the enzymes with DON or cPRPP was larger, indicating that the GPATase–DON–cPRPP structure was more compact than others (Supplementary Figure S1B). These results demonstrated that the GPATase with two substrate analogues was more stable than those with DON or cPRPP only.

FIGURE 2

FIGURE 2. The free energy landscape for Rg and RMSD of (A) GPATase, (B) GPATase–DON, (C) GPATase–cPRPP, and (D) GPATase–DON–cPRPP. (E) Relative frequency of RMSD of the four systems, (F) relative frequency of Rg of the four systems.

Root-mean-square fluctuation (RMSF) was used to compute the fluctuations of a single residue. RMSF values can display the extent of freedom for atomic motion, indicating the flexibility of a protein region. RMSF values were calculated at equilibration state. The glutamine active site domain (R73-S79) and the PRTase active site domain (K326-L350) are highlighted in Figures 3A,B. The RMSF values fluctuated, and typical conformations are also displayed in the cartoon. The corresponding RMSD values during the 200 ns cMD simulation of R73-S79 and K326-L350 are shown in Figures 4A,B. The solvent accessible surface area (SASA) contribution of the active binding pocket (i.e., the glutamine and the PRTase active sites) of the GPATase–DON–cPRPP complex was the greatest among the four systems, indicating that this complex could provide a decent hydrophilic environment for ligand binding (Figures 4C,D). Changes in the secondary structures of the complexes, primarily those in Q339 to K349 (which are a part of the PRTase active site domain), are plotted in Figure 5A. The α-helix (Q339 to K349) of the GPATase bound to DON and cPRPP was almost contained, whereas that of the GPATase with cPRPP completely disappeared (Figure 5B). In the two other systems, the α-helix partly disappeared. The corresponding cartoon structures during the 0 and 100 ns cMD simulations of the four systems are exhibited in Figure 5C. Results of cMD analysis revealed mutual stabilization between DON and cPRPP, and they exerted a synergistic effect on the stability of the enzymes.

FIGURE 3

FIGURE 3. RMSF values of the backbone atoms of the four systems over time during the simulations. Comparison of the RMSF plots of proteins in (A) GPATase, GPATase–DON, GPATase–DON- cPRPP complexes; (B) GPATase, GPATase–cPRPP, GPATase–DON- cPRPP complexes.

FIGURE 4

FIGURE 4. Stability analysis of DON binding site and cPRPP binding site for systems during the 200 ns simulations. (A) The RMSD values of DON binding site, (B) RMSD plot of cPRPP binding site of the four systems during 200 cMD simulations. (C) SASA plot for DON binding pocket. (D) SASA plot for cPRPP binding pocket.

FIGURE 5

FIGURE 5. Analysis of changes in secondary structures. (A) DSSP results of the four systems, the color bar represented different secondary structures. (B) The probability of α-helix (Gln339-Lys349) during 200 cMD simulations. (C) The corresponding cartoon structure of the 0 and 100 ns cMD simulations for the four systems.

Cross Correlations and PCA Calculations

DCCA values were computed using cMD trajectories at the equilibrium state to explore the internal dynamics of GPATase. Cross-correlation maps of correlated intermolecular motions between the remote regions of the proteins in different complexes are presented in Figures 6A–D. GPATase with DON or cPRPP displayed larger fluctuations than the two other systems. The GPATase–DON–cPRPP complex manifested the lightest color in the cross-correlation matrix maps, indicating that it experienced the weakest flexibility and it was the most stable during the simulations. Furthermore, the correlation analysis revealed negative correlation between the RTase flexible loop (K326 to L350) and the glutamine loop (R73 to E84) (the relevant regions are highlighted by black boxes). The binding of DON and cPRPP to the GPATase–DON–cPRPP complex weakened the negative correlation between the two loop domains.

FIGURE 6

FIGURE 6. The residues cross-correlation maps for four systems (A) GPATase, (B) GPATase–DON, (C) GPATase–cPRPP, and (D) GPATase–DON- cPRPP. The positions of glutamine loop and flexible loop were labeled by black boxes.

The dominant conformational changes were identified via PCA analysis of the Ca atoms in the cMD trajectories. The overall motion was determined via PCA through the eigenvectors of the covariance matrix. The two-dimensional projections of the total cMD conformational space on the first two principal components, namely, PC1 and PC2, are given in Figures 7A–D. The sum percentage of the first two principal components (PC1 and PC2) in the four systems accounted for about 50%, which ensured the dependability of the observed motions (Table 2). The conformations of the lowest and the sub-lowest energies shown in Figure 7 indicated that the two active sites had different movement trends. The low-energy regions of the GPATase–DON–cPRPP complex were larger and more centralized than those of the GPATase complexes with only DON or cPRPP. This phenomenon also explained why the latter reached equilibrium after a long time during the MD dynamic simulations. This result was consistent with that of RMSD analysis. These observations implied that the presence of different ligands produced varying effects on the dynamic behavior of GPATase.

FIGURE 7

FIGURE 7. The two-dimensional projection of total cMD conformational space on the first two principal components (PC1 and PC2). FEL maps and PC1 and PC2 structures for (A) GPATase, (B) GPATase–DON, (C) GPATase–cPRPP, and (D) GPATase–DON–cPRPP systems. The conformation of GPATase proteins with one global minimum is marked in red. The depth of the energy landscape indicates the value of the minimum free energy.

TABLE 2

TABLE 2. Principle component probability during cMD simulations.

Subnetwork Analysis of Protein-Ligand Interaction

The subnetworks of the interfaces between the ligands and the proteins were extracted to understand differences in the binding affinity of DON and cPRPP to the complexes (Figures 8A–D). The interfaces suggested that the synergism of DON and cPRPP may strengthen the interactions between the ligands and the protein. Thus, the number of interactions in the GPATase–DON–cPRPP complex was greater than that in the DON or cPRPP binding complexes. This result was further verified by analyzing the hydrogen bonds between the ligands (DON or cPRPP) and GPATase in the three complexes during the cMD simulations. The probabilities of hydrogen bond formation are listed in Table 3. The probability of hydrogen bond formation was higher between the DON/cPRPP and GPATase in the DON–cPRPP–GPATase complex than that in the DON- or cPRPP-bound complexes. These results indicated that the combination of cPRPP and DON may induce more stable interaction between the ligands and GPATase, as well as an inseparable communication between the nodes in the glutamine and the PRTase active sites.

FIGURE 8

FIGURE 8. The subnetwork interaction analysis between ligands and protein. The subnetwork between GPATase and DON in the (A) GPATase–DON complex and (B) GPATase–DON- cPRPP complex. The subnetwork between GPATase and cPRPP in the (C) GPATase–cPRPP complex and (D) GPATase–DON- cPRPP complexes.

TABLE 3

TABLE 3. Probability of hydrogen bond formation between protein (GPATase) and ligands (DON and cPRPP) for the GPATase–DON, GPATase–cPRPP, and GPATase–DON-cPRPP structures during the 200 ns cMD simulations.

The previous studies indicated that NH₃ can transfer between two active sites without external potential in GPATase by performing LES/PMF simulations (Xiang et al., 2009), so two domains mutually effecting each other induced by substrate binding was likely to be the main factor to facilitate ammonia travel. We performed the above cMD simulations to further decode atomic level mechanism relating to structure changes of GPATase. The cMD analysis results provide a basis for further performing aMD simulations.

The Stability Analysis of R73-DON Salt Bridge during aMD Simulations

The conformational states and structural dynamics of GPATase in each system were evaluated via 400 ns aMD simulations. The RMSD values were clearly convergent for the aMD simulations (Supplementary Figure S2). R73 was optimally positioned for binding with glutamine, and the PRTase sites were connected by a 20 Å NH3 channel. R73 formed a salt bridge with the carboxyl group of glutamine analogues, which was improperly positioned for this interaction (Bera et al., 1999). A salt bridge formed between the R73 and DON occurred only in the active conformer of the enzyme. The representative conformations of the GPATase–DON and GPATase–DON–cPRPP are shown in Figures 9A,B, respectively. The ability of cPRPP to bind or not to bind to GPATase and its influence on the formation of salt bridge were explored. The distance between R73 and DON during the 400 ns aMD simulations is given in Figure 9C. A distance of less than 4 Å denotes the formation of a salt bridge (Kumar and Nussinov, 2002). The distance had larger fluctuations for GPATase bound to DON only than the GPATase–DON–cPRPP complex (Figure 9C). Moreover, the peak value for the GPATase–DON complex was 3.2 Å, whereas that for the GPATase–DON–cPRPP complex was 4.8 Å (Figure 9D). According to the results of the subnetwork analysis summarized in Figure 8 and Table 3, the probability of hydrogen bond formation between R73 and DON was higher for the GPATase–DON–cPRPP complex than for the complex without cPRPP. Hence, the hydrogen bond interactions between them moved the side chain of R73 closer to the DON carboxylate group. This movement may be useful to the formation of a salt bridge. Furthermore, salt bridges between R73 and DON carboxylate formed only in the active form of the enzyme, indicating that binding cPRPP to the PRTase flexible loop was important for reconstructing the glutamine loop and inducing GPATase in its closed form (i.e., active state).

FIGURE 9

FIGURE 9. Comparison of R73–DON salt bridges between the GPATase–DON and GPATase–DON–cPRPP complexes during aMD simulations. The salt bridges are shown as green dashed lines. The glutamine loop is marked in light blue. R73 and DON are presented as sticks. Representative structures of the (A) GPATase–DON and (B) GPATase–DON–cPRPP complexes. (C) Variations in the distance between R73 and DON of the two systems. (D) Relative frequency of distance of the two systems.

Conformational Free Energy Landscape of PRTase Active Site Domain

In each system, 400 ns aMD simulations were carried out to examine the conformational changes and atomic-level molecular mechanism of GPATase. Free energy profile maps of GPATase were explored according to the extension and changes in the shape of the PRTase flexible loop, both of which were important for GPATase activation. The free energy distribution map of each system was projected along with two collective variables (abscissa and ordinate). CV1 reflected the deviation in the PRTase flexible loop from the crystallographic structure (PDB code 1ECC) (Krahn et al., 1997), which was able to describe the opening-closing state of the PRTase flexible loop. The residues moving outward formed an open conformation of the PRTase flexible loop, which was in its inactive state because the binding cavity in the shape was not complementary. The negative value represents the closed conformation status of the flexible loop, whereas the positive value indicates the open status of the loop. CV2 was the distance between I335 and R342, which were the residues in the PRTase flexible loop and had interactions with cPRPP. The values of CV2 characterized the shape of the PRTase active site domain, in others words, the values helped determine whether there was sufficient space for cPRPP to bind.

The free energy profile maps showed remarkable differences in conformational changes in the PRTase flexible loops of each system (Figures 10A–D). To prove the number of aMD samples was sufficient, cMD free energy landscape maps were drawn as a comparison (Supplementary Figure S3), it is obvious that aMD can sample more phase spaces than cMD simulations. As the flexible loop for GPATase with cPRPP binding deviated, the CV1 value sampled in this system mainly distributed between 3 and 5 Å (Figure 10C). The outward deviation suggested that the flexible loop was fully extended and in an open state. A representative conformation is presented in Figure 10G (blue). By contrast, the three other systems were in a closed state (Figures 10E,F,H) (also marked in blue). CV2 fluctuated between 12.5 and 13.5 Å for the GPATase–DON–cPRPP complex. Moreover, the CV2 value of this complex was considerably higher than that of the three other systems, indicating that it had sufficient space to bind cPRPP and it contributed to the enzyme catalytic reaction. The corresponding representative conformations are marked in red in Figures 10E–H. Notably, although the loop status of GPATase bound to DON only was closed, the shape of the active site domain obstructed cPRPP binding because of the relatively small amount of space. By contrast, the PRTase domain of GPATase with DON and cPRPP provided sufficient space for cPRPP binding.

FIGURE 10

FIGURE 10. The conformation change of PRTase flexible loops during 400 ns aMD simulations. Free energy profile maps of flexible loops for the (A) GPATase; (B) GPATase–DON; (C) GPATase–cPRPP; and (D) GPATase–DON–cPRPP complexes as a function of CV1 and CV2 in Å; representative structures of flexible loops for the (E) GPATase; (F) GPATase–DON; (G) GPATase–cPRPP; and (H) GPATase–DON-cPRPP complexes. The PRTase domain active site is shown in cartoon. Cartoons marked in blue and red represent different perspectives.

The crystal structures determined by Krahn et al. revealed that the most dramatic structural differences between the inactive and active form of GPATase was its open/close state (Krahn et al., 1997). But dynamics information related to opening-closing conformational transition for the flexible loop remained vague. The above analysis suggested that although the ordering of the PRTase flexible loop was induced by binding cPRPP, the restructuring of the glutamine loop also had a direct steric effect on the PRTase flexible loop that kept it closed and had enough space to assist cPRPP binding to GPATase in the DON–cPRPP–GPATase complex.

The Relative Position Variation Analysis between I335 and Y74

Since the aromatic nucleus of Y74 and methyl group of I335 can form one wall of the ammonia transfer channel in GPATase, the relative position variation analyses between two residues during aMD simulations were performed (Bera et al., 1999). The relative position variation analysis between I335 and Y74 were performed. The distances between Y74 and I335 of the four systems were distinctly different (Figures 11A–D). GPATase bonded to both cPRPP and DON had the shortest distance. I335 and Y74 formed one wall of the NH3 channel because the structure of the NH3 channel depends on hydrophobic amino acid side chains. Changes in the distance between Y74 and I335 during 400 ns aMD simulations are provided in Figure 11E, and the corresponding average distance and standard deviation are presented in Figure 11F. The distance of GPATase with DON and cPRPP remained stable at 6.2 Å, whereas that of GPATase with DON or cPRPP reached about 8 Å. The distance of the latter was not stable relative to that of the former. Therefore, DON and cPRPP played a role together to induce GPATase catalytic action.

FIGURE 11

FIGURE 11. Comparison of the relative position between the glutamine loop and flexible loop of the four systems. Both loops are displayed as cartoons. DON and cPRPP are shown as spheres. Y74 and I335 are depicted as sticks. Representative conformation of the (A) GPATase, (B) GPATase–DON, (C) GPATase–cPRPP, and (D) GPATase–DON–cPRPP complexes during aMD simulations. (E) Distance between Y74 and I335 in the four systems. (F) The average distance between Y74 and I335 in the four systems; standard deviations are labeled in the histogram.

The work of Kim et al. revealed that Y74 is an important residue for coupling between the Gln domain and PRTase domain by analysis of the x-ray structure model and the mutant enzyme (Kim et al., 1996). Subsequently, research from Bera et al. utilized fluorescence monitoring and determination of a K_d experiment to prove that Y74 and I335 were key residues for the interdomain signaling in GPATase (Bera et al., 1999). However, dynamic behavior changes of two residues induced by ligand binding remained unclear. Our dynamic analysis results demonstrated that the sole presence of DON or cPRPP in GPATase may disturb the interaction of these two key residues as the distance between them become larger, and only binding both DON and cPRPP maintains a close enough distance between them to have more interaction contacts. Thus it can be seen that DON and cPRPP functioned together for GPATase catalytic action.

The results of aMD simulations suggested that binding DON to GPATase affected the status and space of the PRTase flexible loop, and binding cPRPP had a great influence on the formation of a salt bridge in the glutamine loop of the enzyme. The salt bridge allowed contact between the glutamine domain and the PRTase domain and helped NH3 entry into the channel. These results were consistent with those of cMD simulations. The mechanism by which GPATase maintains stable states was clarified.

Conclusion

GPATase catalyzed the synthesis of PRA from PRPP and glutamine at separate catalytic sites in different domains. In general, differences in conformations between active and inactive forms of GPATase occurred in flexible loops and extended to core domains. Moreover, binding PRPP to the PRTase domain depended on activating the reaction at the glutaminase domain, and the C-terminal reaction remained stable to keep the N-terminal in active state. Although binding cPRPP to the active site of PRTase could organize the flexible loop, only DON bonded to GPATase could retain the stable state of the loop and keep it closed. The glutamine domain and the PRTase domain approached each another when DON and cPRPP co-existed, a condition that is beneficial for enzymes to play catalytic roles. This study provides important dynamical evidence of conformational changes in DPEase for designing effective inhibitors that can target DPEase.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

CL wrote and revised this paper. SC prepared the tables and figures. TH and HC made the Supplementary Materials. FZ and JY provided some revision advice. WL and WH provided the ideas and modified the papers.

Funding

This work was supported by the Overseas Cooperation Project of Jilin Province (20200801069GH) and the Graduate Innovation Fund of Jilin University (101832020CX097). This work was performed at the High-performance Computing Center of Jilin University.

Conflict of Interest

Author HC was employed by the company Jilin Province TeyiFood Biotechnology Company Limited.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a shared affiliation, though no other collaboration, with several of the authors (CL, SC, TH, FZ, JY, WL, and WH).

Acknowledgments

We would like to thank D.A. Case for providing free access to the AMBER software.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.640994/full#supplementary-material.

References

Atilgan, A. R., Akan, P., and Baysal, C. (2004). Small-world communication of residues and significance for protein dynamics. Biophys. J. 86 (), 85–91. doi:10.1016/s0006-3495(04)74086-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bera, A. K., Chen, S., Smith, J. L., and Zalkin, H. (1999). Interdomain signaling in glutamine phosphoribosylpyrophosphate amidotransferase. J. Biol. Chem. 274 (51), 36498–36504. doi:10.1074/jbc.274.51.36498

PubMed Abstract | CrossRef Full Text | Google Scholar

Berendsen, H. J. C., Postma, J. P. M., Van Gunsteren, W. F., Dinola, A., and Haak, J. R. (1998). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81 (8), 3684–3690. doi:10.1063/1.448118

CrossRef Full Text | Google Scholar

Brannigan, J. A., Dodson, G., Duggleby, H. J., Moody, P. C., Smith, J. L., Tomchick, D. R., et al. (1995). A protein catalytic framework with an N-terminal nucleophile is capable of self-activation. Nature 378 (6555), 416–419. doi:10.1038/378416a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Brinda, K. V., and Vishveshwara, S. (2005). A network representation of protein structures: implications for protein stability. Biophys. J. 89 (6), 4159–4170. doi:10.1529/biophysj.105.064485

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, X. L., Du, B. W., Han, F. J., Zhou, Y., Ren, J. H., Wang, W. H., et al. (2020). Crystal structure of the chloroplastic glutamine phosphoribosylpyrophosphate amidotransferase GPRAT2 from Arabidopsis thaliana. Front. Plant Sci. 11, 157. doi:10.3389/Fpls.2020.00157

PubMed Abstract | CrossRef Full Text | Google Scholar

Chittur, S. V., Klem, T. J., Shafer, C. M., and Davisson, V. J. (2001). Mechanism for acivicin inactivation of triad glutamine amidotransferases. Biochemistry 40 (4), 876–887. doi:10.1021/Bi0014047

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornell, W. D., Cieplak, P., Bayly, C. I., and Kollman, P. A. (1993). Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J. Am. Chem. Soc. 115 (21), 9620–9631. doi:10.1021/ja00074a030

CrossRef Full Text | Google Scholar

Cornell, W. D., Cieplak, P., Bayly, C. I., Gould, I. R., Merz, K. M., Ferguson, D. M., et al. (1995). A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117, 5179–5197. doi:10.1080/07391102.2005.10507028

CrossRef Full Text | Google Scholar

Tom, D., York, D., Pedersen, L. (1993). Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 98 (12), 10089–10092. doi:10.1063/1.464397

CrossRef Full Text | Google Scholar

Delano, W. L. (2002). The PyMOL molecular graphics system. San Carlos, CA: DeLano Scientific.

Doncheva, N. T., Klein, K., Domingues, F. S., and Albrecht, M. (2011). Analyzing and visualizing residue networks of protein structures. Trends Biochem. Sci. 36 (4), 179–182. doi:10.1016/j.tibs.2011.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, T. W. M., Bruntz, R. C., Yang, Y., Song, H., Chernyavskaya, Y., Deng, P., et al. (2019). De novo synthesis of serine and glycine fuels purine nucleotide biosynthesis in human lung cancer tissues. J. Biol. Chem. 294 (36), 13464–13477. doi:10.1074/jbc.RA119.008743

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, Y., Lund, L., Shao, Q., Gao, Y. Q., and Raushel, F. M. (2009). A combined theoretical and experimental study of the ammonia tunnel in carbamoyl phosphate synthetase. J. Am. Chem. Soc. 131 (29), 10211. doi:10.1021/ja902557r

PubMed Abstract | CrossRef Full Text | Google Scholar

Götz, A. W., Williamson, M. J., Xu, D., Poole, D., Le Grand, S., and Walker, R. C. (2012). Routine microsecond molecular dynamics simulations with AMBER on GPUs. 1. Generalized born. J. Chem. Theor. Comput. 8 (5), 1542–1555. doi:10.1021/ct200909j

CrossRef Full Text | Google Scholar

Gasper, P. M., Fuglestad, B., Komives, E. A., Markwick, P. R., and Mccammon, J. A. (2012). Allosteric networks in thrombin distinguish procoagulant vs. anticoagulant activities. Proc. Natl. Acad. Sci. U.S.A. 109 (52), 21216–21222. doi:10.1073/pnas.1218414109

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Segura, L., Witte, J. F., Mcclard, R. W., and Hurley, T. D. (2007). Ternary complex formation and induced asymmetry in orotate phosphoribosyltransferase. Biochemistry 46 (49), 14075–14086. doi:10.1021/bi701023z

PubMed Abstract | CrossRef Full Text | Google Scholar

Grant, B. J., Rodrigues, A. P., ElSawy, K. M., McCammon, J. A., and Caves, L. S. (2006). Bio3d: an R package for the comparative analysis of protein structures. Bioinformatics 22 (21), 2695–2696. doi:10.1093/bioinformatics/btl461

PubMed Abstract | CrossRef Full Text | Google Scholar

Grubmeyer, C., Hansen, M. R., Fedorov, A. A., and Almo, S. C. (2012). Structure of Salmonella typhimurium OMP synthase in a complete substrate complex. Biochemistry 51 (22), 4397–4405. doi:10.1021/bi300083p

PubMed Abstract | CrossRef Full Text | Google Scholar

Hünenberger, P. H., Mark, A. E., and van Gunsteren, W. F. (1995). Fluctuation and cross-correlation analysis of protein motions observed in nanosecond molecular dynamics simulations. J. Mol. Biol. 252 (4), 492–503. doi:10.1006/jmbi.1995.0514

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamelberg, D., De Oliveira, C. A., and Mccammon, J. A. (2007). Sampling of slow diffusive conformational transitions with accelerated molecular dynamics. J. Chem. Phys. 127 (15), 155102. doi:10.1063/1.2789432

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamelberg, D., Mongan, J., and McCammon, J. A. (2004). Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J. Chem. Phys. 120 (24), 11919–11929. doi:10.1063/1.1755656

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. 14 (1), 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Iida, S., Mashimo, T., Kurosawa, T., Hojo, H., Muta, H., Goto, Y., et al. (2016). Variation of free-energy landscape of the p53 C-terminal domain induced by acetylation: enhanced conformational sampling. J. Comput. Chem. 37 (31), 2687–2700. doi:10.1002/jcc.24494

PubMed Abstract | CrossRef Full Text | Google Scholar

Jie, Q., and Liu, D. (2003). Modified conjugate gradient method for diagonalizing large matrices. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 68 (5 Pt 2), 056706. doi:10.1103/PhysRevE.68.056706

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, J., Kuroyanagi, S., Akisada, T., Hagiwara, Y., and Tateno, M. (2012). Unidirectional mechanistic valved mechanisms for ammonia transport in GatCAB. J. Chem. Theor. Comput. 8 (2), 649–660. doi:10.1021/ct200387u

CrossRef Full Text | Google Scholar

Kappel, K., Miao, Y., and Mccammon, J. A. (2015). Accelerated molecular dynamics simulations of ligand binding to a muscarinic G-protein-coupled receptor. Q. Rev. Biophys. 48 (4), 479–487. doi:10.1017/S0033583515000153

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. H., Krahn, J. M., Tomchick, D. R., Smith, J. L., and Zalkin, H. (1996). Structure and function of the glutamine phosphoribosylpyrophosphate amidotransferase glutamine site and communication with the phosphoribosylpyrophosphate site. J. Biol. Chem. 271 (26), 15549–15557. doi:10.1074/jbc.271.26.15549

PubMed Abstract | CrossRef Full Text | Google Scholar

Kokh, D. B., Amaral, M., Bomke, J., Grädler, U., Musil, D., Buchstaller, H. P., et al. 2018). Estimation of drug-target residence times by τ-random acceleration molecular dynamics simulations. J. Chem. Theor. Comput. 14 (7), 3859–3869. doi:10.1021/acs.jctc.8b00230

CrossRef Full Text | Google Scholar

Krahn, J. M., Kim, J. H., Burns, M. R., Parry, R. J., Zalkin, H., and Smith, J. L. (1997). Coupled formation of an amidotransferase interdomain ammonia channel and a phosphoribosyltransferase active site. Biochemistry 36 (37), 11061–11068. doi:10.1021/bi9714114

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., and Nussinov, R. (2002). Close-range electrostatic interactions in proteins. Chembiochem 3 (7), 604–617. doi:10.1002/1439-7633(20020703)3:7<604::aid-cbic604>3.0.co;2-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., and Merz, K. M. (2016). MCPB.py: a python based metal center parameter builder. J. Chem. Inf. Model. 56 (4), 599–604. doi:10.1021/acs.jcim.5b00674

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., Roberts, B. P., Chakravorty, D. K., and Merz, K. M. (2013). Rational design of particle mesh Ewald compatible Lennard-Jones parameters for +2 metal cations in explicit solvent. J. Chem. Theor. Comput. 9 (6), 2733–2748. doi:10.1021/ct400146w

CrossRef Full Text | Google Scholar

Lund, L., Fan, Y., Shao, Q., Gao, Y. Q., and Raushel, F. M. (2010). Carbamate transport in carbamoyl phosphate synthetase: a theoretical and experimental investigation. J. Am. Chem. Soc. 132 (11), 3870–3878. doi:10.1021/ja910441v

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyu, W. L., Arnesano, F., Carloni, P., Natile, G., and Rossetti, G. (2018). Effect of in vivo post-translational modifications of the HMGB1 protein upon binding to platinated DNA: a molecular simulation study. Nucleic Acids Res. 46 (22), 11687–11697. doi:10.1093/nar/gky1082

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theor. Comput. 11, 3696–3713. doi:10.1021/acs.jctc.5b00255

CrossRef Full Text | Google Scholar

Markwick, P. R., and Mccammon, J. A. (2011). Studying functional dynamics in bio-molecules using accelerated molecular dynamics. Phys. Chem. Chem. Phys. 13 (45), 20053–20065. doi:10.1039/C1CP22100K

PubMed Abstract | CrossRef Full Text | Google Scholar

Musick, W. D. L., and Nyhan, W. L. (1981). Structural features of the phosphoribosyl-transferases and their relationship to the human deficiency disorders of purine and pyrimidine metabolis. Crit. Rev. Biochem. Mol. Biol. 11 (1), 1–34. doi:10.3109/10409238109108698

CrossRef Full Text | Google Scholar

Nicolaï, A., Delarue, P., and Senet, P. (2013). Decipher the mechanisms of protein conformational changes induced by nucleotide binding through free-energy landscape analysis: ATP binding to Hsp70. PLoS Comput. Biol. 9 (12), e1003379. doi:10.1371/journal.pcbi.1003379

PubMed Abstract | CrossRef Full Text | Google Scholar

Parry, R. J., Burns, M. R., Jiralerspong, S., and Alemany, L. (1997). Synthesis of (+)-(1S)-1-Pyrophosphoryl-(2R,3R)-2,3-dihydroxy-(4S)-4-( phosphoryloxymethyl)cyclopentane, a stable, optically active carbocyclic analogue of 5-Phosphoribosyl-1-pyrophosphate (PRPP). Tetrahedron 53 (21), 7077–7088. doi:10.1002/chin.199739259

CrossRef Full Text | Google Scholar

Paschek, D., Day, R., and García, A. E. (2011). Influence of water-protein hydrogen bonding on the stability of Trp-cage miniprotein. A comparison between the TIP3P and TIP4P-Ew water models. Phys. Chem. Chem. Phys. 13 (44), 19840–19847. doi:10.1039/c1cp22110h

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, X., Zhang, Y., Li, Y., Liu, Q., Chu, H., Zhang, D., et al. (2018). Integrating multiple accelerated molecular dynamics to improve accuracy of free energy calculations. J. Chem. Theor. Comput. 14 (3), 1216–1227. doi:10.1021/acs.jctc.7b01211

CrossRef Full Text | Google Scholar

Pierce, L. C., Salomon-Ferrer, R., Augusto F de Oliveira, C., McCammon, J. A., and Walker, R. C. (2012). Routine access to millisecond time scale events with accelerated molecular dynamics. J. Chem. Theor. Comput. 8 (9), 2997–3002. doi:10.1021/ct300284c

CrossRef Full Text | Google Scholar

Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23 (3), 327–341. doi:10.1016/0021-9991(77)90098-5

CrossRef Full Text | Google Scholar

Semmelmann, F., Hupfeld, E., Heizinger, L., Merkl, R., and Sterner, R. (2019). A fold-independent interface residue is crucial for complex formation and allosteric signaling in class I glutamine amidotransferases. Biochemistry 58 (22), 2584–2588. doi:10.1021/acs.biochem.9b00286

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, J., Tanner, S. W., Thompson, N., and Cheatham, T. E. (2007). Clustering molecular dynamics trajectories: 1. Characterizing the performance of different clustering algorithms. J. Chem. Theor. Comput. 3 (6), 2312–2334. doi:10.1021/ct700119m

CrossRef Full Text | Google Scholar

Singh, S. A., and Christendat, D. (2006). Structure of Arabidopsis dehydroquinate dehydratase-shikimate dehydrogenase and implications for metabolic channeling in the shikimate pathway. Biochemistry 45 (25), 7787–7796. doi:10.1021/bi068019n.10.1021/bi060366+

PubMed Abstract | CrossRef Full Text | Google Scholar

Skjærven, L., Yao, X. Q., Scarabelli, G., and Grant, B. J. (2014). Integrating protein structural dynamics and evolutionary analysis with Bio3D. BMC Bioinformatics 15 (1), 399. doi:10.1186/s12859-014-0399-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Slanina, Z. K., Boltalina, O. V., Kobayashi, K., and Nagase, S. (2004). B3LYP/6-31G* computations of C60F36(g) isomers. Fuller. Nanotub. Carbon Nanostruct. 12 (3), 691–695. doi:10.1081/fst-200029529

CrossRef Full Text | Google Scholar

Spindel, E. R. (2012). Muscarinic receptor agonists and antagonists: effects on cancer. Handb. Exp. Pharmacol. (208), 451–468. doi:10.1007/978-3-642-23274-9_19

PubMed Abstract | CrossRef Full Text | Google Scholar

Teplyakov, A., Obmolova, G., Badet-Denisot, M. A., and Badet, B. (1999). The mechanism of sugar phosphate isomerization by glucosamine 6-phosphate synthase. Protein Sci. 8 (3), 596–602. doi:10.1110/ps.8.3.596

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, C. J., Ma, B., and Nussinov, R. (1999). Folding and binding cascades: shifts in energy landscapes. Proc. Natl. Acad. Sci. U.S.A. 96 (18), 9970–9972. doi:10.2307/48671

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., and Miao, Y. (2019). Mechanistic insights into specific G protein interactions with adenosine receptors. J. Phys. Chem. B. 123 (30), 6462–6473. doi:10.1021/acs.jpcb.9b04867

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Shen, D., Xia, G., Shen, W., Macdonald, R. L., Xu, D., et al. (2016). Differential protein structural disturbances and suppression of assembly partners produced by nonsense GABRG2 epilepsy mutations: implications for disease phenotypic heterogeneity. Sci. Rep. 6, 35294. doi:10.1038/srep35294

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Sheridan, R., Sumer, S. O., Schultz, N., Xu, D., and Gao, J. (2018). G2S: a web-service for annotating genomic variants on 3D protein structures. Bioinformatics 34 (11), 1949–1950. doi:10.1093/bioinformatics/bty047

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X. S., Roitberg, A. E., and Richards, N. G. (2009). Computational studies of ammonia channel function in glutamine 5'-phosphoribosylpyrophosphate amidotransferase. Biochemistry 48 (51), 12272–12282. doi:10.1021/bi901521d

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X. S., Roitberg, A. E., and Richards, N. G. (2009). Computational studies of ammonia channel function in glutamine 5'-phosphoribosylpyrophosphate amidotransferase. Biochemistry 48 (51), 12272–12282. doi:10.1021/bi901521d

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolan, D. W., Cheong, C. G., Greasley, S. E., and Wilson, I. A. (2004). Structural insights into the human and avian IMP cyclohydrolase mechanism via crystal structures with the bound XMP inhibitor. Biochemistry 43 (5), 1171–1183. doi:10.1021/bi030162i

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., Li, C., Yang, H., Guo, X., Huang, T., and Han, W. (2020). Computational study on the effect of inactivating/activating mutations on the inhibition of MEK1 by trametinib. Int. J. Mol. Sci. 21 (6), 2167. doi:10.3390/ijms21062167

CrossRef Full Text | Google Scholar

Zhu, Y., Li, T., Ramos da Silva, S., Lee, J. J., Lu, C., Eoh, H., et al. (2017). A critical role of glutamine and asparagine γ-nitrogen in nucleotide biosynthesis in cancer cells hijacked by an oncogenic virus. Mbio 8 (4), e01179. doi:10.1128/mBio.01179-17

PubMed Abstract | CrossRef Full Text | Google Scholar