Design of novel potential aromatase inhibitors via hybrid pharmacophore approach: docking improvement using the QM/MM method

Abstract

Considering the potent cytotoxic activities of hybrid benzofuran–imidazolium and quinazolinone derivatives on the breast cancer cell line MCF-7, novel hybrid derivatives incorporating benzofuran, imidazole and quinazolinone pharmacophores were designed using a molecular hybridization approach. Since aromatase is highly expressed in the MCF-7 cell line, we tried to put these pharmacophores together in such a way that they would arrange themselves in a symmetrical shape, similar to aromatase inhibitors. Subsequently, the binding of these novel hybrid compounds to aromatase have been investigated using a docking procedure applying a combined quantum mechanical/molecular mechanical (QM/MM) method. The QM/MM calculation was performed on the reference structures to obtain atomic charges on the ligand atoms. The results indicated that the hybrid compounds were adopted properly within the aromatase binding site, suggesting that they could be potential inhibitors of aromatase. These novel designed compounds engage in hydrophobic and H-bond interactions with the aromatase binding site, which are in agreement with the basic physicochemical features of known aromatase inhibitors. To obtain more accurate results for the binding energies of the ligands, the structures of the ligands with the best interaction energies, obtained from the docking study, were re-optimized using a three-layer ONIOM method (QM:QM:MM) in which the binding pocket of the enzyme was considered as medium-level. The results demonstrated that when the optimized geometrical structures were subjected to re-docking, a better interaction energy was obtained, which strengthens the ability of these compounds to act as potential inhibitors of aromatase.

---

1. Introduction

The application of informatics to the discovery, design and optimization of biologically active compounds, which is referred to as computer-aided drug design (CADD), has emerged during recent decades as an important tool to help and complement experimental information. This is because computer simulations provide a systematic and economical tool to draw biologically relevant conclusions and also propose new hypotheses based mainly on computer-generated data. A large number of techniques, which differ both in accuracy and speed, are available for the modeling of different phenomena.^1–3

One of the more computationally inexpensive methods for predicting if and how a ligand will bind to a protein binding pocket, followed by an estimation of the strength of the ligand binding affinity, is molecular docking. This aims to achieve a relative orientation of the protein and ligand such that the free energy of the overall system is minimized. Given the biological and pharmaceutical significance of molecular docking, considerable efforts have been directed towards promoting use of this method.^4–6 In the earliest docking approaches, both the ligand and the protein were treated as rigid bodies, while later semi-flexible docking was used, in which the ligand is treated flexibly by allowing bonds to rotate.^7,8 Molecular docking simulation of proteins is a fast and inexpensive method for description of the ligand–protein interactions, but poses some difficulties.⁵ Two of the most important limitations of conventional docking are the assumption of no protein flexibility upon ligand binding and the use of a force field based fixed dielectric charges for both protein and ligand atoms and therefore having false positives and/or false negatives in the energetic quantification of protein–ligand binding.⁹ For the first problem, it should be mentioned that flexible protein docking methods, which treat the protein in a flexible manner, require high computational costs.¹⁰ About the latter, it has been shown that assuming fixed dielectric charges for both protein and ligand atoms leads to low accuracy in protein–ligand docking results.¹¹ Therefore, to increase the accuracy in the docking result, it is reasonable to expend additional effort to improve the quality of the charge. The comprehensive study of polarization and charge transfer requires quantum mechanical methods but their use in biological models requires sophisticated computing systems.^12,13 To limit the computational complexity, combined quantum mechanical and molecular mechanical (QM/MM) methods have been developed as an economical approach, in which a small portion of the protein–ligand system is treated in QM detail.^14–18 It has been shown that by using an ab initio quantum chemical approach via QM/MM methods, a better assumption of the ligand charges, which take polarization into account, is obtained. Therefore this method could be used as a promising start towards the development of more accurate docking methods for lead optimization applications.¹¹

Molecular hybridization as a rational approach for drug design has attracted much attention of researchers to discover new chemical entities with a potential to afford some promising drugs of the future.^19,20 In this method, active compounds and/or pharmacophoric units which recognize and derive from known bioactive molecules are fused to each other directly or with spacer. This method has been employed to develop new anticancer, anti-alzheimer, and antimalarial agents.^21–23

Imidazole, as one of the most important pharmacophores in medical chemistry, has a critical role especially in antifungal, anti-bacterial, anti-cancer and sedative drugs.²⁴ Its biocompatibility provides a scaffold for the preparation of different derivatives to afford new bioactive compounds.²⁴ It has also been reported that imidazolium salts show potent cytotoxic activities towards different cancerous cell lines.^25–27 For example, benzofuran–imidazolium hybrid compounds (1, 2) have good cytotoxic activity towards the MCF-7 cell line (Scheme 1).^26,27 Quinazolinone is another pharmacophore which has been explored for developing pharmaceutically important molecules. Its derivatives have drawn considerable attention due to their profound chemotherapeutic properties including anticancer, antiinflammatory, anticonvulsant, and antidiuretic activities.^28,29

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

