Quantum mechanical/molecular mechanical and docking study of the novel analogues based on hybridization of common pharmacophores as potential anti-breast cancer agents.

In an attempt to identify some new potential leads as anti-breast cancer agents, novel hybrid compounds were designed by molecular hybridization approach. These derivatives were structurally derived from hybrid benzofuran-imidazole and quinazolinone derivatives, which had shown good cytotoxicity against the breast cancer cell line (MCF-7). Since aromatase enzyme (CYP19) is highly expressed in the MCF-7 cell line, the binding of these novel hybrid compounds to aromatase was investigated using the docking method. In this study, due to the positive charge on the imidazole ring of the designed ligands and also, the presence of heme iron in the active site of the enzyme, it was decided to optimize the ligand inside the protein to obtain more realistic atomic charges for it. Quantum mechanical/molecular mechanical (QM/MM) method was used to obtain more accurate atomic charges of ligand for docking calculations by considering the polarization effects of CYP19 on ligands. It was observed that the refitted charge improved the binding energy of the docked compounds. Also, the results showed that these novel hybrid compounds were adopted properly within the aromatase binding site, thereby suggesting that they could be potential inhibitors of aromatase. The main binding modes in these complexes were through hydrophobic and H bond interactions showing agreement with the basic physicochemical features of known anti aromatase compounds. Finally, the complex structures obtained from the docking study were used for single point QM/MM calculations to obtain more accurate electronic interaction energy, considering the electronic polarization of the ligand by its protein environment.

INTRODUCTION

Breast cancer is a life threatening disease in women whose incidence is growing annually (1).

Approximately two-thirds of breast cancer tumors are hormone-dependent, implying that endogenous estrogens are essentially required for their proliferation. Thus one of the most important approaches for the treatment of breast cancer is interfering with endogenous hormone production (23).

Aromatase is always considered as the most promising target for the selective lowering of estrogen levels in patients with estrogen- dependent breast cancer (4).

Joining two or more biologically active pharmacophores in a single molecular framework might result in pharmaceutically important hybrid molecules which could address the active site of different targets and/or offer the possibility of overcoming drug resistance or reducing the unwantedside effects (5). Among the various pharmacophores in medical chemistry, imidazole, benzofuran, and quinazolinone moieties have received substantial attraction because of their critical role in many structures within bioactive compounds (678).

Computational docking of small molecules to macromolecular targets has been widely used in rational drug design and discovery (9). Although molecular docking simulation of proteins is a fast and inexpensive method for the description of ligand-protein interactions, this technique is associated with several limitations (10). Two most important limitations of the conventional docking are: (a) assuming non protein flexibility upon ligand binding, and (b) using force field based fixed dielectric charges for both protein and ligand atoms; therefore false positive or false negative protein–ligand binding energy may be obtained (11). It has been shown that the polarized binding site of proteins could affect the atomic charges on ligand, so the use of fixed dielectric charges for both protein and ligand atoms (especially for metalloproteins with highly polarized binding sites) lead to the low accuracy of the docking results (12). Therefore, to get superior accuracy in the docking studies, it is reasonable to use the corrected charge of the ligand according to the polarized active site environment (12).

Chen, et al. reported that bezofuran- imidazolium hybrids (1) bearing naphthalene moiety showed excellent cytotoxic effects on MCF-7 cell line (7). Additionally, quinazolinone derivatives (2) have also been effective on this cell line (13). According to the bioactivities revealed by quinazolinone and bezofuran-imidazolium analogs, the aim of the present investigation was the design of novel hybrid structures incorporating these moieties into a single molecular scaffold to evaluate their potential additive effects as cytotoxic agents on the MCF-7 cell line (Scheme 1). On the other hand, benzofuran was used in the structure of some potent aromatase inhibitors (14). Since aromatase is overexpressed in the MCF-7 cell line, interaction of designed compounds was studied with aromatase enzyme through the docking simulation. By using the obtained atomic charges through QM/MM calculations, the docking was performed and changes of the free binding energies were evaluated. Finally, via single point QM/MM energy calculation on the complexes obtained from the docking, the interaction energies in which electronic polarization of the ligand by its protein environment was considered were estimated.

