Molecular determinants of thyroid hormone receptor selectivity in a series of phosphonic acid derivatives: 3D-QSAR analysis and molecular docking

Chemico-Biological Interactions 240 (2015) 324e335

Contents lists available at ScienceDirect

Chemico-Biological Interactions journal homepage: www.elsevier.com/locate/chembioint

Molecular determinants of thyroid hormone receptor selectivity in a series of phosphonic acid derivatives: 3D-QSAR analysis and molecular docking Fang-Fang Wang a, Wei Yang b, Yong-Hui Shi a, Guo-Wei Le a, * a b

The State Key Laboratory of Food Science and Technology, School of Food Science and Technology, Jiangnan University, Wuxi 214122, China Department of Biochemistry and Molecular Biology, Faculty of Medicine, Monash University, Melbourne, VIC 3800, Australia

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 May 2015 Received in revised form 16 July 2015 Accepted 3 September 2015 Available online 9 September 2015

A mathematical study was performed on a set of phosphonic acid derivatives that are substrates for thyroid hormone receptor b (TRb) and thyroid hormone receptor a (TRa), three-dimensional quantitative structure-activity relationship (3D-QSAR) models using comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) methods were employed to investigate the structural requirements for this series of compounds with improved activity. Some descriptors were also employed to significantly improve the performance of the derived models. The CoMFA model for TRb exhibited R2cv of 0.612, R2pred of 0.7218, whereas CoMSIA model showed R2cv of 0.621, R2pred of 0.7358; the CoMFA model for TRa displayed R2cv of 0.678, R2pred of 0.6424, and the CoMSIA model had R2cv of 0.671, R2pred of 0.6932, which indicate that the constructed models are statistically significant. The derived contour maps further pointed out the regions where interactive fields may influence the activity. In order to validate the QSAR models and explore the origin of the selectivity at the amino acid level, molecular docking was developed, and the results indicate that Arg282, Arg320, Asn331, Gly332, Thr329 and His435 for TRb, but Ala225, Arg228, Met259, Arg262 and His381 for TRa, respectively are important residues. The information obtained from the QSAR models can be used in the design of more potent TR agonists. © 2015 Elsevier Ireland Ltd. All rights reserved.

Keywords: TRb TRa CoMFA CoMSIA Molecular docking

1. Introduction Thyroid hormones have extensive effects on growth, development, metabolism and homeostasis [1,2], which are synthesized in the thyroid gland as thyroxine (T4) and 3,5,3-triiodo-L-thyronine (T3). T4 is the major form of thyroid hormone whereas T3 is considered the active form produced in target tissues by deiodination of T4 [3]. Most of the effects of T3 are mainly mediated by binding to the thyroid hormone receptors (TRs), which are members of the nuclear receptor superfamily of transcription regulators, and these receptors are involved in a variety of physiological and physiopathological processes [4e6]. TR is composed of some organizations that include a central conserved DNA-binding domain (DBD), a less conserved carboxyterminal ligand-binding domain (LBD), and an amino-terminal

* Corresponding author. E-mail address: [email protected] (G.-W. Le). http://dx.doi.org/10.1016/j.cbi.2015.09.008 0009-2797/© 2015 Elsevier Ireland Ltd. All rights reserved.

domain. The DBD domain anchors the TRs to specific DNA response elements (TREs). The amino-terminal domain serves as transcription regulator to influence the receptor functions. The LBD domain mainly binds the ligand and various corepressors and coactivators. When ligand binds to the LBD domain, the conformation of the LBD would be altered, and then facilitate dissociation of repressors and association of activators, further influence the way TRs bind to DNA [7]. Two different isoforms (TRa and TRb) are produced by alternative splicing of the two genes, TRa and TRb [8]. TRa is expressed during fetal development and is mainly located in the heart, while TRb emerges in development and predominates in the liver, kidney and lung [9]. Similar to the different expression patterns, the two isoforms also function distinct, for instance, TRa mainly controls cardiac output, whereas TRb helps in mediating the cholesterollowering [10e13]. Thus, ligands act on TRs would serve as promising therapeutic agents for the treatment of various diseases. In recent years, a number of TRs ligands have been developed, for example, L-94901 [14], Triac [15], CGS23425 [16] and GC-1 [17],

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

etc. However, selectivity issues in TRs binding exist and originate from the application of nonselective compounds, leading to serious side effects, particularly deleterious cardiovascular, i.e. arrhythmias, tachycardia, and heart failure [18]. Therefore, the identification and development of TRs ligands that bind preferentially to TRb is required. More recently, a series of phosphonic acid derivatives were synthesized and found to be TRb-selective agonists [19]. In fact, in vitro assessment of suspected TRs agonists remains a labor-intensive and time-consuming process. However, an in silico model as a more efficient and economical alternative way, could be employed for the purpose of predicting the binding effects of compounds on TRs. The ease of use, low cost and speed, combination with the known atomic resolution structure of TRs have made the efforts be mainly directed at the development of QSAR models [20e22], as well as the application of a 4D-6D QSAR model to design selective thyroid compounds [23,24]. To quantitatively identify the relationship between the structures and the activities, comparative molecular field analysis (CoMFA) [25] and comparative molecular similarity indices analysis (CoMSIA) [26], which make it possible to guide rational design of future ligands with improved activities and selectivity, have been employed. CoMFA and CoMSIA sample the steric, electrostatic, hydrophobic, hydrogen-bond donor and hydrogen-bond acceptor fields surrounding the aligned ligands. In addition, to reveal the detailed information of the receptoreligand interactions at molecular levels [27,28], molecular docking done by Autodock program was conducted. In conclusion, the developed QSAR models associating with the results of molecular docking can help to understand the binding process and provide additional insights into the potential structural modification and development of potent and selective ligands of TRs. 2. Materials and methods 2.1. Datasets The in vitro biological activity reported as Ki for agonists of TRb and TRa [19] was used for the current study. The Ki (nM) were transformed to pKi values and used as dependent variable in the CoMFA and CoMSIA 3D-QSAR models. Selection of appropriate training and test set is a significant step in determining the quality of the model. In the present work, the 3D-QSAR models were generated using a training set of twenty agonists. The predictive power of the generated models was further validated by a test set of seven molecules. The test set was selected by consideration of the structural diversity and wide range of binding activities in the dataset [29], and all of the ligands with their structures, Ki values, and the division of the training set and test set are listed in Table 1. All the agonists were minimized using the Powell conjugate gradient minimization algorithm and the Tripos force field [30], the energy gradient convergence criterion was set to 0.05 kcal/mol Å to obtain stable conformations, and the maximum iterations for minimization were set to 1000. Gasteiger-Hückel charges were assigned to each agonist. The prepared agonists were used as initial structures for 3D-QSAR model generation and molecular docking. 2.2. Molecular modeling and alignment Generally, the predictive ability of the models and the reliability of the contour maps mainly depend on the alignment procedure [31]. In this work, two different alignment methods were employed: template ligand-based alignment (superimposition I) and receptor-based alignment (superimposition Ⅱ). For superimposition I, the most active compounds (compound 4/TRb, compound 1/TRa) were employed as templates, then the remaining

325

Table 1 Molecular structures of Phosphonic Acid Derivatives and their TR binding affinity values (pKi).

Compound

R1

R2

R3

X

TRb1 pKi a

TRa1 pKi

