Synergistic effect of conformational changes in phosphoglycerate kinase 1 product release

Abstract In the glycolysis pathway, phosphoglycerate kinase 1 (PGK1) transfers one phosphoryl-group from 1,3-diphosphoglycerate (1,3BPG) to ADP to product 3-phosphoglycerate (3PG) and ATP. The catalytic process is accompanied with the conversion between the open conformation and the closed conformation of PGK1. However, the dynamic collaboration mechanism between the PGK1 conformation transition and the products releasing process remains poorly understood. Here using molecular dynamics simulations combined with molecular mechanics generalized born surface area (MM/GBSA) analysis, we demonstrated that PGK1 in the closed conformation first releases the product ATP to reach a semi-open conformation, and releases the product 3PG to achieve the full open conformation, which could accept new substrates ADP and 1,3BPG for the next cycle. It is noteworthy that the phosphorylation of PGK1 at T243 causes the loop region (residues L248-E260) flip outside the protein, and the phosphorylation of Y324 leads PGK1 become looser. Both modifications cause the exposure of the ADP/ATP binding site, which was beneficial for the substrates/products binding/releasing of PGK1. In addition, the other post translational modifications (PTMs) were also able to regulate the ligands binding/releasing with different effects. Our results revealed the dynamic cooperative molecular mechanism of PGK1 conformational transition with products releasing, as well as the influence of PTMs, which would contribute to the understanding of PGK1 substrates/products conversion process and the development of small molecule drugs targeting PGK1. Communicated by Ramaswamy H. Sarma.


Introduction
Tumors are heterogeneous diseases, and the cellular and structural heterogeneity results in the complex metabolic patterns.With the in-depth understanding of tumor biology and the complexity of tumor metabolism (Reinfeld et al., 2021), the metabolic reprogramming represented by energy metabolism reprogramming is found as a hallmark of malignant tumors (Kim & Deberardinis, 2019;Mart� ınez-Reyes & Chandel, 2021).Most cancer cells, similar with other rapidly proliferating cells, prefer a phenomenon known as the warburg effect (Heiden et al., 2009;Koppenol et al., 2011) that to produce energy in the cytoplasm with the glycolytic pathway rather than the oxidative phosphorylation pathway in mitochondria (Yang & Lu, 2013;Yang & Lu, 2015),.Phosphoglycerate kinase 1 (PGK1), as an essential enzyme for ATP generation in the glycolytic pathway, plays a crucial role during the progression of various tumors (S.Sun et al., 2015), including promoting tumor cells proliferation, migration and invasion (Qian et al., 2017;Wilson et al., 2019).PGK1 is composed of two domains, N lobe and C lobe, connected by the linker region which formed with two adjacent helixes as the core structure (Bowler, 2013;Figure 1A), and the ligand binding domain lies between these two lobes (Figure 1B-D) (Pey et al., 2013;Young et al., 2007).The transfer of the phosphoryl-group from 1,3-diphosphoglycerate (1,3BPG) to ADP with the formation of 3-phosphoglycerate (3PG) and ATP (Bernstein et al., 1997;Watson et al., 1982) is the main function of PGK1 (Figure 1E) (Yan et al., 2012;Zhang et al., 2005).Three conformation states are known under different substrate binding states, that the open state with the absence of substrates, the semi-open state after binding to 1,3BPG, and the closed state with ADP and 1,3BPG (Bernstein et al., 1997).However, the dynamic collaboration mechanism between the PGK1 conformation transition and the substrates/products conversion process remains poorly understood.
Recent researches show that protein post translational modifications (PTMs) of PGK1 play an important role in tumor metabolism and growth regulation.The phosphorylation of PGK1 Threonine 243 (pT243) mediated by recombinant phosphoinositide dependent protein kinase 1 (PDPK1) causes the changes in affinity between PGK1 and substrates (Y.Zhang et al., 2018), thus promoting the activity of PGK1 to the glycolytic direction, which is conducive to aerobic glycolysis of human glioblastoma multiforme.Qian and colleagues found the autophosphorylation of PGK1 Tyrosine 324 (pY324) could enhance the rate of glycolysis by increasing the release of ATP, to promise the malignant tumors, and to down-regulate the release by dephosphorylation with phosphatase PTEN (Qian et al., 2019).Moreover, ELP3 family histone acetyltransferase KAT9-mediated acetylation of PGK1 at Lysine 220 (aK220) inhibits protein activity by blocking binding to ADP (Wang et al., 2015), which could be reversed by deacetylation with histone deacetylase 3 (HDAC3).The other reversible acetylation modification at Lysine 323 (aK323) by P300/CBP-associated factor (PCAF) or NAD-dependent protein deacetylase sirtuin-3 (SIRT3) could regulate the PGK1 activity (Hu et al., 2017), and eventually regulate the occurrence and development of liver cancer.Furthermore, the proliferation and glycolysis of colon cancer cells are boosted under the glycosylation of Threonine 255 (gT255), which promotes the glycolysis with the PGK1 activity enhancement, and inhibits the tricarboxylic acid cycle by inducing the translocation of PGK1 to mitochondria (Nie et al., 2020).Nevertheless, how these protein PTMs regulate PGK1 activity remain unclear, which largely restricts the research of PGK1 in tumor.
Using computational simulation to describe the flexibility and dynamics of biological macromolecules in terms of dynamic spatial conformation and thermodynamic properties, which couldn't observed directly via experimental methods, exhibits distinct advantages in revealing the molecular mechanism.Researchers have already revealed the mechanism of PGK1 conformational regulation with substrate binding by molecular dynamics simulations (Palmai et al., 2014).In this study, we investigate the cooperative mechanism of PGK1 conformational transition for substrates releasing, as well as the molecular mechanism of PTMs (T243 and Y324 phosphorylation, K220 and K323 acetylation, T255 glycosylation) in PGK1 activity regulation with molecular dynamics simulations and molecular mechanics generalized born surface area (MM/GBSA) analysis.

