In order to ensure the reliability of activity values and reduce accidental errors, a set of 78 PGAM1 inhibitors were retrieved from different literature sources in terms of the same group (Wang et al., 2018a; Wang et al., 2018b; Huang et al., 2019a; Huang et al., 2019b). The molecular structure and experimental bioactivity of all chemicals are listed in Table 1. First, corresponding IC₅₀ values of experimental bioactivity expressed in nM were converted into negative logarithm (–lgIC₅₀) and acted as the dependent variable for the QSAR modeling. According to the diversity of the molecular structure and activities, all compounds were split into a training set and a test set at a ratio of approximately 4:1. Finally, 62 compounds were selected randomly as the training set and the remaining 16 compounds as the test set. The molecular structure of each compound was determined using ChemDraw 18.0 and then imported to SYBYL 6.9 (SYBYL, XX) to minimize the energy based on the Tripos force field with a convergence criterion of 0.01 kcal/mol. The Gasteiger–Hückel method was employed to calculate the partial atomic charges. Then, the multisearch strategy was performed to obtain the lowest energy conformation, and the lowest energy geometry after being filled with energy was reserved for alignment.

Molecular alignment in terms of the same structure is considered to be one of the most significant elements in the process of built 3D-QSAR modeling. Hence, molecular alignment based on the most active molecule, 35, was employed by atom-by-atom fits. After a common substructure is set, the dominant conformations of the remaining 77 compounds are selected for superimposition.

The 3D-QSAR model for the training set compound was built after alignment by using SYBYL 6.9 software. The CoMFA (Cramer et al., 1988b) and CoMSIA (Cramer et al., 1988b) are the most widely used methods for constructing 3D-QSAR. The CoMFA and CoMSIA descriptors were obtained by placing the superposed compound in a 3D cubic lattice with a grid spacing of 2 Å. Using the SP³ hybrid carbon as the probe atom, the Lennard–Jones and the coulomb potential were applied to obtain the steric field energy and electrostatic field energy of each lattice point. The contributions of the hydrogen bond acceptor field, hydrogen bond donor field, and hydrophobic field were calculated by the probe atom. The partial least squares method (Cramer et al., 1988a) was employed to deal with the linear correlation between the CoMFA and CoMSIA fields and biological activity. The cross-validation correlation coefficient (q²) and optimum number of components (N) were obtained using the leave-one-out method for cross-validation analysis. In addition, the r m 2 (Roy et al., 2013; Cardoso et al., 2016), r p r e d 2 , external standard deviation error of prediction (SDEPext), and applicability domain (Roy et al., 2015; de Assis et al., 2016) were also calculated to evaluate the performance of built models.

The predictive capabilities of built 3D-QSAR models were evaluated via the test set of 16 compounds. After all compounds were superimposed upon compound 35, the pIC₅₀ values of all compounds were estimated through the built CoMFA and CoMSIA models.

To obtain more accurate docking results, the resolution of all crystal structures of PGAM1 in complex with small molecules obtained from the RSCB Protein Data Bank (PDB) was compared, and 5Y35, with the best resolution of 1.99 Å, was preserved as the docking template. Subsequently, the Protein Preparation Wizard module within (Schrödinger, 2015) was utilized to preprocess the crystal structure, including adding hydrogens and side chains, deleting water molecules, and calculating partial charges and protonation states by using the OPLS2005 force field (Jorgensen et al., 1996). Then, a grid box centered at the native ligand with a similar size was produced to determine the binding pocket of PGAM1 by using the Grid Generation module of the Schrödinger package. All molecules were preprocessed using the LigPrep module implemented in the Schrödinger package, and the ionization states were calculated using Epik (Shelley et al., 2007) at pH = 7.0 ± 2.0. Finally, all chemicals were docked into the binding pocket of PGAM1 and evaluated using the standard precision (SP) mode of Glide. The scale factor was set at 0.8, and the partial charge intercept was set at 0.15. The 10,000 poses of each ligand during the initial docking phase were preserved for evaluation.

