Probing the origins of aromatase inhibitory activity of disubstituted coumarins via QSAR and molecular docking.

This study investigated the quantitative structure-activity relationship (QSAR) of imidazole derivatives of 4,7-disubstituted coumarins as inhibitors of aromatase, a potential therapeutic protein target for the treatment of breast cancer. Herein, a series of 3,7- and 4,7-disubstituted coumarin derivatives (1-34) with R1 and R2 substituents bearing aromatase inhibitory activity were modeled as a function of molecular and quantum chemical descriptors derived from low-energy conformer geometrically optimized at B3LYP/6-31G(d) level of theory. Insights on origins of aromatase inhibitory activity was afforded by the computed set of 7 descriptors comprising of F10[N-O], Inflammat-50, Psychotic-80, H-047, BELe1, B10[C-O] and MAXDP. Such significant descriptors were used for QSAR model construction and results indicated that model 4 afforded the best statistical performance. Good predictive performance were achieved as verified from the internal (comprising the training and the leave-one-out cross-validation (LOO-CV) sets) and external sets affording the following statistical parameters: R (2) Tr = 0.9576 and RMSETr = 0.0958 for the training set; Q (2) CV = 0.9239 and RMSECV = 0.1304 for the LOO-CV set as well as Q (2) Ext = 0.7268 and RMSEExt = 0.2927 for the external set. Significant descriptors showed correlation with functional substituents, particularly, R1 in governing high potency as aromatase inhibitor. Molecular docking calculations suggest that key residues interacting with the coumarins were predominantly lipophilic or non-polar while a few were polar and positively-charged. Findings illuminated herein serve as the impetus that can be used to rationally guide the design of new aromatase inhibitors.

Introduction

class="Disease">Breast cancer is commoclass="Chemical">nly fouclass="Chemical">nd iclass="Chemical">n class="Chemical">n class="Species">women and is reported as the second leading cause of women death (Desantis et al., 2011[7]). Estrogen is associated with the development of breast cancer by activating intracellular signaling cascades (Yager and Davidson, 2006[55]). Aromatase is an enzyme involved in the catalysis of androgens (i.e. androstenedione) to estrogens (i.e. estradiol). Therefore, the ability to inhibit such enzyme offers great therapeutic benefits for the treatment of breast cancer (Narashimamurthy et al., 2004[35]). Aromatase (i.e. 19A1) is a member of the cytochrome P450 family and is comprised of 503 amino acids. The proximal ligand of aromatase (i.e. Cys437) is coordinated to the heme prosthetic group and it was previously reported that the presence of the characteristic cysteine proximal ligand is crucial for the catalytic activity observed in cytochrome P450s (Auclair et al., 2001[1]; Yoshioka et al., 2001[56]). Aromatase inhibitors comprising of steroidal and non-steroidal compounds have either been synthesized by introducing or modifying functional groups on the core structure of prototype lead compounds (Ferlin et al., 2013[11]; Nativelle-Serpentini et al., 2004[37]; Neves et al., 2009[38]; Stefanachi et al., 2011[42]; Varela et al., 2012[48]) or extracted from natural sources (Balunas et al., 2008[3]). In addition, computational analysis such as quantitative structure-activity relationship (QSAR), molecular docking and modeling had previously been employed for constructing models and interpreting the interaction between compounds of interest with the aromatase enzyme (Bheemanapalli et al., 2013[4]; Galeazzi and Massaccesi 2012[13]; Ghosh et al., 2012[15]; Nantasenamat et al., 2013[34], [32]; Narayana et al., 2012[36]).

class="Chemical">Coumarin is a class="Chemical">naturally occurriclass="Chemical">ng bioactive class="Chemical">n class="Chemical">benzopyrone found in many plant species (Venugopala et al., 2013[49]). It represents an important structural scaffold in medicinal chemistry as it can afford a wide range of bioactivities such as anti-bacterial (Nagamallu and Kariyappa, 2013[28]), anti-cancer (Wu et al., 2014[54]), anti-inflammatory (Hemshekhar et al., 2013[17]), anti-oxidant (Guinez et al., 2013[16]) as well as neuroprotective (Sun et al., 2013[43]) properties. Furthermore, 7-hydroxycoumarin is the primary metabolite of coumarin in human. A number of 7-oxy substituted coumarin analogs have been shown to exhibit an array of bioactivities (Venugopala et al., 2013[49]). Recently, 3,7- and 4,7-disubstituted coumarin derivatives (1-34, Figure 1(Fig. 1)) have been reported to display aromatase inhibitory potency in the nanomolar range (Stefanachi et al., 2011[42]).

Figure 1

Chemical structures of disubstituted coumarins 1-34

