1. Introduction

The heat shock protein 90 (HSP90) is a chaperone protein responsible for maturation and activity of a variety of key client proteins involved in cellular processes [1] [2] [3]. This protein is overexpressed in tumour cells, and causes uncontrolled proliferation of transformed cells [1] [4] [5]. HSP90 function is dependent on its ability to bind and hydrolyse ATP at the N-terminal domain [1] [4] [5]. An intrinsic ATPase activity is required for the operation of a functional chaperone cycle leading to the stabilization of client proteins [1] [2] [6].

The investigation of HSP90 as anticancer target has arisen from evident progress in clinical trials using the natural inhibitors geldanamicyn (GMD) and radicicol as well as the potent inhibitor 17-allylaminogeldanamycin (17-GMD) [7]. However, several toxicities, extensive metabolism and cell resistance hampered the clinical progression. There is, therefore, a growing interest in design of new therapeutic inhibitors to overcoming these limitations [1] [4] [7] [8]. This spurred development of new classes of pharmacophores using varied pharmacological approaches. Computer-assisted drug development is currently the most used approach as it requires less time and capital. Computer-assisted drug development based on molecular docking, molecular dynamic simulation and quantitative structure activity relationships (QSARs) have been used to predict the antagonist activity of various heterocyclic compounds against the target protein HSP90 [4] [8]. Investigations on compounds containing selenium are of some interest as it has been previously reported that selenium reduced the risk of numerous cancers like prostate, colon, live, lung and mammary cancers [9].

The present work explored the binding ability of GMD selenoderivatives at the N-terminal domain of human HSP90 using molecular docking. HSP90 was derived from the Protein Data Bank (PDB code: 1YET). The docking was performed with the GOLD 5.6 program (Genetic Optimization for Ligand Docking) implemented in the Cambridge Structural Database (CSD). Gold fitness scores and interaction energies were calculated to determine the affinity and the thermodynamic stability of the binding, respectively.

2. Materials and Methods

2.1. Preparation of HSP90

The crystal structures of human HSP90 N-terminal domain with the co-crystallized inhibitor GMD were retrieved from the Protein Data Bank (PDB codes: 1YET, 1.9 Å) and imported into Accelrys Discovery Studio 2019 and visualized with Hermes 1.6 implemented in the Cambridge Structural Database (CSD). Survey of the inhibitor binding site was carried out with the reference of amino acid residues of the ATP-N-terminal domain as previously reported [4]. Hydrogen atoms were added to the protein for correct ionization and tautomeric states of amino acid residues. Co-crystallized water molecules, forming one or more hydrogen bonds to either the protein or the ligand in the protein binding site, were kept in the calculation, whereas the remaining water molecules and inhibitors were extracted before the docking. HSP90 was prepared for docking without minimization of protein energy. The molecular docking was performed after conversion of the HSP90 using GOLD 5.6 program implemented in CSD.

2.2. Preparation of Ligands

Geldanamycin (GMD) was the co-crystallized ligand in the inhibitor site of HSP90. The selenoderivatives of GMD, labelled as GMDSe (i) were designed using Mercury 3.10 program implemented in CSD by replacing an oxygen atom (O) according to Scheme 1 by a selenium atom (Se) in geldanamycin at the position given by the i label $\left(i=1,2,\cdots ,n\right)$.

Scheme 1. Geldanamycin (GMD) et its GMDSe (i) derivatives (i = 1, 2, ∙∙∙, n).

The ligand conformations were generated and optimized after structure conversion from 2D to 3D as illustrated by the GMDSe1, GMDSe6, GMDSe7and GMDSe8 examples in Figure 1 using Mercury 3.10. Only energetically least conformers were selected for further use.

The molecular electrostatic potentials (MEP) of ligands were calculated at the Hartree-Fock theory (HF) implemented in the Gaussian 09 program [10] using the basis sets 6-31G.

2.3. Molecular Ligand-Protein Docking