To obtain the structural basis and significant residues involved in the process of ligand binding, molecular dynamics simulations were employed in terms of the crystal structure of compounds 23 and 49 using Amber16 (Case et al., 2005). The general AMBER force field (GAFF) (Wang et al., 2004) was employed to parameterize the compounds, while the AMBER ff14SB force field (Maier et al., 2015) was employed for the PGAM1 structure. The partial charges of compounds were calculated by using the restrained electrostatic potential fitting procedure (Bayly et al., 1993; Cieplak et al., 1995; Fox and Kollman, 1998) based on the electrostatic potentials calculated using the Hartree–Fock (HF) method with the 6-31G* basis set in the Gaussian 09 package (Frisch et al., 2009). Then, the complex was solvated in a cubic box of TIP3P waters, with the solute 10 Å away from the water box boundary. After adding sodium ions to neutralize each system, the steepest descent method followed by the conjugate-gradient method were employed to minimize the system every 2,500 steps. Subsequently, each system was heated in the NVT ensemble from 0 to 300 K in 50 ps restraint on backbone atoms. The restraint force was gradually decreased from 5 to 0.1 kcal/(mol Ų) within 0.9 ns. Under a periodic boundary condition, 50 ns MD simulations were performed at 300 K and 1 atm without any restraint. The particle mesh Ewald method (Linse and Linse, 2014) was used to calculate the long-range electrostatic interactions, and the SHAKE method (Ryckaert et al., 1977) was employed to constrain all covalent bonds containing hydrogen atoms.

After the MD simulation finished, trajectories were dissected via the Cpptraj module (Roe and Cheatham, 2013) in AmberTools 16. First, the root mean square deviations (RMSDs) value was calculated in terms of the last 10 ns of each MD trajectory. Second, the molecular mechanics/generalized born surface area (MM/GBSA) approach (Massova and Kollman, 2000) was applied to calculate the binding free energy. After withdrawing a total of 2,500 snapshots, the MM/GBSA calculation was executed on each snapshot. The binding free energy (ΔG_bind) was calculated as follows (Hou et al., 2011; Sun et al., 2014): Δ G bind = G complex − ( G protein + G ligand ) where the energy term (G) is estimated as follows: G = E v d w + E e l e + G G B + G G B S U R

In the equations above, the Evdw, Eele, G _GB , and G _GBSUR represent van der Waals, electrostatic energy, the electrostatic contribution to the solvation free energy, and non-polar contribution to the solvation free energy, respectively. The changes of conformational entropy were ignored. Moreover, the total free energy was decomposed to each residue in PGAM1 to obtain the crucial residues contributed to the ligand binding process.

In the present study, a series of 78 PGAM1 inhibitors were obtained. The molecular structures and pIC₅₀ values of all molecules are listed in Table 1. The quality of molecular superposition is considered to be one of the important factors affecting 3D-QSAR prediction accuracy (Cho et al., 1996). On the basis of the structure and bioactivity of PGAM1 inhibitors, the compounds in the training set were aligned to compound 35, which had the highest activity based on the common substructure. It can be seen from Figure 1 that the common skeleton of all molecules is overlapped. However, the side chains of several compounds surround the common skeleton due to the large difference. Then, the 3D-QSAR models of CoMFA and CoSIA were successfully developed.

To examine the predictive ability and reliability of the built model, q² and r² were applied to evaluate the predictive power of the built 3D-QSAR model, r², F, and SEE values were employed to assess the reliability of the model, and r m 2 , r p r e d 2 , and SDEP_ext values were utilized for external validation of the model. Table 2 lists the classical parameter statistics of CoMFA and CoMSIA models. In general, r² > 0.7 and q², r m 2 , r p r e d 2 >0.5 are necessary for a good model (Pratim Roy et al., 2009). As shown in Table 2, the values of q², N, SEE, r², r m 2 , r p r e d 2 , SDEP_ext, and F are 0.81, 6, 0.106, 0.97, 0.78, 0.89, 0.22, and 258.06, respectively. The results show that the built CoMFA model exhibits a good stability and predictive ability. The contribution of the steric field and the electrostatic field is 81 and 19%, respectively, indicating that the biological activity of compounds is more affected by the steric field. In addition, the predicted activity of the new chemical is only valid when the predicted compound falls within the applicability domain of the developed model (Roy et al., 2015). The calculated results show that all compounds are within the application domain of the built CoMFA model, so this prediction result is reliable.

Different field combinations of CoMSIA models were constructed, and it had been proved that CoMSIA-SEHA is the best model. Based on this model, the values of q², N, SEE, r², r m 2 , r p r e d 2 , SDEP_ext, and F are 0.82, 6, 0.11, 0.96, 0.79, 0.89, 0.23, and 228.71, respectively. In this model, the contribution of the steric field is 20%, that of the electrostatic field is 22%, that of the hydrophobic field is 40%, and that of the hydrogen bond acceptor field is 18%, respectively. The results show that the hydrophobic field has a greater effect on the bioactivity of the PGAM1 inhibitors. The calculation results of the application domain show that almost all the compounds are within the application domain of the CoSIA model, except for compound 24 with an S_new of 3.87 and compound 25 with an S_new of 4.06. By analyzing the descriptors in CoMSIA, we found that compounds 24 and 25 have the largest electrostatic field contribution. The experimental and predicted values of the biological activity of the training set and the test set in the established CoMFA and CoMSIA models are shown in Table 3.

