3D-QSAR, Molecular Docking, and MD Simulations of Anthraquinone Derivatives as PGAM1 Inhibitors

PGAM1 is overexpressed in a wide range of cancers, thereby promoting cancer cell proliferation and tumor growth, so it is gradually becoming an attractive target. Recently, a series of inhibitors with various structures targeting PGAM1 have been reported, particularly anthraquinone derivatives. In present study, the structure–activity relationships and binding mode of a series of anthraquinone derivatives were probed using three-dimensional quantitative structure–activity relationships (3D-QSAR), molecular docking, and molecular dynamics (MD) simulations. Comparative molecular field analysis (CoMFA, r² = 0.97, q² = 0.81) and comparative molecular similarity indices analysis (CoMSIA, r² = 0.96, q² = 0.82) techniques were performed to produce 3D-QSAR models, which demonstrated satisfactory results, especially for the good predictive abilities. In addition, molecular dynamics (MD) simulations technology was employed to understand the key residues and the dominated interaction between PGAM1 and inhibitors. The decomposition of binding free energy indicated that the residues of F22, K100, V112, W115, and R116 play a vital role during the ligand binding process. The hydrogen bond analysis showed that R90, W115, and R116 form stable hydrogen bonds with PGAM1 inhibitors. Based on the above results, 7 anthraquinone compounds were designed and exhibited the expected predictive activity. The study explored the structure–activity relationships of anthraquinone compounds through 3D-QSAR and molecular dynamics simulations and provided theoretical guidance for the rational design of new anthraquinone derivatives as PGAM1 inhibitors.

Introduction

Reprogramming energy metabolism has been regarded as one of the 10 essential hallmarks of cancer cells (Hanahan and Weinberg, 2011), which was called the “Warburg effect.” In 1924, Warburg found that cancer cells are more likely to metabolize glucose by means of aerobic glycolysis instead of oxidative phosphorylation as in normal cells (Wang et al., 2018a; Huang et al., 2019b). Cancer metabolic reprogramming is the performance of adapting to the environment during tumor formation or metastasis. More and more scientists are focusing on the pivotal enzymes in the metabolic reprogramming of cancer cells in order to find new cancer treatment targets (Wang et al., 2018b).

Phosphoglycerate mutase 1 (PGAM1) is a key enzyme that catalyzes the invertible conversion of 3-phosphoglycerate (3-PG) and 2-phosphoglycerate (2-PG) during the process of glycolysis (Fothergill-Gilmore and Watson, 1989). Recent studies have proven that once the expression of PGAM1 is upregulated, it will promote tumor cell proliferation and tumor growth in coordination with glycolysis and biosynthesis (Hitosugi et al., 2012). PGAM1 regulates the proliferation of cancer cells in term of biosynthesis regulation, partly by regulating intracellular levels of its product 2-PG and 3-PG (Hitosugi et al., 2012). In the oxidative pentose phosphate pathway (PPP), 3-PG inhibits 6-phosphogluconate dehydrogenase after binding, while 2-PG feedback control of the levels of through activates 3-phosphoglycerate dehydrogenase. In addition, PGAM1 is overexpressed in multiple cancers (Li and Liu, 2020), including ovarian cancer (Zhang et al., 2020), non–small-cell lung cancer (NSCLC) (Li et al., 2020), colorectal cancer (Liu et al., 2008; Lei et al., 2011), pancreatic ductal adenocarcinoma (PDAC) (Liu et al., 2018), prostate cancer (PCa) (Wen et al., 2018), and glioma (Xu et al., 2016). Particularly, high expression of PGAM1 was associated with poor prognosis in NSCLC patients (Sun et al., 2018; Li et al., 2020). Downregulation of the expression of PGAM1 or suppression of its metabolic activity will lead to weakened cell proliferation and tumor growth (Hitosugi et al., 2012; Peng et al., 2016; Liu et al., 2018). Thus, PGAM1 is considered to be an emerging target for cancer treatment.

