Binding mode of triazole derivatives as aromatase inhibitors based on docking, protein ligand interaction fingerprinting, and molecular dynamics simulation studies.

Aromatase inhibitors (AIs) as effective candidates have been used in the treatment of hormone-dependent breast cancer. In this study, we have proposed 300 structures as potential AIs and filtered them by Lipinski's rule of five using DrugLito software. Subsequently, they were subjected to docking simulation studies to select the top 20 compounds based on their Gibbs free energy changes and also to perform more studies on the protein-ligand interaction fingerprint by AuposSOM software. In this stage, anastrozole and letrozole were used as positive control to compare their interaction fingerprint patterns with our proposed structures. Finally, based on the binding energy values, one active structure (ligand 15) was selected for molecular dynamic simulation in order to get information for the binding mode of these ligands within the enzyme cavity. The triazole of ligand 15 pointed to HEM group in aromatase active site and coordinated to Fe of HEM through its N4 atom. In addition, two π-cation interactions was also observed, one interaction between triazole and porphyrin of HEM group, and the other was 4-chloro phenyl moiety of this ligand with Arg115 residue.

INTRODUCTION

Breast cancer is mostly due to environmental and genetic factors such as diet, exercise, cigarette smoking, alcohol consumption, family history, early menarche, and late menopause. It is the most observed malignancy in the women and is also the second mortal cause of death after lung cancer(1). In the recent years, there has been an explosion of life-saving treatment advances against breast cancer to bring new hope and excitement. Instead of only one or two options, today there is an overwhelming menu of treatment choices that fight the complex mix of cells in each individual cancer. There are different types of treatment for patients suffering from breast cancer including surgery, sentinel lymph node biopsy followed by surgery, radiation therapy, and target and hormone therapy(2).

About 80% of breast cancers, once established, rely on supplies of the hormone estrogen to grow. Therefore, they are known as hormone-sensitive or estrogen-receptor-positive (ER+) cancers. Indeed, by binding to their receptors in the tumor, estrogens initiate signals that cause proliferation of immature cells in the neoplastic tissue. The main source of estrogen is the ovaries in premenopausal women, while in post-menopausal women most of the body's estrogen is produced in peripheral tissues including skin, adipose tissue and breast. Thus, one of the main strategies to treat this kind of cancer is hormone therapy by suppressing the production of estrogen in the body(345).

Aromatase, an important membrane HEM protein of the endoplasmic reticulum, is a key enzyme in the biosynthesis of estrogens. It catalyzes the last step of estrogen biosynthesis from androgens. This step includes three successive hydroxylation of the 19-methyl group of androgens, followed by simultaneous elimination of the methyl group as formic acid and aromatization of the A-ring(678),. There are two types of aromatase inhibitors (AIs) approved to treat breast cancer: first, irreversible steroidal inhibitors such as exemestane and Formestane and second non-steroidal inhibitors, such as anastrozole and letrozole, that inhibit the synthesis of estrogen via reversible competition mechanism inside aromatase binding site(39). Although, AIs are currently common and successful in the treatment of postmenopausal ER+ώ breast cancer, designing novel drugs is noteworthy and necessary for some reasons such as avoiding the risk of possible arising resistances towards available drugs, reducing toxicity, and undesirable side effects associated with a prolonged use(10).

Consequently, publishing a high resolution X-ray structure of human aromatase (PDB code 3EQM, resolution 2.90 Å) has led to a revolution towards understanding the structure and function of the human aromatase enzyme.

Furthermore, this structure made it possible to design new potent aromatase inhibitors for prevention and treatment of hormone-dependent breast cancer(11). On the basis of the proceeding consideration and in order to propose a reasonable inhibitor of aromatase, 300 compounds were designed by considering the 1,24-triazole ring of letrozole and anastrozole.

The most important parts of flavones and isoflavones based on the previous studies, were also utilized in molecular library design(12).

The designed structures were subjected to modeling studies. First of all the compounds passed into a drug likeness filter using DruLito software (V. 1.0). Subsequently, molecular docking simulations and protein ligand interaction fingerprints analysis were applied using Autodock4 (V. 4.2) and AuPosSOM (V. 2.1) softwares, respectively. Due to great role of molecular dynamic (MD) simulations in drug design(1314), the best compound was selected for further MD experiment. It should be mentioned that in some studies, QM/MM methods for biological systems were presented(1516).

Here, with the use of these results, a possible binding mode of the compounds was therefore suggested for future studies of drug design for this target(910).

MATERIALS AND METHODS

