Identification of PPA1 inhibitor candidates for potential repurposing in cancer medicine

Inorganic pyrophosphatase 1 (PPA1) is pivotal to cellular metabolism as it facilitates the hydrolysis of PPi—a by‐product of various metabolic processes that influence cell growth and differentiation. Overexpression of PPA1 enzyme has been linked to diminished patient survival and was shown to influence tumor cell dynamics, thereby positioning it as a potential therapy target for a variety of cancers including colorectal cancer, diffuse large B‐cell lymphoma, and lung adenocarcinoma. Despite this therapeutic promise, there are no known inhibitors of PPA1 as of today. In this study, we searched for potential PPA1 inhibitors using a molecular docking screen of 30 470 compounds with a history of clinical trials and/or US Food and Drug Administration approval. We specifically targeted the active pocket that coincides with the established catalytic domain. Our screen identified promising hits, which we further subjected to ADMET (absorption, distribution, metabolism, excretion, and toxicity) filtering. Subsequent molecular dynamics (MD) analyses were conducted on devazepide, quinotolast, and tarazepide—the three substances that successfully navigated all filters. MD analyses reinforced the stability of the protein‐ligand complexes and confirmed ligand binding, as substantiated by our root mean square deviation, radius of gyration and secondary structures of proteins analyses. Furthermore, Molecular Mechanics Poisson‐Boltzmann Surface Area calculations post‐MD identified devazepide and quinotolast as showing higher binding affinities; being supported by principal component analysis, free energy landscape, and dynamic cross‐correlation matrix results. Overall, our study reveals devazepide and quinotolast as potential candidates for PPA1 inhibition which could be considered for repurposing studies that need further experimental validation. These results not only reveal a potential for clinical repurposing for PPA1 inhibition but they also offer valuable insights into the development of future compounds for targeting the crucial PPA1 enzyme.

for PPA1 inhibition which could be considered for repurposing studies that need further experimental validation.These results not only reveal a potential for clinical repurposing for PPA1 inhibition but they also offer valuable insights into the development of future compounds for targeting the crucial PPA1 enzyme.

| INTRODUCTION
Phosphate is one of the most prevalent elements in the human body, accounting for approximately 1% of total body mass. 1 Phosphate plays multiple roles in the regulation of cellular processes, including membrane formation and energy conversion, signal transduction, as well as gene expression. 2,3The maintenance of phosphate is highly important for the body and the dysregulation of phosphate metabolism can result in the pathogenesis of various diseases 4 including cancer. 5In cells, pyrophosphatases play a significant role in converting inorganic PPi, which is generated as a by-product of many metabolic reactions, 6 to inorganic phosphate (Pi) through hydrolysis (Figure 1, lower middle panel). 7Two types of family I pyrophosphatases are found in humans: inorganic pyrophosphatase 1 (PPA1; 289 aa) and −2 (PPA2; 305 aa). 8PPA1 is present in all human tissues and plays a more prominent role than PPA2. 9PA1 exerts a wide range of functions including lipid metabolism, bone formation, collagen synthesis, 10 DNA synthesis, and neurite growth. 11It indirectly regulates pathways involving key cancer-related proteins such as p53, β-catenin, Bcl − 2, and caspase-3 through directly dephosphorylating c-Jun N-terminal kinase 1 (JNK), thereby controlling cell proliferation and apoptosis (Figure 1). 12PPA1 also indirectly activates the PI3K pathway 13 which is heavily emphasized in cancer.In line with these, PPA1 expression has been notably reported to be upregulated in various types of cancer, such as colorectal cancer, 12,14,15 diffuse large B-cell lymphoma, 16 lung adenocarcinoma, 17,18 prostate cancer, 19 hepatocellular carcinoma, 20 breast cancer, 21 gastric cancer, 22,23 intrahepatic cholangiocarcinoma, 24 and ovarian cancer. 18,25Remarkably, survival analyses revealed a negative impact of PPA1 levels on patient survival in nonsmall cell lung cancer, 26 colon adenocarcinoma, 12 and ovarian cancer. 27he intriguing relationship between PPA1, cancer and survival prompted scientists to investigate the tantalizing possibility of a potential therapeutic outcome of PPA1 inhibition.As no such inhibitors are known yet, scientists were able to demonstrate its potential use by employing RNA interference mediated knockdown experiments, which resulted in reduced cell proliferation and increased apoptosis in lung and breast cancer cells. 28imilar outcomes were also seen in diffuse large B-cell lymphoma, with an emphasis on the emerging relationship between PPA1 and p53. 16There were also promising attempts to inhibit PPA1-related pathways and hence several potentially therapeutic agents such as the JNK inhibitor JNK-IN-8, 29 PI3K inhibitor BKM120, 30 and AKT inhibitor capivasertib 31 were introduced.
Despite the significant biological and patho-physiological implications, there are still no known inhibitors for this pivotal enzyme as of today.In line with the emerging need of PPA1 inhibitors in experimental cancer medicine, the aim of this study was to identify potential competitive inhibitors of PPA1 from a pool of currently available small molecules.For potential repurposing, we focused on named molecules available for sale, which either are of biogenic origin or reached clinical trials and/or US Food and Drug Administration (FDA) approval.We followed a computational pipeline as given in Figure 2. Initially, we performed a drug pocket analysis with DogSiteScorer 32 and identified only one potential druggable pocket (not shown) that aligns well with the enzyme's catalytical site as reported by Niu et al. in 2021.Targeting this catalytical site, we employed PyRx molecular docking screen system 33 using 30 470 molecules, as obtained from ZINC15 compund library. 34Small molecules that displayed relatively higher affinity (binding free energy ≤ −9.5 kcal/mol) through this process were subjected to absorption, distribution, metabolism, excretion, and toxicity (ADMET) filters, to determine the inhibitor candidates that could reach clinics as fast as possible within a potential repurposing scheme.Successful candidates were then subjected to molecular dynamics (MD) simulations against the catalytical site of PPA1, followed by detailed calculations (principal component analysis [PCA], free energy landscape [FEL], Dynamic Cross-Correlation Matrix [DCCM], define secondary structures of proteins [DSSP] and MM-PBSA analyses).