The GOLD 5.6 genetic algorithm was used to generate the bioactive binding poses of GMD selenoderivatives in the active site of protein HSP90. Before docking HSP90 with ligands, the reproduction of experimental poses using the N-terminal domain of 1YET and 49 additional HSP90 complexes with co-crystallized ligands derived from PDB (Table 1) was performed to validate the program GOLD 5.6. A good reproduction was observed, if the value of Root Mean Square Deviation (RMSD) was lower or equal to 2 Å [11].

The docking of HSP90 with the GMD selenoderivatives was performed using the default GOLD fitness function ChemPLP by considering the flexibility of the ligand and specific amino acid residues of the protein site. Greater the ChemPLP fitness score, better the binding affinity. Rescoring with Chemscore was performed in order to determine the free enthalpy of ligand binding (ΔG_bind) according to the equation below (Scheme 2).

Scheme 2. The equation describes the free enthalpy of binding [12]. The ΔG₀ (−5.48 kcal∙mol⁻¹) and S (bar) components are experimental and entropic terms of the binding energy, while ΔG_HB, ΔG_Met, ΔG_Lip, ΔG_Rot represent energy terms related to hydrogen-bonding, metal, lipophilic and rotation interactions, respectively.

Each simulation was performed ten times yielding ten docked conformations unless three of the ten poses were within 1.5 Å RMSD of each other. Conformer pose with the lowest energy was considered as the binding conformation in the protein site. The program Accelrys Discovery Studio 2019 was used to model non-bonded polar and hydrophobic contacts in the inhibitor site of HSP90.

3. Results and Discussion

3.1. Validation GOLD 5.6 Potential for Ligand-Protein Docking

The distribution of RMSD values derived from the docking of 50 HSP90 examples with co-crystallized ligands was performed to validate use of the GOLD 5.6 program. In addition, RMSD values calculated with GOLD 5.6 for the docking of 1YET, 2BMS and 2VCI were compared to those previously reported by Lauria et al. [13] using the program IFD (Induced Fit Docking). The superposition of experimental poses of the co-crystallized ligands and those derived from the molecular docking with GOLD 5.6 was also carried out in order to validate the potential of GOLD 5.6 to modelling ligand-protein interactions.

The distribution of RMSD values determined for 50 examples of HSP90 complexes is shown in Figure 2.

As it can be seen from the Figure 2, experimental poses of co-crystallized inhibitors were reproduced by the docking with GOLD 5.6 for 44 HSP90 examples (88%) with RMSD values less than 2 Å, which is recommended by the literature as a limit value for a good pose reproduction [11] [14] [15] [16] [17] [18].

Table 2 reports the RMSD values from docking of 1YET, 2BMS and 2VCI using GOLD 5.6 program and those reported by Lauria et al. [13] using the program IFD (Induced Fit Docking).

The results of docking of 1YET, 2BMS and 2VCI complexes listed in Table 2 reveal that the values of RMSD calculated with GOLD 5.6 were smaller than those with IFD [13] suggesting GOLD 5.6 as a powerful program for ligand-protein docking.

Figure 3 illustrates the experimental pose of co-crystallized GMD from PDB (1YET) and typical superposition of the pose derived from the molecular docking with GOLD 5.6.

The experimental poses and the poses derived from the docking were superposed with RMD values less than the 2 Å cut-off level as illustrated in Figure 3(b). This also confirms the potential of the GOLD 5.6 1 program to modelling ligand-protein interactions.

3.2. HSP90 Docking with GMD Selenoderivatives

Figure 4 shows the contact areas involved in the interaction between the N-terminal domain site of human HSP 90 and the selenoderivatives of GMD as derived from the docking with GOLD 5.6.

Figure 4 demonstrates that hydrogen-bonding and lipophilic interactions were involved in the bonding of the selenoderivatives of GMD to the inhibitor site of HSP 90. The affinity of the ligands to the binding site values was determined based on calculation of GOLD fitness scores and interaction energies. The values of GOLD fitness scores, interaction energies and RMSD and are collected in Table 3.