Data set

As it was mentioned, 300 compounds were designed. The structural features and docking binding energy of the studied compounds are listed in Table 1.

Table 1

Chemical structure and the docking binding energies of the triazole derivatives used in this study.

These structures were designed based on many backbones of flavonoid (flavones, isoflavones, and flavnans) containing imidazole or triazole ring(12) and the structure of compounds was initially synthesized by S. Castellano, et al(17). Here, important pharmacophore parts of these compounds were incorporated together and new structures were designed. These structures are shown in Fig. 1.

Fig. 1

Structure of (A) flavone, (B) isoflavone, and (C) 4-triazolylflavans derivatives were characterized by the presence of a common scaffold containing a central azole moiety. (D) Structure of some imidazole and triazole derivatives as aromatase inhibitors which used in the design of new entities.

Preparation of the ligands

For the ligand preparation, the 3D structures of all proposed triazole compounds were generated by marvin sketch and converted to 3D mol2 using Open Babel 2.3.2(18).

Energy minimization of the structures was operated by means of an in house TCL script using Hyperchem8 with MM+ and AM1 methods(1920). Polak-Ribiere algorithm with RMS gradient value of 0.1 was taken and the maximum number of cycles was set to 32767 in order to obtain convergence with all structures(21).

Drug likeness studies

All the optimized structures were subjected to an open source virtual screening toolbox (DruLiTo). Lipinski's rule of five was calculated in this study. The passed structures were selected for further studies(22).

Preparation of the protein structure