Due to the important role of PGAM1 in the occurrence and development of tumors, many researchers have focused on the discovery and characterization of small molecules that can target and modulate the metabolic activity of PGAM1 (Huang et al., 2019a). MJE3 was first revealed as a covalent PGAM1 inhibitor on Lys 100 by the Cravatt group in 2005 (Evans et al., 2005). (-)-Epigallocatechin-3-gallate (EGCG) is a natural product extracted from green tea, which was first discovered as a non-substrate competitive PGAM1 inhibitor with potent inhibition activity against PGAM1 (Li et al., 2017). Anthraquinone derivatives PGMI-004A (Hitosugi et al., 2012) and xanthone derivatives (Wang et al., 2018b) were identified as allosteric PGAM1 inhibitors by the Zhou group, which exhibited moderate inhibition activity on PGAM1. As another anthraquinone derivative, HKB99 was identified to allosterically obstruct the activation of PGAM1, thereby affecting its catalytic activity and the intermolecular interaction of ACTA2 (Huang et al., 2019c; Liang et al., 2021). Based on the excellent anticancer activity of PGMI-004A and HKB99, new small molecules with the anthraquinone core have been synthesized, which may have similar mechanisms of action and therapeutic potential. Therefore, the design and development of novel small molecules with an anthraquinone core targeting PGAM1 may prove to be an effective strategy for the treatment of cancer cells.

Computer-aided drug design is an effective tool in the drug discovery and design process. It can not only be used to predict the activity of small molecules, explain the action mechanism, and provide guidance for the design of more effective drug molecules but also reduce the consumption of manpower and material resources (Jorgensen, 2004). To elucidate the structure–activity relationships and provide optimization guidance for anthraquinone derivatives, 62 collected compounds were employed to construct 3D-QSAR models using CoMFA and CoMSIA methods. According to the contour maps by 3D-QSAR and the crucial residues by MD simulations, 7 compounds with high predictive activity were designed. This study will provide a valuable theoretical basis for the activity prediction and structural modification of targeted PGAM1 inhibitors containing anthraquinone structures.

Materials and Methods

Data Sets and Preparation

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.

TABLE 1

TABLE 1. Structure and corresponding activity data of reported PGAM1 inhibitors.

Molecular 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.

Construction of CoMFA and CoMSIA Models

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_{pred}^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.

Evaluation of the 3D-QSAR 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.

Molecular Docking

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.

Molecular Dynamics Simulations

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.

Trajectory Analysis

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):

$\Delta G_{\mathrm{bind}}=G_{\mathrm{complex}}-\left(G_{\mathrm{protein}}+G_{\mathrm{ligand}}\right)$

where the energy term (G) is estimated as follows:

$\mathrm{G}=E_{vdw}+E_{ele}+G_{GB}+G_{GBSUR}$

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.

Results and Discussion

CoMFA and CoMSIA Models

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.

FIGURE 1

FIGURE 1. Structural alignment of all the molecules in the training set based on the common substructure of compound 35.

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_{pred}^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_{pred}^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_{pred}^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.

TABLE 2

TABLE 2. Summary of CoMFA and CoMSIA models.

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_{pred}^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.

TABLE 3

TABLE 3. Experimental pIC₅₀ (Exp.), predicted pIC₅₀ (Pred.), and corresponding residuals (Res.) of the anthraquinone derivatives.

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.

FIGURE 2

FIGURE 2. Scatter plot of experimental and predicted bioactivity values (pIC₅₀)of the CoMFA (A) and CoMSIA models (B), respectively.

Contour Maps Analysis of CoMFA and CoMSIA

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.

FIGURE 3

FIGURE 3. Steric contour map (A) and electrostatic contour map (B) of the CoMFA model based on molecule 35. Green regions represent bulky groups that increase anti-PGAM1 activity, while yellow regions represent sterically unfavored regions. Blue regions show where positive groups are beneficial for increasing anti-PGAM1 acitivity, and red regions show where negative groups are favored.

FIGURE 4

FIGURE 4. Steric contour map (A), electrostatic contour map (B), hydrophobic contour map (C), and hydrogen bond acceptor contour map (D) of the CoMSIA model based on molecule 35. Green regions are sterically favored regions, while yellow regions are sterically unfavored regions. Blue regions are where electron-donating groups are favored, and red regions are where electron-withdrawing groups are favored. The cyan regions are where the hydrophobic group is favorable to the activity, while the white regions are where the hydrophilic group is favorable to the activity.

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.

FIGURE 5