Scheme 1

Design of novel analogues based on the hybridization of benzofuran, imidazole, and quinazolinone pharmacophores.

MATERIALS AND METHODS

Docking method

To elucidate the binding mode of aromatase and novel hybrid compounds, molecular docking was performed by AutoDock4 software (15). For this purpose, the atomic coordinate of the protein was obtained from Protein Data Bank (PDB) using PDB ID 3EQM. The protein structure was visualized and all water molecules were removed from the protein structure. Before using the ligands in docking calculations, their structures were optimized using the PM6 semiempirical method and Gaussian 09 quantum chemistry package. AutoDockTools was used to prepare the protein and ligand structures and parameters before submitting them for docking analysis. A grid box size of 60 × 60 × 60 Å points with a grid spacing of 0.375 Å was applied. The center of the grid box was defined as the center of the co-crystallized inhibitor. AutoDock parameter was set and distance-dependent dielectric functions were used for calculating the van der Waals and the electrostatic terms, respectively. After applying the docking protocol, the free energies of binding (ΔGb) and inhibition constants (Ki) were calculated by AutoDock. The best conformers were chosen according to the lower docked free energy and top-ranked cluster then used to perform docking analysis with AutoDockTools and PyMOL. Since simulations of proteins, where ligand binding involved prosthetic groups like heme, posed a great challenge for molecular docking, in this study a set of charges, obtained according the work of Favia, et al. (16) for the heme cofactor, were used in molecular docking.

QM/MM methodology

For QM/MM calculations, we employed the Gaussian 09 quantum chemistry package (17). Our own N-layered integrated molecular orbital and molecular mechanic (ONIOM) was implemented in Gaussian 09 and used for QM/MM calculations (18). The ONIOM scheme is more general in the sense that it can combine any number of molecular orbital methods as well as molecular mechanics method and so, it is considered as the QM/MM method. This method enables different ab initio or semi-empirical methods to be applied to different parts of a system. The interactions between ligand and protein are exclusively non-covalent; therefore, ligand was assumed as the QM region and the protein as the MM region. PM6 semi-empirical method, one of the best semi-empirical methods in quantum mechanics, was used to represent the QM region (ligand) and the universal force field (UFF) was employed for the MM region (aromatase). Therefore, a two layer ONIOM (PM6:UFF) was used for the calculations. The partial atomic charges of the MM region were assigned using the QEq formalism (18) and the polarization of the ligand due to the partial atomic charge of aromatase (charge embedding) was also considered in the calculations. The atomic charges of Fe(II)- heme complex in these QM/MM calculations had been taken from the work by Favia, et al. (16), where the charges were calculated using DFT method and considering B3LYP functional and 6-31G(d) basis set. Therefore, although Fe(II)-heme complex in our ONIOM calculation was in the MM region, its atomic charges were from DFT calculations performed in the mentioned reference (16).

RESULTS

Considering the anticancer activities of naturally occurring quinazolinones as well as the potent cytotoxic activities of synthetic hybrid imidazole-benzofuran derivatives on the MCF-7 cell line (78), novel hybrid compounds bearing both these moieties, were designed by molecular hybridization approach (Scheme 1). Since aromatase is highly expressed in this cell line, the binding mode of these novel hybrid compounds to aromatase (CYP19) enzyme was investigated through the docking method.

The prediction of the binding affinity in a molecular docking tool is estimated by a scoring function which generally needs to be both fast and accurate. However, accurate scoring of protein-ligand interactions remains challenging due to some limitations associated with the docking method. This method has some limitations including the assumption of non protein flexibility upon ligand binding and the use of force field based fixed dielectric charges for protein and ligand atoms, leading to low accuracy in the calculated results. This problem is more important for proteins which involve metal in its active site. In fact, the presence of metal in the active site induces a higher polarization effect and enhances the restriction of docking in the prediction of electronic interactions.

To enhance the accuracy of the results, polarization effects should therefore be considered. QM/MM methods can help to offer a superior estimate of the electronic interactions. Previous studies have shown that the docking program gives better results if the ligand partial charges are refitted with QM/MM (12).