During protein preparation, the pdb file for the crystal structure of human placental aromatase complexing with androstenedione (3EQM.pdb, Resolution 2.9 Å and R-value 0.244) was retrieved from protein data bank (http://www.rcsb.org/pdb/home/home.do). All water molecules and co-crystal ligand were removed, missing hydrogens were added and after determining the Kollman united atom charges, non-polar hydrogens were merged into their corresponding carbons using AutoDock Tools. An in house application (MODELFACE) was used for this purpose(23). Subsequently, the enzymes were converted to PDBQT using MGLTOOLS 1.5.6.(24).

Docking procedure

The docking procedure was performed using an in house batch script (DOCKFACE)(1920) for automatic running of AutoDock 4.2 in parallel mode and using all system resources. The batch script was designed to facilitate the virtual ligand screening stepwise. The procedures include ligands and receptor preparation, grid maps generation, dpf files preparation and performing docking runs. The molecular docking was conducted with a Genetic algorithm (GA) method to find the best pose of each ligand in the active site of the target enzyme. Hundred independent GA runs were considered for each ligand under study. For Lamarckian GA; 27000 maximum generations; 2500000 maximum numbers of energy evaluations, a gene mutation rate of 0.02; and a crossover rate of 0.8 was applied. The grid maps of the protein were calculated using AutoGrid (part of the AutoDock package). The size of grid was set in a way to include not only the active site but also considerable portions of the surrounding surface. For this purpose, a grid of 65 × 65 × 65 Å in x, y, and z directions was built on the center of mass of the co-crystal ligand in the active site of aromatase with a spacing of 0.375 Å. AutoDock Tools was employed to produce both grid and docking parameter files i.e. gpf and dpf. Cluster analysis was performed on the docked results using an RMS tolerance of 2 Å.

Rigid protein-flexible ligand docking protocol were applied. For this purpose, random orientations of the conformations were generated after translating the center of the ligand to a specified position within the receptor active site, and making a series of rotamers. This process was recursively repeated until the desired number of low-energy orientations was obtained. Non-polar hydrogens of compounds were merged and then rotatable bonds were assigned. All visualization of protein ligand complexes were done using VMD software(25).

For docking validation, 20 active ligands and 70 inactive decoys were retrieved from ChEMBL database as SMILES format(26). Iterative runs of Open Babel 2.3.2 through a shell script provided the primary 3D generation of the structures as mol2 format(27). The docking of these compounds was obtained based on the applied docking procedure for our designated ligands. The two metrics of virtual screening including the area under the curve (AUC) for receiver operating characteristic (ROC) plot and the maximum value of enrichment factor (EFmax) were calculated for active ligands and decoys using our application(28).

Protein ligand interaction fingerprint

In order to carry out protein ligand interaction fingerprint (PLIF) studies on docking results, the poses of docking were extracted from dlg files using an in-house VB.Net application (preAuposSOM). The resulted pdbqt's and the receptor were converted to mol2 by means of a batch script using Open Babel 2.3.1. The resulted mol2 files were submitted to AuposSOM 2.1 web server(29). Two training phases with 1000 iterations were set in the self-organizing map settings of AuposSOM conf files. Other parameters of the software were remained as default. The output files were subjected to Dendroscope 3.2.10 for visualization of the results(30).

Molecular dynamics simulation

To perform MD simulation, the best docking pose in complex with aromatase enzyme was subjected to Gromacs 4.5.3. The input file of the protein was prepared using GROMOS96 53a6 force field and the topology file of the ligands was provided by means of PRODRG server. A cubic box was made around the system and the system was subsequently solvated with SPC water. A concentration of 0.15 M NaCl was added to the system to mimic physiological conditions of the body. Steepest descent minimization with maximum 50,000 steps was performed on the system in order to prepare it for further equilibration. Two 1 ns equilibration states were done on the system based on NVT and NPT ensembles with restrains on the ligand structures. LINCS constraint algorithm was used in both states of equilibration. The thermal groups were defined as the protein-ligand and solvent-Na-Cl groups. After equilibration, the main run of MD started using NPT ensemble and leap-frog integrator together with Berendsen thermostat was used for the whole system. The system was simulated for 50 ns on a 24 core intel server running on Linux Ubuntu. The trajectories of the simulation were thereafter subjected to VMD for further analyses and visualization of binding mode.

RESULTS

First, we designed 300 compounds by considering the 1,24-triazole ring of letrozole and anastrozole. The most important parts of flavones and isoflavones based on the previous studies(12) were also utilized in molecular library design.

All of the compounds were subjected to drug-likeness analysis. One hundred structures met the primary criteria of drug likeness. Drug-likeness rules are guidelines for the structural properties of the molecules based on the Lipinski's rule of five. The binding free energy values of those structures with suitable drug likeness properties were calculated using docking studies. The result of Gibbs free energy changes for the 20 best compounds are listed in Table 1. As seen, the compound 15 was showing the best binding energy among all docking poses of this series of compounds. It is important to know that HEM is a cofactor in aromatase enzyme and a metal binding site exist near the binding cavity. Therefore, it was thought that the binding energies of the ligands are a function of the distance between nitrogens and Fe inside HEM. To check out this hypothesis the distance of nitrogen (N4) inside triazole ring of each ligand from Fe was measured in case of all 20 best docking poses. A plot of distances and binding energies is depicted in Fig. 2A. It was seen that the obtained binding energies are co-related to the distance between N4-Fe (R = 0.85). Visualization of the poses in the active site with regard to Fe of HEM is depicted in Fig. 2B. The least distance between N4 and Fe is 2.39 which is suitable for metal coordination. To compare the protein ligand interaction fingerprint of the structures with positive controls (anastrozole and letrozole), all docking poses of the top 20 structures with less negative binding energy values were subjected to AuposSOM software (V. 2.1).

Fig. 2

(A) correlration (linear regression) of distance between Fe of HEM and N4 of triazole ring with binding free energy. The obtained binding energies are correlated to the distance between N4-Fe (R = 0.85). (B) Visualization of seven best protein ligand complexes using VMD software. The least distance between N4 and Fe is 2.39 which is suitable for metal coordination.

As shown in fig. 3, three structures (compounds 8, 9 and 15) were co-clustered with positive controls. This is showing that the interaction pattern of these ligands is similar to that observed for anastrozole and letrozole. Interestingly, all these ligands were those with less binding energy values compared to other structures. Among these three compounds, the structure with the best binding energy was selected for MD simulation. After performing 50 ns MD simulation, a plateau was obtained based on the root mean squared distance (RMSD) of Cα atoms (Fig. 4A). This was showing the system at equilibrated state. To analyze the ligand protein fluctuations during MD, heat map analysis of the system was performed. The most fluctuations occurred at 1–141 (from N-terminal residues) and 271–451 (from C-terminal residues) (Fig. 4B). The structure of the enzyme and ligand in a frame with the most fluctuation is displayed in Fig. 4C. The residues with most fluctuation are colored as red and orange and those with more stability during simulation are colored green. It shows that most residues in contact with the ligand are in stable parts of the enzyme.

Fig. 3

Representation of the protein ligand interaction fingerprint of all docking poses for the top 20 structures with less negative binding energy values using AuposSOM software. Three structures (8, 9 and 15) were co-clustered with positive controls, letrozole and anastrozole.

Fig. 4

(A) root mean squared distance (RMSD) of Cα for the enzyme residues during simulation. After performing 50 ns molecular dynamic (MD) simulation, a plateau was obtained based on the RMSD of Cα atoms.This was showing the system at equilibrated state. (B) heatmap analysis of the trajectories during simulation. The most fluctuations occurred at residues 1–141 (N-terminal) and 271–451 (C-terminal). These regions are terminal parts of the enzyme. (C) the frame with high fluctuation of terminal residues.

DISCUSSION

The binding mode of the ligand in the active site of the enzyme at equilibrated state is displayed in Fig. 5A.

Fig. 5

(A) 3D of ligand-receptor interactions for Lig 15 in aromatase active site. The triazole of ligand 15 pointed to HEM group in aromatese active site and coordinate to Fe of HEM through its N4 atom (dotted cyan line). In addition, two π-cation interaction was also observed, one interaction between triazole and porphyrin of HEM group (dotted green line), and the other one was 4-chloro phenyl moiety of this ligand with Arg115 residue (dotted green line). (B) Receiver operating characteristic (ROC) and enrichment factor (EF) diagrams for aromatase. The more AUC for ROC value means that the docking protocol is more able to discriminate between active ligands and decoys. EF diagram also validate this protocol of docking.

The triazole of ligand 15 pointed to HEM group in aromatase active site and coordinate to Fe of HEM through its N4 atom. In addition, two π-cation interactions was also observed, one interaction between triazole and porphyrin of HEM group, and the other was 4-chloro phenyl moiety of this ligand with Arg115 residue. It was observed that this ligand occupied hydrophobic pocket with the residues Leu477, Trp224, Leu372, Ala306, Ile133, Val373 and Val370. The proposed structures were sent to our synthesis division to verify the obtained results of design. The application of relative operating characteristic curve (ROC) in computational medicinal chemistry as a useful metric tool to evaluate the validity of docking protocols was first reported by Triballeau et al.(28). Nowadays, it is widely used as a validating procedure. First of all, about 100 aromatase inhibitors were retrieved from ChEMBL database as SMILES format(26). The structures based on their experimental activities categorized into two subsets of active ligands and inactive decoys. Twenty five ligands and 75 decoys were generated. Subsequently, by means of Open Babel 2.3.2 through a shell script provided the primary 3D generation of the structures as mol2 format(27). Ionization states at PH = 7 were also calculated for all structures. The shell script was provided by means of batch scripting in windows operating system. The screening method should be able to discriminate between ligands and decoys. ROC value is the AUC for the plot of the true positive rate (TPR or sensitivity) against the false positive rate (FPR or 1 – specificity) at various threshold settings. The ROC curve is thus the sensitivity as a function of 1 – specificity. The AUC for ROC is calculated by trapezoidal integration method as implemented in our in-house application. The more AUC for ROC value means that the docking protocol is more able to discriminate between active ligands and decoys. As it is depicted in Fig. 5B, the AUC of 0.828 showed that our applied docking procedure was a validated not a random protocol.

Here, EFmax as a key parameter was used as another tool to evaluate the efficiency and quality of docking protocol. Its calculations were based on the Li et al. work(31). EFmax factor in comparison to ROC curves, is highly dependent on the number of actives in a data set(28). It means that early enrichment can be easily obtained if the number of active ligands is increasing in a dataset.

To explore the interaction of the compound with receptor based on simulation data interaction of nitrogen in triazole ring and Fe of HEM was considered as the most important structure and was studied during simulation.

CONCLUSION

Aromatase inhibitors exert their activity by inhibiting the enzyme, aromatase (cyp19), which is involved in the conversion of androgens to estrogens. Obtaining a proper binding mode for the newly proposed structures is important in developing more potent and rational compounds for this target. Here, we designed some novel structures based on main parts of flavones, isoflavones, and triazole compounds as AIs. The interaction of aromatase with its inhibitors which were elucidated by docking studies indicated that the binding energies of the ligands are a function of the distance between nitrogens of triazole and Fe inside HEM. PLIF study showed that the interaction pattern of these ligands is similar to that observed for anastrozole and letrozole. The equilibrium structure of aromatase with ligand 15 was obtained after 50 ns of simulation depicting that the distance between N and Fe decreased from the beginning of the simulation to facilitate the observed interaction during simulation. Based on the experiments of this study it was seen that some structural features of the ligands are taking role in designing potent inhibitors against aromatase. These features include the ability to coordinate with Fe and accommodation features of the ligands based on π-π interaction and hydrogen bonding interactions. The knowledge acquired from this study has important implications for the development of new and more selective AIs agents.