Preparation of the protein structures
The closed state structure of human PGK1 (PDB ID: 2X15) binding with different ligands (Magnesium ions (MG), ATP, ADP, 1,3BPG and 3PG) was chosen to conduct the following research with different ligands binding state of substrate/ product (Bowler, 2013).The missing loop (residues D134 to V140) was completed by SWISS MODEL using the crystal structure of (PDB ID: 2X15) as the template (https://swissmodel.expasy.org)(Waterhouse et al., 2018).Then, The complete structure was divided into three types: 1) the structures of PGK1 binding with MG, ATP and 3PG was obtained by removing ADP and 1,3BPG from complete structure; 2) the structures of PGK1 binding with MG, ADP and 1,3BPG was obtained by removing ATP and 3PG from complete structure; 3) the structures of apo PGK1 was obtained by removing all ligands from complete structure.Subsequently, modified the T243 (phosphorylation), Y324 (phosphorylation), K220 (acetylation), K323 (acetylation), and T255 (glycosylation) in PGK1 to obtain the structures of the post-translationally modified PGK1 binding with different ligands.

Molecular dynamics simulations
Atomistic MD simulations of wild type and post-translationally modified PGK1 binding with different substrates/products were performed using AMBER16 software (Case et al., 2016), embracing Amber ff14SB force field for protein (Maier et al., 2015).The parameters of ATP, ADP, phosphorylated tyrosine, phosphorylated threonine, glycosylated threonine and acetylated lysine were obtained from the parameters reported by Bryce (http:// amber.manchester.ac.uk (Homeyer et al., 2006)).The electrostatic potential of 3PG and 1,3BPG were calculated using the Gaussian 09 program with the B3LYP functional under 6-311 G � basis set.The partial charges of the substrate molecules were derived using the RESP charge fitted with the antechamber module in Amber 16 (Case et al., 2016).Each simulation complex was immersed in a hexahedral solvent box filled with TIP3P (Jamsawang, 2018) water molecules accompanied by a distance of 10 Å from the surface of the protein to the edge of the water box, which was neutralized with a number of magnesium and chloride ions.Energy minimization was conducted by imposing a strong restraint on simulations system and was followed by minimizing whole system for a few thousands steps.Then, a NVT simulation (Berendsen et al., 1984) is performed to heat whole system slowly from 0 to 300 K, which was followed by a 5 ns NPT equilibration run.Subsequently, the production run was performed for 1ls with a time step of 4 fs, and the coordinates for all models were saved every 2 ps (each simulation system was repeated for five times).During the production run, the SHAKE algorithm was employed to constrain all bonds associated with hydrogen atoms.Electrostatic interaction was treated by Particle Mesh Ewald (Darden et al., 1993) and the cutoff value of nonbonded interactions is set to 9 Å.

Data analysis
We performed the analysis of the distance between relative atoms, root mean square fluctuation (RMSF), root mean square deviation (RMSD, using the crystal structure as the reference conformation), principal component analysis (PCA), and dynamical cross-correlation mapping (DCCM) using either CPPTRAJ (Roe & Cheatham, 2013) or TCL scripts (Humphrey et al., 1996) according to the last 100 ns trajectories of each simulations, then further processed and plotted using Python.All structures were extracted from the simulation trajectories and processed with Pymol.To investigate the extent of correlation motions caused by post-translational modification in each PGK1-ligand complex and apo PGK1, the cross correlation matrix, C ij , which reflects the fluctuations in the coordinates of the atoms relative to their average positions from the last 100 ns of the simulations, was determined by the following equation (Maisuradze et al., 2010): where Dri/Drj indicated the deviation of the Ca atom of the ith residue from its mean position.The value of C ij fluctuates between À 1 and 1.A positive C ij value indicated the correlation motion between the ith residue and the jth residue, while a negative C ij value represented the inverse correlation motion.The eigenvalues indicated the magnitude of motions along a direction.In general, the first few principal components (PCs) described the most important slow system modes related to the functional motions of a biomolecular system.

The calculation of binding free energy
The binding free energy between PGK1 and substrates/products as well as the energy contribution of each residue to the binding energy were obtained by Molecular Mechanics Generalized Born Surface Area (MM-GBSA) method (Kollman et al., 2000;Xu et al., 2013) base on the 1000 snapshots extracted from the last 100 ns MD trajectory.The total binding energy (DG bind ) was computed according to the following equation: where DG bind represented the binding free energy between PGK1 and substrates/products, calculated according to the discrepancy between the sum of the free energy of the ligand (G ligand ) and the protein (G protein ) and the total free energy of the complex (G complex ).The binding energy was resented as: where E MM was the molecular mechanics energy of the molecule expressed as the sum of the internal energy of van der Waals energies and molecule plus electrostatic.The solvation free energy was expressed as nonpolar and polar contributions to the solvation energy: G nonpolar was calculated with the LCPO algorithm based on solvent-accessible surface area (SASA) (Weiser et al., 1999): where c¼ 0.0072 kcal/mol/Å, and b¼ 0 kcal/mol (H.Sun et al., 2018).

The cooperative mechanism of PGK1 conformational transition with products releasing
To investigate the dynamic mechanism between the PGK1 conformational transition and the products releasing, we performed a series of molecular dynamic simulation with durations of 1 ls for various wild type PGK1 systems at different ligands binding states with products and substrates (ADP/ 1,3BPG, ATP/3PG, ATP, 3PG and apo form), respectively.The root-mean-square deviations (RMSD) of the Ca atom (the crystal structure 2X15 as the reference, Figure S1) and the radius of gyration (Rg) were used to characterize the global dynamic state of the protein.In these systems, RMSD and Rg show a similar variation tendency, that the mean values step up along with the conformational change from closed to open in different systems (Figure S2).However, the most direct manifestation of protein conformational changes during different states was the center of mass distance variation between N lobe and C lobe.The main distribution of this distance increased from 32.8 Å with substrates binding (WT/ADP/1,3BPG with lemon solid line in Figure 2) to 34.6 Å with products (WT/ATP/3PG with orange solid line in Figure 2), which is more suitable for products releasing.
The difference also could be observed in the complex systems of PGK1 with one product (ATP or 3PG respectively), which aims to study the conformational dynamics at the absence of the other product (3PG or ATP) to find the products releasing order of PGK1.The distribution of distance between N/C lobes in PGK1/ATP and PGK1/3PG systems have totally different trends in the simulation trajectories, that in the former system the distance decreases from 34.6 Å to 32.7 Å (Figure S3).While the main distribution increases to around 35.8 Å in PGK1/ 3PG system (WT/3PG with cyan solid line in Figure 2).This distinction indicated that the protein prefers a much open conformation when bound with 3PG rather than ATP, which would benefit for the following release of the second product.Furthermore, removed 3PG from the PGK1/3PG system result in the further increase of the distance between N lobe and C lobe up to 39.4 Å (WT-apo with pink solid line in Figure 2), which is agreement with the conclusion reported in previous studies (Bernstein et al., 1997) that PGK1 reach to the fully open state after releasing all of the products.
To verify of the results, we compared the result conformations with other PGK1 crystal structures (Figure S4).The distance between N lobe and C lobe is 32.4 Å in the closed structure (PDB ID: 2WZC) (Cliff et al., 2010) and is 39.8 Å in the open structure (PDB ID: 1FW8) (Tougard et al., 2002), which are consistent with the distance results from simulations (Figure S4A and C).In addition, the conformational alignment also showed that the simulated structure had a high degree of overlap with the crystal structure (Figure S4B and D).
The above results of distance and similarity could help us to illustrate a preliminary hypothesis of products releasing process related with conformational changes of PGK1, that after the reaction, the protein release ATP first accompanied with conformational change from closed to semi-open, and then release 3PG to reach the full open form (Cartoon structures in Figure 2).
To further confirm the presumption, the principal component analysis (PCA) was used to study the dynamics of PGK1 under different ligand binding states during simulations.From the movement tendencies of PGK1 under different ligand-binding states by PCA analysis, we found that in PGK1products systems (PGK1/ATP/3PG in Figure 3B and PGK1/3PG in Figure 3C) shows a distinct trend against PGK1-substrate system (Figure 3A), that these two lobes in the products systems have strong tendency to move far away from each other during the simulations.Moreover, the movement trend in PGK1/3PG system is significantly larger than the trend in PGK1/ATP/3PG system.This result also proved that the release of ATP further increases the distance between the N lobe and the C lobe, and urges PGK1 to reach the semi-open state.In order to understand the mechanism between products releasing and protein conformational changes, we used the cross correlation matrix analysis to detect the specific protein regions with intense relevant dynamic motions, which would be helpful in clarifying the movement difference of PGK1 under different ligand binding states.The most significant negative correlation motions in the covariance matrix maps of protein from three different simulation systems, PGK1/ADP/1,3BPG (Figure 3D), PGK1/ATP/3PG (Figure 3E) and PGK1/3PG (Figure 3F), are mainly occurred between V20-N180 (N lobe) and K220-R330 (C lobe), which is consistent with the results of PCA.And it's worth noting that after ATP releasing, there is an obvious increasing trend of the negative correlation motion between N lobe and C lobe of PGK1.The results of residues motion trend analysis of PKG1 (Figure 3) are consistent with the results of distance distributions between N lobe and C lobe (Figure 2), which are able to well confirm our previous presumption.

T243 phosphorylation enhances the ADP/ATP binding of PGK1 by causing the loop L248-E260 to flip to the outside of the protein
Several PTMs of PGK1 were reported in previous researches and could affect the occurrence and development of a variety of diseases, but known little with the mechanisms.We performed MD simulations on PGK1 with several PTMs (pT243, pY324, aK220, aK323 and gT255), to discuss the mechanism of PTM regulating the biological function of PGK1.We made some priori assumptions on how these PTMs regulate PGK1: 1) to regulate the conversion efficiency of substrate/product through affecting the conformational transition of PGK1; 2) to change the binding states of substrates/products by affecting the contribution binding free energy of key residues.We analyzed the effect of PTMs on the conformational transition and activity of PGK1 with these assumptions.
From the results of RMSD (Figure S5) and Rg (Figure S6) indicate that compared to the wild type PGK1, the modification systems have opposite effect on the protein conformation.Two modifications aK323 and gT255 could cause PGK1 to more  stable and tightness in structure, whereas in pT243, pY324 and aK220 systems PGK1 structure tend to be more unstable and flexible.The same results could also get from the results of distance between N/C lobes, that the modifications of pT243, pY324 and aK220 could increase the distance significantly (Figure S7), which mean an easier conformation for transiting to the open state.Since K220 located at the substrates binding pocket, it is easy to understand that acetylation of K220 would the property of the pocket and the interaction with ADP and 1,3BPG, to induce the conformational transition towards the open state.The root mean-square fluctuation (RMSF) of PGK1 wild type and PTMs systems during the simulations were monitored to characterize the fluctuation of protein, to investigate the influence of the modifications outside the substrate binding pocket.The analysis of protein residue fluctuation showed that the effect of phosphorylation on T243 is mainly reflecting residues L248-E216 in protein movement (Figure 4A).The increase of fluctuation leads the loop region of residues L248-E260 flip outside the protein, makes the exposing of ATP, which is benefit for the following releasing (Figure 4B-C).From the structural analysis we found the phosphorylated T243 induces a long-range interaction pathway, which eventually results in the exposure of ATP.The phosphorylated T243 could form stable interactions with K246, and cause the flipping of K246 side chain to interact with the main chain of L257 to make the L257 loop move away from ATP (Figure 5).
Since that the change in residues interaction is caused by T243 phosphorylation, which is not related with the products binding; we can audaciously speculate that T243 phosphorylation could also lead to the exposure of the substrates.To further corroborate this conjecture, we analyzed the effect of T243 phosphorylation on the conformational transition of PGK1 protein under different substrate/product binding states.The results of RMSD (Figure S1) and Rg (Figure S2) indicated that compared to wild type PGK1, T243 phosphorylation was able to increase the RMSD and Rg values of PGK1 regardless of whether PGK1 binding substrates or products.The data of the distance distribution between N lobe and C lobe (Figure S8) showed that T243 phosphorylation of PGK1 could increase the distance from N lobe to C lobe under either substrates or products binding.This suggested that T243 phosphorylation could enhance the substrates/products conversion rate of PGK1 by exposing the substrates/products to the surface of PGK1 protein.