Cytochrome P450 enzymes, as a family of isozymes containing heme cofactors, are responsible for many critical enzymatic reactions in living systems. They are, in general, the terminal oxidase enzymes in electron transfer chains.^30–33 Cytochrome P450 19A1, commonly known as aromatase, is an enzyme of the cytochrome P450 superfamily that catalyzes the final and rate-limiting step of the conversion of androgens, testosterone and androstenedione, into estrogens, estradiol and estrone, respectively. It is comprised of a polypeptide chain of 503 amino-acid residues and a prosthetic heme group at its active site. An androgen-specific cleft consisting of hydrophobic and polar residues is situated within the confines of the aromatase binding site. This cleft is specific for androstenedione binding to catalyze conversion of androgen to estrogen via a three-step process. Each step requires one mol of O₂, one mol of NADPH and NADPH cytochrome reductase. This reaction converts androstenedione, testosterone and 16α-hydroxytestosterone to estrone, 17β-estradiol and 17β, 16α-estriol, respectively. The two initial steps are the typical C₁₉-methyl hydroxylation, while aromatization of the steroid A-ring is catalyzed in the final step.³⁴ To block estrogen production, it is necessary to inhibit the enzyme through the use of aromatase inhibitors.³⁵

According to previous studies, benzofuran–imidazole analogs, as well as quinazolinone derivatives (2), have good cytotoxic effects on the MCF-7 cell line (Scheme 1).^26,27,36,37 Therefore, combining quinazolinone and benzofuran–imidazole moieties together as a potential cytotoxic agent towards this cell line would be logical. On the other hand, benzofuran has been used in the structure of some potent aromatase inhibitors.^38–40 For example, Whomsley et al.³⁸ identified substituted l-[(benzofuran-2-yl)-phenylmethyl]-imidazoles (3) as a class of potent aromatase inhibitors with in vitro IC₅₀ values <10 nM, which is 80–1000 times the inhibitory activity of aminoglutethimide (Scheme 1). Since aromatase is overexpressed in the MCF-7 cell line, we want to rationally design a novel structure that incorporates these moieties into a single molecular scaffold which could act as an aromatase inhibitor. According to a previous study,⁴¹ besides the physicochemical properties of the inhibitor, molecular shape is expected to be extremely important for the access and fit within the active site of the aromatase. Most aromatase inhibitors used to build the common features model have a similar shape that is expected to be complementary to the volume of the aromatase active site. For example letrozole (Scheme 1), a quite rigid third generation inhibitor of aromatase, has a high degree of symmetry. To obtain a final symmetrical shape, which consists of three pharmacophores, the best arrangement is to put each of the pharmacophores on the vertices of a triangle. As the basic physicochemical features of known aromatase inhibitors are a high degree of hydrophobicity and the potential to establish hydrogen bonds,^42,43 it is expected that this structure with heterocyclic rings will provide the desired hydrophobicity. On the other hand, hetero atoms in the ring and substitutions located on the structure provide the hydrogen bond potential for this structure.

After designing them, the energetic and structural properties of these compounds were investigated as aromatase inhibitors using molecular docking and QM/MM calculations. 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, QM/MM calculations were performed, in which the partial charges of the ligand were re-fitted according to the polarized active site environment of the enzyme to increase the accuracy of the docking results. Then, the structures with the best interaction energy obtained from the docking study were optimized using a three-layer ONIOM method and finally re-docked to the 3D structure of enzyme to obtain the interaction energy. The results demonstrated that these novel designed compounds have good interaction energy and could be used as potential aromatase inhibitors.

2. Experimental

2.1. Docking method

Molecular docking was performed using the AutoDock4 software to elucidate the binding mode of aromatase with the novel hybrid compounds.⁴⁴ The atomic coordinates of the protein were obtained from the Protein Data Bank using PDB id 3EQM. The protein structure was visualized and all water molecules were eliminated from the protein structure. The structures of the ligands were optimized using the PM6 semi-empirical method in the Gaussian 09 quantum chemistry package. AutoDock Tools was used to prepare the molecules and parameters before submitting it for docking analysis with AutoDock. For this purpose, polar hydrogen atoms were added while non-polar hydrogen atoms were merged, assuming a physiological pH value of 7.0. Then, Gasterier partial atomic charges were assigned to the ligands. All rotatable bonds of the ligands, defined by the default of the program, were allowed to rotate during the automated docking process and the prepared protein and ligand structures were saved in the PDBQT format suitable for calculating energy grid maps. A grid box size of 60 × 60 × 60 Å points with a grid spacing of 0.375 Å was considered, with its center defined as the center of the co-crystallized inhibitor. The Van der Waals and electrostatic terms were calculated using the AutoDock parameter functions. The initial position, orientation and torsions of the ligand molecules were set randomly. Each docked compound was derived from 100 independent docking runs that were set to terminate after a maximum of 2.5 × 10⁶ energy evaluations with a mutation rate of 0.02 and crossover rate of 0.8. The population size was set to use 250 randomly placed individuals. The search for low-energy binding orientations was performed by the Lamarckian genetic algorithm using a translational step of 0.2 Å, a quaternion step of 5 Å and a torsion step of 5 Å. To conclude, the free energies of binding (ΔG_b) and inhibition constants (K_i) were calculated using AutoDock. Several clusters and binding energies were obtained for the docked hybrid compounds, in which the best conformers were selected according to the lower docked free energy and top-ranked cluster to perform docking analysis with AutoDock Tools and PyMOL. To evaluate the validity of the docking process, the substrate of the enzyme, androstenedione, was extracted from the binding cavity and re-docked to the aromatase enzyme. According to the root mean square deviation (RMSD) of 1.47 Å, the orientation of the re-docked substrate was nearly identical to the X-ray crystallographic conformer, which confirmed the validity of the docking. Molecular docking simulations of proteins where ligand binding involves prosthetic groups, like NAD or HEM (in our case the iron atom of the heme cofactor), pose a great challenge. In this study, a set of charges obtain from the work of Favia et al., calculated using the DFT method considering the B3LYP functional and 6-31G(d) basis set, was applied for the heme cofactor and used in the molecular docking.⁴⁵