| Collection of the molecular structures of small molecules
We downloaded 30 470 compounds from the ZINC15 database. 34To effectively manage the extensive libraries within this database, we implemented a focused approach by narrowing the investigation to "named" and "commercially available" compounds that are readily accessible in the market and have biogenic origins, currently in clinical trials, historically used in clinical trials and/or approved by FDA; with the purpose of finding hits that could proceed to clinical repurposing as soon as possible.The molecules that we selected were saved in the *.mol2 file format and imported into the workspace of PyRx docking screen tool. 33Following this, 10 000 steps of energy minimization with the GAFF and conjugate gradient method were used to optimize the atomic positions of the ligands, ensuring their precision for reliable compound screening.

| Preparation and validation of PPA1 structure
To create a viable PPA1 protein structure, the X-ray diffraction crystal structure of PPA1 (Protein Data Bank [PDB]: 6C45) was obtained from PDB 35 at a resolution of 2.39 Å. 36 To complete the missing residues, we performed homology modeling with the MODELLER software 37 making use of other human PPA1 structures found in the PDB database (7CMO, 7BTN).Created models were verified using the Discrete Optimised Protein Energy [DOPE] algorithm implemented in MODELLER. 38The best model determined by the DOPE score was selected for further assessment to comprehensively evaluate its structural integrity and accuracy with the help of advanced tools such as ProSA-web 39 and PROCHECK. 40

| Molecular docking screen and ADMET analysis
To determine targetable pockets in the protein structure, we employed DogSiteScorer, a powerful tool for identifying potential binding sites. 32This revealed only one suitable ligand pocket, which perfectly aligned with the catalytic domain of PPA1 as Niu et al. previously published (2021).The center of mass of the gridbox for this target site was calculated as 19.55, 60.46, and 37.41 on X, Y, and Z-axes, respectively.To ensure that the docking gridbox encompases all essential atoms of the catalytical pocket, box dimensions were set to 20 Å × 20 Å × 20 Å.Molecular docking was serially performed using PyRx, which utilizes AutoDock Vina. 33,41Following the docking process, binding energies were obtained for each candidate molecule.Molecules with a binding free energy ≤ −9.5 kcal/mol were projected into ADMET analysis.
ADMET analysis was conducted using two software tools, SwissADME 42 and AdmetSar. 43,44The SMILES format of each molecule was extracted from the ZINC15 database.SwissADME, a platform featuring rule-based filters from major pharmaceutical companies, was employed for Lipinski (Pfizer) 45 and Muegge (Bayer) 46 filters.Candidates with one or less filter violation and bioavailability scores ≥ 0.55 were included for further consideration as suggested by previous literature. 47,48AdmetSar was utilized for further analysis, human intestinal absorption was classified using the HIA data set and framework to establish its viability. 49Acute oral toxicity was determined based on a model considering toxicity measurements, with lower toxicities (<2.5) indicating promising candidates. 50