The scatter plot of the experimental and predicted values of the studied PGAM1 inhibitor is shown in Figure 2. It can be seen from Figure 2 that the experimental and predicted bioactivity values of all molecules are distributed around the Y = X equation, indicating that the predicted values are in good accord with the experimental values, which further demonstrates that the model has good predictive ability.

The structure–activity relationships between PGAM1 inhibitors and activity can be well demonstrated by using 3D contour maps to display the QSAR equation. The field type Stdev* Coeff was used to generate 3D contour maps. As shown in Figures 3, 4, compound 35 with the best anti-PGAM1 activity was selected as the template compound to dissect the results of CoMFA and CoMSIA models.

The contour map of the steric field of CoMFA is shown in Figure 3A, and the effect of the steric field on the activity is shown in green and yellow. The presence of green regions around the molecule indicates that the group with a large connecting space contributes to increasing the activity of the compound, while the presence of yellow regions indicates that the group with a large connecting space may decrease the activity of the compound. As can be seen from Figure 3A, there is a green region distributed on the R₁ substituent, so the introduction of a slightly larger volume of groups at the R₁ substituent site is conducive to the improvement of the activity of the compound. For example, compound 22 (pIC₅₀ = 5.82) with a benzene ring was significantly higher than compound 19 (pIC₅₀ = 5.27) in bioactivity. The contour map of the electrostatic field of CoMFA is shown in Figure 3B, and the effect of the electrostatic field on the activity is shown in blue and red. The blue regions around the molecule indicate that the connection of the electron-donating group is beneficial to the improvement of the activity of the compound, while the red regions indicate that the connection of the electron-withdrawing group is beneficial to the improvement of the activity of the compound. From Figure 3B, we can see that the connection of electron-withdrawing groups near the R₁ substituent is conducive to improving the activity of the compound, so it can explain how the activity of compound 22 (pIC₅₀ = 5.82) is higher than that of compound 19 (pIC₅₀ = 5.27). There is a blue region around the R₂ substituents of anthraquinone, where the introduction of electron groups is beneficial. For example, the bioactivity of compound 72 (pIC₅₀ = 5.77) with a hydroxyl group was significantly higher than that of compound 8 (pIC₅₀ = 5.22).

The contour map of the steric field (Figure 4A) and the electrostatic field (Figure 4B) of the CoMSIA is very similar to the CoMFA model, so they will not be explained here. The contour map of the hydrophobic field of the CoMSIA model is shown in Figure 4C. The cyan regions represent how the introduction of the hydrophobic group is favorable to the activity, while the white regions represent how the introduction of the hydrophilic group is favorable to the activity. There is a cyan region near the R₁ substituent, indicating that the introduction of the hydrophobic group is very helpful to the improvement of the activity. Therefore, the biological activity of compound 22 (pIC₅₀ = 5.82) is higher than that of compound 19 (pIC₅₀ = 5.27). The contour map of the hydrogen bond receptor field of CoMSIA is shown in Figure 4D. The orange area is where the hydrogen bond acceptor group is conducive to the activity of the compound, and the purple area is where the hydrogen bond donor group is conducive to the activity of the compound. As shown in Figure 4D, there are purple regions with substituents of R₆ and R₂, where hydrogen bond donors can be imported to improve the anti-PGAM1 activity of the chemical. Moreover, a large purple region is near the nitrogen atom on the amino group, suggesting that the group may be a hydrogen bond donor.

Based on the outcome of CoMFA and CoMSIA analysis, we obtained the structure–activity relationship diagram of anthraquinone compounds (see Figure 5). The introduction of hydrogen bond donors in Region A is beneficial to improving the activity of the compounds, such as the carbonyl group. The group with a large space in Region B is conducive to the activity of the compounds, such as biphenyl or p-cyclohexylbenzene (Huang et al., 2019b). The introduction of the hydrophilic group in Region C is conducive to the activity, such as hydroxyl groups (Wang et al., 2018a). The group with a small space in Region D can improve the activity of the compound, such as hydrogen.

The molecular docking method was employed to interpret the 3D-QSAR result and study the structural basis between PGAM1 and inhibitors. First, the reliability of the glide docking algorithm with the SP mode was evaluated by redocking analysis. It can be seen from Figure 6 that the redocking conformations of the molecule are well superimposed with the initial structure in PGAM1 protein. The RMSD value between docking conformation and native conformation is 0.005Å. The results suggest that the glide algorithm exhibits a good performance for the PGAM1 protein, which can reproduce the binding pose of the native ligand. Subsequently, all chemicals were docked into the binding site of PGAM1. However, we discover that the docking scores of these compounds are not correlated with the inhibitory activity, and the r² of pIC₅₀ vs. the docking score is 0.051, which demonstrates the fact that glide docking is not appropriate for all compounds. We speculate that one of the most important reasons is that 3-PG plays an important role in the process of compounds binding to PGAM1, and the glide scoring function currently used is not suitable for this system. In addition, because PGAM1 catalyzes the conversion of 3-PG to 2-PG in the physiological process, the current docking simulation methods cannot completely simulate this process. Therefore, the docking score and activity do not show a correlation.

