Novel aromatase inhibitors selection using induced fit docking and extra precision methods: Potential clinical use in ER-alpha-positive breast cancer

Aromatase (CYP19A1) the key enzyme of estrogen biosynthesis, is often deregulated in breast cancer patients. It catalyzes the conversion of androgen to estrogen, thus responsible for production of estrogen in human body. However, it causes over-production of estrogen which eventually leads to proliferation of breast cancer cells. Identification of new small molecule inhibitors targeted against CYP19A1 therefore, facilitates to increase drug sensitivity of cancer cells. In this scenario, the present study aims to identify new molecules which could block or suppress the activity of aromatase enzyme by molecular docking studies using Schrödinger-Maestro v9.3. In this study we used in silico approach by modeling CYP19A1 protein the strcture was subjected to protein preparation wizard; to add hydrogen and optimize the protonation states of Thr310 and Ser478 and Asp309 residues. Active site of the CYP19A1 protein was identified using SiteMap tool of Scchrodinger package. We further carried out docking studies by means of Glid, with various ligands. Based on glid score, potential ligands were screeened and their interaction with CYP19A1 was identified. The best hits were further screened for Lipinski’s rule for drug-likeliness and bioactivity scoring properties. Thus, we report two rubivivaxin and rhodethrin compounds that have successfully satisfied all in silico parameters, necessitating further in vitro and in vivo studies.


Background:
Breast cancer is the most prevailing malignancy among females, and the second most common cause of cancer-related deaths in women in the United States [1]. An estimated 235,030 new cases of invasive breast cancer are diagnosed among both males and females in the US during 2014 withan estimated death of 40,430 people [2]. Breast cancer is a complex and heterogeneous diseasewith varying clinical outcomes, disease progression, and responses to specific treatmentsattributed by a wide array of elements ranging from tumorintrinsic genetic factorstoextrinsic tumor micro-environmental factors [3]. Gene expression studies using DNA microarrays have identified several distinct breast cancer subtypes based on an intrinsic gene list that includes 496 genes that differentiate breast cancers into separate groups based only on gene expression patterns [4]. Estrogens are important players in breast cancer tumorio-genesis. The estrogen-bound Estrogen Receptor (ER) complex regulates the transcriptome of breast cancer cells by interaction with different transcription factors.Despite the plethora of physiological and pathophysiological functions of estrogen, a large number of ©2016 studies support the role ofestrogens in growth and development of ER positive breast cancer [5]. The production of estrogens in human body is meditated by the enzyme Cytochrome P450 aromatase, commonly known as estrogen synthase [6].
Cytochrome P450 aromatase (CYP19A1; EC 1.14.14.1) catalyzes the rate limiting step of estrogen biosynthesis, by the aromatization of C-19 aliphatic androgens to C-18 aromatic estrogens [7,8]. It is encoded by CYP19A1gene [9] located on chromosome 15q21.1.The 57.9 kDA molecular weight protein,is madeof a single polypeptide chain of 503 amino-acid residues and a prosthetic heme group at its active site [10]. It is widely expressed in many tissues including the gonad of both sexes [11]. In the male gonads, it is generally expressed in the testis and accessory glands wherein it maintains high levels of estradiol needed for normal spermiogenesis, maturation, and motility of sperm [12]. In females, it is mostly expressed in the ovaries of premenopausal women, and in the placenta of pregnant women. Furthermore, it is also expressed in the peripheral adipose tissues of men and postmenopausal women.
Since ER plays a critical role in breast cancer [13,14], several therapies have been developed over the past few decades among which endocrine therapy is the foremost treatment for ERPositive breast cancer patients [15]. The therapy includes (AIs) aromatase inhibitors, (SERMs) selective estrogen receptor modulators and (SERDs) selective estrogen receptor down regulators. These agents function by either reducing circulating estrogen levels or competing with estrogen for binding to its receptor [16]. AIs bind and block the active site of aromatase thus preventing the conversion of androgens to estrogens. Hence they are widely utilized for inhibition of growth and proliferation of breast cancer cells [17] and serve as first-line therapy for metastatic breast cancer.