| Molecular dynamics simulations of PPA1-ligand complexes
Promising small molecule candidates underwent MD simulations to gain a better insight into their mode of binding in dynamic motion.Here, we employed GROMACS 2021.3 software 51 using the AMBER99SB-ILDN all-atom force field, 52 which is specifically compatible for our downstream analyses (including Molecular Mechanics Poisson-Boltzmann Surface Area [MMPBSA]) to calculate the outcome of our MD simulations.Our MD simulations were initiated with *.pdb files obtained as a result of molecular docking.To generate ligand topologies and force field parameters, Acpype (AnteChamber PYthon Parser interfacE) was employed. 53Protein-ligand complexes were solvated in a cubic box (20 × 20 × 20 Å) using TIP3P water molecules.To maintain system neutrality, solvent molecules were replaced with either Cl − or Na + ions.Long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method.Throughout the simulations, a constant temperature of 300°K was maintained, while the pressure was controlled isotopically at 1 bar using Parrinello-Rahman coupling. 54he LINCS algorithm 55 was applied to constrain bond lengths, employing a time step of 2 fs.Before the MD run, structures were relaxed through energy minimisation using the steepest descent algorithm for 50 000 steps, followed by 1 ns of equilibration at NVT (constant number of particles, system volume, and temperature) and NPT (constant number of particles, system pressure, and temperature).Subsequently, independent MD simulations were performed, with atomic coordinates saved every 10 ps.The resulting trajectories were collected for further analysis and characterization.

| Analysis of molecular dynamics trajectories
We performed post MD analyses using the useful tools in the GROMACS kit that includes root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bond count, FEL, and the radius of gyration (Rg).The estimation of secondary structure mobility was calculated using the dictionary of protein secondary structure method (do_dssp). 56,57Results of these analyses were visually represented with Xmgrace.To assess dynamic cross-correlations, the MD-TASK suite 58 was employed.Principal components were calculated using the MODE-TASK PCA script. 59Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach, which has been utilized as a useful tool in recent publications 60,61 was adopted to approximate the free energy of binding (ΔG) for each candidate molecule using the gmx_MMPBSA tool. 62,63Per-residue decomposition analyses during MM-PBSA calculations were performed for all residues located within the proximity of 20 Å from the ligands.

| Preparation and validation of PPA1 structure
PPA1 is predominantly found in a homodimeric ortetrameric form. 36We extracted its homodimeric form from 6C45 PDB structure and performed homology modeling with MODELLER 37 using other human PPA1 structures (PDB:7CMO, 7BTN) to fill the missing aminoacid residues.To generate the best model, we used the Z-dope method, which approximates a protein model's score depending on its composition and size. 38mong 25 structures produced by MODELLER, the one with lowest normalized Z-dope score (−1.36785) was selected for downstream processing and analysis.Once the complete structure of PPA1 dimer was obtained, it was superimposed with 6C45, which resulted in an RMSD of 0.295 Å (Figure 3A).To further verify the validity of our model, we employed ProSA 39 and PROCHECK 40 web-tools.ProSA evaluated the quality of our model by comparing the C α atom positions with similar structures found in the PDB database revealing a Z-score of −8.36 (Figure 3B).We compared our structure with other human PPA1 structures found in the PDB database, including 6C45, 7CMO, and 7BTN, which produced Z-scores of −7.82, −7.02, and −7.66, respectively.We also checked if the catalytical pocket could fit into an appropriate docking gridbox (Figure 3C).For further structural evaluation, we generated a Ramachandran plot with PROCHECK.The majority of residues (91.9%) were found in most favored regions, 7.3% in additional allowed regions, 0.8% in generously allowed regions, and none in disallowed regions (Figure 3D).Taken together, these results indicated a well-modeled structure of PPA1, which was ready for further processes and analyses in this study.