The aim of this study is to explain the origins of class="Gene">aromatase iclass="Chemical">nhibitory activity via QSAR modeliclass="Chemical">ng aclass="Chemical">nd molecular dockiclass="Chemical">ng. QSAR models were built as a fuclass="Chemical">nctioclass="Chemical">n of sigclass="Chemical">nificaclass="Chemical">nt descriptors accouclass="Chemical">nticlass="Chemical">ng for class="Chemical">n class="Gene">aromatase inhibitory activity. Insights into structure-activity relationship are also discussed and it is anticipated that such information could serve as a pertinent guideline for the rational design of novel aromatase inhibitors based on the coumarin chemotype.

Materials and Methods

Data set

A set of class="Chemical">coumarin derivatives with reported class="Chemical">n class="Gene">aromatase inhibitory activity was taken from the work of Stefanachi et al. (2011[42]). Briefly, the aromatase inhibitory activity was determined in vitro using human placental microsomes as the source of aromatase while [1β-H] androstenedione was used as the substrate. IC50 values were logarithmically transformed to pIC50 as summarized by the following equation:

pIC50 = -log10(IC50) [1]

Chemical structures are depicted in Figure 1(Fig. 1) while its corresponding descriptors and bioactivity are shown in Table 1(Tab. 1). Schematic workflow of this study is displayed in Figure 2(Fig. 2).

Table 1

Molecular descriptors and aromatase inhibitory activity of coumarins (1-34)

Figure 2

Schematic workflow of QSAR and molecular docking studies performed herein

Geometrical optimization and descriptors calculation

Molecular structures were drawn using ChemAxon MarvinSketch version 6.2.1 (ChemAxon Ltd., 2014[5]) and converted to the appropriate file format with Babel version 3.3.0 (OpenEye Scientific Software, 2014[39]). An initial geometrical optimization was performed at the semi-empirical Austin Model 1 (class="Species">AM1) level followed by further reficlass="Chemical">nemeclass="Chemical">nt at the declass="Chemical">nsity fuclass="Chemical">nctioclass="Chemical">nal theory (DFT) level usiclass="Chemical">ng the Becke's three parameter hybrid method aclass="Chemical">nd the Lee-Yaclass="Chemical">ng-Parr (B3LYP) fuclass="Chemical">nctioclass="Chemical">nal togclass="Chemical">n class="Chemical">ether with the 6-31g(d) basis set using the Gaussian 09 software package (Frisch et al., 2009[12]). A set of 6 quantum chemical descriptors were obtained from the aforementioned low-energy conformer that included the following: dipole moment (µ), the total energy, the highest occupied molecular orbital energy (HOMO), the lowest unoccupied molecular orbital (LUMO), difference in energy values of HOMO and LUMO (HOMO-LUMO) and the mean absolute atomic charge (Q) (Karelson et al., 1996[19]; Thanikaivelan et al., 2000[46]). Furthermore, the low-energy conformer were also subjected to the generation of additional descriptors from the Dragon software package, version 5.5 (Talete, 2007[45]) to derive a set of 3,224 molecular descriptors comprising 22 categories: 48 Constitutional descriptors, 119 Topological descriptors, 47 Walk and path counts, 33 Connectivity indices, 47 Information indices, 96 2D autocorrelation, 107 Edge adjacency indices, 64 Burden eigenvalues, 21 Topological charge indices, 44 Eigenvalue-based indices, 41 Randic molecular profiles, 74 Geometrical descriptors, 150 RDF descriptors, 160 3D-MoRSE descriptors, 99 WHIM descriptors, 197 GETAWAY descriptors, 154 Functional group counts, 120 Atom-centered fragments, 14 Charge descriptors, 29 Molecular properties, 780 2D binary fingerprints and 780 2D frequency fingerprints.

Descriptors selection

Constant variables from the initial set of 3,224 molecular descriptors obtained from the Dragon software were subjected to removal and subsequently combined with a set of 6 quantum chemical descriptors. Significant descriptors correlated with the n class="Gene">aromatase iclass="Chemical">nhibitory activity were derived from stepwise multiple liclass="Chemical">near regressioclass="Chemical">n (Worachartcheewaclass="Chemical">n et al., 2014[53]) usiclass="Chemical">ng SPSS Statistics versioclass="Chemical">n 18.0 (SPSS Iclass="Chemical">nc., USA). Iclass="Chemical">ntercorrelatioclass="Chemical">n matrix of descriptors was coclass="Chemical">nstructed from Pearsoclass="Chemical">n's correlatioclass="Chemical">n coefficieclass="Chemical">nt values as to deduce the preseclass="Chemical">nce of variable reduclass="Chemical">ndaclass="Chemical">ncy.

Generation of internal and external sets