In this study, owing to the positive charge on the imidazole ring of designed ligands and also, the presence of the heme iron in the active site of the aromatase enzyme, polarization effects played an important role in energy calculations. Thus, to consider polarization effects, after obtaining a top scoring structure from the first standard docking protocol with AutoDock, QM/MM calculation was performed and a new set of charge values for ligand atoms, according to the polarized protein environment, were obtained. For this purpose, ligand was treated via PM6 calculation as the QM level and Mulliken population analysis, which was implemented in the QM software, was used to modify the charges on the ligand atoms. Subsequently, with these new charges, another docking protocol was performed and the binding energy was obtained.

Since top-scoring pose in docking is dependent on the charges of the ligand atoms, these two steps were repeated until change in the charge values became insignificant. Subsequently, by the selection of the best complex between ligand and protein, according to its cluster and binding energy, the main interactions were evaluated. The free energies of binding (ΔGb) and inhibition constants (Ki), as estimated by AutoDock, are summarized in Table 1.

Table 1

Chemical structures and the changes of free binding energy (kcal/mol) during the determination of more accurate charges, as well as the inhibition constants (Ki) of ligands calculated by AutoDock.

The presence of the metal in the active site of metalloproteins (in our case aromatase) induced a higher polarization effect. Hence the electronic interaction appeared to be important in determining the optimal state of the ligand in the active site. Electrostatic energies were calculated to provide qualitative insights into polarization effects on these novel designed ligands in the aromatase active site. For this purpose, ten of the complexes with the lowest energy were selected and subjected to QM/MM energy calculations. Electronic interactions were computed for ligands according to ΔEinter = Ecomplex - Eprotein – Eligand, as shown in Table 2. Finally, these ten ligands were extracted from complexes and used for rigid docking in which all rotatable bonds were to be held constant and ligand charges were replaced with the obtained values from QMMM calculation.

Table 2

Electrostatic energy information obtained from QM/MM energy calculations for ten complexes with their free energies of binding (ΔGb) of rigid docking, ascalculated by AutoDock.

DISCUSSION

According to the results, it was observed that hybrid compounds were accommodated properly within the aromatase binding site, suggesting that they could be potential inhibitors for aromatase (Fig. 1A). The aromatase active site was highly hydrophobic and it was dominated by aliphatic and aromatic amino acid residues such as Met 374, Val 373, Val 370, Ile 305, Ala 306, Ile 133, Trp 224, Leu 372, Leu 477, Phe 134, and Thr 310. Consequently, the hydrophobic inhibitors with alkyl or aromatic groups were supposed to bind with high affinity to this enzyme (19). Given the nature of the quinazolinone, benzofuran, and imidazole, it was not surprising that the main binding modes in these complexes would be through hydrophobic interactions. Examination of the best-ranked docking results revealed that among the three hetrocycles located in the substrate cavity, a relatively quinazolinone nucleus was positioned in the vicinity of the heme, which might be able to be coordinated to the heme iron.

Fig. 1

The binding mode of (A) new hybrid scaffold, (B) compound 2, and (C) compound 12 in the active site of aromatase.

As shown in Fig. 1A, this interaction could be through the N1 of the quinazolinone with heme. Additionally, π-π conjugate interactions were observed between the phenyl ring of the quinazolinone and Arg 115. It was found that Val 369, Phe 221, and Trp 224 were able to interact with the phenyl ring substituted on the benzofuran through van der Waals interactions. Benzofuran ring formed hydrophobic interactions with the Phe 134 and leu372 (Fig. 1A). The Ser 478 was in close proximity with the oxygen atom of the benzofuran group. However, no hydrogenbond between oxygen and Ser 478 was predicted by the docking results. Unlike other aromatase inhibitors bearing imidazol ring, the docked protocol did not find the typical imidzolium scaffold orientation for these hybrids molecules toward the heme iron in the binding site (Fig. 1A).

Due to the presence of the positive charge on the nitrogen atom of imidazole, it seemed logical to be accommodated away from the heme iron.