| In silico compound screen and ADMET analysis
Molecular docking can be iteratively used to screen a vast number of compounds, 64 thereby enabling the discovery of molecules with potentially high binding affinities against a target protein.In this study, a total of 30 470 compounds sourced from the ZINC15 database 34 were subjected to docking onto the active site of the Chain-A of PPA1. 36Chain-A was selected as the better-scored chain in the original crystallized structure.Docking was performed using a suitable gridbox configuration, as illustrated in Figure 3C.We employed PyRx virtual screening software, 33 which employs Autodock Vina. 41he interaction between the ligand and the protein was determined by evaluating the binding free energy (Supporting Information: Table 1).Our virtual screening process was concluded by establishing a threshold of ≤ −9.50 kcal/mol.This allowed reducing the number of molecules, selecting only the significant compounds as demonstrated in Table 1.
Following the identification of 12 promising small molecules through molecular docking, their absorption, distribution, metabolism, elimination, and toxicity (ADMET) properties were evaluated using the widely recognized in silico ADMET profiling tools, Swis-sADME 42 and AdmetSar, 43,44 as presented in Table 2.In SwissADME, the Lipinski, 45 and Muegge 46 filters were applied.These filters are well-established indicators that provide valuable insights into a molecule's features, including molecular weight, lipophilicity and rotatable bonds. 47Additionally, only the molecules that had good Abbott bioavailability scores (≥0.55) were selected as suggested previously. 47,48AdmetSar was utilized to assess crucial factors such as human oral bioavailability, intestinal absorption characteristics and acute oral toxicity scores; where a threshold of ≤2.5 was applied as described before. 50We focused on candidate hits that could reach clinics easily with fewer potential side effects and the aforementioned factors served as eliminative parameters to screen and prioritize candidate molecules. 65Consequently; three ligands-devazepide, quinotolast and tarazepide-successfully met all the applied filters (given in bold characters in Table 2).To better understand the relationship between the promising ligands and PPA1, we visualized their binding modes using LigPlot+, 66 which allows one to observe the hydrogen bonds and the interacting amino acids (Figure 4).Next, the dynamics of these potential interactions between the selected ligands and PPA1 were analyzed using MD simulations.

| Molecular dynamics simulations
Stability and dynamics of the most promising small molecule ligands; devazepide, quinotolast and tarazepide, and their interaction with PPA1 were analyzed with a MD simulation using GROMACS. 51Upon molecular docking, all three ligands were separately subjected to a 40 ns MD simulation with the PPA1 homodimer at 300°K.The stability of each complex was demonstrated using RMSD of backbone-to-backbone trajectories to see possible backbone shifts.The backbone RMSD values for devazepide, quinotolast and tarazepide complexes reached stability at 20 ns and remained stable for the rest of MD simulation.The RMSD values for Devazepide-PPA1 complex averaged at 0.250 nm, with a standard deviation of 0.031.Similarly, tarazepide-PPA1 complex exhibited an average RMSD of 0.222 nm with a standard deviation of 0.036, while quinotolast-PPA1 complex displayed an average RMSD of 0.306 nm with a standard deviation of 0.045.The results pointed out similar RMSD characteristics for all three compounds in a dynamic relationship with PPA1 (Figure 5A).The Rg was then evaluated as a measure of compactness and protein structure size. 67This value is calculated by the distribution of atoms around the protein axis, facilitating the estimation of pressure exerted on a specific location and showing the strength of bonds between two cross-sections.In all inhibitorprotein complexes, we observed similarly stable results (Figure 5B).The Rg for devazepide-PPA1 complex averaged 2.746 nm with 0.015 standard deviation, the Rg for tarazepide-PPA1 complex averaged 2.749 nm with 0.012 standard deviation and the Rg for quinotolast-PPA1 complex averaged 2.750 nm with 0.014.To determine whether PPA1 protein undergoes structural fluctuations, we calculated RMSF.Overall, all protein-ligand complexes exhibited similar behavior with this analysis.The Devazepide-PPA1 complex demonstrated an average RMSF of 0.137 nm, accompanied by a standard deviation of 0.066.In the case of the tarazepide-PPA1 complex, the RMSF averaged at 0.134 nm, with a standard deviation of 0.085.Similarly, the quinotolast-PPA1 complex exhibited an average RMSF of 0.148 nm, with a standard deviation of 0.066.Peak fluctuation was observed between residues 100 and 120 for all ligands, but tarazepide exhibited a notable increase in fluctuation in this interval (Figure 5C).In contrast to PPA1 protein per se, the ligands inside the pocket exhibited a more dynamic behavior.Devazepide and tarazepide underwent a significant change in binding mode at the beginning of the simulation, culminating in a stable state later on.Quinotolast, on the other hand, underwent a considerable change in its binding mode at the 8 ns mark, which was sustained until 28 ns before changing again and remaining at this final state until the end of simulation (Figure 5D).The devazepide complex exhibited an average RMSD of 0.782 nm, with a standard deviation of 0.082.For tarazepide, the average RMSD was 0.504 nm, accompanied by a standard deviation of 0.058.Additionally, the quinotolast complex displayed an average RMSD of 0.435 nm, with a standard deviation of 0.122.The behavior of all ligands were also reflected in hydrogen bond analysis (Figure 5E).Devazepide and tarazepide formed stable hydrogen bonds throughout the simulation, whereas quinotolast drifted in hydrogen bond formation at around 8 and 25 ns marks.Among the compounds examined, quinotolast demonstrated the highest count of hydrogen bonds averaging 1.160 throughout the simulation, while devazepide exhibited a relatively lower average of 0.334, and tarazepide achieved the lowest number of H-bonds averaging at 0.011.

