In silico study on β-aminoketone derivatives as thyroid hormone receptor inhibitors: a combined 3D-QSAR and molecular docking study

In order to explore the structure–activity correlation of a series of β-aminoketone analogs as inhibitors of thyroid hormone receptor (TR), a set of three-dimensional quantitative structure–activity relationship (3D-QSAR) models based on comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA), for the first time, were developed in the present work. The best CoMFA model with steric and electrostatic fields exhibited , for TRβ, and , for TRα. 3D contour maps produced from the optimal models were further analyzed individually, which provide the areas in space where interactive fields would affect the inhibitory activity. In addition, the binding modes of inhibitors at the active site of TRs were examined using molecular docking, the results indicated that this series of inhibitors fit into the active site of TRs by forming hydrogen bonding and electrostatic interactions. The docking studies also revealed that Leu305, Val458 for TRβ, and Asp407 for TRα are showing hydrogen bonds with the most active inhibitors. In any case, the 3D-QSAR models combined with the binding information will serve as a useful approach to explore the chemical space for improving the activity of TRβ and TRα inhibitors.


Introduction
Nuclear hormone receptor (NR) superfamily consists of the largest families of ligand-dependent transcription factors and are important for the regulation of development, growth, and metabolism (Aranda & Pascual, 2001;Ribeiro et al., 1995). The family receptors are composed of estrogen receptors, androgen receptors, retinoid X receptors, peroxisome proliferator-activated receptors, orphan receptors, as well as receptors those do not have endogenous ligands or for which the ligand has yet not to be identified (Gronemeyer, Gustafsson, & Laudet, 2004), which can regulate gene transcription in response to molecule hormones.
As part of the NR superfamily, the thyroid hormone receptor (TR) mediates the actions of thyroid hormone, mainly T3 (3,5,3′-triiodo-L-thyronine), which is secreted by thyroid gland, however, the transcriptional activity can be excited with and without ligands, it is because that they bind to the chromatin constitutively, and form heterodimers with retinoid X receptors.
Two different TR subtypes are expressed by two different genes, TRα and TRβ (Lazar, 1993). TRα 1 , TRα 2 , TRβ 1 , and TRβ 2 are generated through the process of differential RNA splicing (Wagner et al., 2001), and these isoforms mainly differ in the amino-terminal transactivation domain and the carboxyl-terminal region. The TRα 1 is prevalent most in brown fat and skeletal muscle, while TRα 2 is expressed heavily in brain, does not bind thyroid hormones, and acts as a transcriptional repressor. TRβ 1 is expressed in most tissues, such as liver and brain. TRβ 2 is widely distributed in the pituitary and hypothalamus. Similar to the tissue distribution, a number of transgenic and knock-out mouse models indicate that the TRα and TRβ isoforms possess distinct and tissue-specific functions. For instance, TRα regulates heart rate and contractility, while TRβ controls cholesterol metabolism and thyroid-stimulating hormone production (Harvey & Williams, 2002).
As we know, it is time-consuming for the synthesis of compounds with potent inhibitory activity, however, an in silico model as economic and facile alternative to in vitro methods to evaluate the activity would help in designing novel inhibitors prior to experiment (Ai, Wang, Li, Li, & Yang, 2008;Wang, Li, Yang, & Yang, 2005). The main computational methods employed in drug design can be divided into two parts: (I) ligandbased studies, which include quantitative structure-activity relationship analysis (Ortiz et al., 1997) and (II) receptor-based studies, mainly known as molecular docking (Mouchlis, Mavromoustakos, & Kokotos, 2010). Comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA), wellknown three-dimensional quantitative structure-activity relationship (3D-QSAR) method, are powerful tools in drug design (Jayatilleke, Nair, Zauhar, & Welsh, 2000). CoMFA calculates steric and electrostatic fields, while CoMSIA calculates steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields surrounding the aligned ligands, in addition, the QSAR models are constructed by connecting the fields with the observed inhibitory activities. Molecular docking is an extensively used method in rational drug design, which was performed to identify the probable conformations and investigate the ligand-receptor interactions. Until recently, the crystal structures of TR in complex with related compounds have been reported Souza et al., 2014), which further promote the research done in the present work.
Our goal was to determine the requisite structure features for a series of β-aminoketone analogs that specially interact with the AF-2 domain of TR. Therefore, theoretical computations, including CoMFA and CoMSIA methods have been employed. In addition, molecular docking was also performed to probe the structural properties and binding modes of the inhibitors at the binding site of TR. All the methods applied in the present work could help in understanding the molecular basis of receptorligand interactions, and further promoting the development of novel and potent TR inhibitors.
2. Methods and materials 2.1. Data-sets and biological activity β-aminoketone and its analogs along with their inhibitory activities taken from (Arnold, Kosinski, Estébanez-Perpiñá, & Guy, 2007;Hwang et al., 2009) were used for the 3D-QSAR studies. The IC 50 values, in μM, were converted into pIC 50 (−log IC 50 ) values, which were further used as dependent variables in developing the 3D-QSAR models. The data-set was split into a training set (42 for TRβ, 23 for TRα) and a test set (14 for TRβ, 7 for TRα) by considering the fact that the inhibitors in the test set represent structural diversity and a range of inhibitory activities similar to those of the training set. Then, the predictive ability of the derived models was evaluated by the test set inhibitors. The structures and inhibitory activities of these inhibitors are shown in Supplementary Table S1.