The data set was divided into 2 major subsets in which an external set was obtained by randomly selecting 15 % of the data samples from the full data set while the remaining 85 % served as the internal set (Nantasenamat et al., 2013[34]). The internal set was used to internally assess the predictive performance of QSAR models by using it as training and leave-one-out cross-validated (LOO-CV) sets. LOO-CV refers to the leaving out of one data sample as the testing set while using the remaining N-1 samples as the training set. The process was iteratively repeated until all samples were used as the testing set. As the name implies, the external set was used to externally assess the predictive performance of the QSAR model.

Multiple linear regression

Multiple linear regression (MLR) ws employed for constructing QSAR models as conceptually summarized by the following equation:

Y = m1x1 + m2x2 + ... + mnxn + b [2]

where Y is the class="Gene">aromatase iclass="Chemical">nhibitory activity (pIC50), m is the regressioclass="Chemical">n coefficieclass="Chemical">nt value of descriptors, x is the descriptor aclass="Chemical">nd b represeclass="Chemical">nts the y-iclass="Chemical">ntercept value. MLR calculatioclass="Chemical">ns were used as implemeclass="Chemical">nted by the Waikato Eclass="Chemical">nvclass="Chemical">n class="Chemical">ironment for Knowledge Analysis (WEKA) software, version 3.4.5 (Witten et al., 2011[50]).

Statistical analysis

Evaluation of the predictive performance of QSAR models was performed using squared correlation coefficient for the training (R2Tr) and cross-validated sets (Q2CV) as well as the root mean squared error for the training (RMSETr) and cross-validated sets (RMSECV). Furthermore, the Fisher (F) ratio as well as the difference of n class="Gene">R2 and Q2 (R2−Q2) was used to assess the predictive quality of coclass="Chemical">nstructed models. Moreover, squared correlatioclass="Chemical">n coefficieclass="Chemical">nt (Q2Ext) aclass="Chemical">nd root meaclass="Chemical">n squared error (RMSEExt) of exterclass="Chemical">nal set were used as iclass="Chemical">ndepeclass="Chemical">ndeclass="Chemical">nt validatioclass="Chemical">n of the coclass="Chemical">nstructed QSAR model (Naclass="Chemical">ntaseclass="Chemical">namat et al., 2013[34]).

Molecular docking

The set of 34 class="Chemical">coumarin compouclass="Chemical">nds was used as ligaclass="Chemical">nds for dockiclass="Chemical">ng to the X-ray crystallographic structure of class="Chemical">n class="Gene">aromatase (PDB id 3EQM). Prior to docking, both protein and ligand structures were subjected to a series of pre-processing. A snapshot of the aromatase structure was obtained from 100 ns of molecular dynamics simulation using the AMBER03 force field under YASARA 10.11.28 (Krieger et al., 2002[21], 2003[22]) as reported in our previous investigation (Suvannang et al., 2011[44]). Subsequently, atomic charges of the heme prosthetic group were parameterized according to density functional calculations reported by Favia et al. (2006[9]). Low-energy conformer of ligands as obtained from previously mentioned geometrical optimization. Kollman and Gasteiger charges were finally assigned to the protein and ligands, respectively, using AutoDockTools (Morris et al., 2009[27]).

Docking calculatons were performed using AutoDock version 4.2 (Morris et al., 2009[27]). A grid box was constructed using AutoGrid and centered at the binding cavity using x,y,z coordinates of 84.0236, 49.6568, 50.0293 where it is encapsulated by a box size of 52 × 66 × 82 Å points with spacing of 0.3759 Å. Each docking pose was obtained from 100 independent runs based on 150 randomly placed individuals in the population. The docking simulation made use of the Lamarckian Genetic Algorithm (Morris et al., 1998[26]) in searching for low-energy binding orientation. Translation step of 2.0 Å, mutation rate of 0.02, crossover rate of 0.8, local search rate of 0.06 and maximum energy evaluation of 250,000,000 were utilized in the present docking simulation. Docked conformations were clustered using an RMSD tolerance of 2.0 Å. Re-docking of the co-crystallized substrate was performed as to evaluate the validity of the docking protocol.

Post-docking analyses were carried out using AutoDockTools and graphical images were prepared using PyMOL version 0.99 (Delano, 2002[6]).

Results and Discussion

Chemical space of investigated coumarins

Stefanachi et al. (2011[42]) reported the synthesis and determination of class="Gene">aromatase iclass="Chemical">nhibitory activity for a set of class="Chemical">n class="Chemical">coumarin derivatives (1-34) as shown in Figure 1(Fig. 1). The coumarin core structure featured various substituents (i.e. halogen, methoxy, aryloxy and imidazole) that conferred its bioactivity against aromatase (Table 1(Tab. 1)). In gaining further insights on the general features of potency afforded by these coumarins, compounds having pIC50 values greater than 7 (< 0.081 nM) were considered active, less than 6 as inactive (> 2.82 μM) and values in between were considered to have intermediate activity. It can be seen that the potency of this set of compounds was relatively robust as deduced from IC50 values in the range of 47 nM to 690 nM for 31 out of the 34 compounds. The remaining three compounds afforded poor IC50 values in the range of 2.820 and 4.010μM.