| Principal motion and cross correlation analyses on PPA1 upon ligand binding
To further reveal the dynamic response of PPA1 structure upon ligand binding, we utilized PCA, a mathematical method that reduces multidimensional sets of variables to smaller dimensions based on covariance matrices.PCA determines the diffusive characteristics of proteins during various folding phases and calculates atomic displacements in each conformation comprising the trajectory. 68Here, devazepide was associated with a relatively stable PPA1 structure throughout the simulation, whereas quinotolast exerted some changes on PPA1 structure after the initial quarter of the simulation (Figure 6A).Tarazepide, on the other hand, causes an obvious change in the protein's conformation at later phases of the simulation.The thermodynamics and kinetics of molecular processes in solution could be inferred by FEL calculations. 69When FEL calculations were projected onto the PCA graph we obtained earlier, devazepide and quinotolast exhibited their lowest energy states towards the end of the simulation (Figure 6B).However, tarazepide began with the lowest energy state and transitioned to a higher energy state by the end.To further characterize the conformational dynamics of PPA1, we conducted DCCM analysis to see the correlated motions among all amino acids of PPA1 upon ligand binding.DCCM calculates cooperative atomic motions by analyzing the covariance and time correlation of positional fluctuations of atoms within a protein-ligand complex. 70Upon analyzing the dimer structure, it was observed that PPA1-devazepide complex exhibited the lowest level of correlation, while PPA1-quinotolast complex displayed a less degree of correlation.Conversely, PPA1tarazepide complex exhibited the highest level of correlation within the structure among all the ligand molecules studied here (Figure 6C).We also analyzed the mobility of secondary structures using the DSSP analysis and we curated structural snapshots from diverse time intervals throughout our MD simulations (0−20−40 ns) and superimposed them onto the original 6C45 structure (Supporting Information: Figure 1).These analyses aimed to understand structural dynamics linked to the interaction between the ligands and PPA1.The results indicated that the presence of any molecule did not result in significant conformational changes in PPA1's secondary structure.This could potentially be linked to the compounds' competitive binding behavior, which hinders substrate binding, rather than being solely attributed to inhibition through conformational changes (Supporting Information: Figure 2).

