Structure-based approach for the study of thyroid hormone receptor binding affinity and subtype selectivity

Thyroid hormone (TH) possesses the ability to lower cholesterol and improve cardiac performance, which have prompted the efforts to design analogs that can utilize the cholesterol-lowering property without adversely affecting heart function. In order to gain insights into the interaction mechanism for agonists at the active site of thyroid hormone receptor β (TRβ), quantitative structure–activity relationship (QSAR) models have been developed on TRβ agonists, significant statistical coefficients were obtained (CoMFA, R2cv, .732), (CoMSIA, R2cv, .853), indicating the internal consistency of the models, the obtained models were further validated using the test set, the acquired R2pred values .7054 and .7129 were in good agreement with the experimental results. The key amino acids affecting ligand binding were identified by molecular docking, and the detailed binding modes of the compounds with different activities were also determined. Furthermore, molecular dynamics (MD) simulations were conducted to assess the reliability of the derived models and the docking results. Moreover, TH exerts significant physiological effects through modulation of the two human thyroid hormone receptor subtypes. Because TRβ and TRα locate in different target cells, selective TR ligands would target specific tissues regulated by one receptor without affecting the other. Thus, the 3D information was analyzed to reveal the most relevant structural features involved in selectivity. The findings serve as the basis for further investigation into selective TRβ/TRα agonists.


