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

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.

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.

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

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

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

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.

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.

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

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.

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.

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.

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

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.

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.

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.

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.