Applying the Lipinski's rule-of-five on this set of compounds revealed its compliance in which there were less than 5 H-bond class="Species">donors, less thaclass="Chemical">n 10 H-boclass="Chemical">nd acceptors, molecular weight (MW) of less thaclass="Chemical">n 500 Da aclass="Chemical">nd class="Chemical">n class="Chemical">Ghose-Crippen octanol-water partition coefficient (ALogP) of less than 5 (Supplementary Information). An exploratory data analysis of the three aforementioned sub-classes revealed that intermediate and active sub-classes were generally larger than their inactive counterpart with MW of 361.246 ± 39.971, 353.146 ± 23.348 and 305.997 ± 98.512 Da, respectively. However, it should be noted that of the three compounds in the inactive sub-class, 32 had significantly higher MW of 419.46 Da than the remaining two with MW of 242.25 and 256.28 Da for 26 and 33, respectively. The number of H-bond donors does not appear to be crucial for its bioactivity as can be seen from the sparseness of their occurrence in the coumarins in which a total of two compounds had one H-bond donor belonging to the inactive (26) and intermediate (25) sub-class. As for the number of H-bond acceptors, the intermediate and active sub-classes had slightly higher mean values of 4.958 ± 1.268 and 4.857 ± 0.900, respectively, when compared to that of the inactive sub-class with corresponding value of 4.333 ± 0.577. Finally, analysis of ALogP revealed that the inactive sub-class afforded higher polarity with a mean of 2.459 ± 1.848 than both the intermediate and active sub-classes, which gave 3.634 ± 0.922 and 3.592 ± 0.537, respectively. A closer glance of the inactive sub-class shows that 32 afforded significantly higher ALogP value of 4.588 than the other two compounds (26 and 33) having values of 1.269 and 1.52.

Data pre-processing and feature selection

Quantum chemical descriptors of compounds from low energy conformer of compounds derived from geometrical optimization at B3LYP/6-31g(d). Furthermore, low energy conformer was then subjected to further calculation using the Dragon software for generating an additional set of molecular descriptors. The usefulness of quantum chemical and molecular descriptors in constructing QSAR/QSPR models of chemical properties (Nantasenamat et al., 2005[33], 2007[29],[31]) and biological activities (Khoshneviszadeh et al., 2012[20]; Martinez-Martinez et al., 2012[24]; Nantasenamat et al., 2010[30], 2013[32]; Uesawa et al., 2011[47]; Worachartcheewan et al., 2011[51], 2013[52]) had previously been demonstrated. Constant and multi-collinear descriptors were removed from the initially large set of 3,224 descriptors from Dragon resulting in a reduced subset of 1,129 molecular descriptors. Subsequently, this set of descriptors was then combined with a set of 6 quantum chemical descriptors to give rise to a total of 1,135 descriptors. Such combined set was then subjected to stepwise MLR for further round of feature selection thereby resulting in 7 significant descriptors consisting of F10[N-O], Inflammat-50, Psychotic-80, H-047, BELe1, B10[n class="Chemical">C-O] aclass="Chemical">nd MAXDP. Tables 1(Tab. 1) aclass="Chemical">nd 2(Tab. 2) display the molecular descriptors coclass="Chemical">nstituticlass="Chemical">ng the data set aclass="Chemical">nd the deficlass="Chemical">nitioclass="Chemical">n of sigclass="Chemical">nificaclass="Chemical">nt descriptors, respectively. Aclass="Chemical">n iclass="Chemical">ntercorrelatioclass="Chemical">n matrix of Pearsoclass="Chemical">n's correlatioclass="Chemical">n coefficieclass="Chemical">nt was coclass="Chemical">nstructed as to verify the iclass="Chemical">ndepeclass="Chemical">ndeclass="Chemical">nce of the coclass="Chemical">nstitueclass="Chemical">nt variables iclass="Chemical">n this set of molecular descriptors. Results iclass="Chemical">ndicated that this set of 7 descriptors were iclass="Chemical">ndepeclass="Chemical">ndeclass="Chemical">nt from oclass="Chemical">ne aclass="Chemical">nother as deduced by the low correlatioclass="Chemical">n coefficieclass="Chemical">nt values (Table 3(Tab. 3)).

Table 2

Definition of important descriptors for QSAR model

Table 3

Intercorrelation matrix of significant descriptors for QSAR model

QSAR model of aromatase inhibitory activity