Molecular docking
Molecular docking was performed to identify the optimal conformation of the inhibitors and to calculate the receptor-ligand interactions. Therefore, AutoDock 4.2.5.1 suite (Morris et al., 1998) as the molecular docking tool was used to perform the molecular docking procedure. The crystal structures of TRβ and TRα (PDB entry code 2PIN and 4LNX) were extracted from Protein Data Bank (http://www.rscb.org/pdb). Prior to molecular docking, the original co-crystalized ligands, water molecules and ions were removed, then polar hydrogen atoms were added, Kollman charges were assigned to the receptors. The binding mode of the inhibitors was obtained by a grid-based docking program. The interaction energy between the inhibitor and the receptor was evaluated using atom affinity potentials calculated on a grid similar to that described by Goodford (1985).
AutoGrid was used for the preparation of the grid map using a grid box with a docking box of 60 × 60 × 60 Å, the box spacing was .375 Å. Finally, AutoDock was run using maximum number of 10,000 retries and 27,000 generations, respectively. The docking possibilities were calculated by the genetic algorithm with local search. For each inhibitor, the simulation was composed of 100 docking runs using the standard Auto-Dock parameters.
Before docking, a redocking process of the cocrystallized ligands LEG and T44 into 2PIN and 4LNX, respectively, were performed to verify the docking parameters used and recover the known complex's structure and interactions. Figure 1 represents the superimposed structures of the crystal and the redocked results, as can be seen that the original structure and the docked structure occupy the same position and possess similar conformation, which suggest that there are no significant differences between the structure extracted from the crystal structure and the docked model of the complex, further indicating that the docking model is rational and valid.

Molecular modeling and molecular alignment
All the three-dimensional structures of the inhibitors were constructed and then minimized using the Powell gradient algorithm with the Tripos force field (Clark, Cramer, & Van Opdenbosch, 1989) and Gasteiger-Hückel charges (Gasteiger & Marsili, 1980). Furthermore, the energy gradient convergence criterion and the maximum iterations were set to .05 kcal/mol Å and 1000, respectively. Finally, the stable and fairly fixed conformations were obtained, which would endow the compounds with a comparable conformation and a similar spatial orientation of the major groups.
In building the 3D-QSAR models, the ligand conformation and molecular alignment are two significant factors (Thaimattam, Daga, Banerjee, & Iqbal, 2005). Actually, the accuracy of the 3D-QSAR models and the reliability of the contour maps mainly depend on the alignment procedure. To derive the statistically significant QSAR models, two patterns of alignment rules were applied. The first superimposition is template ligand-based superimposition (superimposition I), during the procedure, the most potent inhibitors (21 for TRβ, 24 for TRα) in the data-set were selected as templates, and then, the remaining compounds were aligned to the common substructures of the templates (Figures 2(A) and (C), marked in blue). And the superimposition of the data-set is shown in Figures 2(B) and (D). The receptorbased superimposition (superimposition II) as the second alignment is shown in Supplementary Figure S1. According to this strategy, molecular alignment was performed on the conformations (with the lowest binding energy) obtained from molecular docking. However, the conformations for all of the compounds with the lowest binding energy in TRs could not exhibit a significant result. Thus, the optimal conformation of each compound was chosen to construct the QSAR models. However, the QSAR models based on superimposition I are superior to those based on superimposition II, actually, some models based on superimposition II also possess high R 2 cv , but exhibit unfavorable predictive power, or the points distribute unevenly around the regression line, in addition, this result is consistent with the works in references (Wang, Ma, Li, Wang, & Wang, 2012;Wang, Yang, Shi, & Le, 2015), in these papers, different alignment rules were also employed to conduct the QSAR models, and the models based on template ligand-based alignment appear to be superior among all of the models derived. This is because that during the process of molecular docking, compounds with different binding activities would show different orientation in the ligand binding pocket, leading to the results of alignment are worse than the template ligand-based alignment. Therefore, discussions are mainly located on the template ligand-based models in the following sections. Figure 1. (A) Structural superposition of the crystal structure (green) and the docked structure (cyan) for TRβ. The projection highlights the structure of the active site with LEG, which is displayed in sticks. (B) Structural superposition of the crystal structure (white) and the docked structure (yellow) for TRα. The projection highlights the structure of the active site with T44, which is displayed in sticks.

CoMFA and CoMSIA analyses
In CoMFA analysis, steric and electrostatic fields were calculated ) based on the alignment results. The aligned inhibitors were placed in a 3D cubic lattice with a grid spacing of 2 Å in x, y, and z directions, and the grid pattern was generated by the CoMFA routine automatically. A sp 3 -hybridized carbon atom with a van der Waals radius of 1.52 Å as a steric probe and a charge of +1.0 as an electrostatic probe was used to calculate the steric (Lennard-Jones potential) filed energies and the electrostatic (Coulomb potential) field energies. The steric and electrostatic fields generated automatically were scaled by the CoMFA standard method with default cut-off energy of 30.0 kcal/mol. For CoMSIA analysis, five similarity indices including steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields were calculated using a sp 3 carbon atom probe with radius 1.0 Å, charge +1.0, hydrophobicity +1, hydrogen bond donor +1, hydrogen bond acceptor +1, and attenuation factor α = .3. Individually, the steric field was related to the third power of the atomic radii, the electrostatic field was derived from atomic partial charges, the hydrophobic field was originated from atom-based parameters, and the hydrogen bond donor and acceptor fields were generated from a rule-based method based on the experimental values. In addition, a Gaussian equation was used to evaluate the distance between the probe atom and each atoms of the molecule. The similarity indices, A F, K between the compounds of interest were computed by placing a probe atom at the intersections of the lattice points and using Equation (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 (steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor), W ik is the actual value of physicochemical property k of atom i, and W probe, k is the value of the probe atom.

3D-QSAR models generation and validation
Partial least squares analysis was employed to correlate the CoMFA and CoMSIA fields with the inhibitory activities (Cramer, Bunce, Patterson, & Frank, 1988), in which fields were used as independent variables and the TR-coactivator inhibitors as the target variables. The performance of the derived 3D-QSAR models was evaluated by leave-one-out (LOO) method in which one inhibitor was removed from the data-set, then its inhibitory activity was predicted using the model derived from the rest of the inhibitors. The cross-validation analysis producing the highest cross-validated correlated correlation coefficient R 2 cv , the optimum number of components (ONC), and the lowest standard error of prediction SEP was considered for further analysis. Then the generated ONC was applied to the non-cross-validation procedure, and the non-cross-validated correlation coefficient R 2 ncv , standard error of estimate (SEE), and F value were generated.
In order to predict the robustness and the predictive capability of the CoMFA and CoMSIA models, the inhibitory activities of the test set inhibitors were further predicted. The predictive correlation coefficient ðR 2 pred Þ was computed using Equation (2): where PRESS represents the sum of the squared deviations between the predicted and the actual pIC 50 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.  Table S1). Thus, in view of the statistical values, the CoMFA model is better than the CoMSIA one. Therefore, we mainly focus on the CoMFA model in the following discussion. The results based on the CoMFA model indicate that the derived model possesses good internal predictive power. The steric contribution dominates with a fraction of 57.8% over the electrostatic contribution with a fraction of 42.2%.

