Repurposing FDA-Approved Agents to Develop a Prototype Helicobacter Pylori Shikimate Kinase (HPSK) inhibitor: A Computational Approach Using Virtual Screening, MM-GBSA Calculation, MD Simulations, and DFT Analysis

1. Introduction

Helicobacter pylori (H. pylori), a gram-negative spiral-shaped bacterium, colonizes the human stomach, leading to chronic infection [1,2]. It is estimated that 50% of the world's population is infected by H. pylori. Although 80-90% of infected individuals show no symptoms, the chronic infection can result in gastritis, peptic ulcers, or advance to cause gastric cancers [3]. Since 1994, the International Agency for Research on Cancer (IARC) has classified H. pylori as a Group 1 carcinogen, owing to its definitively established gastric carcinogenesis, as evidenced by epidemiologic studies [4]. Gastric cancer is a highly lethal malignant tumor that begins at the stomach mucosa level, disseminates to the other layers of the stomach wall, and evades adjacent and distant body parts in the most advanced stages [5]. In 2020, gastric cancer was the fifth most common cancer worldwide, with approximately 1.1 million new cases, and the fourth leading cause of cancer-related deaths, with around 800,000 deaths [6].

For more than two decades, H. pylori has been treated using triple therapy. In the first line, the therapy includes a proton pump inhibitor coupled with amoxicillin and clarithromycin. However, metronidazole-containing triple therapy or bismuth-containing quadruple therapy are used when clarithromycin resistance is present [7,8]. Currently, the major obstacle to eradicating H. pylori is the emergence of drug resistance, which deteriorates the efficacy of eradication regimens [9]. According to the World Health Organization (WHO), clarithromycin-resistant H. pylori is listed as a high-priority bacterium for the discovery, research, and development of new antibiotics [10]. Additionally, other factors such as treatment complexity, side effects, and costs necessitate the discovery of new antibiotics to combat H. pylori [11].

Advancing research and developing new antibiotics with unique mechanisms of action that provide high efficacy and low toxicity profiles are crucial for reducing the risk of drug resistance and expanding the antimicrobial arsenal [12,13]. In this context, the shikimate pathway (SP) represents a goldmine for finding new antimicrobials. The SP is a central metabolic pathway consisting of seven metabolic steps catalyzed by seven enzymes: 3-deoxy-D-arabino-heptulosonate 7-phosphate synthase (DAHPS), 3-dehydroquinate synthase (DHQS), DHQ dehydratase (DHQD), shikimate dehydrogenase (SDH), 5-enolpyruvylshikimate 3-phosphate synthase (EPSPS), and chorismate synthase (CS). Following the sequence of the aforementioned enzymes results in the transformation of phosphoenolpyruvate and erythrose 4-phosphate to chorismate – a precursor of essential metabolites such as aromatic amino acids, folic acid, and ubiquinone (Figure 1) [14,15]. The SP is present in microbes and plants but absent in animals, making it a promising target for the discovery of selective antimicrobials and herbicides [16]. Extensive efforts have been dedicated to uncovering compounds that interfere with the SP. Various natural, semi-synthetic, and synthetic inhibitors targeting distinct enzymes within the pathway have been identified. However, none of these inhibitors have successfully obtained approval as antimicrobials or herbicides, with the exception of glyphosate, a well-known herbicide that exerts its effect by inhibiting EPSPS [17].

With an intriguing drug target that meets the criteria for novelty and selectivity to address resistance, there remains a notable gap in the discovery of approved drugs that specifically target it. Intending to participate in filling this gap, in the current study, we targeted one of the key enzymes in the SP: shikimate kinase (SK), which catalyzes the phosphorylation of shikimate to form shikimate 3-phosphate in the fifth step of the metabolic sequence [18]. SK is classified as a P-loop kinase, sharing a common alpha-beta-alpha fold structure with conserved residues that create its active site, typically occupied by ATP and shikimate [19]. When the gene encoding SK (aroK) was deleted in Mycobacterium tuberculosis, it resulted in the disruption of cell viability, further establishing the attractiveness of targeting SP in general, and SK in particular [20].