1

I

I

O

9.1675

2

Me

Me

O

9.3279

9.2924

3

Me

Me

CH2

9.4949a

8.9626a

4

Cl

iPr

O

9.6198

8.1085a

Compound

R1

R2

X

Y

TRb1 pKi

TRa1 pKi

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

I Me Me Me Me Me Me OMe Me Me Me Me Cl Br Br I I Br I Me Cl Br Me

I iPr Me Et iPr Ph s-Bu iPr iPr iPr iPr iPr iPr iPr iPr iPr I iPr iPr iPr iPr iPr iPr

O O O CH2 CH2 CH2 CH2 CH2 CH2 CH2 CH2 O O O O O O O O O O O O

CH2CH(NH2) NHCO OCH2 OCH2 OCH2 OCH2 OCH2 OCH2 OCH2CH2 CH2CH2 NHCH2 CH2 CH2 CH2 CH2CH2 CH2 CH2 OCH2 OCH2 OCH2 NHCH2 NHCH2 SCH2

6.5670 7.4437a 7.0472 8.0123 8.5302 6.9747 8.1778a 6.8356 7.0605 7.5045 8.0410 6.8928a 7.2815 7.8097 8.6021 9.0757 7.5850a 8.6003 9.0132 7.5952a 7.5157 8.8210 7.3757

5.8489 6.5452a 6.0472 7.2104 7.4535 5.9337 7.5157a 5.4820 6.2262 6.4225 7.5467 5.7545a 6.5186 6.6289 7.3757 8.7352 7.7959a 7.7305 8.4271 6.7375a 6.7305 7.9872 6.8447

a

9.4815

Molecules belonged to the test set.

compounds were aligned to the common substructure of the template, the common substructures (marked in blue) and the aligned results are shown in Fig. 1. During the procedure of superimposition Ⅱ, each agonist was docked into the binding pocket of the receptors (TRb and TRa), then the structure of each agonist with highest docking binding affinities was chosen and aligned automatically inside the binding sites of TRb/TRa, the results are depicted in Fig. S1.

2.3. Molecular docking AutoDock [32], a grid-based docking program, was used for automatic placement of TR agonists in the binding pocket of TRb (PDB code: 1NAX) and TRa (PDB code: 1NAV). The two crystal

326

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

Fig. 1. (A) Compound 4 is used as TRb template for alignment. (B) Compound 1 is used as TRa template for alignment. (C) Present the alignments from superimposition I for TRb. (D) Present the alignments from superimposition I for TRa.

structures extracted from protein data bank (http://www.rscb.org/ pdb) are recombinant proteins expressed in Escherichia coli BL21(DE3), which belong to the species of Homo sapiens. During the process of docking, water molecules and related ligands (IH5) were removed from the crystal structures, then polar hydrogens and united atom Kollman charges were assigned for the receptors. All of the agonists were docked into the well prepared receptors. In order to validate the docking protocol, the interactions of IH5 against 1NAX and 1NAV were employed as standard models. Prior to molecular docking, grid maps were generated by the auxiliary program AutoGrid, the dimension of the grid box was set to 60
60
60 Å, which could enclose the reference ligand IH5, the box spacing was defined as 0.375 Å, and the grid center was designated at dimensions (xyz): 4.457
20.909
32.09 and 47.992
18.291
19.62, respectively for TRb and TRa. In conclusion, the current total grid points per map are 226981 for the two systems. The docking runs for each agonist were set to 100 with the standard AutoDock parameters. Finally, the maximum number of 10,000 retries and 27,000 generations were employed, respectively. The genetic algorithm with local search (GALS) was used to calculate the docking possibilities. Eventually, the conformations with the lowest binding energy were extracted and then aligned together for generating receptor-based 3D-QSAR models. 2.4. Calculation and selection of dragon descriptors Dragon descriptor generation routines have been successfully used for developing various QSAR models in our previous works [33,34]. 1644 molecular descriptors can be calculated, these descriptors are divided into 20 groups: constitutional descriptors, topological descriptors, walk and path counts, connectivity indices, information indices, 2D autocorrelations, edge adjacency indices, Burden eigenvalue descriptors, topological charge indices, eigenvalue-based indices, Randic molecular profiles, geometrical descriptors, RDF descriptors, 3D-MoRSE descriptors, WHIM descriptors, GETAWAY descriptors, functional group counts, atomcentered fragments, charge descriptors and molecular properties. Then feature selection was applied to reduce the dimensionality of data. Firstly, descriptors with constant values and containing 95% zero values were removed. Secondly, a stepwise linear regression method as the variable selection in R software (www.r-project.org) was further employed to reduce the descriptor space [35]. The Edge adjacency indices ESpm15r and the RDF descriptors RDF075v were

incorporated to enhance the robustness and generalization of the 3D-QSAR modes with R2 ¼ 0.651 and 0.630 for TRb and TRa, respectively.

2.5. 3D-QSAR studies 2.5.1. CoMFA interaction energy calculations For CoMFA analysis, the aligned agonists were placed in a grid box with a grid spacing of 2 Å in x, y and z directions. Steric and electrostatic fields were calculated using the Lennard-Jones potential and Coulomb potential, respectively [36e38]. A sp3 carbon probe atom with a van der Waals radius of 1.52 Å and a charge of þ1.0 was employed to compute the steric and electrostatic energies between the probe and the agonists. The energy values of these fields were truncated at 30 kcal/mol. Furthermore, in order to reduce noise, the column filtering was set to 2.0 kcal/mol, and those lattice points whose energy variation less than this value were excluded from the analysis.

2.5.2. CoMSIA interaction energy calculations In CoMSIA analysis, steric (S), electrostatic (E), hydrophobic (H), hydrogen bond donor (D) and hydrogen bond acceptor (A) fields were calculated using the same lattice box as that used in CoMFA analysis. It possesses some superiority over CoMFA, such as greater robustness regarding both region shifts and small shifts within the alignments, and more intuitively interpretable contour maps [39]. Default settings were used for the five fields: a probe atom with radius 1.0 Å, grid spacing 2 Å, charge þ1.0, hydrophobicity þ1, Hbond donor þ1, H-bond acceptor þ1 and attenuation factor a 0.3. A Gaussian method was applied to evaluate the distance between the probe atom and each agonist atom. And the similarity indices between the agonists of interest and the probe atom were calculated according to Equation (1):

AF;K q ðjÞ ¼ 

X

uprobe;k uik ear iq2

(1)

where q represents a grid point, i is the summation index over all atoms of the molecule j under computation, k represents five physicochemical properties (S, E, H, D and A), uik is the actual value of physicochemical property k of atom i, and uprobe;k is the value of the probe atom.

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

2.5.3. 3D-QSAR models generation and statistical validation 3D-QSAR models were generated by running partial least squares (PLS) analysis [40], which was used to linearly correlate the CoMFA and CoMSIA descriptors with the binding activities. The initial cross-validation was performed using leave-one out (LOO) method [41,42] in which one agonist was removed from the dataset and its binding activity was further evaluated by the model developed from the remaining agonists in the dataset, this method is useful to quantitatively evaluate the internal predictive ability of the model. The optimum number of components (ONC), the crossvalidated correlation coefficient (R2cv ) and the lowest standard error of prediction (SEP) were produced by PLS in conjunction with cross-validation option. Then the final model was developed by non-cross validated analysis, with the non-cross-validated correlation coefficient (R2ncv ), standard error of estimate (SEE) and F values generated. Additionally, the predictive power of the derived 3D-QSAR models was evaluated by predicting the biological activities of the test set which was not included in the training set. The predictive correlation coefficient R2pred was calculated by Equation (2):

R2pred ¼ ðSD  PRESSÞ=SD

(2)

where PRESS represents the sum of the squared deviations between the predicted and the actual pKi values for the test set compounds, SD represents the sum of the squared deviations between the biological activities of the test set molecules and the mean activity of the training set compounds. 3. Results and discussion 3.1. QSAR model In the present work, two different alignment rules were applied to investigate the structural and chemical features affecting the biological activities of the TRs ligands. To obtain superlative 3D-QSAR models, all 31 possible combinations of the fields were calculated, the statistical quality parameters of different interaction fields are shown in Table S1eS4. Additionally, the optimum 3D-QSAR models were listed in Table 2, which indicate that the models depend on superimposition Ⅰ are superior to those based on superimposition Ⅱ, this phenomenon might originate from the fact that during the procedure of molecular docking, compounds with different binding activities would show different orientation in the ligand binding pocket, and in the present work, all of the compounds in the binding site were employed to construct the QSAR models, leading to the results of the receptor-based alignment are worse than the template ligand-based alignment. Therefore, further discussions are mainly focused on superimposition Ⅰ models. 3.1.1. For TRb The CoMFA model gives a cross-validated correlation coefficient R2cv of 0.612 with an optimal number of principal components of 2, suggesting that the model is reliable and useful in predicting the binding activities of related ligands. The non-cross-validated analysis produces a non-cross-validated correlation coefficient R2ncv of 0.750, F value of 25.516, and a SEE value of 0.484. As obvious from Table 2, the percentage of the variance explained by steric and electrostatic is 13.0% and 13.9%, respectively, indicating that the two fields contribute equally to the binding activity. Furthermore, two dragon descriptors ESpm15r and RDF075v were incorporated to enhance the fitting degree of the model, which make 46.1% and 27.1% contributions to this model, respectively. The molecular descriptor types and definition of the symbols are shown in Table 3.

327

Table 2 Statistical data of optimal QSAR models. Parameters

R2cv R2ncv SEE F R2pred SEP ONC Field contribution S E H D A ESpm15r RDF075v

TRb

TRa

CoMFA

CoMSIA

CoMFA

CoMSIA

0.612 0.750 0.484 25.516 0.7218 0.602 2

0.621 0.763 0.471 27.315 0.7358 0.596 2

0.678 0.797 0.546 33.317 0.6424 0.688 2

0.671 0.797 0.546 33.410 0.6932 0.695 2

0.130 0.139 e e e 0.461 0.271

0.068 e 0.091 0.126 e 0.455 0.259

0.133 0.139 e e e 0.466 0.263

0.075 e 0.089 0.119 e 0.454 0.263

R2cv ¼ cross-validated correlation coefficient using the leave-one-out methods. R2ncv ¼ Non-cross-validated correlation coefficient; SEE¼Standard error of estimate; F ¼ Ratio of R2ncv explained to unexplained ¼ R2ncv /1  R2ncv ). R2pred ¼ Predicted correlation coefficient for the test set of compounds; SEP ¼ Standard error of prediction; ONC¼Optimal number of principal components; S ¼ steric, E ¼ electrostatic, H ¼ hydrophobic, D ¼ H-bond donor, A ¼ H-bond acceptor. Superimposition Ⅰ, template ligand-based alignment; Superimposition Ⅱ, receptorbased alignment.