In order to further analyze the atomic details of the interaction between small molecules and PGAM1, molecular dynamics simulations were employed based on the co-crystal complex of compounds 23 (PDB ID: 5Y35) and 49 (PDB ID: 6ISN) using Amber 16, respectively. 50 ns simulation was performed for each complex. The RMSD plots of Cα, residues within the range of ligand 5Å, ligand, and 3-PG for complexes were shown in Figure 7. By monitoring the fluctuation of RMSDs, it can be found that the RMSD fluctuation of each system after 20 ns are all within the range of 2Å. Moreover, the fluctuation of binding free energy over time was also monitored. As shown in Supplementary Figure S1, binding free energy of each system fluctuates around 30 kcal/mol after 35 ns. In summary, these results indicate that the two systems finally reached a stable state.

During the process of small molecules binding to PGAM1, the hydrogen bond plays an important role as one of the most important non-bonding interactions. In order to explore the interaction between small molecules and PGAM1, the changes of the hydrogen bond between each residue of PGAM1 and the inhibitor were also monitored. The fraction of the hydrogen bond is greater than 10% as listed in Table 4. The results show that two hydrogen bonds formed between compounds 23 and 49 and Arg116, and the total occupancies are 180.12% and 38.48%, respectively. The results indicate that the hydrogen bonds formed between Arg116 of PGAM1 and inhibitors play a remarkable role in the binding of molecules. Besides, another hydrogen bond is also formed between compound 23 and Arg90 with the occupancy of 12.14%. It is precisely because the small molecules form hydrogen bonds with Arg116 and Arg90 to fix the anthraquinone skeleton of the compounds that compounds 23 and 49 are stably binding with PGAM1.

The binding free energy is used as a reference standard for evaluating the activity of molecules. It is generally believed that the lower the binding value, the more stable the complex formed by the protein and the small molecule. To evaluate the binding affinity of each complex, the MM/GBSA method was performed to calculate the binding free energy of inhibitors. It can be seen from Table 5 that the binding free energy of compounds 23 and 49 are −27.40 kcal/mol and −27.85 kcal/mol, respectively, which are consistent with their biological activities. Among them, van der Waals energies (ΔE_vdw) are −38.68 kcal/mol and −41.63 kcal/mol, respectively, and their values are much lower than other energy terms, indicating that hydrophobic interaction is the major contributor to the ligand binding process. In addition, electrostatic energy (ΔE_ele) also contributes significantly to the binding free energy, which indicates that electrostatic interaction also plays a vital role in ligand binding. It is worth noting that the polar contribution (ΔG_GB) is not conducive to ligand binding, which may be attributed to the large size of the binding pocket and the exposure of the hydrophobic ligand to the solvent.

To further confirm the key residues referred to in the ligand binding process, MM/GBSA calculation was performed to decompose the binding free energy into inhibitor–residue pairs. It can be seen from Figure 8 that the primary residues with binding free energy less than −1 kcal/mol contributing to the ligand binding are F22, K100, V112, W115, and R116. In order to further observe the orientation of compounds and the position of the key residues, we extracted the average structure (see Figure 9). It can be seen from Figure 9 that compounds 23 and 49 adopt a similar binding pose, which is surrounded by those critical residues. Compound 23 forms three hydrogen bonds with R90, W115, and R116. Among the three of them, R90 and R116 show higher fraction in hydrogen bond analysis, while the bond length of W115 is 3.4 Å due to weak potency. For compound 49, there is no hydrogen bond formed between compound 49 and key residues, which may be due to the low occupancy.

According to the structure–activity relationships obtained from CoMFA and CoMSIA models, seven molecules with the anthraquinone skeleton were designed as potential PGAM1 inhibitors by introducing new substituents at different positions of compound 35 (see Table 6). Compounds 79 and 80 were designed by adding the hydrogen bond donor in the R₆ position to form the key hydrogen bond. Compounds 81, 82, and 83 were designed by introducing the substituent in the R₁ position to increase volume. Based on the contribution of the steric and hydrogen bond donor, compounds 84 and 85 were designed. The pIC₅₀ values of designed compounds were predicted by built CoMFA and CoMSIA models. As shown in Table 6, all of the designed compounds exhibit better inhibitory activity targeting PGAM1 than compound 35, and the predictive values are in accordance with the summarized structure–activity relationships.