Molecular dynamics simulation reveals how phosphorylation of tyrosine 26 of phosphoglycerate mutase 1 upregulates glycolysis and promotes tumor growth

Phosphoglycerate mutase 1 (PGAM1) catalyzes the eighth step of glycolysis and is often found upregulated in cancer cells. To test the hypothesis that the phosphorylation of tyrosine 26 residue of PGAM1 greatly enhances its activity, we performed both conventional and steered molecular dynamics simulations on the binding and unbinding of PGAM1 to its substrates, with tyrosine 26 either phosphorylated or not. We analyzed the simulated data in terms of structural stability, hydrogen bond formation, binding free energy, etc. We found that tyrosine 26 phosphorylation enhances the binding of PGAM1 to its substrates through generating electrostatic environment and structural features that are advantageous to the binding. Our results may provide valuable insights into computer-aided design of drugs that specifically target cancer cells with PGAM1 tyrosine 26 phosphorylated.


SUPPLEMENTARY DATA Tunnel analysis
Tunnel analysis was carried out by using the CAVER 3.0 software package [1]. CAVER is a software tool widely used for identification and visualization of tunnels and channels of a specified protein structure. Based on the trajectory generated by MD simulations, CAVER 3.0 was capable of finding the possible egress tunnels from the binding side to the surface of the enzyme, enabling automatic analysis of large ensembles of protein conformations. Several pathways with nearly identical axes may be classified into one cluster. The lowest cost pathway is selected and all the other pathways in the cluster are discarded. A trajectory of MD simulation serves as the input, while the detailed characteristics of individual transport pathways and their time evolution are the outputs. The mass centers of His11, Glu89 and Arg90 were chosen as the starting point for tunnel searching. For each trajectory, 1000 snapshots were extracted from the last 100 ns simulation. The probe radius and the clustering threshold were set to 1.0 and 8.5, respectively. For the other parameters, their default values were taken throughout the calculations. The obtained tunnels were visualized with VMD (Supplementary Figure 1).

Generation of 2pg and 3pg complexes
The crystal structure 3FDZ [3] recorded the accurate atomic coordinates of both bacterial PGAM1 and 2,3-BPG. The crystal structure 4GPZ) [2] recorded only human PGAM1. To create human PGAM1 in complex with 2,3-BPG, one has to dock the 2,3-PBG molecule into the active site of 4GPZ. Here, the superimposition method was used to obtain the complex structure. The sequence identity between 3FDZ [3] and 4GPZ was 58.1% and the result of sequence alignment was shown in Supplementary Figure 2.
Then the two crystal structures were superimposed. After the superimposition, the 3FDZ protein structure was removed, leaving a complex formed by 2,3-BPG and 4GPZ, namely the desired human 2,3-BPG:PGAM1 complex. The 3PG:PGAM1 complex was obtained by the superposition of chain B of 3FDZ and chain A of 4GPZ (The 3FDZ crystal structure contained two chains: chain A was in complex with 2,3-BPG and chain B was in complex with 3PG). In contrast, the 2,3-BPG:PGAM1 complex was obtained by the superposition of chain A of 3FDZ and chain A of 4GPZ. In the protein data bank, we did not find the crystal structure of 2PG:PGAM1. Therefore, the 2PG molecule should be docked into the active site of PGAM1 in order to create the 2PG:PGAM1 complex structure. Here, AutoDock4.2 software package [4] was employed for the docking. The parameter set for docking calculations is summarized as follow: The 2PG molecule was docked into the active site of PGAM1 employing the Auto Dock tools (ADT). The nonpolar hydrogen atoms were removed and only the polar hydrogens were kept. The Gasteiger charges were added to the PGAM1 and 2PG molecules. A box size of 40×40×40 Å 3 with grid spacing 0.375 Å was defined around the binding site of PGAM1 so that it contained all the residues that are critical for interacting with 2PG. Here, Arg10, His11, Ser23 and Lys100 were selected as the active site residues interacting with 2PG. The grid map around the binding site of PGAM1 was generated by the probe atoms employing the Auto Grid program. Each grid in the map represents the potential energy of a probe atom in the presence of all the atoms of the receptor molecule. The Lamarkian Genetic Algorithm (LGA) was used for docking study. One hundred runs with 15000000 maximum evaluations and 270000 generations were used for docking simulation. The docking pose with the lowest binding energy −4.47 kcal mol −1 was chosen as the starting structure for the subsequent MD simulation.
To obtain reasonable structures before generating force field parameters, the three ligands 3PG, 2,3-BPG, and 2PG were optimized by employing Gaussian 09 software package [5] using DFT method (B3LYP) [6] under 6-31G basis set. The force field parameters for these three ligands were generated by Antechamber module in AMBER14 [7].

Binding free energy calculations
Based on the last 200ns simulation of each system, binding free energies were calculated by using MM/GBSA method. For each system, 160 frames were extracted from the last 200ns trajectory files. For all the 160 frames of each system, the binding free energy calculations all took the entropy into consideration. The entropy was calculated by using Normal Mode Analysis (nmode). The results are presented in Supplementary Figure 3.
As Supplementary Figure 3 shows, the binding free energies of all the wild type systems went through larger fluctuations than Y26-phospho systems. The binding free energies of Y26-phospho systems were quite stable during the last 200ns simulation. Besides, the binding free energies of Y26-phospho systems were all lower than those of wild type systems, which implied that all the three molecules (3PG, 2,3-BPG and 2PG) bound more tightly with PGAM1 phos than with PGAM1wt.