In silico identification of potential Hsp90 inhibitors via ensemble docking, DFT and molecular dynamics simulations

Abstract The molecular chaperone heat shock protein 90 (Hsp90) has emerged as one of the most exciting targets for anticancer drug development and Hsp90 inhibitors are potentially useful chemotherapeutic agents in cancer. Within the current study, Hsp90 inhibitors that entered different phases of clinical trials were subjected to Zinc15 structure query to find similar compounds (≥ 78%). Obtained small molecules (1-29) with defined similarity cut-off were docked into ensemble of Hsp90-α NTDs. Docked complexes were ranked on the basis of binding modes and Gibbs free energies as Hsp90 binders (cut-off point; ΔGb ≤ −12 kcal/mol). Top-ranked compounds were subjected to energy decomposition analysis per residue of binding pocket via density functional theory (DFT) calculations in B3LYP level of theory. Subsequent MD simulations of the top-ranked complexes were performed for 100 ns to explore the stable binding modes during a reasonable period in explicit water. Results of molecular docking and intermolecular binding analysis indicated that H-bond, hydrophobic and salt bridge interactions were determinant forces in complex formation. Compounds 19 and 20 were well accommodated in binding pocket of Hsp90 via relatively varied conformations. It was revealed that Asn51 and Phe138 were key residues that interacted stably to 19 and 20. Although primary mechanism of action for proposed molecules are unknown and yet to be explored, results of the present study revealed key structural features for future structure-guided optimization toward potent inhibitors of Hsp90-α NTD. Highlights Hsp90 inhibitors that entered different phases of clinical trials were subjected to Zinc15 based structure query to afford potential enzyme inhibitors 19 and 20. Quantum chemical calculations confirmed docking results and verified pivotal role of a conserved residues (Asn51, Leu103, Phe138 and Tyr139) in making effective hydrogen bonds. MD simulations of top-ranked docked derivatives revealed the achievement of stable binding modes with less conformational variation of 20 than 19 in the active site of Hsp90-α NTD. H-bond, hydrophobic contacts and salt bridge interactions were determinant forces in binding interactions of in silico hits. Resorcinol and isoxazole were important structural motifs of in silico hits in binding to the active site of Hsp90-α NTD. Communicated by Ramaswamy H. Sarma