Y324 phosphorylation leading to PGK1 become more loose
Similarly, we analyzed the effect of Y324 phosphorylation against the protein conformation and ligand binding.As we can see, the shift of protein fluctuations (Figure S9A) caused by Y243 phosphorylation are concentrated in several regions (G237-D259, K291-A296 and A309-A33).The phosphorylation of Y324 also leads an increase in the main distribution value of protein solvent accessible surface area (SASA, Figure S9B) from 16500 Å 2 to 17600 Å 2 , suggesting that the pocket of PGK1 becomes much open after Y324 phosphorylation.Superimposed the simulation result with WT structure showed that the cavity in Y324 phosphorylation-modified PGK1 was larger than that in wild type PGK1(Figure S9C), which correlated with the position changes of the regions where residues G237-D259, K291-A296 and A309-A331 were located.The simulation data shows that Y324 phosphorylation can lead to PGK1 become looser by significantly enhancing the fluctuations of residues G237-D259, K291-A296 and A309-A331.Further analysis revealed that the reason for those changes was that the introduction of a phosphate group could increase the steric hindrance effect of pY324, resulting in significant changes in the spatial positions of residues V340, G238, P283, F286 and C318 around pY324 (Figure S9C-E).The change of the position of these nonpolar amino acids was closely related to the change of SASA value.Similarly, we can speculate that changes in this fluctuation causing by Y324 phosphorylation may greatly affect the conformational transition of PGK1 protein.Therefore, we also analyzed the effect of Y324 phosphorylation on the conformational transition of PGK1 protein.The results of RMSD (Figure S10) and Rg (Figure S11) demonstrated that Y324 phosphorylation led to PGK1 become unstable and fluctuation regardless of whether PGK1 binding substrates or products.As expected, the data of the distance distribution between N lobe and C lobe (Figure S12) indicated that Y324 phosphorylation can greatly increase the distance between N and C lobes in PGK1 regardless of the presence of substrates/ products.Previous study (Qian et al., 2019).reported that PGK1, as a protein kinase, was able to enhance the efficiency of PGK1 to generate ATP through autophosphorylation at the Y324, thereby enhancing the glycolysis capability of tumor cells and promoting tumor growth.Our simulation results well revealed the microscopic molecular mechanism of Y324 phosphorylation enhancing PGK1 activity: Y324 phosphorylation can enhance PGK1 from close conformation to open conformation, which was consistent with the results of biological experiments (Qian et al., 2019).