2.2. QM/MM methodology

For the QM/MM calculations, we employed the Gaussian 09 quantum chemistry package.⁴⁶ “Our own N-layered integrated molecular orbital and molecular mechanic” (abbreviated as ONIOM) method implemented in Gaussian 09 was used for the QM/MM calculations.⁴⁷ 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 is considered as a 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 the ligand and protein are exclusively non-covalent; therefore, the ligand was assumed as the QM region and the protein as the MM region. To obtain the partial atomic charge of the ligand atoms in the active site of aromatase, the 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 used for the MM region (aromatase). Therefore, a two layer ONIOM calculation (PM6:UFF) was used for the calculations. The partial atomic charges of the MM region were assigned using the QEq formalism.⁴⁷ The partial atomic charges of the ligands obtained through QM/MM calculations were used to increase the accuracy of the docking calculations.

To further increase the accuracy of the docking results, three layer ONIOM (QM:QM:MM) calculations were performed on the structures with the best interaction energies obtained from the previous docking study. For this purpose, the geometry of the ligand along with the porphyrin structure of the heme-iron and several residues, including Met 374, Val 373, Val 370, Ile 305, Ala 306, Ile 133, Trp 224, Leu 372, Leu 477, Phe 134 and Thr 310, were optimized in the field of the protein. The density functional theory (DFT) method, employing the B3LYP/LANL2MB/6-31G(d) basis set, was used for the high-level part of system (ligand). The PM6 semiempirical method was used for the medium-level part of system (heme-iron and the residues) and UFF was employed for the rest of protein.

3. Results and discussion

As mentioned above, the benzofuran–imidazole and quinazolinone analogs have shown good cytotoxic effects on the MCF-7 cell line. In this study, by means of molecular hybridization method, novel hybrid compounds bearing imidazole, benzofuran and quinazolinones were designed as potential cytotoxic agents (Scheme 1). It was expected that, if the spatial orientation of these hybrid compounds was symmetrical, similar to azole-type aromatase inhibitors,⁴¹ they could interact with the aromatase enzyme. The binding modes of these novel compounds to aromatase were investigated through the docking method.

The prediction of the binding affinity in a molecular docking tool is estimated using a scoring function, which generally needs to be both fast and accurate. The electronic interaction is one of the important components of the energy model. So, assuming fixed dielectric charge for the protein and ligand atoms, when considered in the docking procedure, leads to low accuracy of the docking results. This problem is more important for proteins with a metal ion in their active site. In fact, the presence of a metal ion induces a higher polarization effect and enhances restriction of docking in the prediction of electronic interactions. In this study, owing to the positive charge on the imidazole ring of the designed compounds and also the presence of the heme iron in the active site, polarization effects are more important in the energy calculations. So, to enhance the accuracy of the results, improvement of the charge model should be considered. QM/MM methods can help to offer a superior estimate of the electronic interactions. Previous studies have shown that a docking program gives better results if the ligand partial charges are refitted with QM/MM.¹¹ In the present study, we performed the QM/MM calculation on the reference structures to obtain atomic charges on the ligand atoms. These calculated charges are presumed to represent with reasonable accuracy the ligand atom charges that are polarized by the surrounding atoms of the receptor molecule within the binding sites. Then, the ligands with improved charges were employed for re-docking into the aromatase enzyme. Since the top-scoring pose in docking is dependent on the charges of the ligand atoms, these two steps were repeated until the change in the charge values became insignificant. AutoDock4 was used to dock the new hybrid compounds into the aromatase enzyme and record the top scoring structure. Next, QM/MM calculation was performed in which the ligand was treated with PM6 calculation as the QM level, and new atomic charges on the ligand atoms were obtained via Mulliken population analysis. Once the charge values in the ligand files were substituted with these new charge values, another AutoDock4 run was performed and the top scoring structure was recorded again. After repeating these operations three times, the charges were almost constant and were used for the final docking process. Subsequently, with the selection of the best complex between the ligand and protein according to its cluster and binding energy, important interactions were investigated. The free energies of binding (ΔG_b) and inhibition constants (K_i) as calculated using AutoDock are summarized in Table 1.

