Identification of novel PHGDH inhibitors based on computational investigation: an all-in-one combination strategy to develop potential anti-cancer candidates

Objective: Biological studies have elucidated that phosphoglycerate dehydrogenase (PHGDH) is the rate-limiting enzyme in the serine synthesis pathway in humans that is abnormally expressed in numerous cancers. Inhibition of the PHGDH activity is thought to be an attractive approach for novel anti-cancer therapy. The development of structurally diverse novel PHGDH inhibitors with high efficiency and low toxicity is a promising drug discovery strategy.

Methods: A ligand-based 3D-QSAR pharmacophore model was developed using the HypoGen algorithm methodology of Discovery Studio. The selected pharmacophore model was further validated by test set validation, cost analysis, and Fischer randomization validation and was then used as a 3D query to screen compound libraries with various chemical scaffolds. The estimated activity, drug-likeness, molecular docking, growing scaffold, and molecular dynamics simulation processes were applied in combination to reduce the number of virtual hits.

Results: The potential candidates against PHGDH were screened based on estimated activity, docking scores, predictive absorption, distribution, metabolism, excretion, and toxicity (ADME/T) properties, and molecular dynamics simulation.

Conclusion: Finally, an all-in-one combination was employed successfully to design and develop three potential anti-cancer candidates.

1 Introduction

As the primary source of a single carbon unit and the third-largest metabolite, serine is taken up by cancer cells to fuel their unregulated and rapid proliferation. The glycolytic intermediate 3-phosphoglycerate (3-PG), as the substrate for three consecutive enzymatic reactions, makes up an intracellular material for de novo serine production. Phosphoglycerate dehydrogenase (PHGDH) catalyzes the first stage of the reactions that convert 3-PG to 3-phosphohydroxypyruvate (3-PPyr) (Figure 1) (El-Hattab, 2016). PHGDH is abnormally expressed in various malignant tumor cells, including hepatocellular carcinoma, breast cancer, melanoma, lung cancer, glioma, colon carcinoma, Ewing sarcoma, pancreatic cancer, leukemia, thyroid carcinoma, multiple myeloma, lymphoma, and gastric and bladder cancer (Zhao et al., 2021). PHGDH dysregulation is a prevalent feature in a significant portion of malignancies with respect to their generation, proliferation, differentiation, and metastatic progress, which suggests that targeting PHGDH is a very promising direction in anti-cancer drug discovery (Rohde et al., 2018; Zhang et al., 2022; Gao et al., 2023).

Figure 1

Figure 1. The serine synthesis pathway.

Over the past decade, numerous potential inhibitors have been reported (Figure 2). Weinstabl et al. (2019) reported a highly potent and selective PHGDH inhibitor, BI-4924, which achieved ∼100 μM hits and was optimized to single-digit nanomolar potency. Prior to this, a series of pyrazole-5-carboxamides analogs representing compound 2 were identified (Zhao et al., 2021). CBR-5884, reported by Mullarky et al. (2016), contained functional groups that targeted the sulfhydryl moiety and could decrease de novo serine synthesis by reacting with the PHGDH cysteine residue. NCT-503 selectively affected the PHGDH oligomerization state in PHGDH-dependent cell lines and xenograft tumors (Rohde et al., 2018). Natural PHGDH inhibitors, such as Azacoccone E and Withangulatin A, were isolated from food and microorganisms (Guo et al., 2019; Zheng et al., 2019). Even though these emerging PHGDH inhibitors have not yet been employed in clinical trials, the rapid discovery of structurally novel potential PHGDH inhibitors with safety profiles is particularly urgent. Therefore, 3D-QSAR-pharmacophore-based virtual screening was performed in this study to identify novel PHGDH inhibitors.

Figure 2

Figure 2. Structures of representative potential PHGDH inhibitors that have been reported.

2 Materials and methods

2.1 3D-QSAR pharmacophore model generation

We performed a literature search and collected 31 compounds with corresponding IC₅₀ values ranging from 0.002 µM to 68 µM (Supplementary Figures S1, S2) (Spillier and Frederick, 2021; Zhao et al., 2021; Zhang et al., 2022; Gao et al., 2023). The Generate Training and Test Data of Discovery Studio (version 4.0, Biovea Inc., Omaha, NE, United States) protocol offered a random way to split a data set into 22 training set compounds and nine test set compounds by setting the training set percentage to 80. The training set was used to create a 3D-QSAR pharmacophore model, and the test set was reserved to assess the quality of the model. The Feature Mapping protocol computed all possible pharmacophore features, including hydrophobic (HYB), hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), negative ionizable (NI), and ring aromatic (RA). The uncertainty value of all compounds was fixed at 1.5, and IC₅₀ values of individual compounds were selected as an active property before hypothesis generation. Test set and Fisher validations with a confidence level of 95% were applied to validate the models. Other parameters were kept as default. The 3D-QSAR pharmacophore model was generated by running the 3D-QSAR Pharmacophore Generation protocol of Discovery Studio, which requires collecting chemically diverse training set compounds with the same bioassay value (Guner et al., 2004). The activities of the input training set (0.003–68 µM) and test set (0.002–39.6 µM) spanned five orders of magnitude. The best model among the 10 generated pharmacophores was selected by combining parameter correlation coefficient (R), root mean square deviation (RMSD), null cost, total cost, error, and fit values (Zhang et al., 2021).