The set of 7 significant descriptors were then used in the development of QSAR models using MLR to deduce a linear equation as described by Eq. [2]. Thus, such multivariate analysis generated a total of 4 models consisting of MLR equations and statistical analysis as summarized in Table 4(Tab. 4). Inherent outliers present in each of the MLR models were identified using absolute standardized residual of 2 as the cutoff.

Table 4

Summary of predictive performance of QSAR model

The initial model 1 produced the following MLR equation and statistical parameters:

pIC50 = - 0.5797(F10[N-O]) + 1.0993(Inflammat-50) + 0.4107(Psychotic-80) + 0.1123(H-047) - 13.2991(BELe1) - 0.6687 (B10[n class="Chemical">C-O]) + 0.4336(MAXDP) + 29.0288 [3]

N = 29, R = 0.9036, RMSE = 0.1412, Q = 0.7693, RMSE = 0.2238, F ratio = 10.00, critical F value = 2.488 and R - Q= 0.1343

The model showed compounds 27 and 29 as outliers, which were removed from the model and the resulting data set was re-calculated to generate Model 2 as shown below.

pIC50 = - 0.5407(F10[N-O]) + 1.0870(Inflammat-50) + 0.4025(Psychotic-80) + 0.1051(H-047) - 15.5871(BELe1) + 0.4285(MAXDP) + 32.8689 [4]

N = 27, R = 0.9347, RMSE = 0.1179, Q = 0.8821, RMSE = 0.1585, F ratio = 24.94, critical F value = 2.599 and R - Q= 0.0526

It was observed that B10[class="Chemical">C-O] was removed from the MLR equatioclass="Chemical">n owiclass="Chemical">ng to the fact that the remaiclass="Chemical">niclass="Chemical">ng 27 compouclass="Chemical">nds had coclass="Chemical">nstaclass="Chemical">nt value of 1 for the B10[class="Chemical">n class="Chemical">C-O] descriptor. Likewise, compound 20 was detected as the outlying compound and following its removal yielded model 3 as described below.

pIC50= - 0.6008(F10[N-O]) + 1.0371( Inflammat-50) + 0.3941(Psychotic-80) + 0.1042(H-047) - 15.7039(BELe1) + 0.4959( MAXDP) + 32.7809 [5]

N = 26, R = 0.9496, RMSE = 0.1082, Q = 0.9111, RMSE = 0.1423, F ratio = 32.44, critical F value = 2.628 and R - Q= 0.0385

Subsequently, compound 18 was identified as the outlying compound and its removal resulted in model 4:

pIC50 = - 0.6067(F10[N-O]) + 1.0108(Inflammat-50) + 0.4262(Psychotic-80) + 0.0884(H-047) - 14.5495(BELe1) + 0.4563(MAXDP) + 30.9432 [6]

N = 25, R = 0.9579, RMSE = 0.0958, Q = 0.9239, RMSE = 0.1304, F ratio = 36.42, critical F value = 2.661 and R - Q= 0.0340

This final model 4 was shown to afford robust predictive performance as verified from its statistical parameters for both the training and LOO-CV sets. Scatter plots of experimental versus predicted pIC50 values for the 4 models of the training and LOO-CV sets are presented in Figures 3a-d(Fig. 3) as white and black squares, respectively. Interestingly, the statistical performance increased from Eq. 3 to Eq. 6 as deduced from the increasing Q2 and F ratio. A metric proposed by Eriksson and Johansson (1996[8]) for evaluating the reliability of predictive models was obtained by calculating the difference of n class="Gene">R2 and Q2. A value iclass="Chemical">n excess of 0.3 is iclass="Chemical">ndicative of chaclass="Chemical">nce correlatioclass="Chemical">n or the preseclass="Chemical">nce of outliers. Thus, the R2−Q2 from the traiclass="Chemical">niclass="Chemical">ng aclass="Chemical">nd LOO-CV sets were fouclass="Chemical">nd to be < 0.3 thereby iclass="Chemical">ndicaticlass="Chemical">ng its statistical sigclass="Chemical">nificaclass="Chemical">nce. The experimeclass="Chemical">ntal aclass="Chemical">nd predicted pIC50 values aloclass="Chemical">ng with their respective residuals are provided iclass="Chemical">n Table 5(Tab. 5).

Figure 3

Plot of experimental versus predicted aromatase inhibitory (pIC50) activity for the training set (white squares; regression line is represented as solid line), leave-one-out cross-validation set (black squares; regression line is represented as dotted line) and external test set (blue triangles) for model 1 (a), model 2 (b), model 3 (c) and model 4 (d)

Table 5

Experimental and predicted aromatase inhibitory activities (pIC50) of coumarin analogs (1-34) obtained from model 4