Table 1 Chemical structures and the changes in free binding energy (kcal mol⁻¹) during determination of more accurate charges according to the protein active site, obtained from single point QM/MM calculations, as well as inhibition constants (K_i) of ligands after refitting the charges calculated using AutoDock

No.	X	Y	Z	E	Cluster	Binding energy^a	Binding energy^b	Binding energy^c	Ki
a The first binding energy calculated with AutoDock.b The binding energy calculated after refitting charge with the values obtained from QM/MM calculation.c The binding energy calculated after fixed change values.
1	H	H	Methyl	H	61	−7.77	−7.83	−7.81	2.00 μM
2	Cl	H	Methyl	H	60	−7.83	−7.97	−7.95	1.50 μM
3	Cl	Cl	Methyl	H	49	−8.21	−8.35	−8.34	632.00 nM
4	OH	H	Methyl	H	52	−7.77	−8.80	−8.82	409.00 nM
5	OH	Methyl	Methyl	H	50	−8.08	−8.32	−8.35	610.00 nM
6	OH	OH	Methyl	H	73	−9.14	−9.32	−9.34	68.00 nM
7	Methoxy	H	Methyl	H	46	−8.14	−8.30	−8.33	640.00 nM
8	Methoxy	Br	Methyl	H	40	−7.61	−7.74	−7.73	2.29 μM
9	Methoxy	H	Propyl	H	32	−6.27	−6.41	−6.43	17.50 μM
10	Br	H	Methyl	H	34	−7.29	−7.39	−7.41	4.25 μM
11	OH	H	Methyl	H	74	−8.18	−8.38	−8.37	1.07 μM
12	OH	Cl	Propyl	H	20	−7.77	−8.01	−8.07	1.20 μM
13	Br	OH	Propyl	H	12	−6.84	−6.54	−6.57	15.00 μM
14	Methoxy	OH	Propyl	H	16	−6.48	−6.37	−6.34	22.00 μM
15	OH	Methoxy	Methyl	H	38	−7.92	−8.36	−8.33	639.00 μM
16	OH	Cl	Methyl	H	54	−8.71	−8.80	−8.86	411.00 nM
17	OH	Br	Methyl	H	48	−8.76	−8.86	−8.85	423.00 nM
18	OH	Methoxy	Propyl	H	12	−6.80	−6.63	−6.66	17.03 μM
19	OH	OH	Propyl	H	36	−7.12	−7.28	−7.21	4.51 μM
20	OH	H	Propyl	H	24	−7.54	−7.09	−7.12	6.04 μM
21	OH	Methyl	Propyl	H	8	−6.19	−6.54	−6.56	15.35 μM
22	OH	Cl	Methyl	Thiophene	20	−6.13	−5.80	−5.83	48.81 μM
23	OH	Br	Methyl	Thiophene	20	−6.56	−6.08	−6.10	25.00 μM
24	OH	Methoxy	Methyl	Thiophene	36	−6.27	−5.98	−6.01	39.00 μM
25	OH	OH	Methyl	Thiophene	42	−6.34	−7.01	−7.08	5.72 μM
26	OH	H	Methyl	Thiophene	66	−5.65	−5.40	−5.44	106.00 μM
27	OH	Methyl	Methyl	Thiophene	32	−5.72	−5.32	−5.35	110.00 μM
28	OH	Cl	Propyl	Thiophene	14	−5.41	−5.31	−5.28	113.00 μM
29	OH	Br	Propyl	Thiophene	10	−5.74	−5.31	−5.38	112.00 μM
30	OH	Methoxy	Propyl	Thiophene	8	−5.62	−5.46	−5.50	213.00 μM
31	Methoxy	H	Methyl	H	36	−8.09	−8.31	−8.34	825.00 nM
32	Methoxy	Methyl	Methyl	H	50	−7.8	−8.01	−8.08	1.18 μM
33	Methoxy	OH	Methyl	H	42	−7.92	−7.97	−8.02	1.19 μM
34	Methoxy	Methoxy	Methyl	H	40	−8.15	−8.55	−8.56	530.00 nM
35	Methoxy	Cl	Methyl	H	58	−7.61	−7.98	−8.03	1.19 μM
36	Methoxy	Br	Methyl	H	54	−7.82	−8.02	−8.13	1.04 μM
37	Methoxy	Methoxy	Propyl	H	30	−6.68	−6.34	−6.33	22.00 μM
38	Methoxy	OH	Propyl	H	79	−6.74	−6.8	−6.79	11.40 μM
39	Methoxy	H	Propyl	H	87	−6.17	−6.30	−6.34	9.26 μM
40	Methoxy	Methyl	Propyl	H	10	−6.03	−5.83	−5.86	31.86 μM
41	Methoxy	Cl	Methyl	Thiophene	41	−6.01	−5.97	−5.96	55.74 nM
42	Methoxy	Br	Methyl	Thiophene	20	−5.78	−5.60	−5.59	70.00 μM
43	Methoxy	Methoxy	Methyl	Thiophene	50	−6.04	−5.93	−5.91	44.60 μM
44	Methoxy	OH	Methyl	Thiophene	36	−6.19	−6.09	−6.06	36.33 μM
45	Methoxy	H	Methyl	Thiophene	50	−6.01	−5.97	−5.98	44.45 μM
46	Cl	H	Methyl	H	30	−7.73	−7.88	−7.85	1.61 μM
47	Cl	Methyl	Methyl	H	60	−8.15	−8.34	−8.36	825.00 nM
48	Cl	OH	Methyl	H	54	−8.30	−8.44	−8.45	642.00 nM
49	Cl	Methoxy	Methyl	H	54	−8.00	−8.65	−8.59	531.00 nM
50	Cl	Cl	Methyl	H	62	−7.78	−8.02	−8.04	1.21 μM
51	Cl	Br	Methyl	H	58	−7.64	−7.81	−7.79	1.69 μM
52	Cl	Methoxy	Propyl	H	14	−5.59	−6.01	−5.93	44.58 μM
53	Cl	OH	Propyl	H	12	−6.02	−6.21	−6.17	30.89 μM
54	Cl	H	Propyl	H	22	−6.60	−6.13	−6.20	26.68 μM
55	Cl	Methyl	Propyl	H	18	−6.16	−5.99	−5.95	5.68 μM
56	Methyl	H	Methyl	H	30	−6.89	−6.97	−6.95	8.78 μM
57	Methyl	Methyl	Methyl	H	40	−7.02	−7.62	−7.65	14.74 μM
58	Methyl	OH	Methyl	H	16	−7.95	−8.45	−8.42	644.00 nM
59	Methyl	Methoxy	Methyl	H	34	−7.35	−7.85	−7.81	13.69 μM
60	Methyl	Cl	Methyl	H	30	−7.97	−8.35	−8.38	820.00 μM
61	Methyl	Methoxy	Propyl	H	36	−5.32	−5.65	−5.62	67.00 μM
62	Methyl	OH	Propyl	H	24	−6.24	−6.01	−6.05	38.12 μM
63	Methyl	H	Propyl	H	52	−6.18	−6.45	−6.44	18.74 μM
64	Methyl	Methyl	Propyl	H	58	−6.57	−6.76	−6.71	11.04 μM