FIGURE 5. Structure–activity relationship diagram of anthraquinone PGAM1 inhibitors.

Molecular Docking Analysis

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.

FIGURE 6

FIGURE 6. Surface of PGAM1 and docking pose of the native ligand based on the alignment. The yellow and cyan carbon atoms represent the native ligand and the docking pose, respectively.

Molecular Dynamics Simulations

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.

FIGURE 7

FIGURE 7. Fluctuation of RMSD values for two complexes during 50 ns MD simulation.

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.

TABLE 4

TABLE 4. Changes of the hydrogen bond over the MD simulations.

Binding Free Energy Calculation

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.

TABLE 5

TABLE 5. Calculated binding energy (kcal/mol) of inhibitor binding to PGAM1.

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.

FIGURE 8

FIGURE 8. Binding free energy decomposition plots for the two systems.

FIGURE 9

FIGURE 9. Average structures of PGAM1 with compounds 23 (A) and 49 (B). The bonds of residues and ligands are represented in stick, and the carbon atoms of compound 23, compound 49, and residues are represented in yellow, cyan, and white, respectively. The red dotted line represents the hydrogen bond.

Design New PGAM1 Inhibitors

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.

TABLE 6

TABLE 6. Newly designed PGAM1 inhibitors and the corresponding predicted activity value.

Conclusion

In the present study, a combined strategy of 3D-QSAR, molecular docking, and molecular dynamics simulations was applied to explore the structure–activity relationships of anthraquinone analogs. The built CoMFA (q² = 0.81, r² = 0.97, $r_m^2$ = 0.78, $r_{pred}^2$ = 0.89) and CoMSIA (q² = 0.82, r² = 0.96, $r_m^2$ = 0.79, $r_{pred}^2$ = 0.89) models have achieved satisfactory results in terms of the statistical results. The results show that the built models have good internal and external predictive power. The acquired contour maps elaborate the structure–activity relationships of anthraquinone derivatives and successfully predict the activity of the test set. According to the results of contour maps, the introduction of hydrogen bond donors in Region A, the group with a large space in Region B, the hydrophilic group in Region C, and the group with a small space in Region D could improve the activity of the compounds. The calculated results of binding free energy suggest that van der Waals interaction is the major contributor to the ligand binding process. The decomposition binding free energy and hydrogen bond show that small molecules with the anthraquinone core mainly interact with F22, R90, K100, V112, W115, and R116 of PGAM1. Based on these findings, 7 new compounds with the anthraquinone core were designed, and the predicted results show that all of the designed compounds exhibit great inhibitory activity against PGAM1. The constructed 3D-QSAR model will provide theoretical guidance for improving the activity of anthraquinone derivatives and help to develop inhibitors with potent anti-PGAM1 activity.

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

YW, HG, and XY designed this study. YG, SQ, and ZL carried out computational modeling and data analysis and wrote the manuscript. RJ, YT, and EL revised this manuscript. All authors have read and approved the final manuscript.

Funding

This study was supported by the Natural Science Foundation of Shaanxi Province (Project No. 2021JQ-734), the National Natural Science Foundation of China (Project No. 82003653), Shaanxi University of Chinese Medicine (Project No. 2020XG01), and the Subject Innovation Team of Shaanxi University of Chinese Medicine (Project No. 2019-PY02).

Conflict of Interest

The 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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