Due to significant costs, slow progress, and high rates of failure in traditional drug discovery, repurposing existing drugs for the treatment of various diseases, both common and rare is gaining popularity. This approach utilizes already approved and de-risked compounds, potentially reducing overall expenses and shortening development duration [21]. In the present study, we utilized computer-aided drug design (CADD) approaches, namely molecular docking, Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) calculations, molecular dynamics (MD) simulations, and Density Functional Theory (DFT) calculations, aiming to investigate the repurposing possibility of US FDA-approved drugs to be used as H. pylori shikimate kinase (HpSK) inhibitors.

2. Results and Discussion

2.1. Molecular Docking and MM-GBSA Calculations

The crystal structure of H. Pylori Shikimate Kinase in its unbound conformation was used for the docking studies (PDB: 1ZUH). As the selected crystal structure of H. Pylori Shikimate kinase lacked a co-crystallized ligand, it was crucial to employ a reliable and effective approach to identify the enzyme’s active pocket where the inhibitor would bind. The active site was predicted with SiteMap and was validated by comparing it with a H. Pylori Shikimate Kinase crystal structure that is bound to a Shikimate ligand (PDB: 3MUF). The SiteMap predicted active site was found to match the shikimate bound active site, thereby, the active site for the studied protein is the competitive orthosteric site and on its wild unbound conformation. Molecular docking occupies an essential place in the process of computational drug discovery. This stems from the fact that it is used to examine the interactions between a target molecule and a ligand through both bonding and non-bonding interactions. Despite its limitation in neglecting the flexibility of biomolecules, it remains a reliable approach for generating preliminary insights on molecular interactions [39,40]. In the Glide module of docking, the HTVS method was applied to screen a total of 1,615 US-FDA-approved agents following their preparation. Subsequently, the top 111 candidates from the initial pool underwent extra-precision (XP) docking. The XP docking mode is known for its superior accuracy and precision compared to other modes, albeit at a higher computational cost. This mode was utilized for initial compound screening and further refinement based on XP docking scores. Glide’s scoring function integrates empirical and force-field-based parameters to calculate binding energy, allowing for the identification of the most favourable docking poses [3,4]. As a reference control, shikimate molecule was prepared as well and docked against the target protein, as the candidate drugs are aimed to be tested as competitive inhibitors of shikimate kinase. The docking score of shikimate was -5.867 kcal/mol and all the of other molecules were compared against it. Worthy of note is that there were 31 compounds that presented docking scores that were better than that of the reference ligand. To account for additional critical factors influencing molecular binding beyond molecular docking, MM-GBSA calculations were conducted on the chosen 31 compounds, along with shikimate.

The MM-GBSA analysis enhances the accuracy of binding energy calculations compared to molecular docking energies, as it follows a more theoretically robust approach rather than relying on empirical scoring functions commonly used in molecular docking [41]. The 31 compounds exhibited binding free energy (dG bind energy) values ranging from 3.61 to -55.92 kcal/mol, while shikimate showed a dG bind energy of -34.24 kcal/mol. Among the 31 compounds, 2 compounds had dG bind energies higher than -10 kcal/mol, 12 compounds fell between -10 and -30 kcal/mol, and 17 compounds had values more than -30 kcal/mol (Table S1). These 17 compounds were further refined based on their identity and indication of use, resulting in the identification of seven promising candidates detailed in Table 1 and their chemical structures in Figure 2.