Results and discussion
Furthermore, the predictive ability of the optimum CoMFA model was determined from the test set inhibitors, the predictive correlation coefficient ðR 2 pred Þ of the CoMFA model is .5580. In addition, there is a good correlation between the actual and the predicted inhibitory activities of the training set and test set based on the best CoMFA model (Figure 3(A)), further suggesting that the CoMFA model is reliable and statistically significant.

For TRα
Based on the criteria that R 2 cv [ :50 is the prerequisite of constructing a reliably predictive model, the results of the CoMSIA models are abandoned. The CoMFA model descripting the steric and electrostatic fields plays significant roles in predicting the inhibitory activities of TRα. Statistical CoMFA model gives R 2 cv of .526 and R 2 ncv of .864 with a lower SEE of .187 and a higher F value of 40.291 (at three components). Corresponding field contributions of SE are 41.4 and 58.6%, respectively, suggesting that the electrostatic field is higher correlated with the inhibitory activity.
The CoMFA model is further validated by the test set that are not included in the training set with an R 2 pred of .6983. The experimental value vs. the predicted value for the training and the test set is shown in Figure 3(B), the predicted values are close to the experimental ones with the points located uniformly around the regression line. Thus, the derived model is able to predict successfully the inhibitory activities of inhibitors that are not used in constructing the model.

3D-QSAR contour maps
The coefficients from the CoMFA models are plotted to generate contour maps that display regions in 3D space around the inhibitors where changes in the respective physiochemical properties are predicted to increase or decrease the binding activity. Analysis of the contour maps with respect to the chemical environment of the TR receptors would reveal the quality of the derived models (Mao, Li, Hao, Zhang, & Ai, 2012). The CoMFA contour maps for TRβ are shown in Figure 4, and those for TRα are listed in Figure 5. To aid in visualization, compound 21 (for TRβ) and compound 24 (for TRα), as the most active TR inhibitors in the series are superimposed with the contour maps.

Steric contour map for TRβ
In the CoMFA steric contour map (Figure 4(A)), regions where steric bulk groups would increase the activity are represented by green polyhedrons, while regions where steric bulk groups would decrease the activity are represented by yellow polyhedrons. In the figure, we can note that the steric contour maps are mainly located at Region A, Region B, Region C, Region D, and the 3-positon of ring A (as shown in Figure 6(A)).
With respect to Region B, a steric map indicates a yellow area around this position, suggesting the importance of the presence of a minor group in this region for the inhibitory activity, thus, compound 15 with -S atom at this positon exhibits higher activity than compounds 17 (-CONH) and 18 (-NHCO) Take compounds 38-43 as an example, a green contour map around Region A shows bulky groups are favored of the inhibitory activity, this is in agreement with the experimental data and the order of the inhibitory activity is: 43 (seven carbon atom) > 42 (six carbon atom) > 40 (five carbon atom) > 39 (four carbon atom) > 38 (three carbon atom). However, a minor yellow contour map is adjacent to the green contour map, resulting in the activity of compound 44 with substituents touching the yellow map is lower than compound 43, which is coincide with the summary in the literature , the most active compounds were those bearing an n-heptyl or n-hexyl para substituent.
Considering the contours at Region C, a sterically unfavorable region (yellow) is presented at this positon, which would increase the inhibitory activity with minor substituents. For instance, compound 52 containing group -NH has lower activity than compound 53 which possesses oxygen atom at the same position.
Concerning the contour maps at the -CO group of Region D, the steric map presents a green region at this area, indicating a steric favorable region, as an example, compound 24 (IC 50 = 1.6) shows a substantial decrease in the inhibitory activity compared to compound 23 (IC 50 = 1.5). Possibly due to increased steric hindrance, replacing the group in compound 56 leads to a further decreased activity in compound 24 (-Cl). And the inhibitory activity of compound 32 ( ) is higher than compound 33 ( ). Furthermore, there are some steric unfavorable yellow contours covering the -NH group at Region D, suggesting that the size of this substituent should not be too bulky, otherwise, it might produce a steric restriction on the activity. With respect to the 3-position of ring A, some yellow contours are positioned, indicating that a large group in this area is disadvantageous to the inhibitory activity. This can be supported by the fact that compound 19 with a -Cl substituent shows less potency than its counterpart 21 with -H group. Another example is that compound 20 with a group -Cl towards this spatial distribution is less active than compound 22 with a group -H.

Electrostatic contour map for TRβ
Electrostatic contour maps are displayed in Figure 4(B), blue and red regions indicate electron-donating or electron-attracting groups would be favored or disfavored, respectively. A large red contour map (Figure 4(B)) at Region B indicates that electronegative groups at this positon would enhance the inhibitory activity. This can be demonstrated by comparing compound 14 (-C) with compounds 15 (-S) and 16 (-SO 2 ). In addition, compounds 19-22 with electronegative -SO 2 in this area exhibit higher potency than compounds 23-56 (with electropositive -C at this position). A blue contour map near the 2-position of ring A indicates that the presence of electropositive substituents is significant for increasing the inhibitory activity. However, in the present work, the most potent compound 21 possesses electronegative substituent -Cl located at this site, thus, modifications can be made at this position to improve the activity. There is a red electronegative region flanking the 5-position of ring A indicates that substitution of an electronegative group would enhance the inhibitory activity. For example, inhibitors 21 and 22, which have -Cl substituents exhibit higher activity than inhibitors 19 and 20 with -C atom at this position. Some red polyhedrons surrounding the substituents at Region C depict their favor for the electronegative substituents. The fact that all inhibitors contain electronegative groups, such as -CO and -NH, occupying this position prove the finding. Moreover, at Region D, some small red contour maps are located, indicating that negatively charged substituents are favored at this area. Compound 45 possessing group at this position exhibits higher activity than that of its counterpart compound 38 ( ). And this can also be illustrated by the case in which compound 51 possesses a larger inhibitory activity than compound 50 for the reason that ( ) occupies Region D for compound 51, while for compound 50, its ( ), this phenomenon is also consistent with the order of toxic activity described in reference (Hwang et al., 2009), compound 51 was nontoxic to HepG2 cells at high concentrations, while compound 50 showed moderate toxicity.

Steric contour map for TRα
Figure 5(A) shows the distribution of steric fields generated using the CoMFA model, the green and yellow contours indicate the favorable and unfavorable steric interactions, respectively. A green plot found at Region A ( Figure 6(B)) indicates the favorable region for the presence of bulky groups for the inhibitory activity, which is further supported by the order of the activity of 39 (four carbon atom) < 40 (five carbon atom) < 42 (six carbon atom) < 43 (seven carbon atom) < 44 (eight carbon atom). Some yellow plots are located in the proximity of substituents at Region B, indicating that the steric bulk groups are detrimental to the inhibitory activity. This observation is consistent with the order of the inhibitory activity in compounds 40 and 41, such that compound 40 (IC 50 = 60.9) with small group is more potent than bulky substituted compound 41 (IC 50 = 69.6). Another yellow contour map occupying Region C suggests a favorable region for the presence of minor groups for the activity, for example, compound 52 (-N) with an increased bulk at this position has lower activity than compound 53 (-O). In addition, some yellow contours around Region D indicate that a small group at this position might be favored, this finding explains why compound with substituent (e.g. compound 32) is less potent than compound 29 with substituent as is the case with the activity order of compound 23 compound 24

Electrostatic contour map for TRα
The electrostatic contour map of the CoMFA model is shown in Figure 5(B). Blue contours indicate regions where electropositive groups would increase the activity, and red contours indicate regions where electronegative groups would enhance the activity. Two blue contour maps around Region B suggest that substitution with electropositive groups at this position would lead to more active inhibitors. Analysis reveals that most of the inhibitors possessing higher inhibitory activity have electropositive substituents in this site. Around the groups at Region C, a large and a small red contour maps are situated. This implies that substitution of electronegative groups at this region would enhance the activity. Compound 53 containing a more electronegative group near the red contours is more active than compound 54 ( ) . Moreover, at Region D, there is a large blue map, which indicates that electropositive groups would be favorable for the activity; however, most of the inhibitors applied in this work possess electronegative groups in this area, suggesting that modifications can be made at this position to improve the inhibitory activity.

Docking analysis and comparison with 3D contour maps
Molecular docking is performed to study the binding environment for the inhibitors interacting within the receptors TRβ and TRα ( Figure 7 describes the docking results of compound 21 in the binding pocket of TRβ, it is seen that some key amino acid residues (Thr277,Ile280,Thr281,Val284,Ile302,Leu305,Lys306,Cys309,Leu454,Glu457,and Val458) interact with the inhibitor at the binding site (Figure 7(B)). As shown in Figure 7(C), the inhibitor is anchored in the binding pocket via hydrogen bonds involving residues Leu305 and Val458. The -NH group at Region D is having hydrogen bond with the -CO group of Leu305 (-O···HN, 2.72 Å, 115.2°) (H-1). Another hydrogen bond is formed between the -CO at Region C and the side chains of Val458 (-O···HN, 2.99 Å, 101.5°) (H-2). This indicates the necessity of the hydrogen bond acceptor groups in the inhibitor for the inhibitory activity. Furthermore, the docked model also indicates that substituents at Region A almost extend outside the ligand binding pocket, suggesting that in this region the steric interaction might be favorable. This result is similar to the CoMFA green contour map at this position (Figure 4(A)). In addition, some amino acid residues Lys306, Leu454, and Val458 are adjacent to the groups at Region C, indicating that too large substituents are detrimental to the inhibitory activity. This result is consistent with the previous contour map (yellow contour at this position). Furthermore, the -CO group at Region D is fell into a large pocket, suggesting that bulky groups are beneficial to the ligand-receptor interactions, this is in agreement with the green contour map in this site. Docking results also demonstrate that the -NH group at Region D is hydrogen bonded with Leu305, and also enclosed by amino acids Cys309 and Val458, thus moderate substituents are preferred at this site, which is in accord with the yellow contour map in this area. Additionally, substituents at Region C adjoin some neutral amino acid residues, such as Leu454, Val458 and face the side chains of electropositive residue Lys306, which would form electrostatic effect with electronegative groups at inhibitors. This conclusion is also inferred from the red contour map of the CoMFA model (Figure 4(B)). A ligand binding pocket is formed among Leu305, Lys306, Cys309, and Val458, therefore, substituents with negative charges at Region D are favorable to the inhibitory activity, and this is in agreement with the CoMFA electrostatic contour maps.

Figures 8(A) and (B)
show the binding orientation of compound 24 in the binding site of TRα. The compound is secured by Thr227, Val230, Asp231, Lys234, Ile248, Leu251, Lys252, Leu400, Glu403, Val404, and Asp407, which are the main contributors to the interactions between the inhibitor and the receptor. The inhibitor is anchored in the binding pocket via hydrogen bond with amino acid residue Asp407 (Figure 8(C)), the -OH at the terminal of Region D acts as a donor to form a hydrogen bond with the side chain of Asp407 (-O···HO, 2.01 Å, 163.3°) (H-1).
It is also noticed that substituents at Region A expose outside the binding pocket, indicating that the size of the group at this position is unrestricted, this is also supported by the steric green contour map (Figure 5(A)). The sterically unfavorable yellow contour map surrounding the groups at Region C is placed in a small pocket composed by Ile248, Lys252, and Val404. Therefore, groups with too large chain at this position are unfavorable to the inhibitory activity, probably due to the crashing with these residues. Besides, at the terminal of Region D, there are some residues (such as Glu403 and Asp407) blocking the substituents in this position, which is supported by the yellow contour map in the CoMFA model.
It is observed that the binding pocket around Region C is surrounded by neutral amino acids Ile248, Val404, and electropositive residue Lys252, thus, electronegative substituents are favorable for improving the inhibitory activity, which is shown by the red contour map at the electrostatic contour maps (Figure 5(B)). In addition, the electropositive group favored blue contour map surrounding Region D is placed near electronegative residues Glu403, Asp407, and neutral residue Val404, which suggest that electropositive groups in this region would increase the binding activity.

Conclusions
In the present work, we have utilized 3D-QSAR method to explore the structural factors influencing the receptor binding activity of TR inhibitors. The CoMFA model performs excellent in terms of some evaluation criteria, such as R 2 cv and R 2 pred . Furthermore, the contour maps bring to light significant structural features that would be responsible for the different inhibitory activities. In addition, the binding modes responsible for the activity are investigated by molecular docking; the contour maps also correlate well with the specific interactions between the inhibitors and the receptors, TRβ and TRα. The obtained information helps us to better interpret the structure-activity relationship of the inhibitors and provides some crucial information for further optimization (Figures 9(A) and (B)): (1) for TRβ, bulky group at Region A, small and electronegative group at Region B, electropositive substituents at the 2-position of ring A, electronegative substituents at the 5-positon of ring A, minor groups at the 3-positon of ring A, electronegative groups at Region C, small and electronegative groups at -NH of Region D, bulky and electronegative substituents at -CO of Region D are likely to increase the inhibitory activity. And the key amino acid residues for the ligandreceptor interactions have been found, i.e. Leu305 and Val458, which form hydrogen bond interactions between the inhibitor and the receptor TRβ; (2) for TRα, a bulky group at Region A is beneficial for the activity; substituents with small and electropositive groups at Region B are favored; minor and electronegative groups at Region C are required for the increase in potency; and the introduction of small and electropositive substituents at Region D would enhance the inhibitor-receptor interactions; moreover, some significant amino acids such as Asp407 are identified.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2015.1124806.

Disclosure statement
No potential conflict of interest was reported by the authors.