| Binding free energy analyses with MM-PBSA and per-residue decomposition
To determine the average binding affinities of stabilized PPA1 protein-ligand complexes, relatively stable final parts of all trajectories covering the last 10 ns of MD simulations were analyzed (Table 3).This analysis was conducted using the MMPBSA calculations. 62,63Devazepide exhibited the highest binding affinity among all tested compounds, indicating a strong interaction with the target protein.Quinotolast also displayed a similarly high binding affinity, implicating another strong interaction with PPA1.In contrast, tarazepide displayed a lower binding affinity, suggesting a relatively weaker interaction with PPA1 in dynamic motion.
We also conducted decomposition analyses to determine the contribution of individual amino acid residues to ligand binding during the MD simulation.Per-residue analysis revealed many interacting residues among which 8 residues exhibited a notable change in free energy (ΔG < −0.5 kcal/mol) for devazepide whereas 7 residues appeared for quinotolast and tarazepide (Figure 7).Although a large number of residues were affected, the decomposition results for significant residues (ΔG < −0.5 kcal/mol) for all candidates generally agreed well with the corresponding residues in the molecular docking poses presented in Figure 4 earlier.Among all the ligands, two residues, LYS153 and TRP187, were found to be shared.These residues have varying degrees of impact on different ligands, with some showing a more significant effect.Devazepide exhibited the highest number of residues that positively influenced the binding, indicating a strong and favorable interaction with PPA1.Compared to devazepide, there were fewer notable residues for quinotolast and tarazepide in this analysis; though with significantly higher binding affinities for certain amino acids such as ASP70 for quinotolast and LYS153 for tarazepide.Moreover, ligand binding to the catalytical domain of PPA1 seems to have minor effects on the stability of this protein when the trajectories are compared to that of ligand-free MD simulations (Supporting Information: Figure 3).Given the critical involvement of the PPA1 enzyme in cancer and the absence of its known inhibitors, our objective was to explore small molecule compounds with the potential to interfere with PPA1 activity.We particularly focused our screen on existing compounds that at least reached clinical trials, with a potential for repurposing strategy that could reach clinics more easily.Our in silico investigation revealed three potential repurposing ligands: devazepide, quinotolast, and tarazepide.Devazapide and tarazepide were successful in preclinical analyses and reached clinical stage, 71,72 while quinotolast progressed to regulatory approval for clinical  use in some countries. 73,74Our MD analyses revealed that the PPA1-ligand complexes exhibited stable behavior with these ligands, as indicated by consistent RMSD backbone and Rg values.The RMSF results exhibited similar fluctuations for all ligands, with the exception of tarazepide displaying greater fluctuation between residues 100−120.These results were in line with PCA and FEL evaluations as well as with MMPBSA calculations followed by decomposition (per-residue) analyses.In particular, devazepide and quinotolast demonstrated greater stability throughout MD and ultimately reached the minimum energy state at the end of MD simulations.Tarazepide, on the other hand, achieved a minimum energy state at the beginning of the simulation followed by a continuous decline as specifically noted in PCA and FEL graphs.DCCM analysis revealed slight differences in the motion of protein-ligand complexes with a notable difference for tarazepide in line with our other results.These observations suggest a compromised binding for tarazepide towards PPA1 in a dynamic simulation as opposed to the earlier docking result.This was also consequently observed in MMPBSA calculations.Devazepide and quinotolast consistently progressed towards the minimum energy complex, maintaining their stability in PPA1-bound form and thereby preserving their high binding affinities.Despite having similar docking poses, interacting amino acids, and molecular structures, devazepide and tarazepide exhibited a substantial difference in their binding to PPA1 during MD simulations.To gain deeper insight into the underlying causation, an intricate analysis of the pharmacophore attributes of these small molecules was conducted 75 as presented in Supporting Information: Figures 4−6.These attributes offer valuable information about the spatial arrangement of functional groups within the molecules, thereby influencing their biological activities. 76Furthermore, a comprehensive pharmacophore model was generated to illustrate the conserved features among these ligands, as depicted in Supporting Information: Figure 7.
The pharmacophore features of devazepide and tarazepide were generally similar, agreeing well with our molecular docking results.However, the difference between the binding stabilities of these ligands could be attributed to the presence of an extra H-bond donor in devazepide.This observation aligns well with our other findings, particularly in terms of hydrogen bond formation during MD simulations.Given that hydrogen bonds are regarded as the primary factors contributing to stronger protein-ligand interaction, 77 the presence of the extra H-bond donor in devazepide likely contributes to its favorable stability within PPA1's target binding pocket, thereby enabling it to have a relatively stable affinity.
Similar to devazepide, quinotolast also exhibited this Hbond donor around the same position and overall exhibiting better H-Bond formation capabilities compared to other ligands (Figure 5E).Intriguingly, quinotolast, which is the most active ligand within the binding pocket (Figure 5B), seems to contradict with the hydrogen bond formation.However, the introduction of hydrogen bond donors or acceptors to bolster proteinligand interactions may not invariably yield heightened binding affinity, 77 and might not necessarily augment the stability of the ligand.Thus, devazepide and quinotolast stand out as inhibitor candidates against PPA1.Devazepide is a cholecystokinin antagonist that targets the neuroactive ligand-receptor interaction, Insulin secretion and pancreatic secretion pathways. 78It is also worth mentioning that devazepide has primarily been tried as an anti-anxiety medication.Quinotolast is an orally administered mast cell stabilizer and exerts a cytoprotective influence on the gastric mucosa, serving as an agent with antiallergic and antiulcer properties. 79Considering all, it becomes crucial to seek further clarification in laboratory experiments to ascertain the safety of these drugs in the context of cancer treatment since computations inherently involve certain approximations and simplifications.Ensuring that these calculations correctly reflected the real life applications and, the anticipated and common side effects, including physiological depression or suppression in the case of quinotolast 80 and gallstone toxicity associated with devazepide, 81 should be effectively addressed before attempting to use them for repurposing.
In conclusion, our study aimed to identify potential competitive inhibitors for PPA1, a pivotal protein involved in phosphate metabolism and cancer cell proliferation.Our primary objective was to identify candidate repurposing molecules that have the potential to bypass the lengthy process of de novo drug development.However, the compounds identified here need further experimental validation.Following our initial analyses, three ligands; namely devazepide, quinotolast, and tarazepide, emerged as promising candidates.Each of these compounds did not cause a significant change in the structural conformation of PPA1.Upon conducting a more comprehensive evaluation, tarazepide did not a demonstrate satisfactory performance.Therefore, preclinical testing is especially warranted for devazepide and quinotolast.To overcome the limitations of our study, further preclinical toxicity/drug interaction analysis would also be invaluable before considering these candidates for PPA1 inhibition.Our MD simulations formed stable complexes within the 40 ns time-frame presented in this study.Still, there is no consensus in the literature on drug discovery in terms of MD simulation times [82][83][84][85] and simulations covering longer time periods would provide deeper insights into the structural dynamics of PPA1 bound by drug candidates.