Thereinto, ESpm15r [43] is an edge adjacency index weighted by resonance integrals, in the present work, the values of this descriptor for all ligands are resembled. In addition, RDF075v [44] belongs to the RDF descriptors which is based on a radial distribution function, the descriptor provides information about atomic van der Waals volumes. When the molecular size is enlarged, the descriptor increased. Further analysis reveals that compounds with less RDF075v values possess higher binding activity, suggesting that small molecules are advantageous to the activity. For CoMSIA model, the model with five fields is not the optimum model in the probability, it is because that the dependencies of each filed often decrease the signal-to-noise ratio in the data [45], and might reduce the statistical significance of the model. Actually, in this work, the optimal model is generated based on the steric, hydrophobic and hydrogen bond donor fields (SHD), statistically significant CoMSIA model gives a cross-validated coefficient R2cv of 0.621 with an optimal component of 2, R2ncv of 0.763, SEE of 0.471, and F value of 27.315. And the corresponding field contributions are 6.8%, 9.1%, 12.6%, 45.5% and 25.9%, for steric field, hydrophobic field, hydrogen bond donor field, ESpm15r descriptor and RDF075v descriptor, respectively, indicating a greater influence of the hydrogen bond donor field. In addition, the hydrophobic field has almost the same influence to the hydrogen bond donor field, indicating that the hydrophobic and hydrogen bond interactions of the molecules with the receptor would be two important factors for TRb agonistic activity. In order to validate the predictive power of the developed models, the binding activity of the test set is predicted, the predictive correlation coefficient R2pred is 0.7218 and 0.7358 for CoMFA and CoMSIA models, respectively, indicating excellent external predictive ability of the derived models. The correlation for the actual pKi versus the predicted pKi values for the whole dataset is shown in Fig. 2, the distribution pattern of the points further validates the robustness of the developed models. 3.1.2. For TRa The best CoMFA model applying steric and electrostatic fields has R2cv of 0.678 using 2 components, this model explains 79.7% of

328

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

Table 3 Descriptors used in model construction. Symbol

Class

Meaning

ESpm15r RDF075v

Edge adjacency indices RDF descriptors

Spectral moment 15 from edge adj. matrix weighted by resonance integrals Radial Distribution Function 7.5/weighted by atomic van der Waals volumes

Fig. 2. Graphs of the predicted versus the experimental pKi values of the optimal models for TRb. (A) CoMFA model. (B) CoMSIA model.

the variance, together with the non-cross-validated correlation coefficient R2ncv of 0.797, SEE of 0.546, and F value of 33.317. And the corresponding field and descriptor contributions of the parameters (steric, electrostatic, ESpm15r and RDF075v) are 13.3%, 13.9%, 46.6% and 26.3%, respectively. The above results indicate that the CoMFA model is statistically significant. Different combinations of the five fields were calculated, and the best predictions were obtained for the CoMSIA-SHD model (R2cv ¼ 0.671, R2ncv ¼ 0.797, SEE ¼ 0.546, F ¼ 33.410). Additionally, the relative contributions of steric, hydrophobic, hydrogen bond donor fields, ESpm15r and RDF075v descriptors are found to be 7.5%, 8.9%, 11.9%, 45.4% and 26.3%, respectively. This suggests that the model is reliable and capable of predicting novel ligands with improved activities. The CoMFA model and the CoMSIA model were also used to predict the binding activities of the test set, the values imply that the two models are able to describe the test set variance with R2pred of 0.6424 and 0.6932, respectively. Fig. 3 shows a relationship between the actual values and the predicted values. Linearity of the plots suggests that the proposed models are able to predict successfully the activities of ligands those are not used in the training process. 3.2. 3D-QSAR contour maps According to the optimal CoMFA and CoMSIA models, the 3D colored contour maps were built. The maps display regions where