PTMs regulate the catalytic activity of PGK1 by affecting the binding of PGK1 to ADP/ATP
The second mode we mentioned previously that protein post translational modifications regulate PGK1 activity may be that PTMs changes the binding state of PGK1 and substrates/products.Therefore, we analyzed the effect of PTMs on the binding free energy between PGK1 and substrates/ products.The MM-GBSA results (Table 1) demonstrated that K323 acetylation led to the binding free energy between PGK1 and ADP increase from À 73.37 kcal/mol to À 89.65 kcal/mol.Similarly, T255 glycosylation of PKG1 caused the binding free energy between PGK1 and ADP achieving to À 87.57kcal/mol.Both of two modifications enhanced the contribution of residues N337 and G374 to the binding energy between PGK1 and ADP (Table S1).The acetylation at K220 site reduced the binding free energy between PGK1 and ADP to À 39.83 kcal/mol, which mainly due to the weakening of the energy contributions of residues G214, A215, K216, K220 and G374 (Table S1).This result was consistent with the previous experimental result that K220 acetylation inhibited PGK1 activity by blocking its binding to ADP (Wang et al., 2015).Interestingly, T243 and Y324 phosphorylation did not alter the binding free energy between PGK1 and ADP.Those results illustrated that K220 acetylation, K323 acetylation and T255 glycosylation regulated the binding of PGK1 and ADP by affecting the contribution of key residues to the binding free energy, thereby realizing the regulation of PGK1 activity.Our results well revealed the microscopic molecular mechanism of the acetylation modification of K323 and glycosylation modification of T255 promoted the proliferation of tumor cells by enhancing the activity of PGK1 enzyme, improving the glucose absorption and cell metabolism rate (Hu et al., 2017;Nie et al., 2020).Since that K220 located at the substrates binding pocket (Figure1), it was easy to understand that acetylation breaks the interaction of K220 with ADP and 3PG, which reduced the activity of PGK1.However, PTMs of the residues (K323 acetylation and T255 glycosylation) did not affect the binding of 3PG with PGK1, which could increase the activity of PGK1 by largely enhancing the binding of PGK1 to ADP.
The T255 glycosylation (gT255/ADP/1,3BPG) reduced the relative occurrence frequencies of the RMSD form 2.1 Å (WT/ADP/1,3BPG) to 1.8 Å (Figure S1), indicating PGK1 become more stable.The results of residues fluctuation analysis (Figure S13) further confirmed residues in T255 glycosylation PGK1 exhibited lower fluctuation compared to residues in wild-type PGK1, which was consistent with the RMSD results.It was worth noting that the fluctuation of residues near T255 decreased more significantly after glycosylation.This change in fluctuation broken the interaction of F258 with V217 and A218 by allowing gT255 to form a stable interaction with F258 (Figure 6 A, B, C, F), leading to a new interaction between L233 and N337 (Figure 6A, B, D, G).The alteration of this interaction resulted in overturning of the side chain of N337 to form a stable interaction with the phosphate group of ADP (Figure 6A, B, E, H), which may be the fundamental reason for the enhancement of the binding between PGK1 and ADP caused by T255 glycosylation.Similarly, the K323 acetylation (aK323/ADP/1,3BPG) led to the relative occurrence frequencies of the RMSD decrease form 2.2 Å (WT/ADP/1,3BPG) to 1.8 Å (Figure S1) as well as reduced the fluctuation of PGK1 protein (Figure S14), indicating that K323 acetylation could promote the stability of PGK1 protein to enhance the binding to ADP.Further analysis results suggested that K323 glycosylation led to the change of protein fluctuation and stability, allowing aK323 to form stable interactions with P283 and K229 (Figure 7), resulting in the PGK1 more compact, which may be the underlying reason for the enhancement of the binding between PGK1 and ADP caused by K323 acetylation.
Given that T243 phosphorylation and Y324 phosphorylation hardly affected the binding of PGK1 to the substrates, did these two post-translational modifications affect the release of PGK1 to the products?Therefore, we calculated the effect of T243 phosphorylation and Y324 phosphorylation on the binding free energy between PGK1 and ATP.It was expected that compared with wild type PGK1, the binding free energy between T243 phosphorylated PGK1 and ATP decreased from À 34.11 kcal/mol tayto À 23.67 kcal/mol (Table S2).Remarkably, T243 phosphorylation regulated the release of PGK1 to the substrate mainly by weakening the contribution of residues K216 and K220 to the binding free energy of PGK1 and ATP (Table S3).In addition, our previous study also confirmed that K216 turnover caused by T243 glycosylation also reduced the binding free energy of PGK1 and 3PG(Y.Zhang et al., 2018).Similarly, phosphorylation of Y324 also resulted in a decrease in the binding free energy between PGK1 and ATP to À 22.78 kcal/mol (Table S2), mainly by weakening the contribution of residues K220 and N337 to the free energy of binding (Table S3).It was worth noting that the binding free energy between PGK1 and substrate (ADP) was significantly stronger than that between PGK1 and the product (ATP), which will be more conducive to the release of the products.The MM-GBSA results indicated that the post translational modification was able to regulate the substrates/products binding (K220 acetylation, K323 acetylation and T255 glycosylation) and products releasing (T243 phosphorylation and Y324 phosphorylation).