Introduction
Cancer is described as a large group of diseases characterized by the quick and uncontrolled creation, growth and spread of abnormal cells beyond their usual boundaries. On the basis of WHO reports, cancer is the second leading cause of death globally, accounting for an estimated 9.6 million deaths, or one in six deaths, in 2018 (Bray et al., 2018).
Several factors, such as population growth and ageing as well as the varying patterns of cancer risk factors linked to social and economic development cause the recurrent increase of certain cancers (Wu et al., 2018). In addition to the aforementioned issues, lifestyle dependent behavioral factors such as diet, physical activity, alcohol consumption and generative patterns have also been important in cancer risk and burden (Haggar & Boushey, 2009;Kabat et al., 2015). Such growing up trend exerts tremendous physical, emotional and financial strain on individuals, families, communities and health systems and therefore, there is an urgent requirement for novel drugs with increased selectivity and efficacy to treat or control the progression of different cancers.
Anticancer molecules or chemotherapeutic drugs, have received a great attention in current efforts to treat cancer (Bhayye et al., 2020). These agents act via numerous biological mechanisms that target cells at different phases of the cell cycle. Recent progression in novel drug design techniques and possibility of predicting ligand-target interactions at molecular level, has led to more scientific efforts for discovering and design of new chemicals entities as privileged medicinal structures against various cancer types (Cui et al., 2020).
One of the major bottlenecks to develop anticancer agents is that within a typical chemotherapy, normal cells are affected along with the cancer cells, and this causes serious adverse effects. According to this issue, identification of novel and specific molecular targets which are crucial for cancer cell production and progression is of great importance for cancer therapy. In this context, the molecular chaperone heat-shock protein 90 (Hsp90) has gained interest as a promising anticancer drug target due to its important role in the stability and function of many oncogenic mutant proteins (Yun et al., 2019). Hsp90 is member of Hsps that include a large family of molecular chaperones classified by their molecular weights. Hsp90 is the most abundant chaperone protein in eukaryotes and includes about 1 to 2% of cytosolic proteins (Wainberg et al., 2013). Hsp90 promotes accurate folding of newly synthesized proteins and also facilitates refolding of denatured proteins under a variety of intracellular and extracellular stressful conditions (Sauvage et al., 2017).
Besides its determinant biological roles in maintaining cellular homeostasis, tumor survival and proliferation, Hsp90 has been regarded as an attractive anticancer target due to its overexpression in different types of cancer (Calderwood et al., 2006). Such abnormal expression leads to the promotion of carcinogenesis through regulating correct folding, stability, and function of oncogenic proteins (Calderwood & Gong, 2016). It has also been reported that the ATPase activity of Hsp90 is raised up to 50-fold in cancer cells (Taiyab et al., 2009). The unique biological characteristics of Hsp90 has inspired its exploration as an attractive anticancer target.
The discovery and development of new drugs is considered to be very expensive and time-consuming. In this respect, computer aided drug design (CADD) techniques could be constructive for facilitated performing of different tasks including protein-interaction network analysis, drugtarget prediction, binding site prediction, pharmacophore development, virtual screening, and many others. Consequently, CADD approaches are commonly used to increase the efficiency of the drug discovery and development. Within the current project, Hsp90 inhibitors in different phases of clinical trials were subjected to Zinc15 structure query to identify similar compounds (!78% similarity). Small molecules with defined similarity cut-off (1-29) were docked into ensemble Hsp90-a NTDs (40 PDB retrieved 3 D structures). Binding modes and affinities of the docked ligand-protein complexes were analyzed and top-ranked Hsp90 binders (Cut-off point; DG b À12 kcal/ mol) were subjected to energy decomposition analysis via DFT method in B3LYP level of theory with Def2-TZVPP split basis set. 100-ns MD simulations of the superior complexes (4NH8) were performed through Gromos 43a1 forcefield to determine all bonding and nonbonding interactions and the stability of complex in explicit water. Obtained data were used to propose in silico hits with hypothetical binding modes (Figure 1).