differences in fields related to differences in binding activity. The generated contour maps can help identifying the significant features affecting the interactions between the ligands and the active site of TRs [46]. To aid in visualization, the representative compounds 4 and 9 for TRb, compounds 1 and 15 for TRa are labeled and displayed in the maps. 3.2.1. CoMFA contour maps for TRb The contour maps of CoMFA indicate the areas in space where the compounds would favorably and unfavorably interact with the receptor. As shown in Fig. S2A and S2B, the green contour maps represent regions of high steric bulk tolerance, while the yellow contour maps represent regions of low steric bulk tolerance. Some green irregular contour maps at the -meta position of ring A (designated in Fig. 1A) show that substituents at this position have favorable steric interactions. This is consistent with the experimental results. Compounds 9 and 20 have -iPr groups at this area exhibit higher activities than compounds 8 (-Et) and 21 (-I). A large green contour map is presented near the -ortho position of ring B, suggesting that favorable steric interaction might occur if bulky substituents are introduced in this region. For example, compound 2 possesses -Me at this position, covering the green map and hence accounts for greater activity than compound 1 (-I). However, for compounds 7e15, due to the linker between ring A and ring B (-CH2) is different from the other compounds (-O), thus the conformation of ring B for these compounds turns to an opposite direction, leading to the substituents at the -ortho position of ring B

Fig. 3. Graphs of the predicted versus the experimental pKi values of the optimal models for TRa (A) CoMFA model. (B) CoMSIA model.

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

towards a yellow contour map (as shown in Fig. S2B), which is unfavorable for bulkier substituents. Therefore, compound 12 with a larger substituent -OMe towards this contour distribution is less potent than compounds 7e11 with a relatively small substituent -Me. Additionally, the large green contour is also surrounding the groups at Region A (Fig. S2A), indicating that bulky substituents would be favorable in this region. In fact, compounds 19 (-CH2CH2), 26 (-NHCH2) and 27 (-SCH2) have bulky substituents in this region and possess higher activities than compounds 16 (-CH2) and 18 (-CH2). Conversely, as shown in Fig. S2B, the group at Region A falls into a sterically unfavorable yellow contour, indicating that the steric bulk at this position is detrimental to the activity, for instance, compounds 13e15 have an order for the activity of 13 < 14 < 15, with the corresponding substituents eOCH2CH2, eCH2CH2 and eNHCH2, respectively. For electrostatic contour maps (shown in Fig. S2C and S2D), the areas indicated by blue contours correspond to regions where electrostatic occupancy with electropositive groups enhance the activity, while the areas highlighted in red contours correspond to regions where electrostatic occupancy with electronegative groups would increase the activity. As can be seen from Fig. S2C, some blue contour maps embedded in the substituents at Region A suggest that an increase in electropositive substitution in this part might lead to enhanced TRb binding activity. This analyses give insight into the possible cause of why compound 19 (-CH2CH2) has higher activity in comparison to compound 22 (-OCH2). A large red colored favorable region is presented at Region B, which indicates that the electronegative substituents in this area are favorable. In fact, all of the compounds used in the present work with electronegative groups such as eCOOH, eH3PO3 touching this contour map, which indicates that the electronegative feature is crucial to the binding activity. Conversely, a red contour map is appearing at Region A for compound 9 (shown in Fig. S2D), which suggests that the high electron-density (negative charge) imparted by suitable groups in this part might lead to enhanced activity, for example, compound 9 and compound 15 with electronegative eOCH2 and eNHCH2 groups extended into the red-colored map have higher activity than compound 14 possessing eCH2CH2 at this position. Furthermore, the appearance of a blue and a red contour map in Region B indicate that careful selection should be made at this position. 3.2.2. CoMSIA contour maps for TRb Furthermore, CoMSIA model reveals two other contour maps that are not presented in CoMFA model. Similar results are obtained for the CoMSIA steric contour map (shown in Fig. S3A and S3B) as those of the CoMFA model. Thus, analyses are mainly focused on the hydrophobic and hydrogen bond donor contour maps (Fig. S3C and S3D), the yellow regions indicate areas where hydrophobic substituents are preferred and the white regions represent areas where hydrophilic groups are favored. A large white polyhedron showing an unfavorable hydrophobic interaction region is presented at the -meta position of ring A. The most potent compound 4 has a hydrophobic group in this area, thus further modifications can be made at this position. Another unfavorable hydrophobic region (white contour) is appeared near Region A. This is in agreement with the binding activity in compounds 19 (-CH2CH2) and 26 (-NHCH2). For compounds 7e15, a similar large white contour surrounding Region A depicts their favor for the hydrophilic substituents, which explains the lower activity of compound 14 (-CH2CH2) compared to that of compound 15 (-NHCH2). Moreover, at the -para position of ring B, a yellow contour map is located, suggesting that hydrophobic groups are favored in this site, in this work, all of the compounds bear hydrophobic substituents in this area. Additionally, the appearance of a large white contour at the linker between ring A and ring B indicates that substitution by

329

hydrophilic groups at this location is advantageous to improving the binding activity. The hydrogen bond donor contour maps are represented by cyan and purple colors (Fig. S3E and S3F), in the cyan regions of contour plot, the hydrogen bond donor substituents may enhance the activity, while the hydrogen bond donor groups in the purple area may decrease the activity. A purple contour around Region B (Fig. S3E) represents the enhanced activity of compounds if they possess a hydrogen bond acceptor at this positon, which is in accordance with the observations that the most active compound (compound 4) in the data set possesses a hydrogen bond acceptor substituent, such as eCOOH, indicating that hydrogen bond acceptor group at this location is necessary for the enhanced binding activity. The CoMSIA hydrogen bond donor contours for compounds 7e15 are shown in Fig. S3F, around the -ortho positon of ring B, the model tends to have prominent hydrogen bond acceptor favored contours, thus modifications can be made at this position to improve the activity. Some purple contour maps near the -para positon of ring B show that hydrogen bond donor groups are disfavored in this region. 3.2.3. CoMFA contour maps for TRa Similar to the contour maps of TRb, the poses for compounds 7e15 are significantly different from other compounds. Therefore, compounds 1 and 15 representing two distinct orientations are shown superimposed with the CoMFA and CoMSIA contour maps in Fig. S4. For CoMFA model, there are several green and yellow contour maps located in the active site (Fig. S4A and S4B). As can be seen in Fig. S4A, a large green contour map falls at the -ortho positon of ring B, indicating the occupancy of this contour by bulky groups would be favorable for the activity. The same green contour map touching the substituents at Region A indicates that a bulky group at this position would enhance the activity. A comparison of compounds 18 and 19 shows that group change from eCH2 to eCH2CH2 in this part might increase the activity. The relatively lower activity of compounds 2 (-NHCO), 3 (-OCH2) and 4 (-CH2) as compared to compound 1 (-CH2CH(NH2)) might due partly to the steric potency. In addition, compound 23 possessing eOCH2 group has higher potency than compound 18 (-CH2) possibly due to the green contour map. The CoMFA contour maps for compound 15 is shown in Fig. S4B, a yellow contour map located at the -ortho position of ring B indicates that minor groups would increase the potency. Compound 12 with an ortho-OMe substituent is less potent than compound 9 with a corresponding meta-Me substituent at this position. The yellow contour maps approaching Region A imply that bulky groups at this area would decrease the activity. When compared compound 13 with compound 9, we find that the eOCH2 group is slightly more favorable than the eOCH2CH2 group (compound 13) in this part. Another unfavorable steric region (yellow contour) is presented near Region B, which indicates that bulky substituents in this area may decrease the activity. The electrostatic fields are analyzed by the presence of blue and red contour maps. The seldom contour maps are mainly located at Region A and Region B. The red contour region around Region A and Region B suggests that the electronegative field might give a positive effect on the activity. Some potent compounds such as 22 (SCH2) and 27 (OCH2) have electronegative groups at this region. Nevertheless, due to only electropositive substituents at this region, thus the binding activity of several compounds (16-CH2, 18CH2, and 19CH2CH2) is lowered. Take the other conformation (compound 15) as an example, a blue contour map around Region A shows that the electropositive substituents are favored for the binding activity. This is in accordance with the experimental results. For instance, the order of the activity of some compounds is:

330

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

14 (CH2CH2) > 13 (OCH2CH2). Another blue contour map around Region B shows that occupancy of the electropositive favorable map would lead to positive effect on activity. Compounds 7e15 possess electron withdrawing groups at this position, which might enhance the replacement of electronegative groups with electropositive groups to improve the binding activity. 3.2.4. CoMSIA contour maps for TRa The CoMSIA steric contour maps (Fig. S5A and S5B) are found to be identical to the corresponding CoMFA contour maps. Therefore, the following discussions are mainly focused on the CoMSIA hydrophobic as well as the hydrogen bond donor contours (Fig. S5CeS5F). In the CoMSIA hydrophobic contour maps (Fig. S5C and S5D), the presence of a white contour at the -meta position of ring A suggests that hydrophilic substituents at this region would be responsible for higher activity, careful analyses suggest that modification is required to design enhanced novel compounds. Two small yellow contour maps are found to cover the substituents located at the -ortho position of ring B indicating hydrophobic groups here are beneficial for improving the binding activity. A large white contour map surrounding Region A indicates the significance of hydrophilic groups for biological activity. For example, compound 1, which has eCH2CH(NH2) substituent exhibits higher activity than compound 4 with group eCH2 at his position. And this phenomenon can explain the order of the binding activities in those compounds: 26 (NHCH2) > 18 (CH2). In addition, the white contour map is also approaching the groups at Region B, which indicates that hydrophilic groups are preferred at this position. For compound 15 (as shown in Fig. S5D), a white contour (different from compound 1) situated at the -ortho position of ring B illustrates that hydrophobic group is disfavored in this area, which further indicate that modifications should be made in this area to improve the activity. The hydrogen bond donor contour maps are shown in Fig. S5E and S5F. A cyan contour map surrounding Region A indicates that the presence of any substituents containing hydrogen bond donor groups would increase the activity. It can display the fact that the activity of compound 26 (NHCH2) is higher than compound 18 (-CH2) and 19 (CH2CH2). Moreover, because of the presence of eNHCH2 in compound 25 instead of eCH2 in compound 17 at Region A, thus the activity of compound 25 is more active than compound 17. Another cyan contour surrounding Region B indicates that hydrogen bond donor groups are beneficial to the binding activity. Conversely, take compound 15 for example (Fig. S5F), a purple polyhedron surrounded Region B, illustrating that hydrogen bond acceptor groups are favorable in the active binding pocket. This observation implies that further design of more potent compounds is required. 3.3. Molecular docking and comparison with 3D contour maps Molecular docking as an efficient in silico method is playing an increasing role in drug design [47,48]. Therefore, in the present work, it was employed to understand the action mechanism between the agonists and the receptors [49], and validate the reliability and accuracy of the derived QSAR models. All of the compounds were docked into the active sites of TRb and TRa, and the optimal conformation of each compound was aligned together as input files for 3D-QSAR analysis. In the present section, some key agonist structures that might demonstrate the most characteristic ligandereceptor interactions will be discussed. 3.3.1. TRb For this class of receptor, the attention will be focused especially

on compound 4 and compound 9. Fig. 4A and Fig. 4B display the binding orientation of compound 4 in the binding site of TRb, characterized by five hydrogen bonds involving Arg282, Arg320, Asn331 and His435 (Fig. 4B): the eOH group at the -para position of ring A is interacting with the nitrogen atom of His435 (-N$$$HO, 2.19 Å, 155.9 ) (H-1); the eCOOH at Region B is having hydrogen bond with the eNH of Asn331 (-O$$$HN, 3.21 Å, 129.9 ) (H-2), interaction is satisfied by the purple contour map appeared at this position (Fig. S3E); one hydrogen bond is also formed between Arg320 and the eCOOH group at Region B (-O$$$HN, 2.87 Å, 165.2 ) (H-3); The eCOOH group at Region B perhaps acts as an hydrogen bond acceptor to form hydrogen bond with amino acid Asn331 (-O$$$HN, 2.46 Å, 104.0 ) (H-4), one purple area near this position indicates that hydrogen bond acceptor substituents at ligand may enhance the binding activity; Furthermore, another hydrogen bond is formed between the eCOOH group and amino acid Arg282 (-O$$$HN, 1.97 Å, 145.8 ) (H-5), this hydrogen bond interaction in this part may favor the activity which is shown by the purple contour map at the hydrogen bond donor contour maps. Closer analysis illustrates that the -iPr substituent at the -meta position of ring A is trapped in a large pocket formed by Phe269, Phe272, Gly344, Gly345, Leu346 and Met442 (Fig. 4A), hence, bulky substitutions at this position will have a favorable steric interaction with TRb as depicted by the green contour map in this area (Fig. S3A and S4A). As shown in Fig. 4A, the groups at the -ortho position of ring B are placed in a large space surrounded by a few amino acid residues, suggesting that bulky groups are beneficial for ligandereceptor interactions, which is in agreement with the results obtained by QSAR studies represented by a large green polyhedra. For Region A, the orientation of eCH2 group bends 90 corresponding to the conformation of ring B, thus it locates into a large cavity, which suggests that large substituents at this position would increase the activity, this result is fully supported by the green contour map in this area. As illustrated in Fig. 4A, the groups at Region B interact with some electropositive amino acid residues, such as Arg282, Arg316 and Arg320, this is consist with the red contour map shown in Fig. S2C, in which the electronegative substituents are preferred. Furthermore, a wihte contour region (Fig.S3C) is in proximity to the substituent at the -meta position of ring A, which is hydrogen bonded with hydrophilic residue His435, indicating that hydrophilic groups are needed to increase the binding acitivity. The docked conformation of the other representitive conformation (compound 9) is shown in Fig. 4C and Fig. 4D, The ligand is anchored into the binding site by forming numerous hydrogen bonds with the side chains of Arg316, Thr329, Asn331, Gly332 and His435. The eOH at the -para position of ring A acts as a hydrogen bond donor and forms hydrogen bond with His435 (-N$$$HO, 2.05 Å, 165.9 ) (H-1); The oxygen atom at Region A forms hydrogen bond with Asn331 (-O$$$HN, 2.24 Å, 150.2 ) (H-2); The H2PO3 group at Region B severs as hydrogen bond donor forming two hydrogen bonds with several residues in the allosteric site: compound 9 … Arg316 (-O$$$HO, 2.00 Å, 105.1 ) (H-3) and compound 9 … Thr329 (-O$$$HO, 1.83 Å, 158.8 ) (H-4); The H2PO3 substituent at Region B also acts as hydrogen bond acceptor to form hydrogen bond with Asn331 (-O$$$HN, 2.13 Å, 98.3 ) (H-5) and Gly332 (-O$$$HN, 1.96 Å, 146.4 ) (H-6). The hydrogen bond interactions (H2, H-5 and H-6) confirm the observation from the CoMSIA hydrogen bond donor contour maps with purple contours occupying these areas. It should be noted that the side chain of Asn331 is appeared at Region A, and forms hydrogen bond with the ligand, therefore, it can be concluded that too large groups might have clash with the receptor, which is consistent with the yellow contour map (Fig. S2B and Fig. S3B). Substituents at Region A are also surrounded by neutral and electropositive amino acids, indicating that

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