Discussion and perspective
Compared with normal cells, tumor cells are characterized by abnormal energy metabolism with a higher rate of glycolysis, called the Warburg effect.Phosphoglycerate kinase PGK1 is a key enzyme in the process of glycolysis, which can catalyze the production of ATP and plays an important role in tumor energy metabolism.In addition to its catalytic activity, PGK1 also participates in the expression and regulation of a variety of oncoproteins, to affect the signal pathways related to proliferation and metastasis in tumors, and significantly promotes the proliferation, growth and metastasis of tumor cells.Hence, PGK1 is a target for tumor therapy and a hotspot in tumor therapy research.PGK1 is consists of a C lobe, N lobe and intermediate linker domains, of which the N lobe can bind the 1,3BPG/3PG, the C lobe can bind ADP/ATP.The catalysis of substrates by PGK1 is achieved by switching between open and closed states (Bernstein et al., 1997), which is related to the change of the distance between N lobe and C lobe.Previous studies suggested that the unbound PGK1 under the inactive open state turned into the semi-open state after binding to 1,3BPG, and further combined with ADP causes PGK1 to reach an active closed state.Then, PGK1 catalyze the transfer of phosphoryl-group from 1,3BPG to ADP with the formation of 3PG and ATP.Our study indicated that after complete the catalytic process, PGK1 releases the product ATP first to reach the semi-open conformation, followed by releasing the product 3PG to achieve the full open conformation (Figure 8).However, how to regulate the glycolysis process of tumor cells by modificating the conversion rate of PGK1 to substrate products is still a problem that needs to be further explored.
Recent research (Qian et al., 2019;Y. Zhang et al., 2018) shows that protein post-translational modifications (phosphorylation, acetylation and glycosylation) could affect the occurrence and development of tumors by regulating the binding of PGK1 with substrates/products and downstream proteins in signaling pathways (Hu et al., 2017;Wang et al., 2015), which would provide a new research direction for exploring the role of PGK1 in tumor development (Nie et al., 2020).However, the microscopic molecular mechanism of protein post-translational modification regulating PGK1 activity is still unclear, which largely restricts the research of PGK1 in tumor progression.Our research revealed PTMs were able to regulate the conversion efficiency of PGK1 to substrate/ product by affecting the conformational transition of PGK1 and change the binding state of PGK1 and substrates/products by affecting the contribution of key residues to the binding free energy (Figure 8).With the increase of protein post-translational modification sites of PGK1 discovered, more ways of protein post-translational modification to regulate PGK1 will be discovered.In addition, how to develop small molecule drugs targeting PGK1 to treat tumors based on the known PGK1 catalytic process and protein post-translational modification mechanism still face many challenges.