In a subsequent filtration process, three drugs—Dolutegravir, Cangrelor, and Isavuconazonium—were selected for further analysis. These compounds were prioritized based on their drug nature, favourable clinical safety profiles, and superior MM-GBSA binding free energy values compared to the reference ligand. Dolutegravir, known for its use in antiretroviral (HIV-1) therapy, it has a well-documented clinical safety profile and excellent pharmacokinetics [42]. Meanwhile, Cangrelor is the sole intravenous P2Y12 receptor inhibitor currently available, known for its strong, consistent, and quickly reversible antiplatelet activity [43]. The third option, Isavuconazonium, is a well-established antifungal prodrug recognized for its demonstrated safety and effectiveness [44].

The ligand interaction diagram generated by Maestro revealed distinct interaction patterns for the control inhibitor, shikimate, and the docked compounds (Dlutegravir, Cangrelor, and Isavuconazonium) with the active binding site residues of the target protein. These differences in binding interactions provide insights into the structural features influencing ligand affinity and specificity.

Shikimate, serving as the control ligand, demonstrated robust hydrogen bonding with key amino acid residues at the active binding site. It formed hydrogen bonds with Gly11 (1.90 Å), Ser15n (1.73 Å), Asp31 (2.34 Å), and Asp33 (1.94 Å). These interactions suggest a strong stabilizing effect of shikimate within the binding pocket, likely contributing to its inhibitory efficacy.

As can be shown in Figure 3, Dolutegravir exhibited a distinct binding profile characterized by hydrogen bonding with Gly80 (1.87 Å) and salt bridges with Lys14 and Asp33. These interactions imply a combination of polar and electrostatic contributions to its binding affinity.

Cangrelor displayed an extensive interaction network within the active site. It formed hydrogen bonds with Gly11 (2.12 Å), Gly13 (1.88 Å), Ser15 (1.72 Å), Ser16 (2.52 Å), Asp31 (2.23 Å), Asp33 (2.55 Å), and Arg107 (2.17 Å). Additionally, two salt bridges with Arg107 were observed, emphasizing the strong electrostatic stabilization provided by this compound.

Isavuconazonium exhibited a unique interaction profile, including salt bridges with Met10, Asp31, and Asp33. It also formed π-π stacking interactions with Phe48 and a π-cation interaction with Arg132. These non-covalent interactions, particularly the π-π stacking and π-cation interactions, likely enhance the stabilization of Isavuconazonium within the binding site.

Comparatively, Shikimate's interaction profile was dominated by hydrogen bonds, contributing to its stable binding. Dlutegravir showed fewer hydrogen bonds but incorporated salt bridges, reflecting a different binding strategy. Cangrelor stood out with the highest number of hydrogen bonds and additional salt bridges, suggesting a potentially stronger and more specific binding affinity. In contrast, Isavuconazonium leveraged a combination of salt bridges, π-π stacking, and π-cation interactions, which may provide unique advantages in terms of binding stabilization and target selectivity. These findings highlight the distinct binding mechanisms of the tested compounds compared to the control inhibitor and underscore the potential of Dlutegravir, Cangrelor and Isavuconazonium as promising candidates for further development. This goal can be achieved via testing the stability of those compounds in molecular dynamics simulation studies.

2.2. Molecular Dynamics (MD) Simulation

Based on the initial docking results, Cangrelor and Isavuconazonium were chosen for further analysis by MD perturbation. Molecular dynamic simulation studies were essential to validate the stability and strength of the binding interaction over a period of 100 nanoseconds. Cangrelor and Isavuconazonium exhibited a robust and stable binding pattern within the binding pocket of the HpSK protein. Root mean square deviation (RMSD) plot analysis was used to measure the average displacement of atoms of a particular frame with respect to a reference frame. The analysis revealed a stable binding profile, as indicated by an average RMSD fluctuation range within 4 Å of the ligand with reference to the protein backbone throughout the simulation period. Figure 4 shows the RMSD plots for the top 2 selected compounds expressing the favorable binding profile with key residues in the binding pocket and most suited for drug repurposing: Cangrelor and Isavuconazonium. The interaction appears to stabilize within 20 nsec of the simulation period and the RMSD fluctuations remain within a narrow 2 Å range.