The results of docking in Table 3 show that the interactions were generated with values of RMSD less than 2 Å. Values of GOLD fitness scores and free enthalpies of ligand binding to the HSP90 site (ΔG_Bind) suggest that GMD as well as selenoderivatives interacted with high affinity to the N-terminal domain of HSP90 leading to thermodynamically stable complexes. According to the GOLD fitness scores, the binding affinity of HSP90 to the ligands decrease in the sequence GMDSe6 > GMDSe7 > GMD > GMDSe8 > GMDSe5 > GMDSe2 > GMDSe4 > GMDSe9 > GMDSe3 > GMDSe1, while the decrease sequence is GMDSe6 > GMDSe7 > GMD > GMDSe8 > GMDSe5 > GMDSe2 > GMDSe4 > GMDSe3 > GMDSe9 > GMDSe1 according to free enthalpies of binding (ΔG_Bind). GOLD fitness scores as well as ΔG_Bind suggest GMDSe6 and GMDSe7 as the best ranked ligands, while GMDSe1 as the last ranked. The values of the ΔG_HB and ΔG_Lip component terms related to hydrogen-bonds and lipophilic interactions, respectively, reveal that these interactions were the main forces stabilizing ligand-protein complexes.

The ΔG_Lip component term related to lipophilic interactions seems to be enhanced for the binding of GMDSe7, GMDSe6 and GMDSe2 to HSP90 suggesting that the replacement of oxygen atoms by divalent selenium atoms increases their hydrophobic contact areas. One could also expect that covalently double-bonded selenium atoms in GMDSe7, GMDSe6, GMDSe3 and GMDSe2 interact with the O and N nucleophile atoms of amino acid residues as divalent selenium atom is a sigma-hole donor [19]. However, sigma-hole interactions, well known as short, weakly attractive and linear contacts between regions of lower electronic density (σ-hole) of donor atoms and electron lone pairs of nucleophile atoms acting as sigma-hole acceptors [20], were not found in the inhibitor site of HSP90.

It is noteworthy that ΔG_Bind markedly decreases when replacing the oxygen atom in the O¹=C or O⁹=C group of GMD by monovalent selenium atom in the Se¹=C or Se⁹=C group to obtain GMDSe1 and GMDSe9, respectively. Indeed, the value of ΔG_Bind for the binding of GMD decreases in comparison with the values of ΔG_Bind for the binding of GMDSe1 and GMDSe9 from −39.65 kcal∙mol⁻¹ to 31.91 kcal∙mol⁻¹ (19.5%) and −35.80 kcal∙mol⁻¹ (9.7%), respectively. Likewise, the value of the ΔG_HB component term decreases from −30.19 kcal∙mol⁻¹ for the interaction of HSP90 with GMD to −19.79 kcal∙mol⁻¹ (34.4%) and −26.86 kcal∙mol⁻¹ (11%) for the interactions with GMDSe1 and GMDSe9, respectively. Bibelayi et al. [21] also reported that monovalent selenium atom is less hydrogen-bond donor than oxygen atom. The authors demonstrated based on crystallographic data derived the CSD and quantum mechanical calculations that the ability of the monovalent Se to forming hydrogen bond stems from resonance-induced C^d+=Se^d− dipoles as illustrated for Se in Scheme 3 for [C,N]C=Se, contrarily to hydrogen-bonding ability of oxygen atom, which is ascribed to C^d+=O^d− permanent dipoles.

Scheme 3. Electron delocalization in (N, C) C=Se system.

The hydrogen-bonding can be ascribed to electrostatic interactions between hydrogen-bond donors and acceptors of polar amino acid (AA) residues present in the HSP90 binding site and those of GMD and its selenoderivatives, while lipophilic interactions are probably due to the contacts between hydrophobic amino acid residues of HSP9 and non-polar areas of the ligand molecules. Contact areas of hydrogen-bonding and lipophilic interactions in GMD and its selenoderivatives were evidenced by molecular electrostatic potentials (MEP) as illustrated in Figure 5.