The purine's imidazole ring had hydrophobic interaction with the Thr 310 and Asp 309 residues. Hydrophobic interaction observed for cationic imidazole ring was interesting. In this cation the charge was delocolized and this delocolization led to the distribution of a positive charge on the five atoms of the ring, such that carbon and hydrogen were still able to participate in hydrophobic interactions. Due to the resonance effect in the imidazole ring, both imidazolium nitrogens were positively charged and therefor, none of them could make hydrogen bonds with the residues in the aromatase active site.

It was expected that apart from hydrophobic interactions, substitution on quinazolinone, benzofuran and imidazol groups might create highly complementary packing and better interactions with the enzyme. For this purpose, efforts were made to substitute various alkyl, aryl, halogen, hydroxyl, and methoxy groups on the ligand structures and investigate their effect on the binding mode of the hybrid compounds. The main interactions in these complexes could be described as follows:

The OH group might play an auxiliary role in stabilizing the interaction between the ligand and the aromatase binding site. The substituted OH on the benzofuran ring (compound 2) could make hydrogen bond with the amino acid residue Met 374 of aromatase (Fig. 1B).

The side chain of the Ser 478 formed a hydrogen bond (Fig. 1C) with the OH group on the quinazolinone ring (compound 12). It is interesting to mention that in this structure (Fig. 1C),benzofuran was accommodated in close proximity to the heme iron which allowed the proper orientation of the oxygen carbonyl group of the quinazolinone ring to establish an additional H-bond with Leu 477. Methoxy group, as an electron-donating substituent, could form hydrogen bond with hydrogen bond donors in the active site of the enzyme and increase the aromatase inhibitory effect.

The Met 374 backbone NH and the methoxy group on the benzofuran ring (compound 14) were in close proximity (Fig. 2A), with favorable hydrogen bonding interactions. The NH residue of Arg 115 formed a hydrogen bond (2.0 Å) with the methoxy group on the quinazolinone ring of the compound 13 (Fig. 2B). In the next step, the effect of chlorine substitution was also studied. Compared to the non-substituted structure, chlorine replacement on the quinazolinone and benzofuran rings enhanced aromatase inhibitiory potency through increased hydrophobic interactions with Val 370, Leu 372, Phe 134, and Val 373. Similarly, methyl substituted analogues formed additional hydrophobic interactions, as compared to the unsubstituted analogue with Leu 372, Phe 134, and Met 374 for methylated benzofuran and with Val 370, Ile 133, and Arg 115 for methylated quinazolinone, respectively.

Fig. 2

Binding modes of (A) compound 14, (B) compounds 13 in the active site of aromatase.

Additionally, phenyl and thiophen groups were also substituted on the positions 2 of the quinazolinone ring for possible π-π stacking interactions. According to the obtained results, by introducing these two aromatic groups on the quinazolinone binding affinity of the ligands could be significantly reduced; this could be attributed to the role of ligand steric hindrance in accessing the aromatase binding sit (Table 1).

The calculated electrostatic energies demonstrated that the novel design ligands had a good electrostatic interaction with its protein environment and redocking of these rigid ligands to aromatase showed that all runs were extended to the creation of one cluster; also, the calculated ΔG was increased as much as one to two kcal/mol for the ligands (Table 2). This indicates that if the structure of ligand were optimized in the active sit of enzyme and then the same optimized ligand were de- docked to protein; the lower binding energy could be obtained.

According to the aforementioned results, it was found that in addition to the van der Waals interactions, hydrogen bond may play an important role in the ligand–receptor interactions for aromatase inhibition. Among these designed analogues, compounds 4 and 6 yielded the highest ΔGb, -11.47 and -10.94 kcal/mol, and showed the best performance of Ki, 9 nM and 36 nM, respectively.

CONCLUSION

In the present work, novel potential aromatase inhibitors were designed by the molecular hybridization approach. To identify the interactions of these ligands with the aromatase enzyme, docking simulation was carried out in which charge values for the ligand atoms were refitted according to the polarized protein environment. The results showed that these novel hybrid compounds were adopted properly within the binding site and illustrated good van der Waals and electrostatic interactions with the aromatase enzyme. Finally, when these optimized ligands were used for some rigid docking via charges obtained from QM/MM calculation, some improvement in the calculated ΔGb was observed, thereby suggesting that QM/MM calculations could help to improve the docking result.