Conclusion
In this study, molecular dynamics simulations combined with MM/GBSA analysis were prepared to illustrate the dynamic cooperative molecular mechanism of PGK1 conformational transition with products release and the influence of protein post translational modifications on this molecular mechanism.The simulations results indicated when the products (ATP and 3PG) bound to the PGK1, PGK1 in the close conformation first released the product ATP to reach the semiopen conformation, followed by releasing the product 3PG to achieve the full open conformation which can accept the new substrates ADP and 1,3BPG for the next cycle.T243 phosphorylation resulted in the loop region where residues L248-E260 located to flipping to the outside of the protein by forming a stable interaction with K246 and Y324 phosphorylation were able to lead to PGK1 become more fluctuation by significantly enhancing the fluctuations of residues G237-D259, K291-A296 as well as A309-A331, both exposing the ADP/ATP to the cavity between N lobe and C lobe of PGK1.In addition, the MM/GBSA results indicated that K220 acetylation, K323 acetylation and T255 glycosylation enhanced binding free energy between PGK1 and ADP, while T243 phosphorylation and Y324 phosphorylation reduced binding free energy between PGK1 and ATP.Those changes were achieved by post translational modifications regulating the contribution of key residues to the binding energy.It suggested that PTMs regulate the substrates/products conversion process of PGK1 by affecting the substrates binding (K220 acetylation, K323 acetylation and T255 glycosylation) and products releasing (T243 phosphorylation and Y324 phosphorylation).Our results illustrated an elegant mechanism in which relatively products release induced changes combine to induce a major conformational change that serves precisely to conversion process of PGK1 to substrates and products.