Ligand dataset
Ligand dataset compromised the molecules with !78% similar structures to Hsp90 inhibitors that have entered into different phases of clinical trials (https://www.clinicaltrials.gov and https://www.selleckchem.com). Zinc15 database was used for identification of structurally similar compounds (!78%).
Physicochemical parameters of studied ligands were estimated through web server SwissADME (Daina et al., 2017).

Molecular docking
Docking validation or checking the potential of docking simulation to regenerate the binding pose of a co-crystallographic ligand, was done on 200 Hsp90-a NTDs which afforded 40 validated macromolecular structures (refer to the previous section). Lamarckian genetic algorithm (LGA) of AutoDock4.2 package was used to run ligand flexible docking simulations on candidate ligands (1-29) (Morris et al., 2009;Sanner, 1999). Genetic algorithm run was set 100 with 5 Â 10 6 maximum number of energy evaluations. Autogrid module was used to define 40 Â 40 Â 40 Å 3 or 60 Â 60 Â 60 Å 3 grid boxes with regard to the size of docked ligands.

DFT calculations
Density functional analysis of intermolecular binding energy components was performed on a system including key interacted Hsp90-a residues and a docked molecule. Polar hydrogens of the ligands and residues were optimized using B3LYP/Def2-TZVPP method through heavy atom fixing (HAF) approximation (Fogarasi et al., 1992). All the ligand-residue binding energies were estimated by the same method and basis set. The whole calculations were performed with the ORCA ab initio quantum chemistry package (Neese, 2012).

Molecular dynamics
The all-atom MD simulations were done by GROMACS 5.1.1 computational package (Spoel et al., 2005). Topology files and force field parameters for ligands were generated using PRODRG server (Schuttelkopf & Aalten, 2004). Water molecules were considered as a simple point charge (SPC216) model. Total charge neutrality of the complex system (Hsp90, ligand, and solvent molecules) was assured via incorporating 7 Na þ ions. GROMOS96 43a1 force field was applied to determine all bonding and nonbonding interactions. Steepest descent followed by conjugate gradient algorithms were used to minimize energy of the complex. Subsequent to energy minimization process, position restraint procedure was conducted in association with NVT and NPT ensembles. An NVT ensemble was adopted at constant temperature of 300 K with a coupling constant of 0.1 ps and time duration of 500 ps. Subsequent to temperature stabilization, NPT ensemble was performed in a way that a constant pressure of 1 bar was employed with a coupling constant of 5.0 ps and time duration of 1000 ps. Pressure was maintained at 1 bar with the Parrinello-Rahman barostat (Parrinello & Rahman, 1981). The lincs algorithm for covalent bond constraints was applied (Hess et al., 1997). The particle-mesh ewald (PME) and cut-off methods were used to treat the long-range electrostatic and van der Waals interactions, respectively (Darden et al., 1993). All other bonding parameters were set due to our previous report (Razzaghi-Asl et al., 2018) and MD simulation was performed with monitoring of equilibration by examining the stability of the energy, temperature, and the density of the system as well as the root-mean-squared deviations (RMSDs) of the backbone atoms during 100 ns.

Ensemble docking
Physicochemical parameters of compounds with !78% similarity to clinical trial derived Hsp90 inhibitors (1-29) are  summarized in Table 2. Chemical structures of clinical candidates are achievable via supplementary file A. One of the advantages of ensemble docking is that target IFMs represent realistic conformational trajectories of the protein and may enhance the quality of docking results if several conformational ensembles are taken into acount. With the aim of considering the effect of macromolecule conformations on docking outputs, compounds 1-29 were docked into the conformational ensembles of Hsp90-a NTDs (IFMs) (Carlson, 2002). Top-ranked ligands associated with highest free binding energy (DG b ) per docked Hsp90-a IFM (Figure 2) were determined and further subjected to DFT analysis. It was observed that compounds 15 and 19 were the most frequent binders per docked IFM. Besides binding frequency of top-ranked compounds, it was decided to set up a cut-off point (DG b À12 kcal/mol) and hence pick up compounds that meet the defined requirement (19; 5-(2,4-dihydroxy-5-isopropyl phenyl)-N-ethyl-4-(4-(piperidin-1-yl methyl) phenyl) isoxazole-3-carboxamide) ( Table 3). It should be notified that compound 20 with the chemical name of 5-(2,4-dihydroxy-5-isopropyl phenyl)-N-ethyl-4-(3-(morpholino methyl) phenyl) isoxazole-3-carboxamide was also selected for DFT studies due to similarity to clinical trial derived Hsp90 inhibitors (! 78%) and also higher binding affinities to several docked proteins in spite of less binding frequency as top-ranked ligand per docked Hsp90-a IFM (Figure 2).
In continuation, binding characteristics of the best interacted complexes i.e. 19-4NH8 and 20-4NH8 is represented in Figure 3 and Table 4. All of the selected compounds included isoxazole ring within their central core. Hsp90 inhibitors containing isoxazole scaffold are attended in treatment of cancer as tumor suppressers and previous studies could also propose new isoxazole-based scaffolds through 3 D-QSAR and MD simulations (Abbasi et al., 2018).

Amino acid decomposition analysis and SBR
In order to determine the relative contribution of each interacted Hsp90 residue in top-ranked complexes (Figure 4), intermolecular binding energies of the best conformational pose of compounds 19 and 20 with key residues of the active site were calculated. The binding energies between intended ligands and each residue of the active site were estimated through using B3LYP/Def2-TZVPP basis set and method. It should be notified that not all residues, but the residue for which the energy calculation was performed, entered in each model under evaluation. All the ligand-residue binding interaction energies were calculated through equation 1: E LR stands for residue ligand interaction energy, E R and E L are indicative of electronic energies for residues and ligand.

Complex of compound 19 and Hsp90-a (4NH8)
Intermolecular binding energies of compound 19 with each residue of docked Hsp90-a NTD (4NH8) was estimated in B3LYP level of theory with Def2-TZVPP basis set ( Figure 4).

Complex of compound 20 and Hsp90-a (4NH8)
Intermolecular binding energies of compound 20 with each residue of docked Hsp90-a NTD (4NH8) was estimated in B3LYP level of theory with Def2-TZVPP basis set (Figure 4). On the basis of decomposed binding energies, less repulsive forces were estimated for compound 20 in comparison with 19. For instance, some repulsive amino acids (Phe138 and Asn51) sitting in the adjacency of 19 were supported by attractive binding forces when interacted to 20. Phe138 was associated with À1.43 kcal/mol binding energy and enhanced binding pose might be attributed to the presence of additional hydrophobic contacts with carbon atoms of resorcinol in 20. Spatial distortion of meta-morpholino methyl substituent of compound 20 avoided the formation of polar interaction with Asp93 in 20 when compared to similar interaction in compound 19 (-9.70 kcal/mol) ( Figure 5). Asp93 was also reported to interact with isoxazole based Hsp90 inhibitors before (Abbasi et al., 2017). Leu103 and Tyr139 participated in attractive binding forces with À8.31 kcal/mol and  À5.45 kcal/mol, respectively. In a similar manner to the complex of 30 with Hsp90-a NTD (4NH8), Tyr139 made tight interactions to 20 through H-bond between O30-H of resorcinol ring and side chain hydroxyl oxygen of the residue (-9.95 kcal/mol). Tyr139 is an important amino acid in binding interactions to the Hsp90-a active site and amino acid decomposition analysis prioritized the binding ability of topranked derivatives to this residue as 20 > 19. From the ligand-protein interaction patterns one possible interpretation would be the orientation of isoxazole ring of 19 and in particular 20 (Figure 6) with regard to Tyr139 side chain and desired hydrophobic contacts between aromatic rings. Asn51 is another key residue that participated in effective H-bond interaction with 20 (-9.95 kcal/mol).

Molecular dynamics
Docking models would be enhanced in quality and reality via MD simulations since flexibility of both ligand and macromolecule are taken into consideration during a reasonable time (Durrant & McCammon, 2011). The stability of topranked docked complexes was elucidated during 100-ns MD simulations in the presence explicit water model. Root mean square deviation (RMSD), root mean square fluctuations (RMSF) and radius of gyration (Rg) were considered as the stability criteria of the evaluated complexes.
RMSD plot vs simulation time indicated that approximate leveling off in deviation from primary Hsp90 structure was gradually achieved after about 50 ns (Figure 6a). Relatively low RMSD average value (2.68 Å) during 100 ns MD simulations revealed the stable binding of 19 to Hsp90 and indicated the convergence of the protein to equilibrium structure. In addition to protein backbone variations, ligand fluctuations were also analyzed in terms of RMSD plot (Figure 6a). RMSD profile of the interacted ligand was relatively stable up to about 30 ns and after a 15-ns descent, a stable condition was re-achieved and continued to 75 ns. The last descent started from 75 ns and persisted up to the end of simulation time (Average RMSD 2.13 Å).
In the case of compound 20, flattening of the RMSD in deviation from primary Hsp90 structure was achieved at longer time (45 ns) (Figure 6b) but the average RMSD value of the complex (2.03 Å) was better than compound 19. Moreover; in the case of compound 20, leveling off the   RMSD during 100 ns (average RMSD 1.41 Å) confirmed that the ligand could preserve its initial binding mode and location within the Hsp90-a active site. In this case, the steady state persisted well for the whole MD period and just at 13530, 29180, 71280 and 96590 ps, maximum RMSD values of about 2.00 Å could be observed.
The flexibility of individual residues during MD simulations can be estimated via average atomic mobility of the Hsp90 Ca atoms through RMSF calculation (Figure 7a). It was recognized that main residual fluctuations (>3 Å) occurred to the amino acids far from binding site with average RMSF values of 1.32 and 1.37 Å for complexes of 19 and 20, respectively. To explain more, very few points exhibited RMSF >3 Å with Gly125 (19; RMSF 5.03 Å), Lys224 (20; RMSF 4.75 Å), Gln85 (20; RMSF 3.87 Å) Glu16 (20; RMSF 3.86Å), Gln85 (19; RMSF 3.71Å), 177 (19; RMSF 3.53 Å), 175 (19; RMSF 3.32Å) and Lys224 (19; RMSF 3.24 Å) showing maximum values, respectively. The narrow range of estimated RMSFs for binding site residues confirmed the RMSD results and demonstrated the capability of three docked molecules and especially 20 (due to their longer RMSD steady state) in forming stable interactions with Hsp90 during MD simulation. In the case of 20, docked conformation seemed to be stable during the MD run while for 19, a relatively new binding mode could be envisaged.
The compactness of a protein during MD simulation can be analyzed by radius of gyration (Rg) which is defined as the RMSD from each atom of the protein to their centroid. It was found that residual backbones and folding of Hsp90 were steadily stable after binding to corresponding ligands. Average RMSF values were 1.63 and 1.62 for Hsp90 in complex to 19 and 20, respectively. The Rg of the backbone atoms of Hsp90 was decreased from about 1.65 nm at the beginning to about 1.62 nm after 5 ns in the presence of ligand for all of the systems (Figure 7b). Maximum value of Rgs were estimated to be about 1.67 nm for the inhibitor 19 (green curve) between 69 and 73 ns.
To explore the stable binding interactions of simulated compounds, corresponding binding trajectories were retrieved after 100 ns MD simulations (Figure 8). Compound 20 was well accommodated via relatively stable binding mode in the active site of the target. Similar binding trajectories were not retrieved for compound 19 due to the shorter leveling off in RMSD plot (refer to Figure 7a, blue curve).
In confirmation to the RMSD plots (Figure 7b), conformational variations of compound 20 during MD simulation is depicted in Figure 9. A characteristic feature was the orientation of morpholine ring to permit a new salt bridge formation with carboxyl moiety of Asp102. Resorcinol OHs turned toward Tyr142 and Leu107 to make hydrogen bonds with these residues. Central isoxazole ring did not participate in any interactions. Some residues that participated in H-bond interactions (Leu103) and hydrophobic contacts (Val150 and Trp162) within docked conformation were not involved in any interactions with the stable binding mode of the ligand.

Conclusion
Cooperative application of molecular simulation methods may be effective in recognizing the role of biological molecules in human disease, and moreover making predictions on binding affinity and selectivity of various molecules. Structure-based query of experimentally validated Hsp90-a NTD inhibitors afforded a few compounds that were subjected to ensemble docking toward different IFMs of Hsp90a NTDs. Compounds 19 and 20 were found to have tighter binding interactions to induced fit models of target. Within the Identified in silico hits, compound 20 interacted more stably during longer MD time in the binding pocket of Hsp90-a NTD. Conserved interactions to some residues such as Asn51, Leu103, Phe138 and Tyr139 in top-ranked compounds regardless of interaction type revealed the key role of these amino acids in complex formation. Analysis of binding trajectories indicated that compound 20 showed more persisted conformational stability with regard to 19 during 100 ns MD simulations and may be an appropriate candidate for further evaluations. Although primary mechanism of action for proposed molecules are unknown and yet to be explored, but analysis of assessed derivatives in complex with Hsp90-a revealed key structural factors for future structure-guided optimization toward potent binders of Hsp90-a NTD as a validated and selective target against cancer.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by the Ardabil University of Medical Sciences under Grant IR.ARUMS.REC.1398.112.