Validation of QSAR model using external set

The predictivity of QSAR models 1-4 was verified by external validation using 15 % of samples from the original data set (i.e. comprising of compounds 4, 6, 12, 17 and 26). It should be noted that the external set was randomly selected and their distribution are shown in Figure 4(Fig. 4) as black squares along with the internal set shown as white squares. The internal and external predictivity of models 1-4 were determined and results suggested good predictive performance especially by model 4, which provided Q= 0.9239, RMSE = 0.1304 and Q = 0.7268. Model 4 was selected for further investigation owing to its absence of outlying compounds and its ability to perform well on both the internal and external sets. Moreover, R2−Q2 calculated from training and external sets afforded a value of < 0.3, which suggests the model's reliability.

Figure 4

Plot of molecular weight versus aromatase inhibitory activity (pIC50) used in selection of external (black squares) and internal (white squares) sets

Plots of experimental and predicted pIC50 values for the 4 models of the external set are presented in Figures 3a-d(Fig. 3) as blue triangles. As previously mentioned for the internal set, the corresponding values of the experimental and predicted pIC50 values and their respective residuals are shown in Table 5(Tab. 5).

Elucidating structure-activity relationship

class="Gene">Aromatase iclass="Chemical">nhibitory poteclass="Chemical">ncy (IC50 value) of class="Chemical">n class="Chemical">coumarin analogs was transformed to pIC50 using [Eq. 1] and was used as the dependent variable. Feature selection performed on a set of molecular and quantum chemical descriptors yielded 7 significant descriptors as follows (listed in order of relative importance as deduced from the MLR regression coefficient values): BELe1 (Burden eigenvalues) > Inflammat-50 (molecular properties) > B10[C-O] (2D binary fingerprints) > F10[N-O] (2D frequency fingerprints) > MAXDP (Topological descriptors) > Psychotic-80 (molecular properties) > H-047 (Atom-centered fragments).

Disubstituted class="Chemical">coumarin pharmacophores beariclass="Chemical">ng 7-class="Chemical">n class="Chemical">ether linkage and 4-methylimidazole moieties displayed different potency in their aromatase inhibitory effect. The QSAR study revealed many significant descriptors plausibly governing its interaction with aromatase. Such descriptors provided information on the types of functional groups related to their lipophilicity, polarity and isomeric effects as well as appropriate distance of functional groups.

In the series of class="Chemical">7-benzyloxy courmarin (1-13, R1 = CH2class="Chemical">n class="Chemical">C6H5, R2 = H), BELe1 values weighted by atomic Sanderson electronegativities were shown to be in the range of 1.920-1.923 where the phenyl groups (R1) were substituted with m- and p-Me, F, Cl, OMe, OCF3 and NO2. Interestingly, m-substituents such as F (3), Cl (4), OCF3 (6) exhibited higher aromatase inhibitory activity than their corresponding p-substituent (8, 9, 11) counterparts. The m-F compound 3 was shown to be the strongest aromatase inhibitor but compound 11 bearing the p-OCF3 moiety displayed the least activity. A closer look revealed that the compounds had the same value of BELel of 1.920 while affording different value for the maximal electrotopological positive variation (MAXDP) with corresponding values of 5.234, 4.937 and 4.914 for m-F (3), m-OCF3 (6) and m-Cl (4), respectively. In comparison with the corresponding p-substituted phenyl compounds (8, 9, 11), their MAXDP values (4.947, 4.930, 4.910, respectively) were relatively less than those of the m-substituted compounds (3, 4, 6). It can thus be presumed that an unsymmetrical m-substituted phenyl (R1) compounds (3, 4, 6) displayed higher MAXDP values as compared to a symmetrical p-substituted compounds (8, 9, 11). Similarly, the unsymmetrical 3,4-di F phenyl compound (13) had higher MAXDP value (5.271) than the symmetrical 3,5-di F analog (12) with MAXDP value of 4.870.

It was noted that for a series of class="Chemical">7-phenoxy coumarins (14-24, R1 = class="Chemical">n class="Chemical">C6H5, R2 = H), the unsymmetrical 3',4'-di F phenyl (24) had higher MAXDP value as compared to the symmetrical 3',5'-di F analog (23). As for the series of R1 (CH3 and C6H5) containing chiral center at the C-4 position of the coumarin ring (28-32), the 3',4'-di F phenyl analog 30 bearing phenyl group at the C-4 methine carbon atom had higher MAXDP value (5.503) when compared to the most potent aromatase inhibitor 24 (without chiral center) having MAXDP = 5.299. So far, the 3',4'-di F phenyl analog 24 displayed 6.74 folds higher aromatase inhibitory activity than analog 30. In addition, analog 24 exhibited higher activity (1.43 folds) than analog 29 containing phenyl moiety (R2) at the C-4 methine carbon but without 3',4'-di F substituted at the phenyl ring (R1). The results indicated that 3',4'-di F phenyl group (R1) predominantly governed the aromatase inhibitory activity as compared to the phenyl ring (R2). This could be attributed to the fact that compounds 29 and 30, which contain steric effect from phenyl group at the C-4 methine carbon, could prevent its binding or interaction with the target site of action. The 7-phenoxy compound 24 (without phenyl group at the methine carbon) was reported to bind the aromatase active site (heme iron) via coordination with the lone pair electron from N atom of the imidazole ring in such a way that the imidazole ring at the C-4 position is perpendicular to the coumarin ring (Stefanachi et al., 2011[42]).