I wish to thank “Xternal Validation Plus” (https://dtclab.webs.com/software-tools) and “application domain” (https://dtclab.webs.com/software-tools) tools for calculating the external validation index.

Supplementary Material

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

References

Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: the RESP Model. J. Phys. Chem. 97 (40), 10269–10280. doi:10.1021/j100142a004

CrossRef Full Text | Google Scholar

Cardoso, G. G., Maria, D. A. T., Assis Letícia, C., Castro, R. T., and Ferreira, D. C. E. F. (2016). Quantitative Structure-Activity Relationship Studies for Potential Rho-Associated Protein Kinase Inhibitors. J. Chem. 2016, 9198582. doi:10.1155/2016/9198582

CrossRef Full Text | Google Scholar

Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., et al. (2005). The Amber Biomolecular Simulation Programs. J. Comput. Chem. 26 (16), 1668–1688. doi:10.1002/jcc.20290

CrossRef Full Text | Google Scholar

Cho, S. J., Garsia, M. L., Bier, J., and Tropsha, A. (1996). Structure-Based Alignment and Comparative Molecular Field Analysis of Acetylcholinesterase Inhibitors. J. Med. Chem. 39 (26), 5064–5071. doi:10.1021/jm950771r

CrossRef Full Text | Google Scholar

Cieplak, P., Cornell, W. D., Bayly, C., and Kollman, P. A. (1995). Application of the Multimolecule and Multiconformational RESP Methodology to Biopolymers: Charge Derivation for DNA, RNA, and Proteins. J. Comput. Chem. 16 (11), 1357–1377. doi:10.1002/jcc.540161106

CrossRef Full Text | Google Scholar

Cramer, R. D., Patterson, D. E., and Bunce, J. D. (1988b). Comparative Molecular Field Analysis (CoMFA). 1. Effect of Shape on Binding of Steroids to Carrier Proteins. J. Am. Chem. Soc. 110 (18), 5959–5967. doi:10.1021/ja00226a005

CrossRef Full Text | Google Scholar

Cramer, R. D., Bunce, J. D., Patterson, D. E., and Frank, I. E. (1988a). Crossvalidation, Bootstrapping, and Partial Least Squares Compared with Multiple Regression in Conventional QSAR Studies. Quant. Struct.-Act. Relat. 7 (1), 18–25. doi:10.1002/qsar.19880070105

CrossRef Full Text | Google Scholar

de Assis, T. M., Gajo, G. C., de Assis, L. C., Garcia, L. S., Silva, D. R., Ramalho, T. C., et al. (2016). QSAR Models Guided by Molecular Dynamics Applied to Human Glucokinase Activators. Chem. Biol. Drug Des. 87 (3), 455–466. doi:10.1111/cbdd.12683

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, M. J., Saghatelian, A., Sorensen, E. J., and Cravatt, B. F. (2005). Target Discovery in Small-Molecule Cell-Based Screens by In Situ Proteome Reactivity Profiling. Nat. Biotechnol. 23 (10), 1303–1307. doi:10.1038/nbt1149

PubMed Abstract | CrossRef Full Text | Google Scholar

Fothergill-Gilmore, L. A., and Watson, H. C. (1989). The Phosphoglycerate Mutases. Adv. Enzymol. Relat. Areas Mol. Biol. 62, 227–313. doi:10.1002/9780470123089.ch6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, T., and Kollman, P. A. (1998). Application of the RESP Methodology in the Parametrization of Organic Solvents. J. Phys. Chem. B 102 (41), 8070–8079. doi:10.1021/jp9717655

CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., and Fox, D. J. (2009). Gaussian 09. Wallingford: Gaussian, Inc.

Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: The Next Generation. Cell 144 (5), 646–674. doi:10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hitosugi, T., Zhou, L., Elf, S., Fan, J., Kang, H. B., Seo, J. H., et al. (2012). Phosphoglycerate Mutase 1 Coordinates Glycolysis and Biosynthesis to Promote Tumor Growth. Cancer Cell 22 (5), 585–600. doi:10.1016/j.ccr.2012.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, T., Wang, J., Li, Y., and Wang, W. (2011). Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations. J. Chem. Inf. Model. 51 (1), 69–82. doi:10.1021/ci100275a

CrossRef Full Text | Google Scholar

Huang, K., Jiang, L., Li, H., Ye, D., and Zhou, L. (2019a). Development of Anthraquinone Analogues as Phosphoglycerate Mutase 1 Inhibitors. Molecules 24 (5), 845. doi:10.3390/molecules24050845

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, K., Jiang, L., Liang, R., Li, H., Ruan, X., Shan, C., et al. (2019b). Synthesis and Biological Evaluation of Anthraquinone Derivatives as Allosteric Phosphoglycerate Mutase 1 Inhibitors for Cancer Treatment. Eur. J. Med. Chem. 168, 45–57. doi:10.1016/j.ejmech.2019.01.085

CrossRef Full Text | Google Scholar

Huang, K., Liang, Q., Zhou, Y., Jiang, L. L., Gu, W. M., Luo, M. Y., et al. (2019c). A Novel Allosteric Inhibitor of Phosphoglycerate Mutase 1 Suppresses Growth and Metastasis of Non-small-cell Lung Cancer. Cell Metab 30 (6), 1107–e1108. doi:10.1016/j.cmet.2019.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L. (2004). The many Roles of Computation in Drug Discovery. Science 303 (5665), 1813–1818. doi:10.1126/science.1096361

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996). Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 118 (45), 11225–11236. doi:10.1021/ja9621760