---

In addition to the physicochemical properties of the inhibitor, molecular shape is also expected to be extremely important for accessing and fitting within the active site of the aromatase. Docking analyses revealed that the novel hybrid compounds obtained from this investigation could fit well within the binding site cavity (Fig. 1) to form a three-branched shape similar to most third generation aromatase inhibitors.⁴¹ As shown previously, a high degree of hydrophobicity and the potential to establish hydrogen bonds with the aromatase enzyme are the basic physicochemical features of known aromatase inhibitors.⁴² These properties are related to the non-polar binding pocket of the enzyme which is dominated by aliphatic 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.⁴² The newly designed compounds could well make the required hydrophobic interactions through their quinazolinone, benzofuran and imidazole moieties. Examination of the best-ranked docking reveals that among the three heteroaromatic rings located in the substrate cavity, relatively, the quinazolinone nucleus was in close proximity with the heme iron. This feature reflects the binding mechanism found for this type of molecule, which explains binding through heterocyclic aromatic coordination to the heme iron of the P450 active site. However, no H-bond was predicted for N₃ of quinazoline from the docking results. Additionally, π–π conjugate interactions are formed among Phe 221, Trp 224, Ile 133 and the phenyl ring of quinazolinone. The benzofuran ring of the hybrid compounds is found to bind through hydrophobic interactions with Phe 134, Val 370, Leu 372 and Met 374. In some cases the oxygen atom of the benzofuran ring formed a hydrogen bond with Asp 309, which is in close proximity to the benzofuran. The imidazole ring was oriented to make hydrophobic interactions with the Trp 224, Val 369 and Thr 310 residues. The hydrophobic interaction observed for the cationic imidazole ring is interesting. The reason that can be stated is that this cation isn’t a simple point charge that can only be involved in electrostatic interactions but is a delocalized charge on the imidazolium ring. This delocalization leads to distribution of a positive charge on the five atoms of the ring such that carbon and hydrogen are still able to participate in hydrophobic interactions, as seen in the figure drawn using Ligplot. It should be mentioned that with this scaffold, due to the presence of a positive charge on the nitrogen atom of imidazole, it seems logical for it to be accommodated away from the heme iron. Also, because of the resonance effect, none of the nitrogens of the imidazole could make hydrogen bonds with the residues in the aromatase active site.

Fig. 1 The binding mode of the new hybrid scaffold in the active site of aromatase, obtained from AutoDock4: (a and b) 3D structure, (c) 2D structure.