The hit compound Cangrelor exhibits one of the most favorable binding profiles, interacting with a number of residues of the HpSK binding pocket as shown in Figure 5. Notably HpSK binding residues Lys 14 and Ser 15 and form stable hydrogen bonds with the ligand for over 90% and 80% of the simulation period respectively. In addition to a stable hydrogen bond interaction formed with the HpSK residue Gly 13 which was stable for 96% of the simulation period. All significant interactions are represented in Figure 5A and B which highlight the residues with the strongest ligand interactions stable throughout the entire simulation period. Figure 5 C depicts all significant interactions displayed by the ligand and interacting residues occurring for over 30% of the simulation period. Cangrelor is an FDA approved anti-platelet medication prescribed for the prevention of myocardial infarction, it acts by inhibiting ADP mediated platelet aggregation through the selective and reversible inhibition of the P2Y12 receptor. The approval status and safety profile of this medication renders it particular amenable to subsequent repurposing studies.

Figure 6 shows the binding profile of hit compound Isavuconazonium with residues in the HpSK binding pocket. Isavuconazonium forms strong and moderate interactions with the key binding residues Met 10, Phe 56 and Phe 48 of the receptor binding pocket: hydrophobic and water bridge interactions with Met 10 and hydrophobic interactions with both Phe 56 and Phe The complex hydrophobic and water bridge interaction observed between Isavuconazonium and the HpSK residue Met 10 was observed for over 40 % of the simulation period. Isavuconazonium is a potent antifungal medication which has been approved by the FDA for the treatment of invasive aspergillosis. Like Cangrelor its approval status and safety profile provides an attractive foundation to launch the development of subsequent precursors that can be modified towards selective HpSK inhibition.

2.3. Density Functional Theory (DFT) Calculations

DFT serves as a powerful approach for assessing the energy levels of orbitals within compounds, enabling the anticipation of ligand stability and charge transfers at the interaction site. This method proves invaluable in computing the electronic structure of a material, providing beneficial insights into the behavior and reactivity of chemical compounds [45]. The electronic structure of the three selected compounds in addition to shikimate was calculated; High Occupied Molecular Orbitals (HOMO), Low Unoccupied Molecular Orbitals (LUMO), and the energy gap (HLG) between them were determined and presented in Figure 7.

For Dlutegravir, HOMO and LUMO orbitals located on the (3S,7R)-13-carbamoyl-12-hydroxy-7-methyl-9-oxo-4-oxa-5,8-diazatricyclo[8.4.0.0^{3,8}]tetradeca-1(14),10,12-trien-1-ylium-11-olate moiety with energies of -0.218 and -0.080 eV, respectively. The docking results revealed that this moiety interacted with the residues LYS14, ASP33, GLY80 of HpSK.

Regarding Cangrelor, the HOMO orbitals relies on the peripheral phosphonate group, with energy of -0.203 eV, while LUMO orbitals relies on purine-6-amine moiety with energy of -0.046 eV. According to the molecular docking, the phosphonate group interacted with ARG107, while the purine-6-amine moiety showed no interaction.

Lastly Isavuconazonium, its HOMO orbitals situated on the 4-(1,3-thiazol-4-yl)benzonitrile moiety with energy of -0.230 eV. The LUMO orbitals situated on 1,2,4-triazol moiety with energy of -0.068 eV. These two moieties interacted with ASP33 and ARG132 of the HpSK protein, respectively.

On the other hand, regarding the shikimate, HOMO and LUMO orbitals disseminated throughout the small molecule with energies of -0.263 and -0.049 eV, respectively. Shikimate exhibited interactions with GLY11, SER15, ASP31, and ASP33 of HpSK.