331

Fig. 4. (A) The TRb active site amino acid residues around compound 4. (B) The enlargement for compound 4 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted red lines, and the nonpolar hydrogens were removed for clarity. (C) The TRb active site amino acid residues around compound 9. (D) The enlargement for compound 9 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted red lines, and the nonpolar hydrogens were removed for clarity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

electronegative groups in this area are favorable for enhancing the binding activity. This is in agreement with the electrostatic contour map, where the eOCH2 group is positioned near a red region. An additional interaction is also observed between the eOH of Region B and Arg316, which coincides with the red contour explained in the electrostatic contour map. Moreover, the other groups at Region B are in close contact with Thr329, Asn331 and Gly332, indicating that electropositive substitutions at this position may enhance the activity, this is fully supported by the electrostatic blue contour map. Besides, a hydrophobic cavity is formed among Ile275 and Met313, which accommodates the groups at the -ortho position of ring B, and this phenomenon can be inferred from the hydrophobic white contours (Fig. S3D). In addition, the hydrophobic contour map also indicates that the substituents at Region A is located near a large white contour map, further illustrating that hydrophilic substituents are favorable for the interactions with the amino acids in TRb, such as Arg282, Asn331 and Gly332. 3.3.2. TRa In Fig. 5A and Fig. 5B, the binding orientation of compound 1 in TRa binding pocket is illustrated. First of all, two hydrogen bonds are formed between compound 1 and some residues in receptor. It can be seen that the eOH at the -para position of ring A acts as a donor to form hydrogen bond with the backbone of His381 (-N$$$HO, 2.05 Å, 162.9 ) (H-1); as a hydrogen bond acceptor, the eCOOH at Region B also forms hydrogen bond with Arg262 (-O$$$HN, 2.39 Å, 173.5 ) (H-2). The active conformation in the binding site further support the generated QSAR contour maps, the groups at the -ortho position of ring B are almost extending outside the binding pocket, which suggests that bulky substituents may

increase the binding activity. This observation can be reflected by the green contour map located at this position (Fig. S4A and S5A). Another large hydrophilic cavity is composed of amino acid residues Ala263, Thr275 and Gly278, analysis reveals that substituents at Region A occupy this pocket, which corresponds to the green and white regions in the steric (Fig. S4A and S5A) and hydrophobic contour maps (Fig. S5C). In addition, the group at Region B is positioned in a hydrophilic pocket created by electropositive residues Arg228, Arg262 and neutral residue Gly278, thus, electronegative groups and hydrophilic interactions are correlated with enhanced activity. In the case of compound 15, it locates in the binding site of TRa by forming a network of hydrogen bonds with Ala225, Arg228, Met259 and His381: 1) the eOH at the -para position of ring A hydrogen bonding interaction with His381 (-N$$$HO, 2.12 Å, 163.5 ) (H-1); 2) hydrogen bonding interactions between the backbone of Met259 and the eNH group at Region A is observed (-O$$$HN, 2.34 Å, 153.1 ) (H-2); 3) the eH2PO3 group at Region B forms hydrogen bond with Ala225 (-O$$$HO, 2.12 Å, 154.2 ) (H-3); 4) another two hydrogen bonds are formed between the eH2PO3 group at Region B and Arg228 (-O$$$HN, 1.95 Å, 145.7 ) (H-4), (-O$$$HN, 3.16 Å, 63.7 ) (H-5), this can be indicated by the hydrogen bond donor purple contour map in the CoMSIA map (Fig. S5F). It is well-known that the eNH group at Region A is hydrogen bonding with the backbone eCO of Met259, thus, we suspect that small and electropositive substituents are suggested to provide favorable interactions, which is consist with the steric (yellow) and electrostatic (blue) contour maps. Another yellow contour map covering the groups at Region B indicates that this position should be occupied by less crowed substituents, it is

332

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

Fig. 5. (A) The TRa active site amino acid residues around compound 1. (B) The enlargement for compound 1 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted red lines, and the nonpolar hydrogens were removed for clarity. (C) The TRa active site amino acid residues around compound 15. (D) The enlargement for compound 15 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted red lines, and the nonpolar hydrogens were removed for clarity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

because that this group falls into a narrow pocket comprised by residues Ala225, Ile226, Arg228 and Arg262. The results of molecular docking and QSAR contour maps can validate each other, and offer useful suggestions for further modification. 3.4. Comparison for TRb and TRa In order to illustrate the relationship between the activities of agonists employed in the present work and the two isoform receptors, TRb and TRa, a mathematical model was developed, the plot (Fig. 6A) reveals that the activities of the compounds to bind the receptors TRb and TRa are obviously linear correlated with a correlation coefficient of R2 ¼ 0.8604 (R2 ¼ coefficient of determination). The correlation between them illustrates that activity for TRb is in agreement with that for TRa, which further suggests that generally the TRb agonists can also bind to TRa, and they share a similar binding mode in the active site, which can also be validated by the docked conformations in the binding pocket of TRb and TRa (Fig. 6B), for TRa, the spatial distribution of favorable and unfavorable contour maps (Fig. S4 and Fig. S5) is nearly similar to TRb (Fig. S2 and Fig. S3), which is corresponding to high degree of conservation of the amino acids in the binding site of TRb and TRa (Fig. 6C). Conversely, differences are also presented in the contour maps when compared TRb with TRa, further analysis indicates that one significant amino acid residue varies between the two receptors (Asn331-Ser277) (Fig. 6C). In addition, differences are also observed for several compounds in the same receptor derived from

the linker between ring A and ring B, the replacement of the oxygen atom between ring A and ring B with a eCH2 group causes an anticlockwise shift of the groups at ring B, this can also be indicated by the contour map analyses. Several interesting phenomenon is discovered by analyzing the two representative conformations (TRb: 4 and 9; TRa: 1 and 15). For TRb-4 and TRa-1, the hydrogen bond interactions in Region A and Region B, to some extents, influence the selectivity for TRb over TRa; For TRb-9 and TRa-15, the charge properties in Region A can account for the selectivity. In a word, differences in the contour maps combined with the molecular docking results indicate that a rational design of agonists selective for TRb over TRa should be possible by the application of the QSAR models.