Introduction
The thyroid hormone (TH), a major physiological regulator, is responsible for the regulation of gene transcription in intestinal, skeletal, cardiac muscles, liver, and the central nervous system, then, influence the overall metabolic rate, cholesterol and triglyceride levels, heart rate, and affect mood and overall sense of well-being (Chiellini, Nguyen, Yoshihara, & Scanlan, 2000;Lazar, 1993;Nguyen et al., 2002). TH is produced exclusively in the thyroid gland as thyroxine (T4) and 3, 5, 3-triiodo-Lthyronine (T3). T4 is the predominant secreted form in human and would be deiodinated by deiodinases to the major active form of the hormone (T3) in thyroid target tissues. Thyroid hormone T3 induces transcriptional responses primarily mediated through its thyroid hormone receptors (TRs), which belong to the superfamily of nuclear receptors (NRs) and would regulate the transcription of target genes positively or negatively in response to hormone binding, thus eliciting a multitude of physiologic effects on growth, metabolism, and homeostasis (Nguyen, Apriletti, Baxter, & Scanlan, 2005). Thus, TR has been targeted, with much success, as a point of intervention in a number of disease states.
There are two major subtypes of the TR, TRα and TRβ, expressed predominantly as the alternatively spliced variants. Mouse knockout studies (Johansson, Vennström, & Thorén, 1998), tissue distribution studies (Lazar, 1993), and studies on evaluation of patients with resistance to TH syndrome (Kahaly, Matthews, Mohr-Kahaly, Richards, & Chatterjee, 2002) have suggested that although the two isoforms are widely expressed, there are distinct patterns of expression with tissue and developmental stage. Therefore, TRα mainly locates in cardiac muscles and regulates cardiac functions, while TRβ predominates in liver and regulates cholesterol homeostasis and thyroid-stimulating hormone production (Schwartz, Strait, Ling, & Oppenheimer, 1992).
Similar to other NRs, TRs have three functional domains: (1) an N-terminal ligand independent transcription activation domain (AF-1); (2) DNA-binding domain; (3) a carboxy-terminal ligand-binding domain (LBD) containing a TH inducible coactivator binding site (AF-2) (Mangelsdorf et al., 1995). TRs normally act in concert with retinoid X receptor, then bind to short DNA sequences (TH response elements or TREs) located in the regulatory regions of target genes (Xu & Li, 2003), further regulating gene expression. The LBDs of TRs have been solved by X-ray crystallography (Darimont et al., 1998;Ribeiro et al., 1998;Wagner et al., 1995), which indicate that the TRs LBDs comprise a three-layer α-helical fold (a hydrophobic core), and the ligand is located in the center of the domain. When ligand is absent, TRs associate with corepressors such as NCoR (nuclear recptor corepressor) acting as strong constitutive repressors, which in turn bind to TREs in positively regulated genes (Privalsky, 2004). Agonist binding induces a conformational change in the receptor, which further causes corepressors release and coactivators complexes recruit in a sequential manner. The helix-12 in the LBD plays a central role in controlling the ligand-induced allosteric change, for instance, agonist binding induces helix-12 tight packing against the receptor, leading to the formation of a hydrophobic cleft accessible for coactivator binding (Weatherman, Fletterick, & Scanlan, 1999;Webb et al., 2002). Conversely, antagonist binding induces a conformational change that is inactive by preventing proper packing of helix-12 (Schapira, Raaka, Samuels, & Abagyan, 2000). Currently, although certain compounds have been discovered as a targeted therapy, it is restricted because of some TRs ligands trigger a mixture of favorable and deleterious actions, for instance, THs lead to weight loss by increasing metabolism and improving lipid balance, however, it also causes heart failure, muscle wasting, and osteoporosis (DeGroot, Larsen, & Hennemann, 1996;Greenspan, 2001;Utiger, 1995). Therefore, it is desirable to develop selective TRs agonists that would induce beneficial effects but do not have unfavorable effects. In addition, previous researches have suggested that TRβ specific agonists would have some desirable properties Wagner et al., 2001).
L-94901 as the first selective analog has been developed, it is approximately 50% as strongly bound to hepatic TRs as T3, while 1.3% as effective in binding with cardiac TRs (Underwood et al., 1986). Triac has an affinity for TRβ that is about three-times greater than T3 (Schueler, Schwartz, Strait, Mariash, & Oppenheimer, 1990). More recently, an analog, named as CGS23425 (Taylor, Stephan, Steele, & Wong, 1997), has a lower threshold for the activation of TRβ than TRα, the concentration required for half-maximal stimulation of TRβ is 2 × 10 −12 M and 10 −10 for TRα. Another compound, referred to as GC-1, showing 10 times preference for binding TRβ (Chiellini et al., 1998). According to the discovered ligands selective for TRβ, quantitative structure-activity relationship (QSAR) methods as an effective tool for drug design, have been established to aid in discovering more selective TRβ ligands, a 4D-6D QSAR model has been developed using Quasar and Raptor software to study the selective thyroid ligands (Vedani, Dobler, & Lill, 2006;Vedani, Zumstein, Lill, & Ernst, 2007). 2D-QSAR models have been constructed for both TRβ and TRα (Li et al., 2010;Ren, Liu, Li, Yao, & Liu, 2007). Additionally, three-dimensional quantitative structure-activity relationship (3D-QSAR) models applying molecular docking have been constructed on a series of agonists of TRβ (Du, Qin, Liu, & Yao, 2008). The structural requirements for more potent affinity would be identified by these models. However, 3D-QSAR models have not been reported on this class of TRβ agonists (Shiohara et al., 2012(Shiohara et al., , 2013. Moreover, in silico modeling approaches (Hao et al., 2011;Wang et al., 2011; such as 3D-QSAR models and molecular docking have been proven to be useful in further optimization and development of novel compounds. 3D-QSAR connects the steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor structural properties of ligands with the binding activity. Molecular docking can elucidate the interactions between the ligands and the receptor, and the results can be analyzed to support the conclusions of 3D-QSAR models (Kirkpatrick, 2004).
Clearly, it is not easy to predict how different ligand structures might impact its biological profile in vitro. Thus, we report the application of comparative molecular field analysis (CoMFA), comparative molecular similarity indices analysis (CoMSIA), molecular docking, and molecular dynamics (MD) simulations to this series of TRβ agonists. The derived models can well predict the binding affinity of the ligands in the data-set, and define the structural requirements for the binding of compounds to the receptor. Furthermore, molecular docking studies on TRβ-25 and TRα-25 have been performed and compared to address the TR receptor selectivity. In conclusion, the good concordance between 3D-QSAR, molecular docking, and MD results might help in identification of the key features of binding mechanism for the employed agonists and further provide some clues for the development of novel drugs possessing increased binding activity and higher selectivity.

Data-set and biological activity
The TRβ agonists together with their binding affinities were taken from the literatures (Shiohara et al., 2012(Shiohara et al., , 2013. The structures and binding affinities of the agonists were listed in Supplementary Table S1. The experimental determined binding affinities of the agonists are given as K i values, which were converted into pK i values. The data-set was divided into a training set and a test set, containing 24 and 9 agonists, respectively. Both the training set and the test set covering a complete range of binding activity and structural features were employed for the development and validation of the QSAR models, respectively.

Molecular modeling
For the agonists, the energy minimization and conformational search were performed by Tripos force field (Clark, Cramer, Opdenbosch, & Iii, 1989) with the Gasteiger-Hückel charges using the Powell conjugate gradient minimization algorithm. To obtain a stable conformation, the energy gradient convergence criterion and the maximum iterations were set to .05 kcal/mol Å and 1000, respectively.

Molecular alignment
Molecular alignment is an important step in the development of 3D-QSAR models (Ram, Daga, Rahul, & Javed, 2005). It is assumed that these agonists share a common structure, thus each agonist binds into the active site in a similar mode. In the present work, three different alignment rules were employed to develop the most reliable CoMFA and CoMSIA models. The first alignment rule is template ligand-based alignment (superimposition I), during the model generation, the most potent agonist (compound 25) was selected as a template, the remaining agonists were aligned to the common substructure of the template using the 'align database' function, which is shown in Figure 1(A) (marked in red), and the aligned result is depicted in Figure 1(B); the second alignment rule (superimposition II) is receptor-based alignment, the pro-cedure is similar to the superimposition I (Figure 1(C)), the conformation of each compound obtained from molecular docking was aligned together to compound 25. Superimposition III is also the docking-based alignment, the structure of each agonist with the highest total score was aligned automatically inside the ligand binding site of TRβ (Figure 1(D)).

CoMFA and CoMSIA analyses
CoMFA employs the active conformation and superposition rule for the agonists to provide quantitative correlation of the binding activity with the steric and electrostatic fields (Cramer, Bunce, Patterson, & Frank, 1988). The aligned agonists were placed in a 3D cubic lattice with a grid spacing of 2 Å in x, y, and z directions. Then, a sp 3 carbon probe atom with a van der Waals radius of 1.52 Å and a charge of +1.0 was placed at every lattice point to measure the steric (Lennard-Jones potential) and electrostatic (Coulomb potential) fields, and the CoMFA fields generated automatically were scaled by the CoMFA-STD method.
CoMSIA method had several advantages over CoMFA such as greater robustness regarding both region shifts and small shifts within the alignments (Pirhadi & Ghasemi, 2010). In addition to the steric (S) and electrostatic field (E), CoMSIA also calculates hydrophobic (H), hydrogen bond donor (D), and hydrogen bond acceptor (A) fields. The five fields were calculated by the standard settings: a probe atom with radius 1.0 Å, charge +1.0, hydrophobicity +1, hydrogen bond donor +1, and hydrogen bond acceptor +1, attenuation factor α .3. The distance between the probe atom and each molecule atom was evaluated by a Gaussian method. Generally, CoMSIA 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 the five physicochemical properties (S, E, H, D and A), 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 statistical validation
Partial least squares (PLS) analysis (Cramer et al., 1988) was used to derive the 3D-QSAR models. Firstly, the predictive values of the 3D-QSAR models were evaluated using leave-one-out method, in which one compound was removed from the data-set and its activity was predicted using the model derived from the rest of the data-set. The optimum number of components produced by PLS was employed to obtain the highest crossvalidated correlated correlation coefficient R 2 cv and the lowest standard error of prediction SEP. Secondly, the non-cross-validated correlation coefficient R 2 ncv was generated by the non-cross-validated analysis. Additionally, the CoMFA and CoMSIA results were graphically represented by contour maps, where the coefficients were generated using the field type 'StDev*Coeff'. In order to assess the robustness and statistical significance of the CoMFA and CoMSIA models, the test set which were not included in the training set were further predicted. The predictive correlation coefficient (r pred 2 ) was defined as follows: where PRESS represents the sum of the squared deviations between predicted and 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.

Molecular docking
Molecular docking simulations of the agonists into the TRβ binding site were performed by AutoDock (version 4) (Morris et al., 1998), the crystal structure of TRβ (PDB entry code 1NAX) was extracted from Protein Data Bank (http://www.rscb.org/pdb). The binding mode of all agonists was obtained by a grid-based docking program. The interaction energy between ligand and receptor was evaluated using atom affinity potentials calculated on a grid similar to that described by Goodford (1985). All ligands and water molecules were removed from the crystal structure, polar hydrogens and united atom Kollman charges were assigned for the receptor. During the docking process, the relative ligand (compound IH5) in the crystal structure 1NAX was used as a standard docking model. For molecular docking, AutoGrid was used for the preparation of the grid map employing a grid box with a docking box of 60 × 60 × 60 Å, the box spacing was set to .375 Å. And the grid center was designated at dimensions (xyz): 4.457 × 20.909 × 32.090. Finally, AutoDock was run with a maximum number of 10,000 retries and 27,000 generations, respectively. The docking possibilities were calculated by the genetic algorithm with local search (GALS). For each agonist, the simulation was composed of 100 docking runs using the standard Auto-Dock parameters. During the process of docking, some residues Arg282, Arg320, Asn331, and His435 were defined as flexible residues. Then the conformations with the lowest binding energy were extracted and aligned together for further 3D-QSAR analysis.

Calculation and selection of dragon descriptors
In the development of various QSAR models, dragon descriptors have been successfully used for quantitatively representing structural and physicochemical features of a molecule . In general, molecular descriptors can be collected from different sources such as substituent constants, physicochemical properties, quantum chemical calculations, and the theoretical structural parameters derived from one or two dimensional molecular structures (Todeschini & Consonni, 2000). In the present work, Dragon Professional, version 5.0 software (http://www.talete.mi.it/in dex.htm) was employed. It could calculate 1644 molecular descriptors that are divided into 20 groups: 48 constitutional descriptors, 119 topological descriptors, 47 walk and path counts, 33 connectivity indices, 47 information indices, 96 2D autocorrelations, 107 edge adjacency indices, 64 Burden eigenvalue descriptors, 21 topological charge indices, 44 eigenvalue-based indices, 41 Randic molecular profiles, 74 geometrical descriptors, 150 RDF descriptors, 160 3D-MoRSE descriptors, 99 WHIM descriptors, 197 GETAWAY descriptors, 154 functional group counts, 120 atom-centered fragments, 14 charge descriptors, and 29 molecular properties. The number of descriptors is subsequently reduced by eliminating those descriptors that are obviously redundant or unrelated to the biological properties. Then, a stepwise linear regression method as the variable selection in R software (www.r-project.org) was employed to reduce the descriptor space (Wang, Li, Ding, Jiang, & Chang, 2008). The WHIM descriptors L2p and the 2D autocorrelations MATS 4p showed significant influences with R 2 = .65. By employing the dragon descriptors, the robustness and generalization of the system can be enhanced.

MD simulations
The classical MD studies were carried on the most potent ligands. The missing atoms and residues (TRβ: Gly209-His210; TRα: Gln181-Trp185, Pro200-Met204) were built rationally by Modeller 9v11 (Eswar et al., 2007). The PDB2PQR web server was then used to model the corrected protonation states at physiological pH 7.0 (Dolinsky, Nielsen, McCammon, & Baker, 2004). The leap in Amber Tools 13 was employed to generate the input files for the simulations. The protein parts were described by FF12SB, for the ligand 25, the RESP charges (Bayly, Cieplak, Cornell, & Kollman, 1993) were assigned by the same Gaussian Base (HF/6-31G) where amber force field was generated from. General Amber Force Field (GAFF) (Wang, Wolf, Caldwell, Kollman, & Case, 2004) was then used to define the bond constants of ligand 25. Sodium ions were added to neutralize the systems. In the end, the systems were saturated in a TIP3P water rectangular box with 12.0 Angstrom to each facade.
All-atom MD simulations were performed using NAMD 2.9 (Phillips et al., 2005). The following parameters, minimization, and thermalization procedures were adopted for both of the MD systems. Equilibration was performed in three stages: firstly, the potential steric clashes in the initial configuration were relieved with 50,000 steps of energy minimization. Initial velocities for each system were assigned randomly from a Maxwell-Boltzmann distribution. Each system was then gradually heated up to 300 K over .1 ns, with the protein harmonically constrained, under the canonical ensemble (NVT) conditions. Following this, each system was simulated for another .1 ns under the isothermal-isobaric ensemble (NPT) conditions with applied constrains. Thereafter, the equilibrated systems were subjected to a 10 ns MD simulations.
The time step of all classical simulation was set to 2 fs, and the nonbonded cut-off length was set at 1 nm. All simulations were run at constant temperature (300 K) and pressure (1 atm), using Langevin damping coefficient of .5 fs −1 . For each simulated system, periodic boundary conditions together with the Particle-Mesh Ewald method were employed to calculate the electrostatics interactions (Phys, 1995). For each system, three different duplicates were used, differ by their distribution of the initial velocities seeds. The conformations of each system were sampled every 10 ps for subsequent analysis.

Results and discussion
3.1. QSAR model To examine the structural features contributing to the binging activity of the agonists, three different superimpositions were quantitatively analyzed using CoMFA and CoMSIA methods. In the present work, all 31 possible combinations of the descriptors were applied to generate different QSAR models, with a full summation represented in Supplementary Tables S2-S4. The best results derived from the CoMFA and CoMSIA models are summarized in Table 1, which suggest that the template ligand-based 3D-QSAR (superimposition I) models exhibit a better statistical result. Thus, the models based on superimposition I will be mainly discussed in the following sections. In addition, the actual and predicted pK i of these compounds are shown in Supplementary Tables S5 and S6. The CoMFA model gives a relatively high R 2 cv value of .732 and R 2 ncv value of .973 with six components, it also exhibits an F value of 102.781 and an SEE value of .161, the high R 2 cv , R 2 ncv , and F values suggest that the derived model is reasonable and possesses a good predictive ability. The contributions from the steric and electrostatic fields are 30.2% and 37.4%, respectively. Additionally, MATS 4P and L 2P descriptors (Julio, Miguel, & Michael, 2006;Nasser, Mohammad, Araujo, & GalvãO, 2009) also make 19.0% and 13.4% contributions to the model, effectively enhancing the fitting degree of the model. MATS 4P is a 2D autocorrelations descriptor, the structural variables introduced by Moran correspond to dimensional autocorrelations between pairs of atoms in the agonist and are also defined to quantify the contribution of a considered atomic property to the analyzed property. Additionally, L 2P belongs to the WHIM descriptors, which are molecular descriptors based on statistical indices calculated on the projection of the atoms along principal axes. And L 2P is built in such a way as to capture relevant molecular 3D information regarding atom polarizabilities with respect to invariant reference frames. The molecular descriptor types and definition of the symbols are shown in Supplementary Table S7.
The CoMSIA model with S, H, D, and A fields gives a cross-validated R 2 cv of .853 with an optimized component of six, which suggests that this model is statistically significant and robust. The non-cross-validated R 2 ncv = .982, F = 154.502 and SEE = .132. The relative field contributions of S, H, D and A field are 6.2%, 13.1%, 28.3%, and 18.3%, respectively. The two descriptors MATS 4P and L 2P exhibit 21.4% and 17.3% contribution rate, respectively, to the CoMSIA model, which indicate that these descriptors are useful for enhancing the reliability of the model.
External validation can be further conducted to estimate the robustness and forecast power of the QSAR models. In the present work, an external test set including nine agonists were employed to validate the derived models. The predictive correlation coefficient R 2 pred values for the test set from CoMFA and CoMSIA models are .7504 and .7129, respectively. The actual vs. predicted pK i values of the training set and the test set are illustrated in Supplementary Figure S1, all the points are rather uniformly distributed around the regression line, indicating the satisfactory predictive capability of the models.

3D-QSAR contour maps
In order to further explore the interaction feature of an agonist with its receptor, the steric and electrostatic contour maps of the CoMFA model are shown in Figure 2. The steric, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor contour maps of the CoMSIA model are shown in Figure 3. The contour maps show regions where differences in molecular fields are associated with differences in the binding activity (Mao, Li, Hao, Zhang, & Ai, 2012). To aid in visualization, the most potent agonist 25 is applied in the Figures.  194 .122 .193 Notes: R 2 cv = cross-validated correlation coefficient using the leave-one-out methods,

CoMFA contour maps
Considering the CoMFA steric contour map first (Figure 2(A)), the green contours show favorable steric interaction and the yellow contours show the region where the steric group is not favored. The steric-favorable region is in the direction of the R 4 group (shown in Supplementary Figure S2), suggesting that bulky groups in this position are crucial for binding activity, for instance, compounds 2-7 have an order for the activity of 2 < 3 < 4 < 5 < 6 < 7, with the corresponding substituents H, Et, iPr, , , , respectively. Two small steric-unfavorable contour maps are distributed around the green region, indicating that a minor group extended to this contour map is advantageous to the activity, compound 12 with 3′-H extending into the small yellow contours exhibits better activity than compound 20 (3′-Cl) and compound 23 (3′-OH). Therefore, at this site, a careful selection of the substituent length is required. A small yellow contour around the R 1 substituent suggests that a small group at this position would benefit the activity, for compound 4, which possesses -Me group shows less potency than compound 29 which has a relative minor group -Br at this site. Around the R 2 and R 3 substituents, it is noted that there is a yellow colored map, indicating that a large group in this position is disadvantageous to the binding activity, this can be supported by the fact that compound 30 with a larger substituent -CH 2 (CH 2 ) 2 CH 2 toward this contour distribution is less potent than compound 29 with a relatively small substituent -CH 2 CH 2 CH 2 . It should be noted that a large yellow contour is located around R 5 group, showing that bulky substituent is detrimental to the binding activity. By comparing compounds 31 (-NHCO) and 32 (-CONHCH 2 ) as well as compound 33 (-NHSO 2 CH 2 ), it is found that the presence of a larger substituent is harmful to the activity. As depicted in Figure 2(A), compound 25 (with highest activity) orients its -OH group to a yellow contour, indicating that its preference for the moderate and less crowed substituent in this position.
In the electrostatic contour map (Figure 2(B)), the red polyhedra represents areas where negative charge is predicted to an increase in the ligand binding activity, whereas the blue polyhedra represents areas where positive charge or groups with lesser electron density may lead to an enhancement in binding activity. The blue contour map around R 4 group reflects a common placement of this function within the highly active ligands, and all of the agonists employed in the present work possess an electropositive group at this area. Two red contour maps around the blue polyhedra discussed above indicate that electronegative groups would increase the binding potency. For example, compound 13 with an electron-donating -F group extended into the red contour map has higher activity than compound 12 without substituent at this position. Other red contour maps located near the R 5 substituent explain the significance of the electronegative group for potent activity. Most of the agonists from the entire series possess electron-donating groups in this area, in addition, compound 11 with -NHCOCH 2 COOH substituent is more potent than compound 9 with a corresponding -CH 2 CH(NH 2 )COOH group at this position. A large blue colored favorable region is presented around the red contour maps near R 5 , which indicates that the electropositive groups are favorable, in case of compound 33, the presence of an electronegative group -SO 2 at this position corresponded with increments in the binding activity, because the introduced group -SO 2 changes the orientation and the conformation of the agonist, and the related substituents are made toward the red contour map.

CoMSIA contour maps
Similar results are obtained for the steric contour of CoMSIA (shown in Figure 3(A)) as those of the CoMFA model, with a large green contour map located around R 4 , a yellow contour around the green contour map, a yellow contour located near R 2 and R 3 and a big yellow contour around R 5 . However, small variation can also be discovered in the CoMSIA contour map, the yellow contours appearing near R 1 have disappeared, which demonstrates its property is less significant.
The hydrophobic contour map of the CoMSIA model is displayed in Figure 3(B). The yellow and white isopleths indicate that hydrophobic and hydrophilic regions are favorable for the binding potency. A yellow contour around the R 4 position suggests that a hydrophobic group at this position may be favored. This is consistent with the fact that the binding activity of compound 13 bearing a hydrophobic group -F is higher than compound 12 (without hydrophobic groups at this position). Moreover, another yellow map is situated near R 1 substituent, which indicates that hydrophobic groups would be favorable for the activity, for example, most of the agonists applied in this work possess hydrophobic groups (-Br, -Cl and -Me) in this area, which would create a hydrophobic environment near this contour. The white contour near the -OH of ring A (shown in Supplementary Figure S2) indicates that a hydrophilic substituent at this position would benefit the binding activity, however, all compounds involved in the present study possess a hydrophobic substituent -OH here, thus, modifications can be made at this position. The presence of a white contour map around R 2 substituent suggests that its occupancy by hydrophilic groups would favor the activity, in fact, the most potent agonist 25 possesses hydrophobic group -Br around this contour map, and therefore, it is maybe better for the binding activity if the hydrophobic substituent is transformed to other hydrophilic groups. The three white contours are proximity to the R 5 substituent, suggesting that occupancy of these contours by hydrophilic groups would be favorable for the activity as exemplified by all of the compounds.
The CoMSIA hydrogen bond donor field (shown in Figure 3(C)) is represented by cyan and purple contours, in which cyan contours represent regions where hydrogen bond donor group would increase the activity, while the purple contours represent areas where hydrogen bond donor group would decrease the activity. Two cyan contour maps near the hydroxyl group of ring A and the end of R 5 reveal that hydrogen bond donor groups are preferred in these areas. One cyan contour and a purple contour are found near R 4 , suggesting that a careful selection of the substituent in this site is required. This is in agreement with the experimental results that compound 11 with a hydrogen bond donor group (-NH) and a hydrogen bond acceptor atom (such as -N, -O) gives a higher activity than compound 12 without hydrogen bond donor and acceptor substituents. A cyan contour map at the -NH of R 5 represents the increased activity if a hydrogen bond donor group presented at this position. Compared compound 4 (-NH) with 32 (-CO), the activity discrepancies can be explained by this contour map, and this is in accordance with the observations that compound 9 in the data-set lacks hydrogen bond donor group, leading to its decreased activity. One purple contour map surrounding the cyan contour (located at the -NH of R 5 ) also reveals that a hydrogen bond acceptor group would improve the activity. In compound 4, the -CO group touching this contour map (serving as hydrogen bond acceptor) is favorable for the activity and its replacement would lead to lower activity, such as compound 32 with -NH extended to this contour. And this is in accord with compound 9 possessing hydrogen bond donor group -NH 2 at this contour map, which is unfavorable for the activity.
In Figure 3(D), only a red and a magenta contour map presented suggest that the hydrogen bond acceptor field is not very important for the potency. The magenta contour explains the favorable hydrogen bond acceptor field while red explains the unfavorable hydrogen bond acceptor field. A magenta contour map at the end of R 5 reveals that a hydrogen bond acceptor substituent is responsible for an increase in binding activity. For example, compound 4 with hydrogen bond acceptor -CO depicts a larger potency than its counterparts 32 with hydrogen bond donor group -NH in this position. In addition, a red contour at the -NH of R 5 depicts the unfavorable interactions with the hydrogen bond acceptor groups. This could be explained by the order of the activity of compounds 4 and 32, compound 4 (-NH) > compound 32 (-CO).

Docking analysis and comparison with 3D contour maps
To validate the reliability of the generated contour maps, molecular docking was applied for comparison. Molecular docking is a method which helps us to understand the binding mode of the ligand-receptor complexes and determine the amino acid residues involved in the recognition of the agonists (Bhadoriya, Sharma, Jain, Raut, & Rananaware, 2013). All agonists in the data-set were docked into the binding pocket of TRβ. Figure 4 displays the conformation of the most potent agonist 25 in the active site of TRβ.
The active site consists mainly of hydrophobic residues (Phe269, Phe272, Ile275, Ile276, Ala279, Met310, Met313, Ala317, Thr329, Leu330, Leu346, Met442, Phe455 and Phe459) and electropositive residues (Arg282, Arg316 and Arg320) (Figure 4(A)). As shown in Figure 4(B), the ligand is anchored in the ligand binding pocket via several hydrogen bonds involving residues Met313, Arg320, Thr329, Asn331, and His435. The -OH of ringA acts as a hydrogen bond donor to form a hydrogen bond with the nitrogen atom of His435 (-N···HO, 1.80 Å, 139.2°) (H-1), the presence of a cyan contour map at this position denotes the areas where hydrogen bond acceptors might exist in the receptor, and thus a hydrogen bond donating group would be favored; the -NH of R 5 interacts through hydrogen bond formation with the backbone -CO of Met313 (-O···HN, 2.62 Å, 147°) (H-2), which is consistent with the cyan and red plots in the hydrogen bond donor and hydrogen bond acceptor contour maps, respectively, indicating that hydrogen bond donor groups are favorable; the -CO of R 5 forms a hydrogen bond with the -NH of Thr329 (-O···HN, 2.39 Å, 133.3°) (H-3); The -OH located at the end of R 5 acts as a donor to form hydrogen bond with the side chain of Thr329 (-O···HO, 2.03 Å, 146.8°) (H-4), therefore, a hydrogen bond donor substituent is needed in this position, evidenced by the presence of cyan contours around this area by the CoMSIA model. In addition, the -CO at the end of R 5 forms hydrogen bonds with -NH of Arg320 (-O···HN, 2.25 Å, 139.3°) (H-5) and (-O···HN, 1.79 Å, 155.9°) (H-6), the hydrogen bond interactions in this area may favor the activity, which is shown by the magenta contour map (Figure 3(D)).
The docked model also exhibits a comparatively large space P1 (shown in Figure 5(A)) around the R 4 substituent, indicating that steric interaction is favorable for the activity in this area, this is similar to the CoMFA and CoMSIA contour maps discussed above (Figures 2(A) and 3(A)). According to the environment in the ligand binding site, the substituent R 5 fits into a small cavity P2 suggesting that too large groups may introduce severe steric clash with the amino acid residues, which is inferred from the steric yellow region of the CoMFA and CoMSIA models. The space at the -OH of ring A is limited by forming hydrogen bond with His435 in the active pocket of TRβ. Therefore, bulky groups in this place are disadvantageous to the binding activity.
According to the environment in the agonist binding site shown in Figures 4(A) and 5(A), the R 4 substituent adjoins neutral amino acids and forms electrostatic interactions with amino acid residue Gly334, indicating that the increased positive charge of substituent might improve the activity, evidenced by the presence of a blue plot around the R 4 group by the CoMFA model. In addition, R 5 substituent is fitted into a binding pocket composed by electropositive and electro neutral residues (Ala279, Ala317, Met313, Thr329, Arg316, Arg282 and Arg320), which means that electronegative groups in this position can result in an increase in the activity, this is fully supported by the CoMFA electrostatic red contour map.
The R 1 hydrophobic substituent can bind to a hydrophobic cavity, consisting of residues Leu330, Met310, and Leu346 (Figures 4(A) and 5(B)), which can be inferred from the hydrophobic yellow maps of the CoMSIA model. In addition, some amino acids are found around the R 5 group in the molecule, such as Ala279, Arg282, Met313, Arg316, Arg320 and Thr329, indicating that hydrophilic groups are beneficial to the activity, this is consistent with the CoMSIA white contour map. The -OH of ring A faces the imidazole group of His435, thus, a hydrophilic substituent is needed, which is evidenced by the presence of a white contour map in this area.
The comparison of the docking results and the contour maps indicates that the derived QSAR models are reasonable and can offer constructive suggestions for further modifications.

Molecular docking simulations of different compounds-bound TRβ
In this section, the molecular docking results of some agonists in the TRβ active site will be discussed. Compound 25 (the most potent active compound), compound 2 (the lowest active compound), and compound 1 (T3) are buried in the binding site of the receptor (shown in Figure 6). The three compounds (Supplementary  Table S1) are structurally similar, thus it is considered that they might bind to the same binding site and possess similar orientation. However, the comparison of the conformations of the three ligand-receptor complexes shows a change in R 5 substituent. The R 5 substituent of compound 1 translates almost 60°from the location of compound 2, and compound 25 rotates in a clockwise direction up to 30°compared to compound 2 in the presence of its target. In addition, compound 1 shows an independent movement into the inside of the P2 pocket, which is consistent with the structure replacement of (compounds 2 and 25) with (compound 1) and the contour map presenting a yellow contour in this position. The detailed chemical structures and binding modes of the three compounds are displayed in Figure 4(B), Supplementary Figures S3 and S4. Similar to compound 25-TRβ, compound 1 and compound 2 are secured by residues Phe272, Ile275, Ile276, Ala279, Arg282, Met310, Met313, Ser314, Arg316, Arg317, Arg320, Thr329, Leu330, Asn331, Leu346, His435, Met442, Phe455, and Phe459. It can be seen that compound 2 and compound 1 are located at the active site of the receptor by forming some hydrogen bonds. Comparing TRβ-2 with TRβ-25, the hydrogen bond interactions are almost conserved. However, the number and the length of the hydrogen bonds have changed, for example, the distance between the -OH group of ring A and His435, the -NH of R 5 and Met313, the -COOH of R 5 and the backbone of Asn331 and Arg320 has changed from 1.80 to 1.87 Å, 2.62 to 1.96 Å, 2.39 to 1.72 Å, 1.79 to 1.87 Å, and 2.25 to 2.10 Å, in addition, a hydrogen bond is formed based on TRβ-25 between the -COOH of R 5 and the -CO of Thr329, which is yet not presented in TRβ-2, it is interesting to note that a cyan contour map is situated at this position (shown in Figure 3(C)), accounting for the lower binding activity of compound 2. Moreover, the analysis of the binding mode of compound 1 suggests that the presence of hydrogen bonds between the ligand and TRβ with lower frequencies, the main findings are listed as follows: (1) the hydrogen bonds (the -OH of ring A···His435, the -COOH at the terminal of R 5 ···NH 2 of Arg320, and the -COOH at the terminal of R 5 ···CO of Thr329) are preserved in TRβ-1 compared with TRβ-25; (2) the length and angle of the conserved hydrogen bonds have altered; (3) Comparing TRβ-1 with TRβ-25, the ligand-receptor Figure 6. Structural superposition of TRβ-1, TRβ-2, and TRβ-25. Note: The projection highlights the structure of the active site with compound 1 (yellow), 2 (cyan) and 25 (blue), which are displayed in sticks. complex possesses less hydrogen bonds; (4) the hydrogen bonds between compound 25 and Met313 and Asn331 have vanished, this is due to the fact that the movement of the R 5 substituent in compound 1 makes it away from the hydrogen bonds involved amino acids. By comparing compounds 25, 1 and 2, it can be concluded that the replacement of the isopropyl at R 4 group by smaller iodine and hydrogen groups would affect the binding activity, which is evidenced by the CoMFA green contour maps and the docking results (with large binding pocket at this position). Compound 2 possessesing a methyl group at R 1 group shows less potent activity than compound 25 (Br), suggesting that bulkier groups at this position will clash with the residues, such as Leu341, Leu346, and Leu353, which is also corroborated by the steric contour maps.

Molecular docking simulations of subtype selectivity for TRβ/TRα receptors
In order to identify the structural divergence of the active sites between TRβ and TRα, the docking model of TRβ-25 is superimposed on the structure of TRα-25, as shown in Figures 7(A) and 8. The regions surrounded by different amino acid residues in the TR subtypes (TRβ Asn331/Ser277 TRα) are significant to discriminate the two complexes. This is the only instance wherein the residue that is altered exposes its side chain to the P2 pocket, contributing to the differences in selectivity of the two receptors towards various agonists. It is evidenced that the R 5 substituent for TRα-25 turns~90°relative to that of TRβ-25, in mechanistic terms, we favor the view that the conformational discrepancy in R 5 is triggered by the favorable substrate-binding environment, such as TRα Arg266/TRβ Arg320 (Arg266 is far away from compound 25). Therefore, the TRβ-selective electrostatic interactions could be established, this is consistent with the red contour maps in this position. Although the amino acid residues in the ligand binding site are conserved between the two TR subtypes, arene-cation interactions are formed between ring A and residue Leu292 in TRα, which is missing in TRβ-25. Additionally, the possible hydrogen bonds interacting model is depicted in Figure 7B. The -OH of ring A is hydrogen bonded with His381 (-N···HO, 2.15 Å, 162.3°) (H-1), the -NH of R 5 forms a hydrogen bond with the side chain of Met259 (-O···HN, 2.10 Å, 158.8°) (H-2), the -COOH at the end of R 5 is involved in hydrogen bonding with the backbone of Ala225 (-O···HO, 1.94 Å, 135.3°) (H-3). Compared with complex TRβ-25, less hydrogen bonds are presented, and the length and angle of the hydrogen bonds have been altered. This demonstrates that hydrogen bonding interactions are also quite essential for ligand binding and selectivity.

MD simulations
With the objective of further validating the 3D-QSAR models and molecular docking, MD simulations (a total of 10 ns) were conducted on the two docked structures, i.e. 1NAX + compound 25, 1NAV + compound 25.

For TRβ
To explore the dynamic stability of the complex, the root mean square deviation (RMSD) of the trajectory with respect to the docked structure was calculated, as shown in Figure 9(A), the RMSD of the ligand stabilizes at 1.5 Å after 5 ns. Furthermore, after 7.5 ns, the RMSD of the complex reaches about 1.8 Å and maintains this value throughout the simulation. This clearly indicates stable conformation for the docked structure during MD simulation, which can be further illustrated by the superimposed conformation of the average structure of the MD simulation and the initial structure, and the hydrogen bonding interactions. It can be seen that there is little difference between the average structure and the docked Compound 25 in the active site is displayed in sticks (yellow for TRβ, green for TRα), the residues are labeled in blue for TRβ, and purple for TRα. (B) The enlargement for the ligand in the binding site of TRα after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted red lines, and the nonpolar hydrogens were removed for clarity. model, as shown in Figure 9(B). And the RMSD of the ligand within 1 Å is lower than that of the RMSD of the complex (1.5 Å), indicating that the ligand binds tightly to TRβ. In addition, the major hydrogen bond interactions are conserved during MD simulation (Figure 9(C)), for example, the hydrogen bonds between the -OH at the para position of ring A···His435, the -NH of R 5 ···Met313, the -CO of R5···Asn331. However, the length and angle of these hydrogen bonds have altered. Additionally, the MD simulation also reveals that the  previous hydrogen bonds between compound 25 and amino acids Arg320 and Thr329 have changed to Ala317 and Arg316, it is witnessed by the conformational alteration occurred at the end of R 5 after MD simulation, illustrated in Figure 9(B). Although the conformation of compound 25 in the active site has undergone some changes, the binding pocket is still stable, leading credit to the reliability of the docking model.
In addition, detailed analysis of root mean square fluctuation (RMSF) vs. the residue number for the complex is shown in Figure 9(D). It reveals that there are two major flexible regions corresponding to residues 247-262 and residues 447-450, which are mainly located at the solvent-exposed areas. On the contrary, residues around the active site, especially, residues Met313, Ala316, Ala317, and Asn331 show a rigid behavior. This analysis further indicates that the analyzed ligand can bind to the receptor TRβ steadily and strongly.

For TRα
A total of 10 ns of MD simulation was performed to explore the key residues controlling ligand binding. Firstly, the RMSD of compound 25 inside the TRα binding pocket vs. simulation time was analyzed ( Figure 10(A)). Interestingly, the RMSD of compound 25 remains at .8 Å during the first 3.5 ns, and it jumps down to .5 Å, then reaches up to .8 Å, furthermore, the jump is stabilized at .5 Å between 8.2 ns and 10 ns, which suggests that the ligand might switch between two different conformations in the binding pocket, and the difference between them is mainly located in the state of the R 5 substituent, as illustrated in Figure 10(B). Additionally, it is worth noting that the RMSD for all atoms fluctuates dramatically, in-depth analysis shows that the stronger structural flexibility of the complex is attributed to the existence of more flexible domain situated at residues 146-157, thus, we further calculated the RMSD values of residues 158-408, note that the lower values were obtained. This phenomenon can be explained by the structure and sequence superimposition of TRα to TRβ, showing that those residues situated at the N-terminal of TRα are independent of that of TRβ at the same position (shown in the upper left corner of Figure 10(A)).
Basically, there is no significant discrepancy between the structure extracted from MD simulation and the initial docked conformation (Figure 10(B)). Although several movements have occurred during MD simulation, the binding pocket and the conformation of compound 25 are stable. In addition, the stability of the complex is examined by the analysis of the hydrogen bonding interactions during MD simulation, the main findings are listed as follows: (1) the main amino acid residues and the position concerning hydrogen bonds in the ligand are conserved after MD simulation; (2) hydrogen bonds formed between the -OH of ring A and His381, the -NH of R 5 and Met259, the -COOH at the end of R 5 and Ala225 are conserved; (3) the length and angle of the conserved hydrogen bonds have changed; (4) the arene-cation interaction between ring A and Leu292 has vanished after MD simulation, conversely, this interaction is compensated by a novel hydrogen bond formed between the ligand and Ser277.
Furthermore, analyses of RMSF vs. the residue number are illustrated in Figure 10D. Although the fluctuations of residues 146-157 and 199-205 are higher, the RMSF values for the residues responsible for ligandreceptor interactions are lower than 1.5 Å. A closer examination of Figures 9(D) and 10(D) reveals that the fluctuations of the residues in TRα are obviously larger than those of the TRβ complex, which is consistent with the experimental activities.
In a word, the results of molecular docking and MD simulation all confirm the fact that the ligands studied in the present work are located properly inside the receptor binding pocket. Additionally, the stable RMSDs and the conformations of the ligands further confirm the rationality and validity of the molecular docking models.

Conclusions
In the present study, 3D-QSAR models, molecular docking, and MD simulations were carried out not only to construct highly predictive QSAR models, but also to understand the interaction mechanism between the agonists and TRβ. On the basis of the template ligand-based alignment, best 3D-QSAR CoMFA (R 2 cv , .732, R 2 pred .7054) and CoMSIA models (R 2 cv , .853, R 2 pred .7129) were developed, the internal cross-validation test, together with the good correlation between the experi-mental and the predicted activity values for an external set of compounds, proved the good predictive power of the models. Molecular docking reveals the probable bioactive conformations of the agonists and the detailed structures of TRβ with the ligands. In addition, the information derived from MD simulation also validates the reliability of the molecular docking procedure. The properties identified from the 3D contour maps correlate well with the specific amino acid residues recognized by docking analysis, and the main findings are summarized as follows ( Figure 11): (1) there should be a small, hydrophobic group at R 1 substituent; (2) the R 2 substituent should be small and hydrophilic; (3) there should be a small group at R 3 substituent; (4) the substituents at the R 4 position should be bulky, electropositive, hydrophobic, and hydrogen bond donor group; (5) around the contour map at R 4 , groups should be small and electronegative; (6) small, electronegative, hydrophilic groups should be located at the R 5 substituent; (7) small, hydrophilic, and hydrogen bond donor groups should be presented near the -OH of ring A. We also conducted molecular docking analysis on complexes TRβ-compound 25, TRβ-compound 1, and TRβ-compound 2 to highlight the structural preferences responsible for the binding activity differences.
Moreover, agonist selectivity is a bottleneck in contemporary ligand design, thus molecular docking analysis was employed to identify the binding modes of agonists at the binding site of TRβ and TRα, and gain insights into possible factors contributing to the subtype selectivity. The results indicate that the active site of TRβ is more accommodating the agonists applied in this study than does that of TRα, the difference in binding site preferences might contribute in part to the selectivity of these compounds. It is implied that an amino acid change, i.e. from Asn331 in TRβ to Ser277 in TRα in the active site would not be considered as a dramatic variation, however, would have a marked effect on the ligand binding. In addition, the R 5 substituent is adjacent to Arg320 in TRβ, but far away from Arg266 in TRα, which creates the possibility for hydrogen bonding interactions between ligand and TRβ. Taken together, the models, molecular docking analyses, and MD simulations could help further design of agonists with better activity, and provide some rationale for the selectivity in the binding effects of this series of agonists against TRβ and TRα.

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

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

Funding
This work was supported by the 12th five-year plan for science and technology development [grant number 2012BAD33B05].