記住我
ORIGINAL ARTICLE
Year : 2022 | Volume
: 17
| Issue : 6 | Page : 635-656
Identification of new small molecules as dual FoxM1 and Hsp70 inhibitors using computational methods
Zahra Alimardan1, Maryam Abbasi2, Ghadamali Khodarahmi1, Khosrow Kashfi3, Farshid Hasanzadeh4, Aghaei Mahmud5
1 Department of Medicinal Chemistry; Isfahan Pharmaceutical Sciences Research Center, School of Pharmacy and Pharmaceutical Sciences, Isfahan University of Medical Sciences, Isfahan, I.R. Iran., I.R. Iran
2 Department of Medicinal Chemistry, Faculty of Pharmacy, Hormozgan University of Medical Sciences, Bandar Abbas, I.R. Iran., I.R. Iran
3 Department of Molecular, Cellular, and Biomedical Sciences, Sophie Davis School of Biomedical Education, City University of New York School of Medicine; Graduate Program in Biology, City University of New York Graduate Center; Department of Chemistry and Physics, State University of New York at Old Westbury, New York, USA., USA
4 Department of Medicinal Chemistry, School of Pharmacy and Pharmaceutical Sciences, Isfahan University of Medical Sciences, Isfahan, I.R. Iran., I.R. Iran
5 Isfahan Pharmaceutical Sciences Research Center; Department of Biochemistry, School of Pharmacy and Pharmaceutical Sciences, Isfahan University of Medical Sciences, Isfahan, I.R. Iran., I.R. Iran
Correspondence Address:
Ghadamali Khodarahmi
Department of Medicinal Chemistry; Isfahan Pharmaceutical Sciences Research Center, School of Pharmacy and Pharmaceutical Sciences, Isfahan University of Medical Sciences, Isfahan, I.R. Iran.
I.R. Iran
Maryam Abbasi
Department of Medicinal Chemistry, Faculty of Pharmacy, Hormozgan University of Medical Sciences, Bandar Abbas, I.R. Iran.
I.R. Iran
Source of Support: None, Conflict of Interest: None
CheckDOI: 10.4103/1735-5362.359431
Background and purpose: FoxM1 and Hsp70 proteins are highly expressed in many cancers. Thus, their inhibition serves as Bonafede targets in cancer treatment.
Experimental approach: FDI-6, an inhibitor of FoxM1, was selected as a template, and based on its structure, a new library from the ZINC database was obtained. Virtual screening was then performed using the created pharmacophore model. The second virtual screening phase was conducted with molecular docking to get the best inhibitor for both FoxM1 and Hsp70 active sites. In silico, ADMET properties were also calculated. Finally, molecular dynamics simulation was performed on the best ligand, ZINC1152745, for both Hsp70 and FoxM1 proteins during 100 ns.
Findings / Results: The results of this study indicated that ZINC1152745 was stable in the active site of both proteins, Hsp70 and FoxM1. The final scaffold identified by the presented computational approach could offer a hit compound for designing promising anticancer agents targeting both FoxM1 and Hsp70.
Conclusion and implications: Molecular dynamics simulations were performed on ZINC1152745 targeting FoxM1 and Hsp70 active sites. The results of several hydrogen bonds, the radius of gyration, RMSF, RMSD, and free energy during the simulations showed good stability of ZINC1152745 with FoxM1 and Hsp70.
Keywords: Cancer; Dual inhibitor; FoxM1 and Hsp70 inhibitors; Molecular dynamics simulation; Pharmacophore modeling; Virtual screening.
The heat shock proteins (Hsp) 70 family, 68 to 74 kDa, are found in all main cell components, such as cytosol, mitochondria, and endoplasmic reticulum. Hsp70 plays a significant role in the protein activation process; folding, transferring from the membrane, degradation of misfolded proteins, suppression of the apoptosis pathway, and aging [1],[2],[3]. This protein is included in the two following major domains: the substrate-binding domain (SBD), which is responsible for binding and refolding of the substrates, and the nucleotide-binding domain (NBD), which is the ATP binding region [4].
Structural variation and open and closed modes of the Hsp70 protein are regulated by ATP and several co-chaperones in the presence of client proteins. Co-chaperones associated with the human Hsp70 are separated into three categories: (i) J-domain proteins stimulate ATPase activity and include Hsp40; (ii) nucleotide exchange factors are employed in promoting the ADP-ATP exchange that releases protein client from Hsp70, (iii) TPR domain cochaperone Hsp70-Hsp90 organizing protein (HOP), chromatin immunoprecipitation (CHIP) binds to the tetra peptide of human Hsp70, EEVD chaperone motif (residues 638-641) at the terminal carboxyl site [5].
In combination with CHIP and HOP proteins, Hsp70 activates the denaturation of the protein through the proteasome. It regulates the Hsp90 protein folding pathway by binding to the HOP and CHIP co-chaperones. Another Hsp70 co-chaperone is the Bcl-2 associated athano gene 3 (BAG3) protein which regulates the FoxM1 transcription factor signaling pathway inside the cell through binding to the ATPase domain of the Hsp70 protein [6],[7].
Hsp70 plays an important role in cancer development. Instead, Hsp70 had multiple effects on signaling pathway components related to tumor initiation, growth, and metastasis, such as FoxM1, hypoxia-inducible factor 1, and nuclear factor kappa B (NF-κB). Hsp70 regulated signaling networks via association with the co-chaperone BAG3, a scaffold protein with the capacity to interact with multiple key regulators of cell signaling. BAG3 bounds to Hsp70 and this complex controlled signaling network that contains FoxM1 [8][Figure 1]. Hsp70 inhibitors prevent the connection of BAG3 to Hsp70, thus disrupting the signaling pathway transcription of FoxM1 [9].
Figure 1: BAG3 protein binds to Hsp70, and this complex controls multiple pathways in cancer cells by regulating FoxM1. BAG3, Bcl-2 associated athano gene 3; Hsp, heat shock protein.FoxM1 is the Achilles heel of cancer and is more highly expressed in cancer cells. It plays a significant role in the induction of cell cycle progression, initiation and migration, differentiation, proliferation, DNA repair, angiogenesis, and suppression of apoptosis [10].
FoxM1 is a member of the forkhead box (FOX) transcription factor family with 100 amino acids and includes two wings directly contacting DNA. This protein consists of three domains: a DNA binding domain (DBD), a negative-regulatory domain (NRD), and a transactivation domain (TDA). The contact point of FoxM1 to DNA is the DBD domain [11],[12].
The classification of Hsp70 inhibitors has been performed based on the SBD and NBD inhibitors and the Hsp70 co-chaperone inhibitors. 2-Phenylethynesulfonamide (PES) is one of the inhibitors that bind to the C-terminal of Hsp70 and leads to the accumulation of misfolded proteins, instability of the lysosome membrane, and induction of autophagy in cancer cells [13]. In addition, the 2-phenylethynesulfonamide derivatives, 2-(3-chlorophenyl) ethyne sulfonamide and triphenyl(phenyl ethynyl)phosphonium bromide (PET-16), have also been found to exhibit more potent anticancer activities, including increased cytotoxicity of cancer cells by binding to the Hsp70 substrate-binding domain [14].
The NBD inhibitors are separated into ATP competitive inhibitors, which bind to the ATP binding pocket, and allosteric inhibitors [15]. VER-155008 (ATP analog) and apoptozole bind to the Hsp70 NBD domain and compete with ATP to stay in the ATP binding pocket [16]. MKT-077, JG-13, JG-98, YM-01, and YM-08 compounds are allosteric inhibitors of the Hsp70 NBD domain [17]. In addition, the MAL2-11B compound could inhibit co-chaperone Hsp40 [Figure 2] [18].
Figure 2: Structure of some Hsp70 inhibitors. (i) substrate-binding domain inhibitors; (ii) NBD, ATP competitive inhibitors; (iii) NBD, allosteric inhibitors; (iv) cochaperone Hsp40 inhibitor. Hsp, Heat shock protein; NBD, nucleotide-binding domain; PES, 2-phenylethynesulfonamide; PET, triphenyl(phenylethynyl)phosphonium bromide.Inhibition of FoxM1 is one of the most important pathways to treating cancers. FoxM1 inhibitors are separated into direct and indirect inhibitors. Thiostrepton and honokiol are natural compounds that inhibit the FoxM1 directly by binding to the FoxM1 DBD domain [19],[20]. Synthetic compounds such as an inhibitor of the forkhead box protein M1 transcription factor (FDI-6), RCM-1, and TFI-10 directly inhibit FoxM1 activity [Figure 3] [21],[22],[23].
Indirect inhibition of FoxM1 is either through inhibition or activation of its upstream proteins. For example, inhibiting proteins such as Hsp70, Hsp90, proteasome, human epidermal growth factor receptor 2, protein kinase B, c-myelocytomatosis oncogene, estrogen receptor alpha, CDK4/6, and NF-κB decrease the expression and activity of FoxM1. On the other hand, FoxM1 could also be inhibited by increasing the activity of p53, retinoblastoma protein, and FoxO3 [10],[24].
As mentioned above, simultaneous inhibition of Hsp70 and FoxM1 proteins offers a critical approach for curtailing cancer cell growth, which is why in recent years these two proteins have been targeted in drug design and development.
Bioinformatic studies are essential in finding cost-effective compounds with enhanced activities in less time [25],[26].
There are several methods to hit prediction and lead optimization in drug design discovery, such as virtual screening, quantitative structure-activity relationship (QSAR), pharmacophore modeling, and molecular dynamics (MD) simulations [27],[28],[29]. Molecular docking virtual screening and pharmacophore virtual screening are the best methods to obtain hit compounds from an extensive database [30],[31].
In this study, we used molecular docking virtual screening, pharmacophore model virtual screening, and MD simulations to predict new compounds as dual inhibitors of Hsp70 and FoxM1. FDI-6 and VER-155008 were selected as lead inhibitors for FoxM1 and Hsp70, respectively. First, a similarity search was performed using the FDI-6 structure on the ZINC database; this gave us 850 compounds with 80% similarity. Next, two pharmacophore models were created based on the binding of FDI-6 with the FoxM1 and the interaction of VER-155008 with the Hsp70 active sites. Then, virtual pharmacophore screening was performed on the completed library. Also, virtual screening was accomplished on the obtained compounds from virtual pharmacophore screening for both FoxM1 and Hsp70 proteins. These compounds were subjected to filters such as Lipinski's Rule of Five to determine drug-likeness properties. Finally, MD simulation studies were performed on the best-predicted compound to investigate further interactions between the ligands and proteins.
Materials and methodsChemical library preparation
The ZINC database contains more than 230 million chemical compounds and is one of the central databases used for virtual screening studies. In the first step, the FDI-6 compound as a FoxM1 protein inhibitor was chosen as a template for creating a data library from ZINC online database (https://zinc.docking.org). Then similarity search was performed using filtering based on 80% similarity with the FDI-6 compound and Lipinski's Rule of Five.
Pharmacophore model creating and pharmacophore virtual screening
The pharmacophore model describes the arrangement of essential chemical features of a ligand that must effectively bind to a receptor site. Structure-based pharmacophore modeling frequently uses the 3D structure of docked ligands. In this study, using the Pharmit online server (http://pharmit.csb.pitt.edu), two pharmacophore models were created based on main interactions between the known ligands, FDI-6 and VER-155008, with FoxM1 and Hsp70 proteins, respectively [32]. Then, based on these models, which included features of hydrogen bond donors and acceptors, and hydrophobic and aromatic interactions, virtual pharmacophore screening was performed on the library obtained from the ZINC database. Compounds were further refined by root mean square deviation (RMSD) and energy minimization.
Docking validation of FoxM1 and Hsp70 proteins
The crystallography structures of FoxM1 and Hsp70 (FoxM1: PDB ID = 3G73; Hsp70 PDB ID = 4IO8) were obtained from the RCSB protein data bank (PDB; https://www.rcsb.org/) [33]. Then, additional molecules such as water, DNA, and ligands were deleted from the protein PDB files using the Accelrys Discovery Studio Visualizer 4.0.2 [34]. The 3D structure of FDI-6 and VER-155008 were drawn by Marvin Sketch v5.7, ChemAxon software [35]. Other steps of the protein preparation were performed in the Auto Dock Tools package [36]. First, all the hydrogens were added, then Kollman and Gasteiger charges were calculated for proteins and ligands, respectively. Polar hydrogens were merged, and the edited files were saved in the pdbqt format. The center of the grid box was set to the coordinates of His287 residue for FoxM1 and Arg272 residue for Hsp70 proteins. A grid box with 70 × 70 × 70 points in x, y, and z directions and the grid point spacing of 0.345 Å was generated. For all docking procedures, 100 runs of AutoDock search were performed. The results of the dockings were ranked based on the binding free energies. The visualization of the ligand-protein interactions based on the resulting dlg files was carried out by Autodock Tools [11],[37].
Structure-based virtual screening
Structure-based virtual screening is an efficient way to determine the best compounds from a chemical library. The obtained compounds from the previous step, 520 for FoxM1 and 120 for Hsp70, were docked into the active sites of FoxM1 and Hsp70 using AutoDock 4.2. The docking results were then ranked based on the lowest binding free energies and the highest interaction between FoxM1 or Hsp70 proteins and their ligands.
ADMET properties prediction
The five important pharmacokinetic parameters, absorption, distribution, metabolism, excretion, and toxicity (ADMET), are essential for a compound designed as a suitable drug candidate [38],[39]. There is a relationship between chemical structures and physicochemical properties, so chemical descriptors can be used to calculate pharmacokinetic properties. For example, the polar surface area has been used to predict the absorption of drugs. Polar surface area is a part of the surface area of oxygen and nitrogen atoms bound to hydrogen atoms [40],[41]. log BB for brain/blood values (() indicate the rate of molecular crossing to the blood-brain barrier. If a compound can cross this barrier, it will be able to act in the central nervous system [42].
Additionally, when a drug crosses the dermal or gastrointestinal barrier, it is transferred by the blood to the target tissue and distributed throughout the body; it can make bonds with the plasma proteins, primarily albumin. Therefore, the LogKhsa parameter indicates a compound's ability to bind to the plasma proteins and their distribution throughout the body. Also, the percentage of oral absorption of compounds in the gastrointestinal tract, Madin-Darby canine kidney (MDCK) permeability, or LogKp 2-Caco dermal permeability coefficient was examined [43]. The octanol-water partition coefficient (logP) is a physicochemical parameter that indicates the hydrophobicity of compounds. Compounds with higher lipophilic properties have higher absorption rates and metabolism in the liver.
ADME descriptors were calculated using the QikProp v.3.0 program. Also, all compounds were subjected to the toxicity risk assessments such as mutagenicity, tumorigenicity, irritant, and reproductive effects using the Data warrior program (http://www.openmolecules.org/datawarrior/). Finally, drug-relevant values, such as drug-likeness values, were determined using the Osiris program (http://www.organic-chemistry.org/prog/peo/).
Molecular dynamics simulation
Molecular dynamics simulation is a helpful in-silico method that simulates the interactions between ligand and protein in the body's physiological environment. Since in-vitro studies of macromolecule structures are challenging and time-consuming, MD simulation studies could be a beneficial alternative. In this study, MD simulation was performed using the GROMACS-2019.3 software package [44]. AnteChamber Python Parcer InterfacE (ACPYPE) was used to generate the necessary topology files and force field parameters for ligands [45]. The PROPKA 3.1 webserver (http://propka.org) determined the amino acids pKa for each protein [46]. Amber force field and TIP3P water model were employed in the next step to assess protein topology [47],[48]. The ligand and protein complexes were dunked in water in a dodecahedral box with a minimum spacing of 1 nm between the protein's surface and the box. The total protein charge was calculated, and for its neutralization and balancing, the number of sodium or chloride ions was replaced with water molecules. At first, the energy minimization of the system was carried out using the steepest descent method. After energy minimization, the primary simulation was carried out in two stages to equilibrate the system at a constant temperature of 300 K during 100 ps; NVT (constant number of particles, volume, and temperature) and NPT (constant number of particles, pressure of 1.0 bar, and temperature). The linear constraint (LINCS) algorithm was applied for covalent bond constraints. Long-range electrostatic interactions were calculated with the Particle Mesh Ewald method [49].
MM-PBSA free binding energy calculation
Molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) is one of the best simulation methods for calculating binding energy. In this method, the binding energy contributions of the protein-ligand complex were calculated such as van der Waals contribution (ΔEvdW), electrostatic contribution (ΔEele), a polar desolvation term (ΔGGB), and a nonpolar termination desolvation term (ΔGSA). These were summarized as follows [50],[51]:
Library preparation
FDI-6 was selected as the template for similarity search in this virtual screening experiment. The online ZINC database was searched for compounds with 80% similarity to FDI-6, resulting in 850 similar compounds, stored in SDF format. These compounds were subjected to the following classified filtering steps.
Pharmacophore model virtual screening
Pharmacophore is capable of including more detailed information about regions available to the ligand for efficiently binding to its receptor. For example, the FoxM1 pharmacophore model possesses six features: one hydrogen bond donor, one hydrogen bond acceptor, and four hydrophobic and aromatic features, as shown in [Figure 4]A. Also, [Figure 4]B displays that the pharmacophore model of Hsp70 that consists of seven components, three hydrogen bond donors, one hydrogen bond acceptor, and two hydrophobic and aromatic features. The virtual pharmacophore screening was performed using two created models on the obtained ZINC database library. The results were refined with RMSD less than 2 Å and energy minimization better than -7 kcal/mol for FoxM1 and -6 kcal/mol for Hsp70, respectively. Finally, 520 molecules from FoxM1 and 120 molecules from Hsp70 were saved.
Figure 4: The structure-based pharmacophore model of (A) FDI-6 inside the active site of FoxM1 and (B) VER-155008 inside the active site of heat shock protein 70, generated with the Pharmit server. These chemical features are color-coded with green as hydrophobic, purple as aromatic, white as hydrogen bond donor, and orange as hydrogen bond acceptor.Docking validation of FoxM1 and Hsp70
Essential interactions of FDI-6 with FoxM1 protein are shown in [Figure 5]A. The binding energy of FDI-6 was obtained -2.7 kcal/mol. Pi-sulfur electrostatic and Pi-Pi interactions were observed between ring A and His287. The ring D formed Pi-Pi interaction with Trp308 and Pi-Alkyl interaction with Leu259 and Arg286. Finally, a Pi-sulfur interaction between ring A and His287 is a significant interaction between FDI-6 and FoxM1’s active site [52].
Figure 5: (A) Main interactions between FDI-6 and FoxM1; (B) main interactions between VER-155008 with heat shock protein 70. Water hydrogen bond (blue), conventional hydrogen bond (green), carbon-hydrogen bond (gray), Pi-alkyl (pink), and Pi-cation (yellow).Figure 5B depicted the main interactions of VER-155008 with the Hsp70 protein. The binding energy of VER-155008 was calculated as -1.77 kcal/mol and RMSD less than 2 Å. A Pi-Alkyl interaction was observed between ring A and Arg272. Ring B created a Pi-alkyl interaction with Arg272, Pi-cation, Pi-Alkyl interactions with Arg342, and also a Pi-stacked interaction with Gly339. A Pi-cation electrostatic interaction was seen between ring C and Arg272 [16].
Virtual screening analysis of FoxM1 and Hsp70 protein
The docking simulation study investigated the binding modes of all selected compounds in FoxM1 and Hsp70 active sites. The docking scores and primary interaction into the active site of proteins were analyzed to find the best compounds. Five compounds were extracted with the lowest energy and the highest interaction with the Hsp70 and FoxM1 active sites. These were ZINC8448537, ZINC1152745, ZINC199847048, ZINC199676474 and ZINC6110690 [Figure 6].
Figure 6: Chemical structure of ZINC8448537, ZINC1152745, ZINC199847048, ZINC199676474 and ZINC6110690.Binding mode of the potential dual inhibitors in FoxM1 active site
The main interactions of ZINC8448537 with the FoxM1 active site are shown in [Figure 7]A. The binding energy of this compound was -4.07 binding energy of this compound was -4.07 kcal/mol. Ring A created a Pi-cation electrostatic bond with His287. Ring D formed a Pi-alkyl interaction with Arg286 and Leu259 and a Pi-Pi interaction with Trp308. The N6, ring B, F23, (C=O)24, and (NH2)29 made five hydrogen bonds of the conventional type and Pi-donor with His287 (2.180 Å), Asn283 (3.499 Å), Trp308 (2.692 Å), Ser290 (1.961 Å), and Arg286 (2.695 Å), respectively.
Figure 7: Main interactions between compounds and FoxM1. (A) ZINC8448537, (B) ZINC1152745, (C) ZINC199847048, (D) ZINC199676474, and (E) ZINC6110690. Conventional hydrogen bond (green), carbon-hydrogen bond and Pi-donor hydrogen bond (gray), Pi-sulfur (yellow), Pi-Pi (purple), and Pi-alkyl (pink).Figure 7B shows the main interactions of ZINC1152745 in the FoxM1 active site. The binding energy of this compound was -4.27 kcal/mol. Pi-cation, Pi-sulfur electrostatic, and Pi-Pi interactions were formed between ring A and His287. In addition, ring C showed a Pi-sulfur electrostatic bond with His287 and a Pi-Alkyl interaction with Arg286. Also, Arg286 made a Pi-alkyl interaction with ring D.
The binding energy of ZINC199847048 in the active site of FoxM1 was -3.83 kcal/mol. As shown in [Figure 7]C, seven hydrogen bonds, Pi-sulfur, Pi-cation, Pi-alkyl, and Pi-Pi interactions were formed between compound ZINC199847048 and FoxM1 active site.
The interactions of ZINC199676474 are shown in [Figure 7]D. The binding energy of this compound with the active site FoxM1 was -4.43 kcal/mol.
Ring A made Pi-sulfur electrostatic and Pi-Pi interaction with His287, while C and D rings formed Pi-alkyl interactions with Arg286. Also, ring D created Pi-alkyl interaction with Leu256 and Pi-Pi interaction with Trp308. Six hydrogen bonds were made between ZINC199676474 and the active site amino acids of FoxM1 at a distance lower than 2.5 Å.
[Figure 7]E shows interactions of ZINC6110690. The binding energy of this compound with the active site of FoxM1 was equal to -3.97 kcal/mol. Ring A formed Pi-sulfur and Pi-Pi interactions with His287. Ring B created a Pi-donor hydrogen bond with Asn283 and an electrostatic bond of Pi-sulfur with His287. Ring D participated in Pi-alkyl interaction with Arg286. Also, several hydrogen bonds were made between compound ZINC6110690 and FoxM1 active site.
Binding mode of the potential dual inhibitors in Hsp70 active site
The interactions of ZINC8448537 binding energy of -2.87 kcal/mol are shown in [Figure 8]A. Rings A, C, D, and E formed Pi-lone pair bond with Gly339, Pi-alkyl with Arg342 and Arg272, Pi-anion electrostatic bond with Glu268, Pi-cation electrostatic bond, and a Pi-alkyl reaction with Lys56 and Pro39, and Pi-anion with electrostatic bond with Asp366, respectively. S25 made a conventional hydrogen bond with crystallographic water at a distance of 2.254 Å. N6 created a carbon-hydrogen bond with Arg272 at a distance of 2.958 Å, and NH11 formed a conventional hydrogen bond with Tyr15 at a distance of 2.864 Å. Also, four conventional hydrogen bonds were shown between ZINC8448537 and Hsp70 active site.
Figure 8: Main interactions between compounds and Hsp70: (A) ZINC8448537; (B) ZINC1152745; (C) ZINC199847048; (D) ZINC199676474; and (E) ZINC6110690. Hydrogen bond (red), water hydrogen bond (blue), conventional hydrogen bond (green), carbon-hydrogen bond and Pi-donor hydrogen bond (gray), Pi-cation and Pi-anion (yellow), and Pi-alkyl (pink).[Figure 8]B shows the interactions of ZINC1152745 with a -3.291 kcal/mol binding energy. Hydrophobic interactions included ring A: Pi-alkyl with Arg342, Ile343, Lys271, Arg272; ring D: Pi-Alkyl with Lys56, and Cl20 and Cl21, alkyl-halogen with Lys56 and Pro39. In addition, electrostatic interactions are shown in ring B, Pi-cation, ring C, Pi-anion, and ring D, Pi-cation with Lys56. Also, (NH2)19 and S25 formed two conventional hydrogen bonds with Tyr15 and crystallographic water at a distance of 2.301 and 2.109 Å, respectively.
The interactions of ZINC199847048 with a binding energy of -2.91 kcal/mol are shown in [Figure 8]C. Ring A made Pi-alkyl interaction with Ile343, Arg342, Arg272, and Lys271. Rings B and D formed Pi-cation electrostatic bond with Lys271 and Lys56; ring D also made a Pi- Pi-alkyl interaction with Lys56. Ring C made Pi-anion electrostatic bond with Glu268. Also, nine hydrogen bonds were shown with Hsp70 active site.
The interactions of ZINC199676474 with a binding energy of -2.88 kcal/mol are shown in [Figure 8]D. Pi-cation and Pi-anion electrostatic bonds were formed between ring B and Lys271 and Glu268. Rings C and D made Pi-anion and Pi-cation interactions with Glu268 and Lys56, respectively. Ring A and D participated in Pi-alkyl interactions with Lys271, Ile343, Arg342, Arg272, and Lys56. Also, six hydrogen bonds were shown with Hsp70 active site at a distance lower than 3 Å.
[Figure 8]E shows interactions of ZINC006110690 with a binding energy of -3.251 kcal/mol. Rings A and D made a Pi-alkyl interaction with Arg342, Ile343, Lys271, Arg272, and Lys56. Rings B and D formed Pi-cation electrostatic bonds with Lys271 and Lys56 residues. Ring C showed a Pi-anion electrostatic bond with Glu268. Also, five hydrogen bonds were displayed in Hsp70 active site, at a distance lower than 3 Å with compound ZINC006110690.
ADMET prediction analysis
QikProp and data warrior programs were used for ADMET prediction of the extracted compounds. Pharmacokinetic parameters which were considered for ligand selection are as follows: logBB, skin permeability coefficient (log Kp), apparent Caco-2, and MDCK permeability (the higher the value of MDCK cell, the higher the cell permeability), aqueous solubility (log S), number of metabolic reactions, logKhsa for serum protein binding, human oral absorption in the gastrointestinal tract, logP for octanol/water controlling hydrophilicity or lipophilicity properties of compounds. The obtained results are shown in [Table 1.] Five predicted compounds with appropriate pharmacokinetic properties passed the screening of ADME prediction. In addition, drug toxicity risk assessment parameters were obtained within acceptable intervals. Also, these compounds did not show parameters related to the risk of mutagenicity, tumorigenicity, and stimulant or reproductive toxicity.
Table 1: Absorption, distribution, metabolism, excretion, and toxicity properties of five proposed compounds for inhibition of FoxM1 and Hsp70.MD studies
MD simulation is one of the primary methods that simulate ligand-receptor interactions in biological environments. MD simulation study was performed in two steps. In the first step, MD was carried out for the five selected ligands and FDI-6 and VER-155008, as references, on FoxM1 and Hsp70 during 10 ns. The obtained MMPBSA energy values of these five compounds were compared with the MMPBSA energy values of the reference compounds, and finally, ZINC1152745 was chosen as the best compound in both proteins. Then, MD simulation in FoxM1 and Hsp70 proteins was performed for 100 ns on ZINC1152745. The resulting trajectories were analyzed to evaluate the system's stability and structural properties after 100 ns MD simulation. The number of hydrogen bonds, the radius of gyration (Rg), root mean square fluctuation (RMSF), RMSD, and free energy were calculated for ZINC1152745 and compared with the reference compounds for FoxM1 and Hsp70 (FDI-6 and VER-155008).
MM-PBSA free energy calculations
MM-PBSA is a method for free energy calculation between protein and ligand by calculating van der Waals, electrostatic, polar solvent, solvent-accessible surface area (SASA), and binding energy. Firstly, for the five ligands obtained from the previous screening, the free binding energy was calculated. Since the increased free energy of a system indicated decreased stability, the results showed that ZINC1152745 had suitable stability with FoxM1 and Hsp70 [Table 2] and [Table 3]. This ligand was simulated for 100 ns with Hsp70 and FoxM1 proteins, as mentioned before. Its free binding energy was calculated, which is shown in [Table 4.] The free binding energy values of ZINC1152745 were lower than of FDI-6 and VER-155008, the known inhibitors of FoxM1 and Hsp70, respectively.
Table 2: The calculated free binding energy for five ligands and FDI-6 with FoxM1 during 10 ns.
Table 3: The calculated free binding energy for five ligands and VER-155008 with Hsp70 during 10 ns.
Table 4: The calculated free binding energy for ZINC1152745, FDI-6 and VER-155008 with FoxM1and Hsp70 during 100 ns.RMSD
RMSD was evaluated to assess the stability and fluctuations of the protein backbone and ligand in complexes. As shown in [Figure 9] and [Figure 10], ZINC1152745 showed the least fluctuations in RMSD means the most stability during the MD simulations compared to FDI-6 and VER-155008. The minor standard deviation commonly indicated the better stability of the simulated system (fluctuation lower than 0.3 nm is suitable for proteins).
Figure 9: The RMSD profile of (A) FoxM1 backbone in complex with FDI-6 and ZINC1152745; (B) FDI-6 and ZINC1152745 as a function of simulation time. RMSD, Root mean square deviation.According to [Figure 9]A, which was related to the RMSD backbone of FoxM1 protein, it can be observed that this protein reached equilibrium and stability after 25 ns. The RMSD for ZINC1152745 ranged from 0.12 to 0.3 nm and for FDI-6 from 0.1 to 0.38 nm. [Figure 9]B shows the RMSD of FDI-6 and ZINC1152745 compounds with the FoxM1 protein. RMSD was between 0.1 to 0.2 nm for FDI-6 and 0.00 to 0.16 nm for ZINC1152745, respectively.
The RMSD backbone of the Hsp70 protein in complex with ZINC1152745 and VER-155008 is shown in [Figure 10]A. A higher stability and equilibrium in the presence of ZINC1152745 relative to VER-155008 was observed from the first nanoseconds. On the other hand, the equilibrium was achieved after 65 ns for VER-155008. The backbone RMSD of Hsp70 was between 0.1 to 0.15 nm in the complex with ZINC1152745 and 0.1 to 0.2 in the complex with VER-155008, respectively.
Figure 10: The RMSD profile of (A) Hsp70 backbone in complex with VER-155008 and ZINC1152745; (B) VER-155008 and ZINC1152745 as a function of simulation time. RMSD, Root mean square deviation.As shown in [Figure 10]B, the VER-155008 in the complex reached equilibrium after 95 ns, while the ZINC1152745 reached equilibrium in the initial nanoseconds. The fluctuation of ligand RMSD was between 0.02 and 0.4 nm for VER-155008 and between 0.04 to 0.16 nm for ZINC1152745.
RMSF
RMSF measurement is the suitable method to characterize the residue stability during MD simulation. This method measures the flexibility of protein in the simulation and could determine the flexibility of residues in the protein. The higher value of RMSF shows the higher amount of mobility of the alpha carbon atoms of the protein in the MD simulations. In the RMSF diagram, averaging of the position is performed per residue.
As shown in [Figure 11], RMSF diagrams of ZINC1152745 and FDI-6 in complex with FoxM1 indicated that the active site residues of Foxm1, including His287, Arg286, Leu259, Asn283, Ser290, and Trp308 had the lowest RMSF indicating the appropriate accommodation of these compounds in the active site.
Figure 11: The RMSF plot of FoxM1-FDI-6 and FoxM1- ZINC1152745. RMSF, Root mean square fluctuation.[Figure 12] is shown that ZINC1152745 was more stable than VER-155008 in Hsp70 active site. The RMSF values were obtained at 0.04-0.35nm for ZINC1152745 and 0.05-0.34 for VER-155008 in a complex with Hsp70. Also, the residues of the active site of Hsp70, which include Try15, Gly202, Gly23, Lys271, Glu268, Ser275, Arg272, Arg343, and Gly339 had the lowest RMSF indicating the appropriate accommodation of these compounds in the active site.
Figure 12: The RMSF plot of Hsp70-VER-1555008 and Hsp70- ZINC1152745. RMSF, Root mean square fluctuation.Rg
The Rg in MD simulation is an indicator of the compactness of the protein. An absolute radius could not be considered for a non-spherical compound. However, some of them are shaped such that they are collectively known as globular or spherical proteins. For such proteins, an approximate radius, called the Rg, could be calculated. During the simulation, if any factor causes the protein to squeeze, the Rg decreases, resulting in protein folding. If a factor causes the protein to outstretch, such as changing solvent, intramolecular interactions of the protein, or denaturation of the protein, the Rg increases and the stability of the protein decreases, causing the protein to unfold.
[Figure 13] shows the diagram of the radius of gyration for the FoxM1 protein complex with FDI-6 and ZINC1152745. The Rg was between 1.1 and 1.3 nm for the FDI-6 ligand and between 0.9 and 1.2 for ZINC1152745, indicating the stability of both complexes.
Figure 13: The radius of gyration plot of the backbone. FoxM1-FDI-6 and FoxM1- ZINC1152745 complexes.[Figure 14] displays the diagram of the Rg for the Hsp70 protein complex with VER-155008 and ZINC1152745. The variation of the Rg was between 1.65 and 1.95 nm for the VER-155008 ligand and between 1.6 to 2 nm for ZINC1152745.
Figure 14: The radius of gyration plot of the backbone of Hsp70-VER-1555008 and Hsp70- ZINC1152745 complexes.Number of hydrogen bonds
The number of hydrogen bonds is important interaction for the stability of FoxM1 and Hsp70-ligand complex. The formation of the hydrogen bond is one of the signs of strengthening the binding and stability of any chemical structure in an active site. [Figure 15] displays the number of hydrogen bonds during the simulation process in the complex of Foxm1 protein with FDI-6 and ZINC1152745. The average hydrogen bonds in FDI-6 and ZINC1152745 with FoxM1 were 2.15 and 1.99, respectively. [Figure 16] also shows the number of hydrogen bonds in the complex of Hsp70 protein with VER-155008 and ZINC1152745. The average hydrogen bonds for VER-155008 and ZINC1152745 complex with Hsp70 protein were 5.54 and 1.74, respectively.
Figure 15: Numbers of hydrogen bonds formed between (A) FDI-6 and (B) ZINC1152745 with FoxM1 during molecular dynamic simulation.
Figure 16: Number of hydrogen bonds formed between (A) VER-155008 and (B) ZINC1152745 with Hsp70 during molecular dynamic simulation. Discussion The critical residues of the active site of FoxM1 involved Asn283, Arg286, and His287, which create a hydrogen bond to the DNA. The essential interactions with the main residues of the FoxM1 protein have been shown in [Figure 17]. The best interaction with the active site of the FoxM1 protein needs sulfur-containing groups in the R1 position for electrostatic interaction with His287. Also, replacing the electron-withdrawing groups in position R2 gave better electrostatic interaction in position R1. Replacing the phenyl ring with electron-withdrawing groups such as halogen, carboxyl, amide, or sulfonamide in the R3 position caused hydrogen bond and hydrophobic interactions with Trp308 and Arg297.
Figure 17: Essential interactions for virtual screening of compounds with FoxM1.On the other hand, the essential residues in the ATP pocket of Hsp70 are Glu268, Lys271, Arg272, Ser275, Arg342, Thr13, The14, and Tyr15. [Figure 18] shows that replacing the thiophene ring in position R1 caused significant electrostatic and hydrophobic interactions and hydrogen bonds with Arg342, Ile343, Lys271, Arg272, and Glu268 residues. Also, electron-withdrawing groups such as CF3 in the R2 position formed hydrogen bonds with Gly201, Gly202, Gly230, and crystallographic water. Replacement of groups such as Cl, COOH, CONH2, and SO2NH2 caused stability of the compound in the active site of Hsp70 protein by creating hydrogen and hydrophobic bonds with Lys56 and Asp53 residues. The NH2 attached to the thiophene ring also made a hydrogen bond with Try 15.
Figure 18: Essential interactions for virtual screening of compounds with Hsp70.ZINC1152745 compound [Figure 6] was selected as the best FoxM1/ Hsp70 dual inhibitor. The binding energy of this compound was calculated as -4.27 and -3.29 kcal/mol in the FoxM1 and Hsp70 active sites, respectively. This compound has suitable substitutions for interaction with FoxM1 and Hsp70 active sites. ZINC1152745 created electrostatic interaction with His287 and hydrogen binding with Asn283 and Arg286, and hydrophobic interactions with Trp308, Leu259 and Leu289 in FoxM1’s active site. Also, this compound formed hydrophobic interactions with Arg342, Ile343, Lys271, Arg272, Lys56, and Pro39, electrostatic interactions with Lys271, Glu286, and Lys56, and hydrogen bonds with Tyr15 and crystallographic water.
ConclusionBased on bioinformatic screening methods, this study proposed five new small molecules as potential dual inhibitors of FoxM1 and Hsp70. A structure similarity search was performed based on 80% similarity using FDI-6, a potent FoxM1 inhibitor, as the template compound and 850 compounds were obtained.
The best docking pose related to FoxM1 with FDI-6 and Hsp70 with VER155008 ligand were selected for designing the pharmacophore model. Then, virtual pharmacophore screening was performed on the created library compounds for FoxM1 and Hsp70 using a Pharmit web server. The second screening was performed by virtual molecular docking using the Autodock program. According to this screening, 120 compounds for FoxM1 and 25 for Hsp70 were obtained with the best interactions and the lowest energy levels. ZINC8448537, ZINC1152745, ZINC199847048, ZINC199676474, and ZINC6110690 were selected for FoxM1 and Hsp70 proteins from two screening (pharmacophore and docking screening). Some ADMET parameters were calculated for these compounds, and all fell in the acceptable ranges for 95% of known drugs. Also, toxicity studies showed that most selected compounds had a low toxicity risk.
MD simulations were performed for ZINC1152745, the compound in FoxM1 and Hsp70 active sites. The results of several hydrogen bonds, the Rg, RMSF, RMSD, and free energy during the simulations showed better stability of ZINC1152745 with FoxM1 and Hsp70. Future experimental investigations on this compound will confirm the predicted computational achievements of this study.
Acknowledgments
We gratefully acknowledge the financial support from Hormozgan University of Medical Sciences through Grant No. 990681 and the Vice-chancellery of Research of Isfahan University of Medical Sciences through Grant No. 397617.
Conflict of interest statement
The authors declared no conflict of interest in this study.
Authors’ contribution
M. Abbasi and G.A. Khodarahmi conceived and supervised the project and manuscript preparation, revision, and final correction; Z. Alimardan performed the experiments, analyzed and interpreted the data; K. Kashfi contribute
留言 (0)