4. Conclusions An important aspect in the present work is the application of the CoMFA and CoMSIA methods to a same dataset of compounds towards TRb and TRa, the models exhibit high predictive power in both internal and external test. The predictive ability of the derived models is further validated by contour maps, which provide a fruitful insight into the different fields contributing to the activity. In addition, the useful information inferred from the contour maps correlates well with the molecular docking results, validating the robustness and reliability of the models, and the relevant structural conclusions are illustrated in the following (Fig. 7): (1) for TRb-compound 4, bulky and hydrophilic groups at the -meta position of ring A, bulky groups at the -ortho position of ring B, bulky,

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

333

Fig. 6. (A) A correlation plot of activities for TRb and TRa of all compounds. (B) Superimposition of TRb and TRa. (C) Alignments of the sequences of TRb and TRa, light green color regions denote that the amino acid residues in the individual column are identical in the sequence alignment. Black triangles are placed under the identical and key residues in the binding pocket. The black rectangle denote that the amino acid residue is different in the binding site. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Structure-activity relationship revealed by QSAR studies for TRs. (A) Compound 4 for TRb. (B) Compound 1 for TRa. (C) Compound 9 for TRb. (D) Compound 15 for TRa.

electropositive and hydrophilic groups at Region A, electronegative, hydrogen bond acceptor groups at Region B are beneficial to the activity (Fig. 7A); (2) for TRa-compound 1, there should be bulky and hydrophobic groups at the -ortho positon of ring B, bulky, electronegative, hydrophilic and hydrogen bond donor groups at Region A, electronegative, hydrophilic and hydrogen bond donor groups at Region B (Fig. 7B); (3) for TRb-compound 9, hydrophilic groups at the linker between ring A and ring B, minor and hydrogen

bond acceptor groups at the -ortho position of ring B, minor, electronegative, hydrophilic and hydrogen bond acceptor groups at Region A can improve the activity (Fig. 7C); (4) for TRa-compound 15, minor and hydrophilic groups at the -ortho position of ring B, minor and electropositive groups at Region A, minor, electropositive and hydrogen bond acceptor groups at Region B are crucial for the activity (Fig. 7D); (5) The activities of agonists are well correlated, indicating that similar binding mode might exist for this

334

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335

series of compounds when bind to TRb and TRa. And key amino acids influencing the receptoreligand interactions are Arg282, Arg320, Asn331, Gly332, Thr329 and His435 for TRb, but Ala225, Arg228, Met259, Arg262 and His381 for TRa, respectively. The selectivity is useful in rational drug discovery and design of novel compounds with enhanced potency and receptor discrimination, the contour maps and the molecular docking analysis reveal the structural requirements for the observed selectivity of the two isoforms of TRs. The changes in amino acid residues like Asn331 in TRb compared to Ser277 in TRa are possibly related to the structural differences and selectivity. As a conclusion, the developed QSAR models give an insight into the different field contributions around the aligned compounds and the influential factors contributing to the activity. The developed models can help in guiding the design of novel agonists with enhanced activity and improved selectivity.

[16]

[17]

[18] [19]

[20]

[21]

[22]

Acknowledgments The study was supported by the 12th five-year plan for science and technology development (No. 2012BAD33B05). Appendix A. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.cbi.2015.09.008.

[23] [24]

[25]

[26]

Transparency document Transparency document related to this article can be found online at http://dx.doi.org/10.1016/j.cbi.2015.09.008.

[27] [28]

[29]

References [30] [1] P.M. Yen, Physiological and molecular basis of thyroid hormone action, Physiol. Rev. 81 (2001) 1097e1142. [2] R.D. Utiger, The thyroid: physiology, thyrotoxicosis, hypothyroidism, and the painful thyroid, Endocrinol. Metabol. 3 (1995) 435e519. [3] A.C. Bianco, D. Salvatore, B. Gereben, M.J. Berry, P.R. Larsen, Biochemistry, cellular and molecular biology, and physiological roles of the iodothyronine selenodeiodinases, Endocr. Rev. 23 (2002) 38e89. [4] N.H. Nguyen, J.W. Apriletti, J.D. Baxter, T.S. Scanlan, Hammett analysis of selective thyroid hormone receptor modulators reveals structural and electronic requirements for hormone antagonists, J. Am. Chem. Soc. 127 (2005) 4599e4608. [5] P.M. Yen, S. Ando, X. Feng, Y. Liu, P. Maruvada, X. Xia, Thyroid hormone action at the cellular, genomic and target gene levels, Mol. Cell. Endocrinol. 246 (2006) 121e127. [6] J. Zhang, M.A. Lazar, The mechanism of action of thyroid hormones, Annu. Rev. Physiol. 62 (2000) 439e466. [7] R.C. Ribeiro, J.W. Apriletti, R.L. Wagner, W. Feng, P.J. Kushner, S. Nilsson, T.S. Scanlan, B.L. West, R.J. Fletterick, J.D. Baxter, X-ray crystallographic and functional studies of thyroid hormone receptor, J. Steroid Biochem. Mol. Biol. 65 (1998) 133e141. € m, Functions of thyroid hormone receptors in mice, [8] D. Forrest, B. Vennstro Thyroid 10 (2000) 41e52. [9] H. Schwartz, K. Strait, N. Ling, J. Oppenheimer, Quantitation of rat tissue thyroid hormone binding receptor isoforms by immunoprecipitation of nuclear triiodothyronine binding capacity, J. Biol. Chem. 267 (1992) 11794e11799. [10] H. Escriva, S. Bertrand, V. Laudet, The evolution of the nuclear receptor superfamily, Essays Biochem. 40 (2004) 11e26. [11] S.Y. Cheng, Isoform-dependent actions of thyroid hormone nuclear receptors: lessons from knockin mutant mice, Steroids 70 (2005) 450e454. [12] Y. Murata, Multiple isoforms of thyroid hormone receptor: an analysis of their relative contribution in mediating thyroid hormone action, Nagoya J. Med. Sci. 61 (1998) 103e116. [13] F.E. Wondisford, Thyroid hormone action: insight from transgenic mouse models, J. Invest. Med. Offi. Publ. Am. Fed. Clin. Res. 51 (2003) 215e220. [14] A. Underwood, J. Emmett, D. Ellis, S. Flynn, P. Leeson, G. Benson, R. Novelli, N. Pearce, V. Shah, A thyromimetic that decreases plasma cholesterol levels without increasing cardiac activity, Nature 324 (1985) 425e429. [15] P. Schueler, H. Schwartz, K. Strait, C. Mariash, J. Oppenheimer, Binding of 3, 5, 30 -triiodothyronine (T3) and its analogs to the in vitro translational products

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41] [42]