Figure 1 .
Figure 1.Overview of PGK1.(A) Organization of the PGK1 protein.The N lobe and C lobe are colored by cyan and yellow respectively, while The Linker is colored by orange.(B) Three-dimensional structure of PGK1 (PDB ID: 2X15).Each domain is highlighted by cartoon and the color is consistent with that in (A), while the substrates and magnesium ion are represented by sticks and spheres respectively.(C) Amino acid sequence of PGK1.Helices and b-strands are showed as red rectangles and green arrows respectively.(D) Binding pocket for ADP and 1,3BPG in PGK1.(E) Schematic diagram of PGK1 catalytic reaction.

Figure 2 .
Figure 2. The dynamic cooperative molecular mechanism of PGK1 conformational transition with substrates release.The lemon solid line, orange solid line, cyan solid line and pink solid line represents respectively the distribution frequency of the distance between the center masse of N lobe and the center mass of C lobe in PGK1 binding with ADP/1,3BPG, ATP/3PG, 3PG and no products/substrates during the simulations.Correspondingly, the cartoons and sticks represent the conformational transition process of PGK1 from close to fully open, accompanied by the release of the substrates.

Figure 3 .
Figure 3. Dominant domain motions of PGK1 accompanied by the release substrates of during the simulations.(A-C) The first principal components in PCA of PGK1 binding with ADP/1,3BPG (A), ATP/3PG (B) and 3PG (C).The PGK1 is showed by tube and the arrow represents the trend of protein movement.(D-F) Crosscorrelated motions of the Ca atoms of PGK1 binding with ADP/1,3BPG (D), ATP/3PG (E) and 3PG (F).Positive regions marked in red and indicate strongly correlated residue motions, whereas negative regions colored in blue associated with anticorrelated movements.