Methodology:
Despite remarkable success of endocrine therapy, resistance toward these therapies has been detected [18,19]. In particular, long-time exposure of targeted agents leads to development of resistant malignant cells [20]. Moreover aromatase inhibitors increase the risk of cardiovascular diseases and musculoskeletal symptoms like arthralgia, osteoporosis and decreased bone mineral density (BMD) in post-menopausal women [21, 22 & 23]. As a result, there is a need for improved and effective therapies which could increase therate of patient survival and in turn improve the quality of life. An important solution towards solving this puzzle would be discoveryof new inhibitors and assessing the potential of currently utilized drugs for treating diversekinds of cancer. Consequently, the present study was conducted to search for novel biomoleculesin-silico which could suppress or inhibit the activity of aromatase enzyme. In this study, we have evaluated various biomolecules by performing molecular docking using Maestro (Schrödinger) [24] focused on novel therapeutic targets may throw light into more moleculear targeted therapies that could specifically and potentially bring down for aromatase enzyme inhibitors.
A comparative molecular docking approach was followed to trace out the most potent inhibitors of aromatase enzyme. We have performedthe molecular screening to filter out and identify most relevant aromatase inhibitors. The identified and screened inhibitor molecules were then further processed for docking studies. Rhodethrin [32] and rubrivivaxin [33] are produced by purple non sulfar photosynthetic bacterium Rubrivivax benzoatilyticus JA2has biological significance in cancer cell lines. Hence, we chose rhodethrin and rubrivivaxin for the molecular docking studies and comparatively analysed the docking results with those of the most prevalent aromatase inhibitors such as troglitazone, toremifene, testolactoneand danazol. Most of the molecular docking analysis was performed using Maestro (Schrödinger).