2.2 Pharmacophore validation

The three techniques for validating the pharmacophore models were cost analysis, test set analysis, and Fischer randomization tests. Three types of costs are reported: total cost, null cost, and fixed cost in the HypoGen algorithm (Li et al., 2000). Generally, ∆Cost (Null cost − Total cost) is important in assessing the pharmacophore model. A ∆Cost of more than 60 bits suggests a significant correlation. A ∆Cost in the range of 40–60 bits means the model falls within the prediction range of 70%–90%. Predicting correlation likelihood will be difficult if the cost difference is less than 40 bits (John et al., 2011). The Fischer randomization method is essential to establish a structure–activity relationship between the structures and biological activity of the training set. The hypotheses were validated using the Fischer randomization approach, which rearranged the activity values of the training set molecules to yield 19 random spreadsheets with 95% confidence levels (Zhang et al., 2021). The pharmacophore model was also validated by inputting the test set consisting of nine compounds. All test and training set compounds were constructed and minimized using comparable procedures.

2.3 High-throughput virtual screening

The virtual screening compounds were downloaded from the Life Chemicals HTS (https://lifechemicals.com/), which consist of 3.04 million structurally novel molecules designed on the basis of promising drug-like scaffolds, carefully selected building blocks, and advanced cheminformatics approach, 460,160 compounds from the Enamine Hit Locator Library (https://enamine.net/compound-libraries/diversity-libraries), and 47,490 compounds from the ChemDiv 3D-pharmacophore database (https://www.chemdiv.com/catalog/screening-libraries/). The virtual screening based on 3D-QSAR pharmacophore was carried out in the Ligand Pharmacophore Mapping module, which compares a set of ligands to a selected pharmacophore. The relevant ligand-3D-QSAR pharmacophore mappings were exported and aligned to the pharmacophore.

2.4 ADMET property prediction

Absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties were predicted on AdmetSAR2 (http://lmmd.ecust.edu.cn/admetsar2/) (Yang et al., 2019) and Swiss-ADME online web (https://www.swissadme.ch) (Daina et al., 2017). The toxicities were predicted in ProTox-II (https://tox-new.charite.de/protox_II/), a virtual lab for predicting the toxicities of small compounds (Banerjee et al., 2018).

2.5 Molecular docking

Molecular docking was performed by using LibDock in Discovery Studio and AutoDock in AMDock software (Version 1.6.1) based on two different algorithms and scoring functions. The LibDock program developed by Diller and Merz (2001) is based on a binding site comprising lists of polar and apolar hot spots and presents docking results as a LibDock score (Raychaudhury et al., 2023). The protein was prepared by removing atomic clashes, water molecules, and unnecessary atoms, deleting alternate conformations, inserting the missing atoms in incomplete residues, and adding hydrogen. The site sphere was defined from PDB site records by Discovery Studio software. The Prepare Ligands module of Discovery Studio was utilized to generate the three-dimensional structures of anticipated bioactive small molecules. AMDock software (Version 1.6.1), which integrates with Autodock Vina, AutoDock4, and AutoDock4Zn, in which Autodock was based on an empirical free-energy force field and rapid Lamarckian genetic algorithm (Fuhrmann et al., 2010) search method, was used to docking study. The exported binding energy was used to assess the receptor–ligand affinity (Forli et al., 2016). Its graphical tool is simple to use and assists in molecular docking studies. The ligands and proteins were prepared by running the Prepare Input module. The co-crystal structure of PHGDH in complex with BI-4924 (PDB ID: 6RJ6) (Weinstabl et al., 2019) was downloaded from the PDB database (https://www.pdbus.org/). The search space was defined from the Center on Hetero: a box (Center: 18.66, −10.27, −1.74; Size: 24, 24, 24) was placed on the geometric center of an existing ligand. Subsequently, the docking simulations were run using Autodock (Valdes-Tresanco et al., 2020).

2.6 Molecular dynamics (MD) simulation

The binding stability and intermolecular interactions of the target macromolecule with promising hits after the docking process were assessed using Desmond2023-1 (Schrodinger, LLC, New York, NY, 2024) with OPLS_2005 force field tools (Kaminski et al., 2001) from a dynamic point of view. The orthorhombic box with the TIP3P water model (Mark and Nilsson, 2001) was used to predefine a 10 Å × 10 Å × 10 Å buffer region between the complex and box sides. The solvated system was neutralized with Na⁺ and Cl⁻ ions. The system builder panel allows minimization, which in turn allows the system to relax into a local energy minimum. Subsequently, elapsed 100-ns MD simulations were carried out at a periodic boundary condition, and the recording interval was set to 100 ps, while the ensemble class was set to NPT (T = 300 K, 1 atm pressure). Finally, approximately 1,000 frames were obtained, and the stability of ligand–receptor complexes was assessed by calculating several important parameters such as RMSD and root mean square fluctuations (RMSF). Their interactions also were analyzed when the simulation jobs were completed (Lanka et al., 2023).

2.7 Growth scaffold

Hit1, with the lowest estimated IC₅₀ value (0.0016 µM), was selected as the most promising lead compound. We performed reaction-based ligand enumeration within the PHGDH active pocket for lead optimization by using the Grow Scaffold protocol of Discovery Studio. Based on isosterism of the sulfonic acid moiety, the carboxyl was used to replace the sulfonic acid moiety in hit1 based on rational synthesizability. The intermediate could be formed by Claisen–Schmidt condensation using 2,5-dioxopyrrolidin-1-yl succinic acid as the starting material. Beginning with the intermediate position in the binding site of the receptor, the hydroxyl was selected to act as a reaction vector for the enumeration. Then, the amide synthesis and esterification reactions were selected. Other parameters were kept as default.

3 Results and discussion

3.1 Pharmacophore model generations

3D-QSAR pharmacophore models were generated from the features of known compounds that correspond to their activity, which could elucidate the spatial arrangement of chemical features of active compounds and facilitate the quick and effective discovery of promising hit compounds through ligand-based virtual screening (Yu et al., 2015; Desai et al., 2023).

The statistical parameters and chemical features of ten pharmacophore models are shown in Table 1. The total cost of the generated pharmacophore models ranged from 133.247 to 203.149, with a null cost of 564.644 and a fixed cost of 74.9106. Typically, the selected hypothesis has a significant total cost value, a low RMSD, a high correlation coefficient, and the highest cost difference. According to the results shown in Table 1, the △costs (Null cost − total cost) of the 10 pharmacophore models generated by the HypoGen algorithm were all larger than 40, indicating that these models are credible. Thus, Hypo_1 should have been selected as the optimal pharmacophore. However, it showed lower predictive ability in the subsequent validation results using the input test set (q² = 0.707). This is in contrast to Hypo_2, which had a more reliable pharmacophore model efficacy (q² = 0.918). Finally, Hypo_2 was selected as the best pharmacophore model, characterized by the lower total cost value (137.012), RMSD (2.37587), correlation coefficient (0.937476), HBA, HYB, HYB, NI, and RA (Figure 3) and the best predictability.

Table 1

Table 1. Statistical parameters of the ten pharmacophore hypotheses generated by the HypoGen algorithm.

Figure 3

Figure 3. The best pharmacophore model (Hypo_2); hydrogen bond acceptor (HBA, in green); hydrophobic feature (HYD, in blue); ring aromatic (RA, in brown); negative ionizable (NI, in dark blue).

Inhibitors with known activity were collected from published studies and classified into four groups: the most active (++++, IC₅₀ ≤ 0.01 µM); highly active (+++, IC₅₀ 0.01–0.1 µM); moderately active (++, IC₅₀ 0.1–1 µM); inactive (+, IC₅₀ > 1 µM). The experimental and estimated activities of the training set based on Hypo_2 are shown in Table 2. The most active compound in the training set, T1, was found to be well-mapped with the essential features of Hypo_2, whereas the inactive compound, T22, was not well linked with three chemical features (Figure 4).

Table 2

Table 2. Experimental and estimated activity of the training set based on Hypo_2.

Figure 4

Figure 4. Alignment of the most active and inactive ligands of the training set on the Hypo_2 pharmacophore model. (A) compound T1 and (B) compound T22.

3.2 Pharmacophore validation

3.2.1 Test set analysis

The ability of Hypo_2 to predict the biological activities of the test set, including nine chemically diverse compounds with IC₅₀ values ranging from 0.002 µM to 39.6 µM, was used to assess its dependability. As shown in Table 3, the most active compounds, highly active compounds, and inactive compounds were anticipated to be in the same range as their experimental range, except for two moderately active compounds. Simple regression has been used to analyze the correlation between experimental values and estimates of activity. The results illustrated in Figure 5 show the correlation coefficients between the actual and predicted PHGDH inhibitory activity values for both the training set (R = 0.937476) and the test set (R = 0.958384). With Hypo_2, it was apparent that the test set of nine compounds with known activity had been properly mapped.

Table 3

Table 3. Experimental and estimated activity of the test set based on Hypo_2.

Figure 5

Figure 5. The correlation coefficient (R) of the actual and predicted values of the test set and training set based on the Hypo_2 pharmacophore model.

3.2.2 Fischer’s randomization test

The pharmacophore hypothesis was assessed using the Fischer randomization test (Jacquez and Jacquez, 2002), which was designed with a 95% confidence level. Then, 19 spreadsheets were produced at random. Using this strategy, hypotheses were produced by randomly rearranging the bioactivity values of the compounds in the training set. Figure 6 depicts the differences in correlations (Figure 6A) and cost values (Figure 6B) between the HypoGen and Fischer randomizations. None of the randomly produced pharmacophores performed better statistically than HypoGen.

Figure 6

Figure 6. Correlation and total cost values of hypo-PHGDH and 19 random spreadsheets. (A) The correlation values; (B) The total cost value.

3.3 Virtual screening results and analysis

The 3.04 million compounds from the Lifechemicals database, 460,160 compounds from the Enamine database, and 47,490 compounds from ChemDiv were first filtered by Lipinski and Veber rules (Veber et al., 2002; Lipinski, 2004). The resulting compounds were then virtually screened using the constructed Hypo_2 model, and 1,006 compounds were mapped to Hypo_2. Ten top-ranked compounds with IC₅₀ 0.0016 µM–0.0325 µM are shown in Figure 7, and compounds 1–4 were selected as hit compounds with estimated IC₅₀ values of less than 0.01 µM. The estimated activity and predictive ADME properties of 10 top-ranked compounds are shown in Table 4. All parameters were within Lipinski’s rule of five (ROF) cut-off range for the test compounds resulting from the pre-filter. Toxicity restricts the development of specific compound families at any stage of drug development. Toxicity may take the form of Ames toxicity, carcinogenicity, hepatotoxicity, mutagenicity, or cytotoxicity and is the most noteworthy property of any potential drug candidates. Ideal drug-like compounds would be more highly potent yet not harmful (Ferreira and Andricopulo, 2019). In silico prediction of compound cytotoxicity and mutagenicity has drawn great attention from researchers, can assist the early identification of potentially harmful and mutagenic compounds, and can minimize the time and funding involved with hit-to-lead optimization (Kazius et al., 2005; Kramer et al., 2007; Chu et al., 2021). Therefore, we focused on investigating the predictive Ames toxicity, carcinogenicity, hepatotoxicity, cytotoxicity, mutagenicity, and LD₅₀ of the top 10 compounds. As shown in Table 5, compounds 1–4 and compounds 7–8 had good ADME properties without any toxicity.

Figure 7

Figure 7. The structures of the top 10 compounds from virtual screening based on Hypo_2.

Table 4

Table 4. The estimated activity and predictive ADME properties of top 10 screening compounds.

Table 5

Table 5. Predictive toxicities of the top 10 screening compounds.

Molecular docking provides a peripheral vision for investigating ligand and receptor interactions in active pockets. The most active ligand (T23) was chosen for active control with a LibDock score of 134.386 kcal/mol and an AutoDock score of −7.1 kcal/mol. Its interactions with PHGDH were first analyzed. As illustrated in Figure 8, the binding affinity of compound T23 was acquired from two strong hydrogen bonds interacting in the adenine pocket, which was accompanied by a bidentate polar anchor interaction of an -NH group and an -OH group with ASP174. In addition, T23 was stabilized by two other strong hydrogen bond interactions between residues ARG154 and ILE155. The aromatic ring exhibited Pi–Sigma interactions with PRO207. Using hit4 as an example for docking analysis, the PHGDH-hit4 complex was mainly stabilized by six strong hydrogen bond interactions between residues GLY151, GLY153, ARG154, ILE155, GLY156, and HIS205, as well as an attractive charge between HIS205. Moreover, the linker was bound within the hydrophobic cavity consisting of residues PRO207 and LEU192. It was interesting that both T23 and hit4 bound into the active pocket of PHGDH in a reversible “V” shape. The docking scores of hit1 to hit4 bound to PHGDH were lower than T23, especially for hit1, which suggested further structure optimization. Meanwhile, derivatives of hit1 were slightly higher than or comparable to T23 (Table 6). However, docked scores alone are insufficient to determine the potential of a compound.

Figure 8

Figure 8. The interactions of co-crystal structure (PHGDH-T23 complex); docking result diagrams of PHGDH with hit4 (below). (A) Solid ribbon model; (B) active sites of amino acid residues represented in the parent color line model; (C) 2D interaction diagram.

Table 6

Table 6. Docking scores of hit 1 to hit4 and derivatives from hit1.

3.4 MD simulation

MD simulations were employed to further identify the virtual leads. Hit1 to hit4 were analyzed and compared with the most potent inhibitor known, T23, in terms of the stability and steady nature of the empty protein and the respective contacts and mobility in the pocket of PHGDH. Following the alignment of each protein frame on the reference frame backbone, the atom selection was used to compute the RMSD, which could give insights into the structural conformation by monitoring RMSD in the simulation. When the simulation ends, its fluctuations are centered around a thermal average structure if it has reached equilibrium. Changes of the order of 1–3 Å are acceptable. On the other hand, larger fluctuations suggest that the protein is changing significantly during the simulation. The 100-ns empty protein simulation reached an equilibration phase with RMSD fluctuations in the range of 6 Å–8 Å (Figure 9A). In the Figure 9B plot, peaks indicate areas of the protein that fluctuate during the simulation. Typically, the tails (N-and C-terminal) fluctuate more than any other part of the protein.

Figure 9

Figure 9. The RMSD and RMSF plots of the empty protein backbone (PDB ID: 6RJ6). (A) The RMSD plot of the empty protein; (B) The RMSF plot of the empty protein (secondary structure elements: alpha-helical and beta-strand regions are highlighted in red and blue backgrounds, respectively).

In order to demonstrate that the screened novel hit compounds have the potential to be more effective inhibitors of PHGDH than the reported available inhibitors, a comparative MD simulation study was conducted using hit compounds and the most potent known PHGDH inhibitor (T23) as reported by Weinstabl et al. (2019). The PHGDH-T23 complex showed initial light fluctuation in RMSD-P within the first 40 ns, which stabilized in the range of 2.5 Å–4.5 Å (Figure 10A). Except for the N-and C-terminal, the P-RMSF fluctuated less than 3 Å. Comparing Figures 9B, 10B shows that protein chain binding to the respective ligand makes the loop region more stable. As shown in Figure 10C, the PHGDH-T23 interactions are categorized into three types: hydrogen interactions, hydrophobic interactions, and water bridges. During the MD simulation, the residues ASP174 (99%) and ARG154 (72%) of PHGDH exhibited long-time H-bond interaction with the inhibitor T23. Water bridging and hydrophobic interactions can also be observed (Figure 10D). The top panel of Figure 10E shows that an average of 10 residues formed contacts with T23 over the course of the trajectory. Some residues make more than one specific contact with the ligand, which is represented by a darker shade of orange. ARG154, ASP174, and PRO175 are shown with dark orange shades (Figure 10E).

Figure 10

Figure 10. The MD simulation results of the PHGDH-T23 complex (PDB ID: 6RJ6). (A) PHGDH-T23 complex RMSD plot; (B) PHGDH-T23 complex RMSF plot. (C) The main interactions of T23 with active residues during MD simulation. (D) A schematic of detailed T23 interactions with the protein residues. (E) A timeline representation of the PHGDH-T23 complex contacts.

The 3D docking poses and MD simulation of the PHGDH-hit1, PHGDH-hit2, and PHGDH-hit3 complexes are shown in Figure 11. The protein backbone RMSD and ligand fit protein initially revealed a large deviation that suggests irregular interactions during the simulations. Throughout simulations, the RMSD fluctuations of ligand fit proteins did not equilibrate or converge, which indicated their weak ligand–protein complex binding behavior and suggests the need for further structural optimization.

Figure 11

Figure 11. The 3D docking pose and MD simulation results of PHGDH–hit1, PHGDH–hit2, and PHGDH–hit3 complexes. (A–C) The 3D docking pose of the hit1, hit2 and hit3 with PHGDH, respectively; (D–F) The protein RMSD of PHGDH–hit1, PHGDH–hit2, and PHGDH–hit3 complexes, respectively.

A similar MD simulation was carried out for hit4. The RMSD of the PHGDH-hit4 complex was steadily maintained from 10 ns to 100 ns, which was approximately 2.4 Å to 4.8 Å (Figure 12A). The residues on the main protein chain showed low fluctuations with an RMSF value of less than 2.4 Å (Figure 12B). The bar histogram displays the four hydrogen bonds between hit4 and residues ARG154, ILE155, ASP174, and HIS205, as well as hydrophobic contacts with residues PRO175, ILE176, PRO207, and LEU209 (Figure 12C). The residues, including ARG154, HIS105, and TASP174, had a major contribution of hydrogen bond interaction to the hydroxyl and carbonyl groups of hit4 with the percentage of simulation time of 80%, 87%, and 86%, respectively (Figure 12D). Averages of eight residue interactions were maintained until the end of the simulation, as illustrated in the protein–ligand contact plot. Like the PHGDH-T23 complex, the residues ARG154, ASP174, and HIS205 showed more than one specific contact in 100-ns MD simulations (Figure 12E). Hit4 might, therefore, be identified as a candidate PHGDH inhibitor for combating various cancers.

Figure 12

Figure 12. The MD simulation results of the PHGDH–hit4 complex. (A) RMSD plot of the PHGDH–hit4 complex; (B) RMSF plot of the PHGDH–hit4 complex; (C) The main interactions of hit4 with active residues during MD simulation; (D) A schematic of detailed hit4 interactions with the protein residues; (E) A timeline representation of the PHGDH–hit4 complex contacts.

3.5 Newly designed compounds

Hit1, with the lowest estimated activity value, was selected as the lead compound for subsequent structure optimization. The newly designed derivatives were based on the inclusion of 3D-QSAR model insights and molecular docking. A total of 198 ligands were produced and were subsequently mapped to Hypo_02, and 133 compounds were mapped to Hypo_02. The criterion of estimating activity values less than 0.01 µM was used to further narrow down the scope of virtual hits. Two potential hits (Figure 13) were selected from output ligands. In the next MD simulation, two resulting compounds were evaluated individually.

Figure 13

Figure 13. The structures of eliciting hits from hit1 using Grow Scaffold.

The MD simulation results of the PHGDH-hit1-1 and PHGDH-hit1-2 complexes are shown in Figure 14. Through MD simulation, the RMSD of the protein backbone, Lig_fit_on_Prot, was stabilized in the range of 3.6 Å–5.2 Å, 1.8 Å–4.8 Å, 2.4–4.8 Å and 1.6 Å–4.0 Å, respectively, which indicated the existence of a stable binding pose between protein and ligands. The RMSF of the protein backbone showed high fluctuation at the tails (N- and C-terminals) and loop regions of the protein, especially hit1-1. In addition, the other RMSF values interacting with ligands were less than 2.0 Å. The newly designed hit1-1 and hit1-2 showed relatively stable binding poses with PHGDH and good results from the MD simulation studies compared to hit1. Therefore, these two newly designed hits could act as potential candidate inhibitors of PHGDH.

Figure 14

Figure 14. (A, B) The RMSD of the PHGDH–hit1-1 complex and RMSF of the protein plot during a 100-ns simulation; (C, D) The RMSD of the PHGDH–hit1-2 complex and the RMSF of the protein plot during a 100-ns simulation.

4 Conclusion

Since PHGDH became an attractive target for cancer research in the past decade, the development of PHGDH-targeted drug discovery has been limited due to the scarcity of reported inhibitors. In this study, a 3D-QSAR pharmacophore with credible predictive was established, which was validated by test set, cost analysis, and Fischer randomization tests. A 3D-QSAR-based high-throughput virtual screening strategy was comprehensively integrated to identify structurally novel PHGDH inhibitors from a commercial database, including approximately 3.54 million small molecules. The drug-like properties of the top 10 compounds were predicted through an in-silico study. Safety was considered the primary criterion for manual screening. Hit1, with the highest estimated active value, distinguished from the published compounds, was selected for lead optimization, even though unsatisfactory MD simulation results were found. Two hits with an estimated activity of less than 0.01 µM were obtained by lead optimization. Further MD simulation was employed to determine the dynamic binding behavior and binding stability of protein–ligand complexes. The related structural data of identified hit1 to hit4 and derivates from hit1 could be helpful in various cancer-related drug discovery projects.

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

YX: funding acquisition, project administration, software, writing–original draft. ZY: data curation, investigation, project administration, and writing–original draft. JY: investigation, methodology, validation, and writing–review and editing. CG: formal analysis, funding acquisition, project administration, and writing–review and editing. NQ: conceptualization, resources, software, supervision, and writing–review and editing. XW: conceptualization, methodology, software, supervision, validation, visualization, writing–original draft, and writing–review and editing.

Funding

The authors declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by grants from the Science and Technology Development Fund of Tianjin Education Commission for Higher Education (2020KJ204) and the Science and Technology Development Fund of Quzhou (2023K253).

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.

Supplementary material

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

References

Banerjee, P., Eckert, A. O., Schrey, A. K., and Preissner, R. (2018). ProTox-II: a webserver for the prediction of toxicity of chemicals. Nucleic Acids Res. 46, W257–W263. doi:10.1093/nar/gky318

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, C. S. M., Simpson, J. D., O’Neill, P. M., and Berry, N. G. (2021). Machine learning - predicting Ames mutagenicity of small molecules. J. Mol. Graph Model 109, 108011. doi:10.1016/j.jmgm.2021.108011

PubMed Abstract | CrossRef Full Text | Google Scholar

Daina, A., Michielin, O., and Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7, 42717. doi:10.1038/srep42717

PubMed Abstract | CrossRef Full Text | Google Scholar

Desai, S. P., Mohite, S. K., Alobid, S., Saralaya, M. G., Patil, A. S., Das, K., et al. (2023). 3D QSAR study on substituted 1, 2, 4 triazole derivatives as anticancer agents by kNN MFA approach. Saudi Pharm. J. 31, 101836. doi:10.1016/j.jsps.2023.101836

PubMed Abstract | CrossRef Full Text | Google Scholar

Diller, D. J., and Merz, K. M. (2001). High throughput docking for library design and library prioritization. Proteins 43, 113–124. doi:10.1002/1097-0134(20010501)43:2<113::aid-prot1023>3.0.co;2-t

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Hattab, A. W. (2016). Serine biosynthesis and transport defects. Mol. Genet. Metab. 118, 153–159. doi:10.1016/j.ymgme.2016.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, L. L. G., and Andricopulo, A. D. (2019). ADMET modeling approaches in drug discovery. Drug Discov. Today 24, 1157–1165. doi:10.1016/j.drudis.2019.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Forli, S., Huey, R., Pique, M. E., Sanner, M. F., Goodsell, D. S., and Olson, A. J. (2016). Computational protein-ligand docking and virtual drug screening with the AutoDock suite. Nat. Protoc. 11, 905–919. doi:10.1038/nprot.2016.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuhrmann, J., Rurainski, A., Lenhof, H. P., and Neumann, D. (2010). A new Lamarckian genetic algorithm for flexible ligand-receptor docking. J. Comput. Chem. 31 (9), 1911–1918. doi:10.1002/jcc.21478

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, D., Tang, S., Cen, Y., Yuan, L., Lan, X., Li, Q. H., et al. (2023). Discovery of novel drug-like PHGDH inhibitors to disrupt serine biosynthesis for cancer therapy. J. Med. Chem. 66, 285–305. doi:10.1021/acs.jmedchem.2c01202

PubMed Abstract | CrossRef Full Text | Google Scholar

Guner, O., Clement, O., and Kurogi, Y. (2004). Pharmacophore modeling and three dimensional database searching for drug design using catalyst: recent advances. Curr. Med. Chem. 11, 2991–3005. doi:10.2174/0929867043364036

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Gu, X., Zheng, M., Zhang, Y., Chen, L., and Li, H. (2019). Azacoccone E inhibits cancer cell growth by targeting 3-phosphoglycerate dehydrogenase. Bioorg Chem. 87, 16–22. doi:10.1016/j.bioorg.2019.02.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacquez, J. A., and Jacquez, G. M. (2002). Fisher’s randomization test and Darwin’s data -- a footnote to the history of statistics. Math. Biosci. 180, 23–28. doi:10.1016/s0025-5564(02)00123-2

PubMed Abstract | CrossRef Full Text | Google Scholar

John, S., Thangapandian, S., Arooj, M., Hong, J. C., Kim, K. D., and Lee, K. W. (2011). Development, evaluation and application of 3D QSAR Pharmacophore model in the discovery of potential human renin inhibitors. BMC Bioinforma. 12 (Suppl. 14), S4. doi:10.1186/1471-2105-12-S14-S4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, G. A., Friesner, R. A, Tirado-Rives, J, and Jorgensen, W. L. (2001). Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B 105, 6474–6487. doi:10.1021/jp003919d

CrossRef Full Text | Google Scholar

Kazius, J., McGuire, R., and Bursi, R. (2005). Derivation and validation of toxicophores for mutagenicity prediction. J. Med. Chem. 48, 312–320. doi:10.1021/jm040835a

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramer, J. A., Sagartz, J. E., and Morris, D. L. (2007). The application of discovery toxicology and pathology towards the design of safer pharmaceutical lead candidates. Nat. Rev. Drug Discov. 6, 636–649. doi:10.1038/nrd2378

PubMed Abstract | CrossRef Full Text | Google Scholar

Lanka, G., Begum, D., Banerjee, S., Adhikari, N., P, Y., and Ghosh, B. (2023). Pharmacophore-based virtual screening, 3D QSAR, Docking, ADMET, and MD simulation studies: an in silico perspective for the identification of new potential HDAC3 inhibitors. Comput. Biol. Med. 166, 107481. doi:10.1016/j.compbiomed.2023.107481

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Sutter, J., and Hoffmann, R. (2000). “HypoGen: an automated system for generating 3D predictive pharmacophore model,” in Pharmacophore perception, development, and use in drug design. Editor O. F. Güner (La Jolla, CA: International University Line), 171–189.

Google Scholar

Lipinski, C. A. (2004). Lead- and drug-like compounds: the rule-of-five revolution. Drug Discov. Today Technol. 1, 337–341. doi:10.1016/j.ddtec.2004.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Mark, P., and Nilsson, L. (2001). Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J. Phys. Chem. A 105 (43), 9954–9960. doi:10.1021/jp003020w

CrossRef Full Text | Google Scholar

Mullarky, E., Lairson, L. L., Cantley, L. C., and Lyssiotis, C. A. (2016). A novel small-molecule inhibitor of 3-phosphoglycerate dehydrogenase. Mol. Cell Oncol. 3, e1164280. doi:10.1080/23723556.2016.1164280

PubMed Abstract | CrossRef Full Text | Google Scholar

Raychaudhury, C., Srinivasan, S., and Pal, D. (2023). Identification of potential oral cancer drugs as Bcl-2 inhibitors from known anti-neoplastic agents through docking studies. J. Math. Chem. 62, 317–329. doi:10.1007/s10910-023-01537-w

CrossRef Full Text | Google Scholar

Rohde, J. M., Brimacombe, K. R., Liu, L., Pacold, M. E., Yasgar, A., Cheff, D. M., et al. (2018). Discovery and optimization of piperazine-1-thiourea-based human phosphoglycerate dehydrogenase inhibitors. Bioorg Med. Chem. 26, 1727–1739. doi:10.1016/j.bmc.2018.02.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Spillier, Q., and Frederick, R. (2021). Phosphoglycerate dehydrogenase (PHGDH) inhibitors: a comprehensive review 2015-2020. Expert Opin. Ther. Pat. 31, 597–608. doi:10.1080/13543776.2021.1890028

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdes-Tresanco, M. S., Valdes-Tresanco, M. E., Valiente, P. A., and Moreno, E. (2020). AMDock: a versatile graphical tool for assisting molecular docking with Autodock Vina and Autodock4. Biol. Direct 15, 12. doi:10.1186/s13062-020-00267-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Veber, D. F., Johnson, S. R., Cheng, H. Y., Smith, B. R., Ward, K. W., and Kopple, K. D. (2002). Molecular properties that influence the oral bioavailability of drug candidates. J. Med. Chem. 45, 2615–2623. doi:10.1021/jm020017n

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinstabl, H., Treu, M., Rinnenthal, J., Zahn, S. K., Ettmayer, P., Bader, G., et al. (2019). Intracellular trapping of the selective phosphoglycerate dehydrogenase (PHGDH) inhibitor BI-4924 disrupts serine biosynthesis. J. Med. Chem. 62, 7976–7997. doi:10.1021/acs.jmedchem.9b00718

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Lou, C., Sun, L., Li, J., Cai, Y., Wang, Z., et al. (2019). admetSAR 2.0: web-service for prediction and optimization of chemical ADMET properties. Bioinformatics 35, 1067–1069. doi:10.1093/bioinformatics/bty707

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Li, X., Ge, C., Si, H., Cui, L., Gao, H., et al. (2015). 3D-QSAR modeling and molecular docking study on Mer kinase inhibitors of pyridine-substituted pyrimidines. Mol. Divers 19, 135–147. doi:10.1007/s11030-014-9556-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Xiang, J., Xie, Q., Zhao, J., Zhang, H., Huang, E., et al. (2021). Identification of influenza PA(N) endonuclease inhibitors via 3D-QSAR modeling and docking-based virtual screening. Molecules 26, 7129. doi:10.3390/molecules26237129

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, F. M., Yuan, L., Shi, X. W., Feng, K. R., Lan, X., Huang, C., et al. (2022). Discovery of PHGDH inhibitors by virtual screening and preliminary structure-activity relationship study. Bioorg Chem. 121, 105705. doi:10.1016/j.bioorg.2022.105705

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J. Y., Feng, K. R., Wang, F., Zhang, J. W., Cheng, J. F., Lin, G. Q., et al. (2021). A retrospective overview of PHGDH and its inhibitors for regulating cancer metabolism. Eur. J. Med. Chem. 217, 113379. doi:10.1016/j.ejmech.2021.113379

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, M., Guo, J., Xu, J., Yang, K., Tang, R., Gu, X., et al. (2019). Ixocarpalactone A from dietary tomatillo inhibits pancreatic cancer growth by targeting PHGDH. Food Funct. 10, 3386–3395. doi:10.1039/c9fo00394k

PubMed Abstract | CrossRef Full Text | Google Scholar