In the next step, butyl, halogen, hydroxyl and methoxy groups were appropriately substituted on the heteroaromatic rings of the designed compounds and their binding modes were investigated through docking. The new analogs share the same binding mode, similar to the unsubstituted derivatives, with an extra anchoring point which might cause stronger biding to the aromatase active site. The main binding modes in these complexes can be described as follows: introduction of OH on the benzofuran ring caused formation of hydrogen bonds with the amino acid residue Ser 478 of aromatase (Fig. 2), while an OH group on quinazolinone ring exhibited a H-bond with Met 374 (Fig. 3). In the structure with two hydroxyl groups substituted on both quinazolinone and benzofuran, both of the mentioned H bonds are visible (Fig. 4).

Fig. 2 Binding mode and hydrogen bond interaction of hydroxyl group on the benzofuran ring with Ser 478 in the aromatase active site: (a) 3D structure, (b) 2D structure.

Fig. 3 Binding mode and hydrogen bond interaction of hydroxyl group on the quinazolinone ring with Met 374 in the aromatase active site: (a) 3D structure, (b) 2D structure.

Fig. 4 Binding modes and hydrogen bond interactions of hydroxyl groups on both quinazolinone and benzofuran rings with Ser 478 and Met 374 in the aromatase active site, respectively: (a) 3D structure, (b) 2D structure.

Insertion of methoxy groups on the quinazolinone and/or benzofuran systems also caused the formation of H-bonds with hydrogen-bond donors in the active site of the enzyme. The methoxy group on the benzofuran ring and Ser 478 are in close proximity (Fig. 5) with favorable hydrogen bonding interactions, while Met 374 formed a hydrogen bond with the methoxy group on the quinazolinone ring (Fig. 6). Additionally, hydrophobic substituents, such as chloro and methyl groups, on the quinazolinone and benzofuran rings somewhat increase aromatase inhibition potency through increased hydrophobic interactions. Chloro and methyl substituents on the quinazolinone ring are stabilized by van der Waals interactions with the non-polar amino acids in the active site (Ala 306, Thr 310, Trp 224, Val 370, Ile 133, Phe 134, Leu 372 and Val 373), while chloro and methyl on the benzofuran ring are involved in hydrophobic interactions with the hydrophobic residues Ile 133, Phe 134, Leu 372 and Val 369. Finally, substitution of a small alkyl group such as methyl on the imidazolium ring make its access to and accommodation within the active site easier compared to the bulkier butyl group.

Fig. 5 Binding mode and hydrogen bond interaction of a methoxy group on the benzofuran ring with Ser 478 in the aromatase active site: (a) 3D structure, (b) 2D structure.

Fig. 6 Binding mode and hydrogen bond interaction of a methoxy group on the quinazolinone ring with Met 374 in the aromatase active site: (a) 3D structure, (b) 2D structure.

The structures with the best Gibbs interaction energies obtained through the previous docking experiments were re-optimized using the three-layer ONIOM method. In other words, the structure of the ligand, along with the interacting residues and heme-iron, was optimized in the electrostatic field of the rest of the protein, which is frozen. The advantage of this method is that the ligand is optimized in the presence of the protein and its structure is more reliable than before. In addition, the geometrical structures of the residues interacting with the ligand are improved by this optimization. For the high-level layer we used the designed inhibitors inside the aromatase active site, for the medium-level we used the heme-iron and several residues from the binding pocket, including as Met 374, Val 373, Val 370, Ile 305, Ala 306, Ile 133, Trp 224, Leu 372, Leu 477, Phe 134 and Thr 310 and the remainder of the protein is treated as the low-level layer, as depicted in Fig. 7. After the geometric parameters of the molecules are fully optimized, the ligands were extracted from the complexes and were subjected to rigid re-docking. This means that all rotatable bonds were to be held constant and the ligand charges were replaced with obtained values from the QM/MM calculation. The results showed that all runs extended to creation of one cluster and also the calculated ΔG_s increases by as much as 1–2 kcal mol⁻¹ for the ligands (Table 2). In the early QM/MM calculation we only used the obtained charges by this method for improving the docking results, but in the final QM/MM calculation the optimized ligand in the protein environment was used for AutoDock. After that the results showed that if the structure of ligand was optimized in the active sit of the enzyme and then the same optimized ligand was re-docked to protein; the lower binding energy was obtained.

Fig. 7 Inhibitor–enzyme model used for ONIOM calculations: the ball and stick representation is used for the high (large spheres) and medium (small spheres) layer atoms and the lines are used for the low-layer atoms.

Table 2 Comparison of free energies of binding (ΔG_b) in kcal mol⁻¹ and inhibition constants (K_i) calculated using AutoDock before and after optimization via the three-layer ONIOM method