When the HOMO value is high, it signifies the compound's strong electron-donating capability, and conversely, a low HOMO value indicates the opposite [46]. For instance, Cangrelor exhibits superior electron donation potential with an energy level of -0.203 eV, whereas shikimate shows weaker electron donation ability at -0.263 eV. Furthermore, a narrow energy gap (HLG) between the HOMO and LUMO implies smoother electron transitions, translating to increased reactivity, improved photo-stability, and ultimately enhancing drug stability. The compounds can be ranked in the following order: Shikimate (0.213 eV) > Isavuconazonium (0.163 eV) > Cangrelor (0.157 eV) > Dlutegravir (0.138 eV). All the three drug candidates have HLG values less than the shikimate, which means they are more reactive and have more facilitating electronic transition than shikimate.

3. Materials and Methods

3.1. Software

In this study, all the in-silico analyses were conducted using Maestro version 2024-4 (procured via the RDIA Grant 12990-iau-2023-iau-R-3-1-HW: P.O. 6947). The Maestro platform is part of the Schrödinger software suite [22]. (License key: 03cb87b8-723c-4fec-9b8c-8a58137d7a76). The desktop workstation was equipped with Intel® Core™ i7-10700F Processor, Linux Ubuntu 22.10 operating system and a RTX 5000 graphics card.

3.2. Preparation of Ligands

A total number of 1615 US-FDA approved agents listed on the ZINC database were subjected to preparation. The LigPrep tool in the Maestro 12.8 software suite was utilized with default settings to process the collected ligands, aiming to minimize potential computational errors. This procedure yielded energy-minimized and structurally optimized 3D chemical assemblies. Structural inconsistencies within the ligand library were corrected, and plausible ionization states were generated under physiological pH conditions (7±2). Then, MacroModel module of Maestro was exploited for the purpose of preparing shikimate (the control ligand) for the consequent processes of molecular docking. In MacroModel module this goal was achieved by refining both the geometry and conformational flexibility of the ligands before molecular docking [23,24]. This refinement involved energy minimization, where bond lengths, angles, and torsions were optimized to reach stable, low-energy states, as well as conformational sampling to explore different ligand shapes and identify potential binding configurations. Utilizing the OPLS force field, MacroModel effectively captured the physical and chemical characteristics of shikimate, ensuring its stability and suitability for the docking process.

3.3. Protein Target Retrieval and Preparation

From the protein data bank that is available at https://www.rcsb.org_the crystal structure of Helicobacter pylori shikimate kinase (HpSK) enzyme was brought to the Maestro interface due to its proper quality (PDB ID: 1ZUH)[25]. The quality features found in this target protein were its high-resolution structure (1.80 Å), and its X-ray crystallographic source. These metrics provided a clear visualization of the overall protein structure, backbone, and side-chain conformations. In order for the protein to be ready for any handling, it should be subjected to a thorough process of preparation to correct the flaws in the protein data bank crystal structure. Usually, the imperfections in the structure occur as an improper bond order and missing side chains. Using the Schrödinger suite's Protein Preparation Wizard, multiple corrections to the protein under study were made [26]. These adjustments included: detecting disulfide bonds, adding hydrogen atoms, assigning correct bond ordering, and fixing mislabelled components. In the end, this preparation method uncovered and fixed structural problems to guarantee a solid basis for further docking. This was assured by making a protein reliability report [27].

3.4. Site Mapping and Generation of Receptor Gride

The SiteMap module in Maestro was used to predict the Crystal Shikimate Kinase active site [28]. In order to choose targetable sites, SiteMap provides an extensive set of tools for identifying, analyzing, and visualizing possible binding sites. This is achieved via evaluating factors such as shape, size, and hydrophobic or hydrophilic properties and ranking all the possible pockets by means of a site score [29,30]. After this, the coordinates of the pocket with the highest site score were used for receptor grid generation, ensuring precise determination of SK active site.

3.5. Molecular Docking