CrossRef Full Text | Google Scholar

Lei, Y., Huang, K., Gao, C., Lau, Q. C., Pan, H., Xie, K., et al. (2011). Proteomics Identification of ITGB3 as a Key Regulator in Reactive Oxygen Species-Induced Migration and Invasion of Colorectal Cancer Cells. Mol. Cel Proteomics 10 (10), M110–M005397. doi:10.1074/mcp.M110.005397

CrossRef Full Text | Google Scholar

Li, F., Yang, H., Kong, T., Chen, S., Li, P., Chen, L., et al. (2020). PGAM1, Regulated by miR-3614-5p, Functions as an Oncogene by Activating Transforming Growth Factor-β (TGF-β) Signaling in the Progression of Non-small Cell Lung Carcinoma. Cell Death Dis 11 (8), 710. doi:10.1038/s41419-020-02900-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, N., and Liu, X. (2020). Phosphoglycerate Mutase 1: Its Glycolytic and Non-glycolytic Roles in Tumor Malignant Behaviors and Potential Therapeutic Significance. Onco Targets Ther. 13, 1787–1795. doi:10.2147/OTT.S238920

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Tang, S., Wang, Q. Q., Leung, E. L., Jin, H., Huang, Y., et al. (2017). Identification of Epigallocatechin-3- Gallate as an Inhibitor of Phosphoglycerate Mutase 1. Front. Pharmacol. 8, 325. doi:10.1073/pnas.191455711610.3389/fphar.2017.00325

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Q., Gu, W. M., Huang, K., Luo, M. Y., Zou, J. H., Zhuang, G. L., et al. (2021). HKB99, an Allosteric Inhibitor of Phosphoglycerate Mutase 1, Suppresses Invasive Pseudopodia Formation and Upregulates Plasminogen Activator Inhibitor-2 in Erlotinib-Resistant Non-small Cell Lung Cancer Cells. Acta Pharmacol. Sin 42 (1), 115–119. doi:10.1038/s41401-020-0399-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Linse, B., and Linse, P. (2014). Tuning the Smooth Particle Mesh Ewald Sum: Application on Ionic Solutions and Dipolar Fluids. J. Chem. Phys. 141 (18), 184114. doi:10.1063/1.4901119

CrossRef Full Text | Google Scholar

Liu, L., Wang, S., Zhang, Q., and Ding, Y. (2008). Identification of Potential Genes/proteins Regulated by Tiam1 in Colorectal Cancer by Microarray Analysis and Proteome Analysis. Cell Biol Int 32 (10), 1215–1222. doi:10.1016/j.cellbi.2008.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Weng, Y., Liu, P., Sui, Z., Zhou, L., Huang, Y., et al. (2018). Identification of PGAM1 as a Putative Therapeutic Target for Pancreatic Ductal Adenocarcinoma Metastasis Using Quantitative Proteomics. Onco Targets Ther. 11, 3345–3357. doi:10.2147/OTT.S162470

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 (8), 3696–3713. doi:10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Massova, I., and Kollman, P. A. (2000). Combined Molecular Mechanical and Continuum Solvent Approach (MM-PBSA/GBSA) to Predict Ligand Binding. Perspect. Drug Discov. Des. 18 (1), 113–135. doi:10.1023/A:1008763014207

CrossRef Full Text | Google Scholar

Peng, X. C., Gong, F. M., Chen, Y., Qiu, M., Cheng, K., Tang, J., et al. (2016). Proteomics Identification of PGAM1 as a Potential Therapeutic Target for Urothelial Bladder Cancer. J. Proteomics 132, 85–92. doi:10.1016/j.jprot.2015.11.027

CrossRef Full Text | Google Scholar

Pratim Roy, P., Paul, S., Mitra, I., and Roy, K. (2009). On Two Novel Parameters for Validation of Predictive QSAR Models. Molecules 14 (5), 1660–1701. doi:10.3390/molecules14051660