Figure 4 .
Figure 4. T243 phosphorylation expose the ATP to the cavity between N lobe and C lobe of PGK1.(A) Comparison of the root mean square fluctuation (RMSF) between wild type and T243 phosphorylated PGK1.(B, C) Comparison of the structure between wild type and T243 phosphorylated PGK1.PGK1 and ATP are highlighted by surface and stick respectively, while the residues L248-E260 are colored by marine.

Figure 5 .
Figure 5. T243 phosphorylation alters the interaction between residues L257 and K326.(A-C) Comparison of the relative positions of residues L248-E260 between wild type and T243 phosphorylated PGK1.The residues L248-E260 in wild type and T243 phosphorylated PGK1 are highlighted by marine and red respectively (A), while the substrates and residues T243/pT243, K246, L257, M312 and F292 are showed as sticks (B, C).The value on the black dotted line represents the distance between the corresponding atoms.(D-F) Distribution frequency of the distance between L257:O and K246:NZ (D), between T243/pT243:OG1 and K246:NZ (E), between L257:CD1 and ATP:N6 (F) in wild type and T243 phosphorylated PGK1 during the simulations.

Figure 8 .
Figure 8. Schematic diagram of the catalytic cycle of PGK1 and the regulation of the catalytic cycle by protein post-translational modifications.