F I G U R E 3
Evaluation of the quality of modeled PPA1 structure.(A) Structural alignment of the current model with the original PDB structure (6C45).(B) The Z-score plot of the model indicating the model's quality in relation to other structures within its range.(C) The placement of the grid box for molecular docking.(D) Ramachandran plot of the processed model.PPA1, inorganic pyrophosphatase 1.T A B L E 1 The outcome of the docking screen as identified molecules and their corresponding binding affinities.
Note: (Molecule structures were taken from the PubChem database).

F I G U R E 4
Candidate small molecule inhibitors against PPA1.Docking poses of (A) devazepide, (B) quinotolast, (C) tarazepide (left panel) and 2D illustration of the binding with LigPlot + (right panel).Ligand molecules are shown in contact with the amino acid residues through hydrogen bonds (yellow) and Van der Waals interactions (red).PPA1, inorganic pyrophosphatase 1.

F I G U R E 5
Molecular dynamics analyses to examine the behavior of potential inhibitors of PPA1.(A) Calculation of RMSD values for backbone C α atoms of PPA1 with different ligands over time.(B) Measurement of the radius of gyration (Rg) for all three ligands.(C) Calculation of RMSF values for all ligands along the aminoacid residues of PPA1.(D) Calculation of the RMSD values of ligands over time.(E) Intramolecular hydrogen bond formation between PPA1 and; devazepide (black), quinotolast (red), tarazepide (green).PPA1, inorganic pyrophosphatase 1; RMSD, root mean square deviation; RMSF, root mean square fluctuation.

F
I G U R E 6 Principal motion and cross-correlation analyses of MD trajectories.(A) Principal component analysis (PCA) where the first frame was highlighted in purple and the last frame was highlighted in yellow.(B) FEL (Free energy landscape) analysis, the function of PC1 and PC2 illustrated by FEL, with the lowest energy state denoted by the blue color and highest with red.(C) DCCM (Dynamic crosscorrelation matrix) plots showing dynamic cross-correlation of PPA1 aminoacids.Complete correlation is represented with red color, while complete anti-correlation is represented with blue.MD, molecular dynamics; PPA1, inorganic pyrophosphatase 1.

F I G U R E 7
The decomposition analysis of PPA1 residues.(Only the residues with ΔG < −0.5 kcal/mol were shown).PPA1, inorganic pyrophosphatase 1.
Evaluation of the ADMET properties of the molecules identified through molecular docking.
T A B L E 2Note: (Substances indicated in bold characters are those that passed all filters).
T A B L E 3 Binding free energies and important attributes of the ligands as calculated from MD trajecories by employing the MMPBSA method.: Avg, average; MD, molecular dynamics; MMPBSA, Molecular Mechanics Poisson-Boltzmann Surface Area; SEM, standard error of the mean. Abbreviations