No	ΔGb docking after refitting of ligand charge	Ki after refitting of ligand charge	RMSD	ΔGb in rigid docking	Ki in rigid docking	RMSD
4	−8.82	409 nM	0.31	−10.11	37.00 nM	0.14
6	−9.34	68 nM	0.27	−11.47	4.39 nM	0.16
16	−8.86	411 nM	0.25	−9.94	65.00 nM	0.14
17	−8.85	423 nM	0.36	−9.75	66.00 nM	0.11
34	−8.56	530 nM	0.41	−9.95	65.00 nM	0.21
48	−8.45	642 nM	0.28	−10.81	9.69 nM	0.11
49	−8.59	531 nM	0.25	−10.41	12.44 nM	0.16
58	−8.42	644 nM	0.38	−10.23	17.00 nM	0.12

---

According to the above results, the designed ligands formed three-branched structures in the protein environment, which were accommodated well into the active site. We found that in addition to the van der Waals interactions, hydrogen bonds are key factors for the ligand–receptor interactions. Compound 6, which yielded the highest ΔG_b (−11.47) and the best performing K_i (48 nM) value, was assumed to be the best ligand. Since the assay method used for evaluating the anti-aromatase potencies of the compounds was theoretical, the precise assessment of aromatase inhibition as the mechanism of action for these compounds needs further studies. Based on the results of the present research, synthesis and cytotoxic evaluation of the designed derivatives are under way in our research group.

4. Conclusions

In summary, in this study we developed novel potential aromatase inhibitors that contain substituted benzofuran, quinazolinone and imidazole rings. Docking simulations were carried out to identify the interactions of these ligands with the aromatase enzyme. Since polarization effects are involved in the interaction energy, QM/MM calculations were performed in which atomic charges on the ligand atoms were calculated with reasonable accuracy according to the polarization effect of the surrounding atoms of the receptor molecule within the binding site. Using the new set of charge values, the docking was performed again. These two steps were repeated three times until the change in charge reached a negligible value. At the end, the lowest energy structure (taking into account the charge polarization) was selected, which showed that the designed hybrid compounds were adopted properly within the aromatase binding site. On the other hand, re-docking of ligands, which were optimized via the three-layer ONIOM method, resulted in better ΔG_s, confirming the ability of these designed compounds as potential inhibitors of aromatase. According to the above results, compound 6 with K_i = 48 nM and ΔG_b = −11.47 kcal mol⁻¹ has the highest interaction with aromatase.

Conflict of interest

The authors declare no competing financial interest.

Acknowledgements

We gratefully acknowledge the financial support from the Iran National Science Foundation (INSF) with project number of 93012407 and the Research Council of Isfahan University of Medical Sciences.

References