Selection of docking molecules
A total set of 45 drug molecules were identified by Transfacand KEGG databases which target and metabolize aromataseenzyme. Out of 45 identified drugs, 42 were downloaded from drug bank [25] and KEGG Drug [26] databases for carrying out docking studies. The remaining 3 drugs (p-Bromophenol, rhodethrin and rubrivivaxin) were sketched in ChemDrawUltra 13.0 (http://www.cambridgesoft.com/EnsembleforChemistry/ ChemBio3D/ChemBio3DUltra13.0Suite/Default.aspx) and saved in MDL molfiles. Subsequently, all the ligands were prepared using LigPrep (Schrödinger) by modifying the torsions of the ligands and assigning them appropriate protonation states (Glide-Schrödinger), a single sterochemical structure was generated per ligand with possible states at target pH 7.0 ± 2.0 using Ionizer and Epik by adding metal binding states, tautomerized, desalted and optimized by producing lowenergy 3D conformation for the ligand under the Optimized Potentials for Liquid Simulations (OPLS-2005) force field while retaining the specified chiralities of the input ligand.

Preparation of target protein
The molecular structure of human aromatase enzyme (PDB ID: 3S79) was retrieved from Protein Data Bank (PDB) with a resolution of 2.75 Å. Before docking the ligands into the protein's active site, the protein was prepared using Schrödinger's molecular docking software, Maestro9.3. The protein was imported to Maestro, in order to investigate the heme-prosthetic group as binding pocket of the enzyme adopted by a ligand in the binding cleft for depiction of both ligand and recepetor was a key importance of all hetero groups using protein preparation wizard. Hydrogen atoms were added to the protein to define correct ionization and tautomeric states of amino acid residues. Finally energy minimization was performed using default constraint of 0.3 Å of Root Mean Square Deviation (RMSD) and OLPS-2005 force field.

Induced Fit Docking (IFD) and Extra Precision (XP) Protocol
IFD and XP were performed using the module IFDof Schrödinger-Maestro v9.3 [19]. Initially Glide docking was carried out for each ligand. The sample ring conformations of ligands were selected and the side chains were trimmed. The prime side chain prediction and minimization was carried out in which residues were refined within 6.0 Å of ligand poses and side chains were optimized. This led tocreation of a ligand structure and conformation that is IF to each pose of the receptor structure. Finally, Glide XP redocking was carried out as per default conditions. The ligand was rigorously docked into the induced-fit receptor structure (IFRS) and the results yielded by IFD score for each output pose.
Drug-likeness, total drug-score and toxicity risk of compounds Calculation for various drugs properties, such as mutagenic, tumorigenic, irritant nature and an adverse effect of the compounds on the reproductive system was performed using Molinspiration (http://www.molinspiration.com/cgibin/properties) and Osiris Property Explorer (http://www.organic-chemistry.org/prog/peo). The tools also give the drug-likeness and total drug-score of the compounds based on the ''Lipinski rule of five''. The overall drug score of a compound is the drug-likeness, i.e. logarithm of the partition coefficient between 1-octanol and water, for the hydrophilicity of the compound, where low log P refers to high absorption or permeation, value less than 5 (c log P) and a unit stripped logarithm (log S) for the aqueous solubility of a compound in mol l_1. The molecular weight and toxicity risk were calculated by Osiris Property Explorer (OPE) for each compound, thus having an overall potential to qualify for a drug.

Results:
The purification and crystallization of human CYP19A1 protein has greatly facilitated the identification of compounds which could inhibit its expression, thus providing therapeutic treatments for ERpositive breast cancer patients. Discovery of novel inhibitors is primarily based on computational techniques, among which IFD plays a vital role in understanding the ©2016 molecular interaction between ligand and the active site of protein. In view of this, all the ligands were docked into the active site of aromatase enzyme ( Table 1). The active site of protein is occupied by its natural substrate androstenedione (ASD). It forms H-bonds with backbone residues metal Asp309 (1.689 Å) and Met374 (1.888 Å) (Figure 1).

Figure 1: Molecular interaction of natural substrate
Androstenedione (ASD) with CYP19A1 protein.
The side chains of residues Thr310 and Ser478 make polar interactions with the bound substrate while positive charge interaction was found with Arg115. The residues Ile133, Phe134, Phe221, Trp224, Ile305, Ala306, Val370, Leu372, Val373 and Leu 477 make hydrophobic interactions with the substrate.The ranking of ligands was based on the IFD score. More negative IFD score indicates better interaction of inhibitor with the target protein (Table 1). A Comparison of IFD score of drug molecules with AIs suggested that the drug molecules namely troglitazone, toremifene, testolactone, danazol, rubivivaxin and rhodethrin showed fairly potent activity against CYP19A1 with more negative IFD Score value. Thus, these drug molecules could prove to be more potent than the available third-generation aromatase inhibitors (Table 2). To understand the in-depth interaction pattern between the ligands and CYP19A1 protein, Maestro Ligand interaction 2-D diagram was used. The moleculesrubivivaxinand rhodethrin shown the drug-likeness and toxicity risk values are less than with troglitazone, toremifene, testolactone, danazol, it indicates that, the molecules has an overall potential to qualify as lead molecules (Table 3).

Binding mode of Troglitazone with the ligand binding region of Aromatase
Docking results showed that the ligand troglitazone occupied the ligand binding region of CYP19A1 with an IFD Score of -981.60, Glide gScore of -13.23 and the Glide energy is -63.51 Kcal/mol. Hydrogen bond interactions were identified with the backbone amino acid residuesLeu372 and Met374.Twelve hydrophobic interactions with the amino acid residues Ile133, Phe134, Phe221, Trp224, Ile305, Ala306, Ala307, Val369, Val370, Val373, Cys437 andLeu477; one positive charge interaction with Arg115 and two polar interaction with the amino acid residuesThr310 and Ser478 in the ligand binding region of CYP19A1 were observed ( Figure  4).

Binding mode of Toremifene with the ligand binding region of Aromatase
Docking results showed that the ligand toremifene occupied the ligand binding region of CYP19A1 with an IFD Score of -976.97, Glide gScore of -10.97 and the Glide energy is -52.09 Kcal/mol. Trp224 was involved in the π-π stacking interaction with ligand.Sixteen hydrophobic interactions with the amino acid residues Met127, Ile133, Phe134, Phe148, Ile217, Tyr220, Phe221, Met303, Ile305, Ala306, Val369, Val370, Leu372, Val373, Met374 andLeu477;one positive charge interaction with Arg115; one negative charge interaction with Glu302 and two polarinteraction with the amino acid residues Thr310 and Ser478 in the ligand binding region of CYP19A1 were observed ( Figure 5).

Binding mode of Testolactone with the ligand binding region of Aromatase
The ligand testolactone occupied the ligand binding region of CYP19A1 with an IFD Score of -976.54, Glide gScore of -10.85 and the Glide energy is -50.76Kcal/mol. Hydrogen bond interactions were identified with the backbone amino acid residue Met374 and metal Ash309.Ten hydrophobic interactions with the amino acid residues Ile133, Phe134, Phe221, Trp224, Ile305, Ala306, Val370, Leu372, Val373 and Leu477; onepositive charge interaction with Arg115; one negative charge interaction with Glu302 and two polar interaction with the amino acid residues ©2016 Thr310 and Ser478 in the ligand binding region of CYP19A1 were observed ( Figure 6).

Binding mode of danazol with the ligand binding region of Aromatase
Docking results showed that the ligand danazol occupied the ligand binding region of CYP19A1 with an IFD Score of -975.58, Glide gScore of -10.88 and the Glide energy is -60.43 Kcal/mol. Hydrogen bond interactions were identified with the backbone amino-acid residue Met374 and metal Ash309. Phe221 was involved in the π-π stacking interaction with ligand. Ten hydrophobic interactions with the amino acid residues Ile133, Phe134, Trp224, Ile305, Ala306,Val370, Leu372, Val 373, Ile398 and Leu477; one positive charge interaction with Arg115 and two polar interactions with the amino acid residues Thr310 and Ser478 in the ligand binding region of CYP19A1 were observed (Figure 7).

Binding mode of rubrivivaxin with the ligand binding region of Aromatase
The ligand rubrivivaxinoccupied the ligand binding region of CYP19A1 with an IFD Score of -975.58, Glide gScore of -11.04 and the Glide energy is -49.96Kcal/mol. Hydrogen bond interactions were identified with the backbone amino acid residuesAla306,Leu372, Met374 and side chain amino acid residues Arg115 and Thr310. Ten hydrophobic interactions with the amino acid residues Ile133, Phe134, Phe221, Trp224, Ile305, Ala307, Val370, Val373, Cys437and Leu477, and one polar interaction with the amino acid residue Ser478 in the ligand binding region of CYP19A1 were observed (Figure 3).

Binding mode of rhodethrin with the ligand binding region of Aromatase
The ligand rhodethrin occupied the ligand binding region of CYP19A1 with an IFD Score of -972.74, Glide gScore of -9.53 and the Glide energy is -42.34 Kcal/mol.Hydrogen bond interactions were identified with the backbone amino acid residues Leu477 and Hem600 and side chain amino acid residue Thr310. Eleven hydrophobic interactions with the amino acid residues Ile133, Phe134, Phe221, Trp224, Ile305, Ala306, Val370, Leu372, Val373, Met374, Leu477; one positive charge interaction with Arg115; one negative charge interaction with Glu302 and two polar interaction with the amino acid residue Thr310 and Ser478 in the ligand binding region of CYP19A1 were observed (Figure 2).

Discussion:
Molecular docking is considered as one of the best methods to identify the protein-ligandinteractions during in silico phases of drug discovery [27]. Previous studies suggest that aromatase activity can be suppressed in the presence of troglitazone [28], toremifene [29], testolactone [30] and danazol [31]. Hence, we have selected these four different aromatase inhibitors and novel drug targets rubrivivaxinand rhodethrinfor our docking studies to comparatively investigate binding efficiency withCYP19A1. The results of docking were verified by considering some top clusters of conformations in to the best scored one. The Novel rubivivaxin and rhodethrin shown thedrug-likeness and toxicity risk values are less than with troglitazone, toremifene, testolactone, danazol, it indicates that, the novel drugs has an overall potential to qualify for a drug. The docking results have shown that troglitazone has relatively lesser binding energy followed by rubivivaxin and toremifene.Thus, our results suggest that the novel compoundrubivivaxin and rhodethrin which are shown a lesser binding energy in comparison with the other molecules may prove potential in structure based drug design to make a major impact on anticancer chemotherapy.

Conclusion:
Aromatase (CYP19A1) is associated with and rosteroidine human breast cancer. Hence, it is important to develop improved inhibitors for Aromatase enzyme. In this study, we described in silico docking studies to identify inhibitors and performed to analyze the inhibition of aromatase activity.Two biomolecues namely rubivivaxin and rhodethrin showed good binding affinity with Aromatase enzyme inhibition. Consequently, the results obtained from this study may be worthwhile, to carry out further in vitro and in vivo studies to design novel and potential inhibitors against Aromatase (CYP19A1).