of c-erbA protooncogenes: differences in the affinity of the a-and b-forms for the acetic acid analog and failure of the human testis and kidney a-2 products to bind T3, Mol. Endocrinol. 4 (1990) 227e234. A.H. Taylor, Z.F. Stephan, R.E. Steele, N.C. Wong, Beneficial effects of a novel thyromimetic on lipoprotein metabolism, Mol. Pharmacol. 52 (1997) 542e547. G. Chiellini, J.W. Apriletti, H. Al Yoshihara, J.D. Baxter, R.C. Ribeiro, T.S. Scanlan, A high-affinity subtype-selective agonist ligand for the thyroid hormone receptor, Chem. Biol. 5 (1998) 299e306. F.H. Epstein, I. Klein, K. Ojamaa, Thyroid hormone and the cardiovascular system, N. Engl. J. Med. 344 (2001) 501e509. S.H. Boyer, H. Jiang, J.D. Jacintho, M.V. Reddy, H. Li, W. Li, J.L. Godwin, W.G. Schulz, E.E. Cable, J. Hou, Synthesis and biological evaluation of a series of liver-selective phosphonic acid thyroid hormone receptor agonists and their prodrugs, J. Med. Chem. 51 (2008) 7075e7093. J. Du, J. Qin, H. Liu, X. Yao, 3D-QSAR and molecular docking studies of selective agonists for the thyroid hormone receptor b, J. Mol. Graph. Model. 27 (2008) 95e104. Y. Ren, H. Liu, S. Li, X. Yao, M. Liu, Prediction of binding affinities to b 1 isoform of human thyroid hormone receptor by genetic algorithm and projection pursuit regression, Bioorg. Med. Chem. Lett. 17 (2007) 2474e2482. F. Li, Q. Xie, X. Li, N. Li, P. Chi, J. Chen, Z. Wang, C. Hao, Hormone activity of hydroxylated polybrominated diphenyl ethers on human thyroid receptor-b: in vitro and in silico investigations, Environ. Health Perspect. 118 (2010) 602e606. A. Vedani, M. Dobler, M.A. Lill, The challenge of predicting drug toxicity in silico, Basic Clin. Pharmacol. Toxicol. 99 (2006) 195e208. A. Vedani, M. Zumstein, M.A. Lill, B. Ernst, Simulating a/b selectivity at the human thyroid hormone receptor: consensus scoring using multidimensional QSAR, ChemMedChem. 2 (2007) 78e87. H. Hong, H. Fang, Q. Xie, R. Perkins, D.M. Sheehan, W. Tong, Comparative molecular field analysis (CoMFA) model using a large diverse set of natural, synthetic and environmental chemicals for binding to the androgen receptor, SAR QSAR Environ. Res. 14 (2003) 373e388. G. Klebe, U. Abraham, T. Mietzner, Molecular similarity indices in a comparative analysis (CoMSIA) of drug molecules to correlate and predict their biological activity, J. Med. Chem. 37 (1994) 4130e4146. P. Kirkpatrick, Gliding to success, Nat. Rev. Drug Discov. 3 (2004), 299e299. X. Sun, Y. Li, W. Li, Z. Xu, Y. Tang, Computational investigation of interactions between human H 2 receptor and its agonists, J. Mol. Graph. Model. 29 (2011) 693e701. X. Li, L. Ye, X. Wang, W. Shi, H. Liu, X. Qian, Y. Zhu, H. Yu, In silico investigations of anti-androgen activity of polychlorinated biphenyls, Chemosphere 92 (2013) 795e802. M. Clark, R.D. Cramer, N. Van Opdenbosch, Validation of the general purpose Tripos 5.2 force field, J. Comput. Chem. 10 (1989) 982e1012. R. Thaimattam, P.R. Daga, R. Banerjee, J. Iqbal, 3D-QSAR studies on c-Src kinase inhibitors and docking analyses of a potent dual kinase inhibitor of c-Src and c-Abl kinases, Bioorg. Med. Chem. 13 (2005) 4704e4712. G.M. Morris, D.S. Goodsell, R.S. Halliday, R. Huey, W.E. Hart, R.K. Belew, A.J. Olson, Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function, J. Comput. Chem. 19 (1998) 1639e1662. F. Wang, Z. Ma, Y. Li, J. Wang, Y. Wang, Structural requirements of pyrimidine, thienopyridine and ureido thiophene carboxamide-based inhibitors of the checkpoint kinase 1: QSAR, docking, molecular dynamics analysis, J. Mol. Model. 18 (2012) 3227e3242. F. Wang, W. Yang, Y. Shi, G. Le, Structural analysis of selective agonists of thyroid hormone receptor b using 3D-QSAR and molecular docking, J. Taiwan Inst. Chem. Eng. 49 (2015) 1e18. Y. Wang, Y. Li, J. Ding, Z. Jiang, Y. Chang, Estimation of bioconcentration factors using molecular electro-topological state and flexibility, SAR QSAR Environ. Res. 19 (2008) 375e395. € hm, J. Stürzebecher, G. Klebe, Three-dimensional quantitative structureM. Bo activity relationship analyses using comparative molecular field analysis and comparative molecular similarity indices analysis to elucidate selectivity differences of inhibitors binding to trypsin, thrombin, and factor Xa, J. Med. Chem. 42 (1999) 458e477. X.H. Liu, Y.X. Shi, Y. Ma, C.Y. Zhang, W.L. Dong, L. Pan, B.L. Wang, B.J. Li, Z.M. Li, Synthesis, antifungal activities and 3D-QSAR study of N-(5-substituted-1, 3, 4thiadiazol-2-yl) cyclopropanecarboxamides, Eur. J. Med. Chem. 44 (2009) 2782e2786. R.D. Cramer, D.E. Patterson, J.D. Bunce, Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins, J. Am. Chem. Soc. 110 (1988) 5959e5967. S. Pirhadi, J.B. Ghasemi, 3D-QSAR analysis of human immunodeficiency virus entry-1 inhibitors by CoMFA and CoMSIA, Eur. J. Med. Chem. 45 (2010) 4897e4903. S. Wold, A. Ruhe, H. Wold, I. Dunn, WJ, The collinearity problem in linear regression. The partial least squares (PLS) approach to generalized inverses, SIAM J. Sci. Stat. Comput. 5 (1984) 735e743. M. Stone, Cross-validatory choice and assessment of statistical predictions, J. R. Stat. Soc. Ser. B Methodol. (1974) 111e147. R.D. Cramer, J.D. Bunce, D.E. Patterson, I.E. Frank, Crossvalidation, bootstrapping, and partial least squares compared with multiple regression in

F.-F. Wang et al. / Chemico-Biological Interactions 240 (2015) 324e335 conventional QSAR studies, Quant. Structure-Activity Relat. 7 (1988) 18e25. [43] Y. Liang, D. Yuan, Q. Xu, O.M. Kvalheim, Modeling based on subspace orthogonal projections for QSAR and QSPR research, J. Chemom. 22 (2008) 23e35. [44] J. Xu, L. Zhu, D. Fang, L. Liu, Z. Bai, L. Wang, W. Xu, A simple QSPR model for the prediction of the adsorbability of organic compounds onto activated carbon cloth, SAR QSAR Environ. Res. 24 (2013) 47e59. [45] U. Norinder, Recent Progress in CoMFA Methodology and Related Techniques, 3D QSAR in Drug Design, Springer, 1998, pp. 25e39. [46] Y. Mao, Y. Li, M. Hao, S. Zhang, C. Ai, Docking, molecular dynamics and

335

quantitative structure-activity relationship studies for HEPTs and DABOs as HIV-1 reverse transcriptase inhibitors, J. Mol. Model. 18 (2012) 2185e2198. [47] I.D. Kuntz, Structure-based strategies for drug design and discovery, Science 257 (1992) 1078e1082. [48] J. Drews, Drug discovery: a historical perspective, Science 287 (2000) 1960e1964. [49] K.S. Bhadoriya, M.C. Sharma, S.V. Jain, G.S. Raut, J.R. Rananaware, Threedimensional quantitative structureeactivity relationship (3D-QSAR) analysis and molecular docking-based combined in silico rational approach to design potent and novel TRPV1 antagonists, Med. Chem. Res. 22 (2013) 2312e2327.