class="Chemical">Coumarin aclass="Chemical">nalogs (R1 = class="Chemical">n class="Chemical">C6H5) without the 3',4'-di F group such as 31 (R2 = C6H4-Cl-p) and 32 (R2 = C6H4-CN-p) afforded high MAXDP (5.480 and 5.504) and high BELe1 (1.957 and 1.962) values, respectively. Notably, the inductive effect of the p-CN group (32) can provide the resonant ionic charge distribution form (32a, Figure 5(Fig. 5)) accounting for its higher MAXDP and BELe1 values than that of the p-Cl analog 31. Amongst the 4,7-disubstituted coumarins (1-32), analog 32 was shown to be the least potent aromatase inhibitor in which such remarkably low activity of 32 might be due to the ionic charge 32a and the bulky p-CN phenyl group (at chiral center).

Figure 5

Chemical structures of ionically charged resonant forms

On the other hand, class="Chemical">7-hydroxy 26 (R1 = H) aclass="Chemical">nd methoxy 27 (R1 = CH3) class="Chemical">n class="Chemical">coumarins showed low values of MAXDP (4.471 and 4.6000) and BELe1 (1.914 and 1.916). Such low values for these descriptors could be possibly a result from the lack of phenoxy group at the C-7 position (26 and 27) in contributing lone pair electron from the O atom into the phenyl ring (R1) leading to higher topological positive variation as compared to that of the 7-phenoxy compound 14 that can provide positively-charged molecule 14a (Figure 5(Fig. 5)).

Similarly, the electronic effect of the 7-phenylamino group (25) on the n class="Chemical">coumarin core structure could plausibly give rise to positively-charged species 25a (Figure 5(Fig. 5)) with high MAXDP value of 4.890. It should be class="Chemical">noted that all compouclass="Chemical">nds (1-34) share the commoclass="Chemical">n ioclass="Chemical">nic charge distributioclass="Chemical">n form (1A) as showclass="Chemical">n iclass="Chemical">n Figure 5(Fig. 5).

However, the highest MAXDP value (5.972) was noted for the class="Chemical">3-substituted imidazolylcoumarin 34 as compared to the class="Chemical">n class="Chemical">4-substituted imidazolylcoumarin analog 29 (MAXDP = 5.447). This could probably be due to the distance or position of substituents on the coumarin ring that affected its resonant ionic charge formation. The 3-imidazole substituent of coumarin 34 lying in close proximity to the carbonyl lactone that induces higher charge distribution (34a) as deduced from the high MAXDP value as compared to the corresponding value of 4-substituted coumarin 29 showing ionic charge formation (29a) as illustrated in Figure 5(Fig. 5). Apparently, the ionic form 34a provided the oxy anion that could attack the imine moiety of the imidazole ring to form negatively-charged N atom on the imidazole moiety (34b, Figure 5(Fig. 5)). This may be accounted for by the higher MAXDP value of 34 than that of compound 29.Taken together, the lipophilic property of R1 substituents play crucial role in affording appropriately high values of MAXDP and BELe1 that accounted for the potent aromatase inhibitory effect. Particularly, the phenoxy analog 24 bearing 3',4'-di F phenyl (R1) exerted the most potent activity as compared to the corresponding 3',4'-di F benzyl (R1) analog 13. This can be attributed to the correlation of descriptors as well as electronic effect of functional groups in contributing ionic charge resonant forms as well as the proper distance between R1 and the coumarin core structure in which the phenoxy (R1 = phenyl) coumarin analog gave the best fit in interacting with the target site.

Molecular docking of coumarins to aromatase