For the molecular docking analysis, Glide module was used. This module enabled precise and reliable interactions between the ligand and target protein. High Throughput Virtual Screening (HTVS) and Extra Precision (XP) are the three docking methods available in the Glide module [31]. To varied degrees, each approach strikes a balance between speed, accuracy, and scoring functions. The XP mode uses a stricter scoring process with improved sampling, although HTVS and SP use similar scoring mechanisms. In particular, chemicals that are less compatible with the protein's binding site are penalized by the XP approach [32,33,34,35]. For these variable ranges of accuracy and precision, the three modes of docking were performed to assess the binding affinities of the ligands to the Shikimate kinase binding site in a consequent manner. In exact, the approved drugs, processed during ligand preparation, were docked into the active shikimate binding site using the HTVS mode. Molecules that achieved high docking scores in HTVS were subsequently docked using the SP mode, followed by an additional docking in the XP mode for refined accuracy.

3.6. Free Binding Energy Calculation (MM/GBSA)

After completing the docking procedure, the ligands with the highest docking scores for both proteins were selected. The docking poses of the approved drugs in the Shikimate kinase binding pocket were then used to compute the free binding energy using the MM/GBSA method, implemented in the Prime module with its default settings [36]. The binding energy was calculated using the following formula:

E = E_complex – E_protein – E_ligand,

(1)

Whereby, E_complex represents the energy of the protein-ligand complex, E_protein denotes the energy of the protein, and E_ligand ​refers to the energy of the ligand.

3.7. Molecular Dynamic (MD) Simulations

The subsequent high-scoring complex was further evaluated with molecular dynamics (MD). Molecular dynamic (MD) simulation studies were conducted using the Desmond Module on Schrödinger’s Maestro platform as previously described [37]. In brief, the ligand-protein complex in its optimal docked configuration was minimized using the Protein Preparation Wizard and prepared for MD simulation using the system builder application of Desmond. The generated simulation environment was comprised of a water-based solvent system: TIP3P water model. An orthorhombic simulation box with a 10 Å buffer parameter from the protein surface was generated, the entire system was neutralized by calculating and adding the required number of counter ions and 0.15 M NaCl in order to attain isosmotic conditions. MD simulation was carried out at a temperature and atmospheric pressure of 300 K and 1.013 bar, respectively. The simulation was run for a total of 100 nanoseconds (1001 frames were saved in order to compile the trajectory). Analysis was run and results presented using the simulation interaction diagram tool of Desmond.

3.8. Density Functional Theory (DFT) Calculations

In order to determine the molecular reactivity of the successfully docked compounds, DFT calculations were performed. DFT is an effective method for evaluating the energy levels of orbitals in compounds, allowing for further prediction of ligands stability and charge transfer at the interaction site. Technically, DFT calculations were conducted using the Jaguar tool within the Maestro platform. The B3LYP method, paired with a 6–31 G** basis set, was chosen along with the Standard Poisson-Boltzmann (SPB) solvent model [38]. This analysis involved calculating the values of HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital).

4. Conclusions

This study successfully employed different computational experiments to investigate the repurposing potential of FDA-approved drugs as inhibitors of Helicobacter pylori shikimate kinase (HpSK). By screening 1,615 FDA-approved drugs through molecular docking and subsequent extra-precision docking, followed by MM-GBSA calculations, three compounds—Dolutegravir, Cangrelor, and Isavuconazonium—were identified as promising candidates based on their superior binding free energies, favourable drug properties, and safety profiles. Robust and stable binding profiles for both Isavuconazonium and Cangrelor were verified via molecular dynamics simulation. Additionally, Density Functional Theory (DFT) revealed that these compounds exhibited lower energy gaps (HLG) than the reference ligand, shikimate, indicating enhanced reactivity and more efficient electronic transitions. The results indicated that the identified drugs, while not optimal for repurposing, may serve as potential candidates for experimental evaluation against Helicobacter pylori shikimate kinase. The findings also indicated that diflunisal is a promising candidate for therapeutic repurposing and warrants further investigation.