1. P. P. Kore, M. M. Mutha, R. V. Antre, R. J. Oswal and S. S. Kshirsagar, Open J. Med. Chem., 2012, 2, 139–148.
2. G. Sliwoski, S. Kothiwale, J. Meiler and E. W. Lowe, Pharmacol. Rev., 2014, 66, 334–395.
3. W. L. Jorgensen, Science, 2004, 303, 1813–1818.
4. R. D. Taylor, P. J. Jewsbury and J. W. Essex, J. Comput.-Aided Mol. Des., 2002, 16, 151–166.
5. S. F. Sousa, P. A. Fernandes and M. J. Ramos, Proteins: Struct., Funct., Bioinf., 2006, 65, 15–26.
6. P. P. Roy and K. Roy, J. Pharm. Pharmacol., 2010, 62, 1717–1728.
7. R. T. Kroemer, Curr. Protein Pept. Sci., 2007, 8, 312–328.
8. D. S. Goodsell, G. M. Morris and A. J. Olson, J. Mol. Recognit., 1996, 9, 1–5.
9. M. Xu and M. A. Lill, Drug Discovery Today: Technol., 2013, 10, 411–418.
10. M. Totrov and R. Abagyan, Curr. Opin. Struct. Biol., 2008, 18, 178–184.
11. A. E. Cho, V. Guallar, B. Berne and R. Friesner, J. Comput. Chem., 2005, 26, 915–931.
12. K. Raha, M. T. B. Peters, B. Wang, N. Yu, A. M. Wollacott, L. M. Westerhoff and K. M. Merz, Drug Discovery Today, 2007, 12, 725–731.
13. H. Liu, M. Elstner, E. Kaxiras, T. Frauenheim, J. Hermans and W. Yang, Proteins: Struct., Funct., Genet., 2001, 44, 484–489.
14. A. Warshel, Angew. Chem., 2014, 53, 10020–10031.
15. S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 211–228.
16. S. C. L. Kamerlin and A. Warshel, Proteins: Struct., Funct., Bioinf., 2010, 78, 1339–1375.
17. A. Warshel, Computer Modeling of Chemical Reactions in Enzymes and Solutions, Wiley, 1997.
18. M. Repic, R. Vianello, M. Purg, F. Duarte, P. Bauer, S. C. L. Kamerlin and J. Mavri, Proteins: Struct., Funct., Bioinf., 2014, 82, 3347–3355.
19. C. Lazar, A. Kluczyk, T. Kiyota and Y. Konishi, J. Med. Chem., 2004, 47, 6973–6982.
20. R. A. Rane and V. N. Telvekar, Bioorg. Med. Chem. Lett., 2010, 20, 5681–5685.
21. K. V. Sashidhara, A. Kumar, M. Kumar, J. Sarkar and S. Sinha, Bioorg. Med. Chem. Lett., 2010, 20, 7205–7211.
22. S. Rizzo, C. Riviere, L. Piazzi, R. Stefano, B. Alessandra, S. Gobbi, M. Bartolini, V. Andrisano, F. Morroni, A. Tarozzi, J. P. Monti and A. Rampa, J. Med. Chem., 2008, 51, 2883–2886.
23. F. W. Muregi and A. Ishi, Drug Dev. Res., 2010, 71, 20–32.
24. P. Molina, A. Tarraga and F. Oton, Org. Biomol. Chem., 2012, 10, 1711–1724.
25. X. D. Yang, X. H. Zeng, Y. L. Zhang, C. Qing, W. J. Song, L. Li and H. B. Zhang, Bioorg. Med. Chem. Lett., 2009, 19, 1892–1895.
26. X. D. Yang, W. C. Wan, X. Y. Deng, Y. Li, L. J. Yang, L. Li and H. B. Zhang, Bioorg. Med. Chem. Lett., 2012, 22, 2726–2729.
27. W. Chen, X. Y. Deng, Y. Li, L. J. Yang, W. C. Wan, X. Q. Wang, H. B. Zhang and X. D. Yang, Bioorg. Med. Chem. Lett., 2013, 23, 4297–4302.
28. D. Wang and F. Gao, Chem. Cent. J., 2013, 7, 95–110.
29. A. S. ElAzab, M. AlOmar and A. Abdel, Eur. J. Med. Chem., 2010, 45, 4188–4198.
30. P. S. Pallan, C. Wang, L. Lei, F. K. Yoshimoto, R. J. Auchus, M. R. Waterman, F. P. Guengerich and M. Egli, J. Biol. Chem., 2015, 290, 13128–13143.
31. P. S. Pallan, L. D. Nagy, L. Lei, E. Gonzalez, V. M. Kramlinger, C. M. Azumaya, Z. Wawrzak, M. R. Waterman, F. P. Guengerich and M. Egli, Biol. Chem., 2015, 290, 3248–3268.
32. S. Goyal, Y. Xiao, N. A. Porter, L. Xu and F. P. Guengerich, J. Lipid Res., 2014, 55, 1933–1943.
33. G. Chowdhury, N. Shibata, H. Yamazaki and F. Peter Guengerich, Chem. Res. Toxicol., 2014, 27, 147–156.
34. J. C. Hackett, R. W. Brueggemeier and C. M. Hadad, J. Am. Chem. Soc., 2005, 127, 5224–5237.
35. R. W. Brueggemeier, J. C. Hackett and E. S. Diaz-Cruz, Endocr. Rev., 2005, 26, 331–345.
36. G. A. Khodarahmi, M. R. Khajouei, G. H. Hakimelahi, D. Abedi, E. Jafari and F. Hassanzadeh, Res. Pharm. Sci., 2012, 7, 151–158.
37. F. Hassanzadeh, E. Jafari, G. H. Hakimelahi, M. R. Khajouei, M. Jalali and G. A. Khodarahmi, Res. Pharm. Sci., 2012, 7, 87–94.
38. R. Whomsley, I. E. Fernandez, P. J. Nicholls, H. J. Smith, P. Lombardi and V. Pestellin, J. Steroid Biochem. Mol. Biol., 1993, 44, 675–676.
39. G. A. Khodarahmi, C. A. Laughton, H. J. Smith and P. J. Nicholls, J. Enzyme Inhib., 2001, 16, 401–416.
40. M. R. Saberi, T. K. Vinh, S. W. Yee, B. J. N. Griffiths, P. J. Evans and C. Simons, J. Med. Chem., 2006, 49, 1016–1022.
41. M. A. C. Neves, T. C. P. Dinis, G. Colombo and M. L. S. Melo, J. Med. Chem., 2009, 52, 143–150.
42. N. Suvannang, C. Nantasenamat, C. Isarankura-Na-Ayudhya and V. Prachayasittikul, Molecules, 2011, 16, 3597–3617.
43. J. Sgrignani and A. Magistrat, J. Chem. Inf. Model., 2012, 52, 1595–1606.
44. R. Huey and G. M. Morris, Using AutoDock4 with AutoDockTools: A Tutorial, The Scripps Research Institute Molecular Graphics Laboratory, California, 2006.
45. A. D. Favia, A. Cavalli, M. Masetti, A. Carotti and M. Recanatini, Proteins: Struct., Funct., Bioinf., 2006, 62, 1074–1087.
46. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman and G. Scalmani, Gaussian Development Version, Revision B.01, Gaussian. Inc., Wallingford CT, 2009.
47. S. Dapprich, I. Komaromi, K. S. Byun, K. Morokuma and M. J. Frisch, J. Mol. Struct.: THEOCHEM, 1999, 462, 1–21.

---