Insights on the binding modality of class="Chemical">coumarins to class="Chemical">n class="Gene">aromatase were elucidated by means of molecular docking. Particularly, low-energy conformers of coumarins obtained from aforementioned quantum chemical calculations were docked to the previously reported protein structure of aromatase (Suvannang et al., 2011[44]) that had been subjected to 100 ns of molecular dynamics simulation. The azole moiety is central to the interaction of non-steroidal AIs with the iron atom of the metalloporphyrrins (Balding et al., 2008[2]; Maurelli et al., 2011[25]; Pearson et al., 2006[40]). The handling of metal ions is not a trivial task and poses a great challenge in molecular docking and heme-containing proteins are particularly difficult to deal with owing to deficiencies in scoring function or the less stringency imposed by heme co-factors when compared to those of other metal ions (Irwin et al., 2005[18]; Seebeck et al., 2008[41]). Thus, to counter this limitation and ensure reliable interpretations we selected only docking conformations in which the azole moiety is oriented towards the heme co-factor for further analysis.

Binding poses with the lowest docked energy and belonging to the top-ranked cluster was selected as the final model for post-docking analysis with AutoDock Tools and PyMOL. Towards understanding the key residues involved in the protein-ligand interaction, we focused on residues within 4 Å around the docked ligands and this is depicted in Figure 6(Fig. 6) where panel a shows all ligands docked to the class="Gene">aromatase biclass="Chemical">ndiclass="Chemical">ng pocket while paclass="Chemical">nel b shows the most poteclass="Chemical">nt compouclass="Chemical">nd 24 as a represeclass="Chemical">ntative example. Residues fouclass="Chemical">nd iclass="Chemical">nteracticlass="Chemical">ng with most of the iclass="Chemical">nvestigated class="Chemical">n class="Chemical">coumarin chemotype were primarily lipophilic or non-polar residues comprising of Ile305, Ala306, Val370, Leu372, Met374, Leu477 and Ser478. Furthermore, polar and positively-charged residues included Thr310 and Arg115, respectively. Thus, it can be seen that aside from the azole-heme interaction and two polar residues, the protein-ligand interactions under investigation were predominantly lipophilic in nature. The aforementioned residues coincided with previously reported key residues in the binding cavity of aromatase (i.e. Ala306, Thr310, Met374 and Ser478) (Favia et al., 2013[10]; Lo et al., 2013[23]; Suvannang et al., 2011[44]). Moreover, as the crystal structure suggests, the dominance of lipophilic residues could be accounted for by its role in forming an androgen-specific cleft that binds snugly to the androstenedione (Ghosh et al., 2009[14]).

Figure 6

Docking poses of all 34 coumarins (a) and the most potent compound 24 (b) within the confinement of the aromatase binding pocket

A closer analysis of the docking poses was carried out by stratifying this set of 34 class="Chemical">coumarin aclass="Chemical">nalogs to 3 sub-classes (i.e. active, iclass="Chemical">ntermediate aclass="Chemical">nd iclass="Chemical">nactive) oclass="Chemical">n the basis of their pIC50 iclass="Chemical">n which compouclass="Chemical">nds haviclass="Chemical">ng values greater thaclass="Chemical">n 7 were coclass="Chemical">nsidered active, less thaclass="Chemical">n 6 as iclass="Chemical">nactive aclass="Chemical">nd values iclass="Chemical">n betweeclass="Chemical">n were desigclass="Chemical">nated as iclass="Chemical">ntermediate. The geclass="Chemical">neral key residues fouclass="Chemical">nd iclass="Chemical">n both active aclass="Chemical">nd iclass="Chemical">nactive sub-classes coclass="Chemical">nsisted of class="Chemical">n class="Chemical">Ala306, Thr310 and Met374. Particularly, active compounds were found to interact with Phe221, Trp224, Leu228, Ile305, Asp309, Leu477 and Ser478 however these residues were not involved in interacting with the inactives.

Conclusion

This study describes the QSAR study of class="Chemical">imidazole derivatives of class="Chemical">n class="Chemical">3,7- and 4,7-disubstituted coumarins having R1 and R2 substituents as inhibitors of the aromatase enzyme. Significant molecular descriptors were identified to include F10[N-O], Inflammat-50, Psychotic-80, H-047, BELe1, B10[C-O] and MAXDP, which were used in the construction of QSAR models using the MLR method. Multivariate analysis afforded good predictive performance for the cross-validated internal set with Q = 0.9239 and RMSE = 0.1304 while an external validation confirmed its robustness with Q = 0.7268 and RMSE = 0.2927. Insights on the structure-activity relationship of compounds were also discussed in light of the selected set of significant descriptors in concomitant with structural details from substituents R1 and R2. Both QSAR model and molecular docking investigations suggest that the aromatase inhibitory activity of compounds was primarily dependent on lipophilic properties and the position of substituent (R1) on the coumarin core structure. Structural knowledge gained from QSAR models and molecular docking could be used to guide the rational design of novel aromatase inhibitors.

Acknowledgements

This research project is supported by the Goal-Oriented Research Grant from Mahidol University to C.N. Talent Management fellowship to A.W. and research assistantship to N.S. are gratefully acknowledged.