The electrostatic potentials calculated on 0.001 molecular surfaces of GMD and its derivatives in Figure 5 feature marked positive (in blue), negative (in red) and non-polar (in grey) regions corresponding to hydrogen-bond donor, hydrogen-bond acceptor and hydrophobic areas, respectively. Figure 6 reports typical hydrogen-bonding and lipophilic contacts involved in the interaction of HSP90 with GMD.

The activity of several potential HSP90 inhibitors has been predicted based on molecular docking in conjunction with molecular dynamic simulation and quantitative structure activity relationship(QSAR) studies [22] [23]. Indeed, therapeutic properties of HSP90 inhibitors have been ascribed to lipophilic and hydrogen-bonding contacts involving key amino acid residues [4] [13]. The modelling of contact areas like illustrated in Figure 6 shows lipophilic and hydrogen-bonding interactions involving typical hydrophobic and polar amino acid residues like Asp93, which is crucial for ATP activity at the N-terminal domain of HSP90 [4] [13]. The modelling of binding interactions also showed that GMD was bonded to the hydrophobic amino acid residues like Leu48, Met98, Leu107, Phe138, Val150 and Val186 as well as the polar amino acid residues Asp93, Lys58, Lys112 and Phe138 in the inhibitor site. Interestingly, the same amino acid residues were also involved in the lipophilic and hydrogen-bonding interactions of HSP90 with all GMD selenoderivatives, except GMDSe5.

Table 4 reports hydrogen-bond donors and acceptors involved in the interactions of HSP90 with GMD and the GMDSe1, GMDSe2, GMDSe3, GMDSe4, GMDSe6, GMDSe7, GMDSe8 and GMDSe9 selenoderivatives. The interaction distances (δ) and angles (θ) were performed according to Scheme 4. The value vdW(H) = 1.10 Å determined by Rowland et al. (1996) [24] was used for van der Waals radius of hydrogen atom, while the values vdW(O) = 1.52Å for oxygen and vdW(Se) = 1.90 Å for selenium were taken from Bondi et al. [25].

Scheme 4. Geometric parameters of the hydrogen-bonding: interaction distance (𝛿) and hydrogen-bond angle (θ).

Table 4 shows that the GMDSe1, GMDSe2, GMDSe3, GMDSe4, GMDSe6, GMDSe7, GMDSe8 and GMDSe9 derivatives formed hydrogen-bonds with the key amino acid residues of HSP90 like GMD did [13]. The values of interaction distances and angles agree well with the geometric parameters of a typical hydrogen-bond. Indeed, θ values greater than the 120˚ recommended limit [26]. The δ bond distances are less than the r_H + r_A sums of the van der Waals radii of the atoms involved in the hydrogen-bonds (r_H + r_O =2.62 Å and r_H + r_Se = 3.0 Å).

4. Conclusion

Earlier clinical trials demonstrated the opportunities for using inhibition of HSP90 as target protein in the treatment of cancer. The present work was an additional support for recent advances in use of the molecular docking as a powerful tool for design and discovery of inhibitors with improved pharmacological properties. Indeed, the study revealed the ability of GMD selenoderivatives to interact with HSP90. The thermodynamic stability of the binding was mainly ascribed to hydrogen-bonds and lipophilic interactions with amino acid residues in HSP90 inhibitor site. The GMDSe6 and GMDSe7 derivatives were the best ranked GMD derivatives. They were bonded to HSP90 stronger than GMD, while maintaining lipophilic interactions and hydrogen bonds with amino acid residues like Asp93, which are catalytically crucial for therapeutic properties of HSP90 inhibitors. Inhibitors containing Se atoms may provide improved therapeutic properties as the anticancer activity of selenium has been previously reported in the literature. Therefore, a pharmacophore study of the GMDSe6 and GMDSe7 compounds should be useful to explore their therapeutic properties. Molecular dynamic simulation and QSAR studies based on ligand chemical and quantum descriptors as well as experimental investigations are suggested to assess the inhibitor activities.

Acknowledgements

Kilembe Thambwe Jason, Bibelayi Dikima Didi and Lundemba Singa Albert are grateful to the Cambridge Crystallographic Data Centre (CCDC) for funding Master and PhD fellowships, and additional financial and material support.