PubMed Abstract | CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. (2013). PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theor. Comput 9 (7), 3084–3095. doi:10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, K., Chakraborty, P., Mitra, I., Ojha, P. K., Kar, S., and Das, R. N. (2013). Some Case Studies on Application of "r(m)2" Metrics for Judging Quality of Quantitative Structure-Activity Relationship Predictions: Emphasis on Scaling of Response Data. J. Comput. Chem. 34 (12), 1071–1082. doi:10.1002/jcc.23231

CrossRef Full Text | Google Scholar

Roy, K., Kar, S., and Ambure, P. (2015). On a Simple Approach for Determining Applicability Domain of QSAR Models. Chemometrics Intell. Lab. Syst. 145, 22–29. doi:10.1016/j.chemolab.2015.04.013

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

Schrödinger, L. (2015). New York, NY: Schrödinger, LLC.

Shelley, J. C., Cholleti, A., Frye, L. L., Greenwood, J. R., Timlin, M. R., and Uchimaya, M. (2007). Epik: a Software Program for pK( a ) Prediction and Protonation State Generation for Drug-like Molecules. J. Comput. Aided Mol. Des. 21 (12), 681–691. doi:10.1007/s10822-007-9133-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, H., Li, Y., Tian, S., Xu, L., and Hou, T. (2014). Assessing the Performance of MM/PBSA and MM/GBSA Methods. 4. Accuracies of MM/PBSA and MM/GBSA Methodologies Evaluated by Various Simulation Protocols Using PDBbind Data Set. Phys. Chem. Chem. Phys. 16 (31), 16719–16729. doi:10.1039/C4CP01388C

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Q., Li, S., Wang, Y., Peng, H., Zhang, X., Zheng, Y., et al. (2018). Phosphoglyceric Acid Mutase-1 Contributes to Oncogenic mTOR-Mediated Tumor Growth and Confers Non-small Cell Lung Cancer Patients with Poor Prognosis. Cell Death Differ 25 (6), 1160–1173. doi:10.1038/s41418-017-0034-y

PubMed Abstract | CrossRef Full Text | Google Scholar

SYBYL, V. (XX). SYBYL. St. Louis, MO: Tripos Inc..

Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and Testing of a General Amber Force Field. J. Comput. Chem. 25 (9), 1157–1174. doi:10.1002/jcc.20035

CrossRef Full Text | Google Scholar

Wang, P., Jiang, L., Cao, Y., Ye, D., and Zhou, L. (2018a). The Design and Synthesis of N-Xanthone Benzenesulfonamides as Novel Phosphoglycerate Mutase 1 (PGAM1) Inhibitors. Molecules 23 (6), 1396. doi:10.3390/molecules23061396

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Jiang, L., Cao, Y., Zhang, X., Chen, B., Zhang, S., et al. (2018b). Xanthone Derivatives as Phosphoglycerate Mutase 1 Inhibitors: Design, Synthesis, and Biological Evaluation. Bioorg. Med. Chem. 26 (8), 1961–1970. doi:10.1016/j.bmc.2018.02.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, Y. A., Zhou, B. W., Lv, D. J., Shu, F. P., Song, X. L., Huang, B., et al. (2018). Phosphoglycerate Mutase 1 Knockdown Inhibits Prostate Cancer Cell Growth, Migration, and Invasion. Asian J. Androl. 20 (2), 178–183. doi:10.4103/aja.aja_57_17

CrossRef Full Text | Google Scholar

Xu, Z., Gong, J., Wang, C., Wang, Y., Song, Y., Xu, W., et al. (2016). The Diagnostic Value and Functional Roles of Phosphoglycerate Mutase 1 in Glioma. Oncol. Rep. 36 (4), 2236–2244. doi:10.3892/or.2016.5046

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Li, Y., Zhao, W., Liu, G., and Yang, Q. (2020). Circ-PGAM1 Promotes Malignant Progression of Epithelial Ovarian Cancer through Regulation of the miR-542-3p/CDC5L/PEAK1 Pathway. Cancer Med. 9 (10), 3500–3521. doi:10.1002/cam4.2929

PubMed Abstract | CrossRef Full Text | Google Scholar