Identification and Molecular Interaction Studies of Thyroid Hormone

Res. Toxicol. , 2016, 29 (8), pp 1345–1354. DOI: 10.1021/acs.chemrestox.6b00171. Publication Date (Web): July 13, 2016. Copyright © 2016 American C...
0 downloads 9 Views 1MB Size
Subscriber access provided by University of Newcastle, Australia

Article

Identification and Molecular Interaction studies of Thyroid Hormone Receptor Disruptors among Household Dust Contaminants Jin Zhang, Yaozong Li, Arun A Gupta, Kwangho Nam, and Patrik L. Andersson Chem. Res. Toxicol., Just Accepted Manuscript • DOI: 10.1021/acs.chemrestox.6b00171 • Publication Date (Web): 13 Jul 2016 Downloaded from http://pubs.acs.org on July 17, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Chemical Research in Toxicology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

Identification and Molecular Interaction studies of Thyroid Hormone Receptor Disruptors among Household Dust Contaminants Jin Zhang†,1, Yaozong Li†,1, Arun A. Gupta†, Kwangho Nam†,*, Patrik L. Andersson†,* †

Department of Chemistry, Umeå University, SE-901 87 Umeå, Sweden.

Corresponding Authors *Email: [email protected]. Tel: +46-90-786-5266 *Email: [email protected]. Tel: +46-90-786-6570 Author Contributions 1

J.Z. and Y.L. contributed equally to this work.

ACS Paragon Plus Environment

1

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

Abstract graphics

ACS Paragon Plus Environment

2

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

ABSTRACT

Thyroid hormone disrupting chemicals (THDCs), often found abundantly in the environment, interfere with normal thyroid hormone signaling and induce physiological malfunctions, possibly by affecting thyroid hormone receptors (THRs). Indoor dust ingestion is a significant human exposure route of THDCs, raising serious concerns for human health. Here we developed a virtual screening protocol based on an ensemble of X-ray crystallographic structures of human THRβ1 and the generalized Born solvation model to identify potential THDCs targeting the human THRβ1 isoform. The protocol was applied to virtually screen an inhouse indoor dust contaminant inventory, yielding 31 dust contaminants as potential THRβ1 binders. Five predicted binders and one negative control were tested using isothermal titration calorimetry, of which four, i.e., 2,4,5-trichlorophenoxyacetic acid (2,4,5-T), bisphenol A (3-chloro-2hydroxypropyl)

(2,3-dihydroxypropyl)

ether

(BADGE-HCl-H2O),

2,2',4,4'-

tetrahydroxybenzophenone (BP2), and 2,4-dichlorophenoxyacetic acid (2,4-D), were identified as THRβ1 binders with binding affinities ranging between 60 µM and 460 µM. Molecular dynamics (MD) simulations were employed to examine potential binding modes of these binders and provided a rationale for explaining their specific recognition by THRβ1. The combination of in vitro binding affinity measurements and MD simulations allowed identification of four new potential THR-targeting THDCs that have been found in household dust. We suggest using the developed structure-based virtual screening protocol to identify and prioritize testing of potential THDCs.

ACS Paragon Plus Environment

3

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

INTRODUCTION Thyroid hormone (TH) homeostasis is important for macronutrient metabolism, thermogenesis and energy balance, cardiovascular function, and normal brain development in humans and wildlife.1 A broad range of environmental xenobiotics have been reported to disturb normal TH functions by interfering with TH signaling, potentially causing malfunctions in blood circulation, brain development and reproduction.2 These chemicals, referred to collectively as thyroid hormone disrupting chemicals (THDCs), are routinely detected in our environment and some e.g., 2,4,6-tribromophenol, 2,3,4,5,6-pentachlorophenol, and polybrominated diphenylethers (PBDEs) have been found in household dust.3-5 Animal and human studies have indicated that THDCs may pose a threat to human health.6-8 Of particular concern is exposure to these chemicals during critical stages of fetal development, such as during neurological development, which may lead to permanent brain function defects, such as cretinism, permanent hearing losses and structural abnormalities in the brain.6, 7 Indoor dust is a major source for human exposure to THDCs.9 Ingestion of dust containing THDCs, such as bromophenols and perfluoroalkyl and polyfluoroalkyl substances (PFASs), has been associated with alteration of TH levels in humans.5 We have previously identified 485 organic contaminants with different chemical composition in indoor dust.10 These dust contaminants were originated mainly from consumer products or outdoor environments.11 Major groups include flame retardants (e.g., PBDEs), pesticides, plasticizers (e.g., phthalates), PFASs, and classical persistent organic pollutants (e.g., polychlorinated biphenyls (PCBs) and dioxins).9 Among them, PBDEs, PCBs, and PFASs have been reported to cause thyroid-related adverse effects, such as metabolic diseases, reduced fertility and thyroid-related tumors.1 For some compounds, molecular targets of the TH system have been identified. For example, PBDE 181

ACS Paragon Plus Environment

4

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

and 190 and certain PFASs, such as perfluorooctanoic acid (PFOA) and perfluorooctanesulfonic acid (PFOS), have been shown to target transthyretin (TTR) and disturb TH transportation.12-14 Thyroid hormone receptors (THRs), which are activated by triiodothyronine (T3) and thyroxine (T4) under normal physiological condition,15 are potential molecular targets for THDCs.16 For example, bisphenol A and certain hydroxylated PBDEs (OH-PBDEs) (e.g., 4-OHPBDE 69 and 4-OH-PBDE 121) have been shown to target THR and affect THR-mediated gene expressions.16-18 The importance of THRs in TH homeostasis is well-studied in the field of computer-aided drug discovery,19-22 however only a few in silico based studies are published on the interaction with THRs of environmental pollutants.23-25 Improved knowledge of interactions of contaminants with THRs is essential for identification of THDCs targeting THRs and human health related risk assessments of dust contaminants. Virtual screening (VS) by molecular docking is a high-throughput approach for predicting protein-compound binding poses and identifying and prioritizing bioactive compounds for experimental testing.26-29 VS has been successfully applied to identify potential environmental pollutants targeting, e.g., the estrogen receptor alpha and haloalkane dehalogenases, from large libraries of industrial chemicals.30-32 VS has the potential to reduce costs and animal testing in chemical risk assessment processes by identifying the most hazardous chemicals or the most promising candidates in drug development.33, 34 In the present study, we have developed a VS protocol to identify dust contaminants targeting human THRβ1. In the protocol, two additional approaches were introduced for robust screening of compounds with diverse scaffolds and functional groups. First, instead of relying on a single X-ray structure, an ensemble of protein structures, i.e., the same protein with slightly different ligand binding site (LBS) conformations,35-37 was used to mimic protein flexibility.38 Second, a

ACS Paragon Plus Environment

5

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

molecular mechanics based method, i.e., the generalized Born surface area (MM-GBSA) method,39 was used to for better estimation of solvation effects, thereby improve the performance of the VS protocol.40 The developed virtual screening model was then applied to screen dust contaminants targeting human THRβ1. THRβ1 was selected among the four possible isoforms of THR (i.e., THRα1, THRα2, THRβ1 and THRβ2),41 because of its importance to normal liver and kidney function, hearing development, and the feedback loop of the hypothalamic-pituitarythyroid axis.42 Disruption of THRβ1 by THDCs has also been suggested to be implicated in inadequate brain development and liver malfunctions.17, 43 Six representative dust contaminants (five predicted binders and one predicted non-binder) were selected for isothermal titration calorimetry (ITC) to determine their binding affinities to THRβ1. MD simulations were then carried out to examine the possible mechanism of THRβ1 ligand recognition and different dynamical behavior of THRβ1 upon binding of diverse binders. MATERIALS AND METHODS The study was conducted by following the procedure (Figure 1): I) development of the VS model, II) initial screening of 485 dust contaminants, yielding 131 initial hits, III) refinement of the initial hits, reducing the number of hits to 31 by ligand-protein interaction analysis, IV) selection of six representative compounds (five predicted binders and one predicted non-binder) based on their structural diversity and environmental relevance, and V) experimental measurement of their binding affinities to THRβ1 using ITC followed by MD simulations. Details of each of these five steps are given below. THRβ1 structure selection and preparation. Among 17 X-ray crystallographic structures of human THRβ1 in complex with various co-ligands available in the RCSB Protein Data Bank (PDB),44 seven structures were selected based on the structural diversity of their co-ligands and

ACS Paragon Plus Environment

6

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

shape/size of their binding sites of the protein. The PDB ID of each structure is given in Table S1. For each structure, any mutated residue was restored to its corresponding wild-type residue, followed by addition of hydrogen atoms, optimization of hydrogen bonding networks and assignment of the protonation state of each residue to its corresponding protonation state at pH 7.0. The pronated state of each protein residue in or near the LBS was further inspected visually and corrected if needed. The structures were then minimized using the OPLS2005 force field to remove unfavorable high energy contacts.45 The preparation steps were carried out using the Protein Preparation Wizard of the Schrödinger Suite.46 As water molecules can mediate interactions between protein residues and ligands by forming hydrogen bonds with them, their presence can improve virtual screening performance.47-49 However, they can also hinder the identification of ligands by blocking direct interactions with protein.50 In the present work, both possibilities were considered by either including (PDB ID with “-W”) or excluding explicit waters (PDB ID with “-nW”) in any X-ray structure with identified crystal waters in its ligand binding site (LBS). This yielded a total of 12 protein structures (Table S2), which were then used in the development of the VS protocol. Ligand datasets for virtual screening (VS). Two different datasets were used in the study: 1) the THR subset of the Database of Useful Decoys: Enhanced (DUD-E),51 and 2) an indoor dust contaminant inventory we compiled previously.10 The THR subset of the DUD-E comprising 103 THR binders and 7450 in silico generated inactive decoys was used for the development and evaluation of the VS protocol. The dust contaminant inventory contained 485 organic contaminants that humans could potentially be exposed to through household dust. These compounds were preprocessed using the ChemAxon JChem Standardizer.52 For each compound, between 1 and 32 conformers were generated depending on the total number of chiral centers

ACS Paragon Plus Environment

7

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

and their ionization and tautomeric states were determined at pH 7.0 ± 1.0, using the Schrödinger LigPrep53 and the Epik modules. Single structure docking, ensemble docking and MM-GBSA rescoring. The LBS of each THRβ1 structure was described by using a 10 Å cubic grid box centered on its corresponding coligand. The prepared compounds of the DUD-E THR subset were docked individually into each structure. For each conformer, a total of 5,000 poses were generated and ranked using the standard prevision (SP) scoring function of the Schrödinger Glide module.54 The top-ranked 400 poses were selected and subjected to 100 steps of conjugate gradient energy minimizations with a dielectric constant of 2.0. The results of each single structure docking were assessed by determining the area under curve (AUC) value of the receiver operating characteristic (ROC) curve and the enrichment factors (EFs) (Table S2). For each THRβ1 structure, the AUC value of the ROC curve was obtained by plotting the true positive rate (sensitivity) against the false positive rate (1-specficity) (Figure S1). We utilized AUC because its value is insensitive to the ratio of actives versus decoys in the DUD-E dataset.55 Enrichment was characterized by two EFs evaluated at 1% (EF1%) and 10% (EF10%) of the ranked database, respectively.56 Ensemble models were prepared by considering all possible combinations of the single structure docking results, from 2 structures up to all 12 structures. The top 20% of the DUD-E subset ranked by SP docking score was extracted for each single structure docking. Then, for a given combination of structures, the extracted results were combined to determine the overall performance, evaluated by the total number of active compounds identified by that particular combination. In Table S3, the results of the best performing ensemble model for each given number of protein structures are presented. The final ensemble model used for screening dust

ACS Paragon Plus Environment

8

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

contaminants was selected by considering both the number of active compounds identified and total computational cost of docking. To further improve the VS performance of the selected ensemble model, we applied the MMGBSA method and rescored the ligand-THRβ1 docking poses. The calculation was performed using the Schrödinger Prime module with the VSGB 2.0 implicit solvent model and the OPLS2005 force field.57 During the rescoring, the MM-GBSA energy minimizations were performed by allowing the conformation of ligand to relax while the rest of the complex was held fixed. The MM-GBSA scores were determined for the resultant optimized poses. The results of the MM-GBSA rescoring were also assessed by the number of identified active compounds, AUC values and enrichment factors, and compared with other docking results, as shown in Table S4. Compounds-THRβ1 interaction analysis. The majority of the amino acid residues constituting the LBS of THRβ1 are hydrophobic. In contrast, only a few residues are hydrophilic, some of which have been suggested to play critical roles in the specific recognition of endogenous ligands and pharmaceuticals targeting THRβ1.22 Since hydrophilic residues in LBS can mediate electrostatic interactions with ligands, we reasoned that these residues, if any, may contribute toward specific binding of dust contaminants to THRβ1. We analyzed these interactions for all seven selected X-ray structures and counted the number of H-bonds and salt bridges between each ligand and THRβ1 based on the poses proposed by the MM-GBSA rescoring. These data were used to refine the 131 initial hits identified by the MM-GBSA scoring, yielding 31 final hits (Table S5). Isothermal titration calorimetry (ITC). Binding affinities of six dust contaminants and T3 (positive control) to THRβ1 were measured by using ITC. The six tested dust contaminants

ACS Paragon Plus Environment

9

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

included five predicted binders and one predicted non-binder as a negative control (Table 1). Their purity, pKa values, and commercial sources are provided in Table S6. The human THRβ1 ligand binding domain (LBD) (6xHis plus Gly209-Asp461) was overexpressed in Escherichia coli BL21 (DE3) host cells by following the procedure outlined by Ye et al.22 The purified THRβ1 LBD was exhaustively dialyzed in phosphate buffer containing 20 mM sodium phosphate, 125 mM NaCl and 1 mM TCEP at pH 8.0. We prepared stock solutions of T3 (20 mM) and the six dust contaminants (100 mM) in dimethyl sulfoxide (DMSO). Before the ITC experiments, all stock solutions were diluted with phosphate buffer to a final level of 5% DMSO, as is a commonly acceptable DMSO concentration for ITC.58, 59 All experiments were conducted at 25°C by injecting aliquots of each compound (in syringe) into THRβ1-LBD buffer solution (in cell) using MicroCal Auto-iTC200 (Malvern Instruments Ltd, Worcester, UK). Initial “buffer scouting” was performed for each compound to avoid buffer mismatch caused by DMSO concentration differences between the cell and syringe. Raw data were collected, corrected for compound heats of dilution and integrated using MicroCal Origin. The ITC results (with “c” value < 5) were fitted to a single-site binding model and fixed stoichiometry parameter (N = 1), allowing estimation of the equilibrium dissociation constant (Kd) and enthalpy change (∆H).60, 61 Molecular dynamics (MD) simulations. The systems for MD simulations were prepared based on the 1Q4X THRβ1 LBD complex structure,62 which was the top performing single protein structure in docking (Table S2), and the poses from the VS protocol were used for the three dust contaminants (2,4,5-T, BADGE-HCl-H2O and BP2; abbreviations given in Table 1). The resulting complex (2,4,5-T system as an example) was solvated in an 85 Å rhombic dodecahedron water box containing 12475 TIP3P water molecules.63 The size of the water box was chosen to allow at least a 10 Å buffer space from each protein atom to the closest boundary

ACS Paragon Plus Environment

10

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

of the water box. 0.15 M NaCl was added to neutralize each system and mimic biophysical conditions. For example, the resulting system for 2,4,5-T contained 41494 atoms (3981 protein atoms, 18 ligand atoms, 37425 water atoms, 40 Na+ and 30 Cl- ions). In addition, one apo THRβ1 LBD system was prepared as a reference by removing its ligand from the X-ray structure. Each system was initially energy minimized for 10,000 steps under a series of position restraints and constraints, heated to 300 K and then equilibrated under the constant temperature and pressure conditions for 1 ns using the CHARMM program (version 38a1).64 50 ns production MD simulation was performed for each system using the NAMD program (version 2.9).65 The temperature of each simulated system was maintained at its target value by using the Langevin thermostat with 5 ps-1 friction coefficient. The pressure was controlled by the Nosé-Hoover Langevin piston method with 200 ps piston period and 100 ps piston decay time. A 2 fs integration time step was used and any bond involving hydrogens was constrained to their corresponding force field bond lengths. Non-bonded van der Waals energies were evaluated by smoothly turning off their interaction energies using a switching function at a distance between 9 Å and 11 Å, and electrostatic interactions were evaluated using the particle mesh Ewald summation (PME) method.66 In all simulations, the CHARMM36 force field with the CMAP correction was used to present the THRβ1 LBD,64,

67

TIP3P for water molecules, and the

CGenFF36 force field for the three dust contaminants.68 Details of the partial optimization of the dust contaminant force field parameters are provided in the Supporting Information. RESULTS AND DISCUSSION Ensemble docking model with MM-GBSA rescoring. The reliability of the VS using Glide was examined by comparing the predicted poses of actual co-ligands found in the 12 human THRβ1 X-ray structures to their corresponding X-ray structures.54 The results showed that the

ACS Paragon Plus Environment

11

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

docking simulation reproduced the actual binding poses with root-mean-squared-deviation (RMSD) lower than 2 Å for all compounds (Table S2).69 The DUD-E THR subset was then subjected to the Glide docking. The VS performance of (single structure) docking was compared between the 12 structures by examining the AUC value, total number of active compounds identified, and EF1% and EF10% values (Table S2) and 1Q4X-W showed the best performance with AUC of 0.863, good early enrichment (EF1% = 32.410 and EF10% = 7.180), and a total of 92 identified binders (out of 103 binders). The results showed that no single structure was able to identify all 103 binders with high AUC value. Thus, we combined different structures as an ensemble of structural models in order to identify a broader range of compounds with a high true positive versus false positive discrimination ratio. We examined all possible combinations of the ensemble structures from 2 structures up to all 12 structures and compared their performances based on the total number of identified binders (see Materials and Methods). Table S3 presents the best performing ensemble model for each given number of structures. The ensemble model with four structures (i.e., 1N46nW, 1Q4X-nW, 1R6G-nW, and 1XZX-nW) was chosen as our “final” ensemble model to identify dust contaminants targeting THRβ1. This model was chosen not only because it was one of the top models tested, successfully identifying 95 compounds out of the total 103 binders, but also because it showed reasonable computational cost, which is essential for high-throughput VS against large compound library. By comparison, the 1Q4X-W model, which is the top single structure model, only covered 75 binders when considering the top 20% hits (Table S4). Finally, the MM-GBSA method was applied to refine the selected ensemble docking model and its result, which yielded an AUC value that increased to 0.865 from 0.804 and EF10% to 7.478 from

ACS Paragon Plus Environment

12

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

6.109, respectively (Table S4). The ROC curve of the model clearly showed the improvement of EF10% and AUC after the MM-GBSA rescoring (Figure S1). Virtual screening of household dust contaminants. Identification of potential binders of THRβ1 was achieved in two steps (Figure 1). In the first step, the VS protocol described above was applied to a recently compiled inventory of 485 indoor dust contaminants.10 A total of 131 compounds (top 20% hits) were identified using the developed VS protocol. For each dust contaminant, all stereoisomers and protonation states, in the range of pH 7.0±1.0, were considered during the virtual screening. Each stereoisomer and protonation form was docked into the THRβ1 structures and the top scored stereoisomer/protonation form was selected for subsequent MM-GBSA calculation. For instance, BADGE-HCl-H2O has four stereoisomers. Among them, the RR stereoisomer has scored the top docking score (Table S7), and selected for the MM-GBSA rescoring and MD simulation. In selecting protonation form, we have calculated the energetic penalty for each protonation state based on its population estimated using the Epik module available in the Schrödinger program70 and the penalty was considered in the docking score. For example, triclosan has a pKa of 7.9. The compound can exist as both the protonated and deprotonated form at pH 7.0±1.0. We selected the protonated form for later MM-GBSA rescoring, as the form has a more negative docking score and lower energy penalty (Table S8). In the second step, ligand-THRβ1 interactions were analyzed for all ligand-THRβ1 complexes, including the 7 X-ray structures (Table S1), to identify protein residues that form hydrogen bonds and/or salt bridge interactions with the bound ligands (see Materials and Method for details). Our rational for this step was two-fold. First, since most endogenous THR binders are recognized via specific hydrophilic interactions, this additional step would aid in selecting chemicals that can form such specific interactions. Second, because hydrophobic residues

ACS Paragon Plus Environment

13

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

dominate the binding pocket of THRβ1, it is possible that our VS energy function underestimates hydrophilic interaction energies. This would then bias our selection toward certain chemicals that are more hydrophobic. With the aid of the second step, we intended to remedy the bias, thereby to improve selectivity between specific versus non-specific binders. Based on the analysis, we identified four hydrophilic residues, i.e. Arg282, Arg320, Asn331, and His435 in the binding pocket of THRβ1. We refined the 131 initial hits by analyzing their electrostatic interactions with these four residues and 31 dust contaminants displayed at least one electrostatic interaction with any of the hydrophilic residues (Table S5). The 31 potential binders can be classified into four major groups: 1) PFASs (10), 2) hydroxylated and/or halogenated aromatic compounds (14), 3) phenoxyacetic acids (3), and 4) miscellaneous (4) (Table S5). Interestingly, the present two-step screening protocol successfully identified known THR binders, such as bisphenol A (Ki = 200 µM),17 tetrabromobisphenol A (TBBPA, EC50 ≈ 100 µM),71 and four PFASs72 including perfluorooctane sulfonate (IC50 =13 µM)72, 73, perfluorononanoic acid (IC50 = 12 µM), perfluoroundecanoic acid (IC50 = 15 µM), and perfluoroheptanoic acid (IC50 = 9 µM), which indicated the selectivity of the present screening protocol. Measuring binding affinities of THDCs to THRβ1. Based on the VS results, six dust contaminants (five predicted binders and one predicted non-binder) (Table 1) were selected for measurement of binding affinities to THRβ1 by ITC. These represent different chemical structures, applications, and are environmentally relevant contaminants. The derived ITC thermograms (Figure 2) show clearly that four out of the five predicted THDCs bound to the THRβ1 LBD, with Kd values between 60 µM and 463 µM (Table 1), which can be compared with a typical hit rate around 40% for commonly practiced VS campaigns.74, 75 The results also

ACS Paragon Plus Environment

14

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

revealed the lack of response for triclosan (the predicted non-binder), suggesting that the developed VS protocol can distinguish binders from non-binders. The derived thermodynamic data revealed that the binding of the four confirmed hits was largely enthalpy driven (Table S9). The docking poses and ligand-THRβ1 interactions of the potential THDCs were compared with the co-crystal pose of T3 (Figure S2). We also measured the binding affinity of the positive control T3 (Kd = 1.6 µM). Our value for T3 is comparable to the previously reported affinity of T3 to THRα1 (Kd = 1.8 µM), measured using a similar approach as ours.59 The similarity between the measured binding affinities of T3 may be due to only a single residue difference in their ligand-binding pockets, namely Ser277 in THRα1 and Asn331 in THRβ1.76 The present value, however, is lower than the affinity (0.26 nM)22 reported previously based on a radio ligand binding assay, which may be due to the ITC protocol that uses only the monomer and the LBD of the THR.59, 77 In the Supporting Information, a more thorough discussion on the difference of the two measurements is provided. In addition, we compared the ITC measured binding free energy (∆G) of the active compounds to their corresponding docking and MM-GBSA scores in Figure S3. Despite displaying reasonable correlation between the two sets of data, the scoring results were mainly used in the present study to prioritize the dust contaminants for the experimental testing, and not for the prediction of binding free energies. One

compound

was

identified

as

false

positive

by

ITC,

i.e.,

1-((2,4-

dichlorophenyl)carbamoyl)cyclopropanecarboxylic acid (cyclanilide). The false identification of cyclanilide may have been due to the negligence of intramolecular strain energy in the VS. Cyclanilide is structurally very similar to T3, but its cyclopropanyl group may have high strain energy, which prevents the formation of the proposed pose in the LBS of THRβ1. Poor buffer solubility of bisphenol A (3-chloro-2-hydroxypropyl) (2,3-dihydroxypropyl) ether (BADGE-

ACS Paragon Plus Environment

15

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

HCl-H2O) could be due to its incomplete saturation of the THRβ1 LBD, allowing only an approximate estimation of its affinity (Figure 2). Binding pose examination by all-atom MD simulations. Three compounds (2,4,5-T, BADGE-HCl-H2O, and BP2) were selected for MD simulations based on their structural variation and the ITC data (Figure 3). Overall, all three compounds were stably bound near their corresponding poses predicted by the VS docking procedure during the entire MD simulations (Figure S5). Most notably, the MD-averaged structure of 2,4,5-T, which was the highest affinity binder among the six tested chemicals (Table 1), was nearly identical to its original pose (Figure 3B). Three types of interactions seemed to contribute to the stable binding of 2,4,5-T. Firstly, the carboxyl group of 2,4,5-T engaged in salt-bridged interactions with the two arginine residues (Arg282 and Arg316) in the THRβ1 binding pocket,78 holding the ligand in place throughout the entire MD simulations (Figure 3B). Secondly, chlorine atoms at the 2 and 5 positions of the phenoxy-group (here named Cl2 and Cl5, respectively; Figure S6) formed hydrophobic contacts with Ile275, Ala279, Met310, and Asn331, thereby restricting the motion of the ligand during the MD simulations. Thirdly, Cl5 formed a hydrogen-bond like interaction with Ser314 (the detailed interaction is shown in Figure S6). We suggest that this interaction, together with hydrophobic contact of Cl5 with Met310, contributed to the higher affinity binding of 2,4,5-T (Kd=60 µM) than 2,4-dichlorophenoxyacetic acid (2,4-D) (Kd=463 µM), which differ by only a single chlorine atom. For the remaining two compounds (BP2 and BADGE-HCl-H2O), the MD simulations showed ligand orientations that were slightly displaced relative to their predicted poses (Figure 3C and D). Analysis of the MD trajectories suggested that the displacement in both cases was partly caused by rotations around the two bonds at the central carbon atom (Figure S7), which may not

ACS Paragon Plus Environment

16

Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

be captured accurately in the VS docking. Nevertheless, the two ligands still bound in the binding pocket, likely stabilized by hydrophobic interaction between the two benzyl moieties and the hydrophobic residues in the THRβ1 LBD. In the Supporting Information, we provide a detailed analysis of the MD trajectories, including the distribution of the two rotatable bonds and the RMSD of ligands and protein backbone atoms. Environmental and toxicological implications of the identified THDCs. Four new potential THDCs targeting THRβ1 were identified in this study. To the best of our knowledge, none of these contaminants have previously been reported to bind to THR. It can be concluded that they are weak binders (Kd = 60-460 µM) in the range of previously reported THR binders, such as OH-PCBs (Ki = 30-90 µM)79, bisphenol A (Ki = 200 µM),17 and tetrabromobisphenol A (TBBPA, EC50 ≈ 100 µM).71 Exposure to weak THR binders e.g. bisphenol A during pregnancy has been associated to reduced TH levels in pregnant women and decreased thyroid-stimulating hormone in male neonates.80 Thus, these contaminants warrant further analysis as humans may be chronically exposed at doses which are similar or significantly higher than human T3 levels. Typical levels in human serum of T3 is 0.9-2.8 nM81 whereas levels of the contaminants in human urine of the general population have been reported to be 167 nM (37 ppb) for 2,4-D82, 220 nM (54.3 ng/ml) for BP283, and 8.7 nM (3.4 ng/ml) for BADGE-HCl-H2O84. Among the four contaminants, 2,4,5-T showed the highest affinity to THRβ1. This compound has been used as a herbicide for broad-leaved plants and found in dust samples collected from private homes and daycare centers, such as in North Carolina, USA.85 2,4,5-T has been confirmed to bind to transthyretin and disturb TH transportation.10 The compound is known to be associated with the development of hypothyroidism, and prostate and skin cancer.86-88 Exposure to 2,4,5-T has been shown to be linked to increased occurrence of gestational diabetes mellitus in

ACS Paragon Plus Environment

17

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

women during pregnancy.89 2,4,5-T has also been reported as a weak estrogen receptor antagonist in reporter cell lines.90 The compound was a component of the defoliant Agent Orange deployed during the Vietnam War, exposure to which has been shown to be associated with the development of hypothyroidism, autoimmune thyroiditis and other endocrine disorders, however, which compounds in Agent Orange cause these diseases is not well understood.91 With a similar structure and application to 2,4,5-T, 2,4-D is still in use and among the most frequently detected herbicides in the USA.92 Although 2,4-D did not show any interaction with the estrogen, androgen, or aromatase/steroidogenesis pathways as evaluated in several in vitro assays,93 the compound has been reported to increase the prevalence of hypothyroidism in male applicators.88 The level of 2,4-D in their urine was detected to be 410-2500 ppb.82 2,4-D has been detected in dust samples from homes and daycare centers in North Carolina, Northern and Central California and Ohio, USA.94, 95 Both 2,4,5-T and 2,4-D are classified as a “possible human carcinogen” by the International Agency for Research on Cancer (IARC).96 BADGE-HCl-H2O is an additive or starting agent in personal care products, food packages and material coatings. The concentrations of BADGEs in the USA urine samples were 3-4 folds higher than the corresponding concentrations of bisphenol A.84 It was found in indoor dust samples from the USA, China, Korea, and Japan (444 ng/g).97 Nevertheless, we have not found any reports in the literature on thyroid related adverse effects of this compound. BP2 is a high volume production chemical98 and a representative of the benzophenone-type UV absorbers used in cosmetics and plastic products.99, 100 It was also found in household dust samples from the USA (49.8 ng/g) and East Asian countries (20.1 ng/g).101 The compound has been reported to disturb the transport of thyroid hormones by binding to transthyretin,10 and it inhibits human thyroid peroxidase and causes endometriosis.102,

103

In vivo experiments have

ACS Paragon Plus Environment

18

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

shown that the compound induces weak anti-thyroidal activity in both rodent studies and an amphibian metamorphosis assay.104 BP2 has also been reported as an agonist for both human and fish (rainbow trout) estrogen receptors (ERα and ERβ).105 CONCLUSIONS Here, we developed a robust VS protocol to identify and prioritize potential THDCs targeting THRβ1 for further in vitro or in vivo toxicological evaluation of their thyroid-related adverse effects. The VS protocol was applied to screen an indoor dust contaminants inventory and the results suggested 31 dust contaminants including musk compounds, PFASs, and bisphenol A derivatives (Table S5) as potential binders to THRβ1. Five representative potential binders were selected and tested in vitro for binding to THRβ1 by ITC. The experiments confirmed that four compounds, i.e., 2,4,5-T, BADGE-HCl-H2O, BP2, and 2,4-D, are weak binders to THRβ1. Their potential binding modes were also validated and refined using MD simulations. The developed VS protocol showed to be capable to identify potential THDCs and we suggest further experimental testing and analysis of the top ranked dust contaminants including the four novel THR binders. ASSOCIATED CONTENT Supporting Information. Details of the crystal structures used in the study, the VS performance and enrichment of the models using single crystal structure, the ensemble models and ‘best’ ensemble models rescored with MM-GBSA; analysis of interactions between compounds and the hydrophilic residues, the compounds tested using ITC; force field parameters of the THRβ1 binders, comparison between docking poses of THR disrupters and co-crystal pose of T3, superposition of existing crystal structures, figures describing MD simulation results, and 3D-coordinates of ligand-THRβ1

ACS Paragon Plus Environment

19

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

complexes used for the MD simulations can be found in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org. Notes The authors declare no competing financial interest. ACKNOWLEDGMENT The MM-GBSA calculations and MD simulations were conducted using the resources provided by the Swedish National Infrastructure for Computing (SNIC) at High Performance Computing Center North (HPC2N) and National Supercomputing Center (NSC). The THRβ1LBD was expressed by the Protein Expertise Platform, Umeå University. Funding Sources This study was financed by the MiSSE project through grants from the Swedish Research Council for the Environment, Agricultural Sciences and Spatial Planning (Formas) (210-2012131), by Umeå University, by the Kempe foundation, and by the Swedish Research Council (VR) (521-2011-6427 and 2015-04114).

ABBREVIATIONS 2,4,5-T, 2,4,5-trichlorophenoxyacetic acid; 2,4-D, 2,4-dichlorophenoxyacetic acid; AUC, area under curve; BADGE-HCl-H2O, bisphenol A (3-chloro-2-hydroxypropyl) (2,3-dihydroxypropyl) ether;

BP2,

2,2',4,4'-tetrahydroxybenzophenone;

Cyclanilide,

1-((2,4-

dichlorophenyl)carbamoyl)cyclopropanecarboxylic acid; DMSO, dimethyl sulfoxide; DUD-E, database of useful decoys: Enhanced; EC50, half maximal effective concentration; EF,

ACS Paragon Plus Environment

20

Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

enrichment factor; ER, estrogen receptor; IARC, international agency for research on cancer; IC50, half maximal inhibitory concentration; ITC, isothermal titration calorimetry; LBD, ligand binding domain; LBS, ligand binding site; MD, molecular dynamics; MM-GBSA, molecular mechanics-generalized

Born

surface

area;

OH-PBDE,

hydroxylated

polybrominated

diphenylether; PBDE, polybrominated diphenylether; PCB, polychlorinated biphenyl; PDB, protein

data

bank;

PFAS,

perfluoroalkyl

and

polyfluoroalkyl

substance;

PFOA,

perfluorooctanoic acid; PFOS, perfluorooctanesulfonic acid; PME, particle mesh Ewald; RMSD, root mean squared deviation; ROC, receiver operating characteristic; SP, standard precision; T3, triiodothyronine;

T4,

thyroxine;

TBBPA,

tetrabromobisphenol

A;

TCEP,

tris(2-

carboxyethyl)phosphine; TH, thyroid hormone; THDC, thyroid hormone disrupting chemical; THR, thyroid hormone receptor; TTR, transthyretin; VS, virtual screening; ∆G, binding free energy; ∆H, enthalpy change;

ACS Paragon Plus Environment

21

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

REFERENCES (1)

(2) (3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

Murk, A. J., Rijntjes, E., Blaauboer, B. J., Clewell, R., Crofton, K. M., Dingemans, M. M., Furlow, J. D., Kavlock, R., Kohrle, J., Opitz, R., Traas, T., Visser, T. J., Xia, M., and Gutleb, A. C. (2013) Mechanism-based testing strategy using in vitro approaches for identification of thyroid hormone disrupting chemicals. Toxicol. In Vitro 27, 1320-1346. Crofton, K. M. (2008) Thyroid disrupting chemicals: mechanisms and mixtures. Int. J. Androl. 31, 209-223. Suzuki, G., Takigami, H., Watanabe, M., Takahashi, S., Nose, K., Asari, M., and Sakai, S. (2008) Identification of brominated and chlorinated phenols as potential thyroiddisrupting compounds in indoor dusts. Environ. Sci. Technol. 42, 1794-1800. Meeker, J. D., Johnson, P. I., Camann, D., and Hauser, R. (2009) Polybrominated diphenyl ether (PBDE) concentrations in house dust are related to hormone levels in men. Sci. Total Environ. 407, 3425-3429. Johnson, P. I., Stapleton, H. M., Mukherjee, B., Hauser, R., and Meeker, J. D. (2013) Associations between brominated flame retardants in house dust and hormone levels in men. Sci. Total Environ. 445-446, 177-184. Zoeller, T. R., Dowling, A. L., Herzig, C. T., Iannacone, E. A., Gauger, K. J., and Bansal, R. (2002) Thyroid hormone, brain development, and the environment. Environ. Health Perspect. 110 (Suppl. 3), 355-361. Miller, M. D., Crofton, K. M., Rice, D. C., and Zoeller, R. T. (2009) Thyroid-disrupting chemicals: interpreting upstream biomarkers of adverse outcomes. Environ. Health Perspect. 117, 1033-1041. Noyes, P. D., Lema, S. C., Macaulay, L. J., Douglas, N. K., and Stapleton, H. M. (2013) Low level exposure to the flame retardant BDE-209 reduces thyroid hormone levels and disrupts thyroid signaling in fathead minnows. Environ. Sci. Technol. 47, 10012-10021. Mercier, F., Glorennec, P., Thomas, O., and Le Bot, B. (2011) Organic contamination of settled house dust, a review for exposure assessment purposes. Environ. Sci. Technol. 45, 6716-6727. Zhang, J., Kamstra, J. H., Ghorbanzadeh, M., Weiss, J. M., Hamers, T., and Andersson, P. L. (2015) In Silico Approach To Identify Potential Thyroid Hormone Disruptors among Currently Known Dust Contaminants and Their Metabolites. Environ. Sci. Technol. 49, 10099-10107. Lioy, P. J., Freeman, N. C., and Millette, J. R. (2002) Dust: a metric for use in residential and building exposure assessment and source characterization. Environ. Health Perspect. 110, 969-983. Weiss, J. M., Andersson, P. L., Zhang, J., Simon, E., Leonards, P. E., Hamers, T., and Lamoree, M. H. (2015) Tracing thyroid hormone-disrupting compounds: database compilation and structure-activity evaluation for an effect-directed analysis of sediment. Anal. Bioanal. Chem. 407, 5625-5634. Hamers, T., Kamstra, J. H., Sonneveld, E., Murk, A. J., Kester, M. H., Andersson, P. L., Legler, J., and Brouwer, A. (2006) In vitro profiling of the endocrine-disrupting potency of brominated flame retardants. Toxicol. Sci. 92, 157-173. Weiss, J. M., Andersson, P. L., Lamoree, M. H., Leonards, P. E., van Leeuwen, S. P., and Hamers, T. (2009) Competitive binding of poly- and perfluorinated compounds to the thyroid hormone transport protein transthyretin. Toxicol. Sci. 109, 206-216.

ACS Paragon Plus Environment

22

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

(15) (16) (17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26) (27)

(28) (29) (30)

Evans, R. M. (1988) The steroid and thyroid hormone receptor superfamily. Science 240, 889-895. Zoeller, R. T. (2007) Environmental chemicals impacting the thyroid: targets and consequences. Thyroid 17, 811-817. Moriyama, K., Tagami, T., Akamizu, T., Usui, T., Saijo, M., Kanamoto, N., Hataya, Y., Shimatsu, A., Kuzuya, H., and Nakao, K. (2002) Thyroid hormone action is disrupted by bisphenol A as an antagonist. J. Clin. Endocrinol. Metab. 87, 5185-5190. Freitas, J., Cano, P., Craig-Veit, C., Goodson, M. L., Furlow, J. D., and Murk, A. J. (2011) Detection of thyroid hormone receptor disruptors by a novel stable in vitro reporter gene assay. Toxicol. In Vitro 25, 257-266. Zloh, M., Perez-Diaz, N., Tang, L., Patel, P., and Mackenzie, L. S. (2016) Evidence that diclofenac and celecoxib are thyroid hormone receptor beta antagonists. Life Sci. 146, 6672. Schapira, M., Raaka, B. M., Das, S., Fan, L., Totrov, M., Zhou, Z., Wilson, S. R., Abagyan, R., and Samuels, H. H. (2003) Discovery of diverse thyroid hormone receptor antagonists by high-throughput docking. Proc. Natl. Acad. Sci. U. S. A. 100, 7354-7359. Wang, F. F., Yang, W., Shi, Y. H., and Le, G. W. (2015) Molecular determinants of thyroid hormone receptor selectivity in a series of phosphonic acid derivatives: 3DQSAR analysis and molecular docking. Chem. Biol. Interact. 240, 324-335. Ye, L., Li, Y. L., Mellstrom, K., Mellin, C., Bladh, L. G., Koehler, K., Garg, N., Garcia Collazo, A. M., Litten, C., Husman, B., Persson, K., Ljunggren, J., Grover, G., Sleph, P. G., George, R., and Malm, J. (2003) Thyroid receptor ligands. 1. Agonist ligands selective for the thyroid receptor beta1. J. Med. Chem. 46, 1580-1588. Vedani, A., Zumstein, M., Lill, M. A., and Ernst, B. (2007) Simulating alpha/beta selectivity at the human thyroid hormone receptor: consensus scoring using multidimensional QSAR. ChemMedChem 2, 78-87. Politi, R., Rusyn, I., and Tropsha, A. (2014) Prediction of binding affinity and efficacy of thyroid hormone receptor ligands using QSAR and structure-based modeling methods. Toxicol. Appl. Pharmacol. 280, 177-189. Li, F., Xie, Q., Li, X., Li, N., Chi, P., Chen, J., Wang, Z., and Hao, C. (2010) Hormone activity of hydroxylated polybrominated diphenyl ethers on human thyroid receptor-beta: in vitro and in silico investigations. Environ. Health Perspect. 118, 602-606. Taylor, R. D., Jewsbury, P. J., and Essex, J. W. (2002) A review of protein-small molecule docking methods. J. Comput. Aided Mol. Des. 16, 151-166. Kitchen, D. B., Decornez, H., Furr, J. R., and Bajorath, J. (2004) Docking and scoring in virtual screening for drug discovery: methods and applications. Nat. Rev. Drug Discov. 3, 935-949. Shoichet, B. K., McGovern, S. L., Wei, B., and Irwin, J. J. (2002) Lead discovery using molecular docking. Curr. Opin. Chem. Biol. 6, 439-446. Jorgensen, W. L. (2004) The many roles of computation in drug discovery. Science 303, 1813-1818. McRobb, F. M., Kufareva, I., and Abagyan, R. (2014) In silico identification and pharmacological evaluation of novel endocrine disrupting chemicals that act via the ligand-binding domain of the estrogen receptor alpha. Toxicol. Sci. 141, 188-197.

ACS Paragon Plus Environment

23

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31)

(32)

(33)

(34)

(35) (36) (37)

(38)

(39) (40)

(41) (42)

(43)

(44) (45)

(46)

Page 24 of 32

Rabinowitz, J. R., Little, S. B., Laws, S. C., and Goldsmith, M. R. (2009) Molecular modeling for screening environmental chemicals for estrogenicity: use of the toxicanttarget approach. Chem. Res. Toxicol. 22, 1594-1602. Daniel, L., Buryska, T., Prokop, Z., Damborsky, J., and Brezovsky, J. (2015) Mechanism-based discovery of novel substrates of haloalkane dehalogenases using in silico screening. J. Chem. Inf. Model. 55, 54-62. Rabinowitz, J. R., Goldsmith, M. R., Little, S. B., and Pasquinelli, M. A. (2008) Computational molecular modeling for evaluating the toxicity of environmental chemicals: prioritizing bioassay requirements. Environ. Health Perspect. 116, 573-577. Muralidharan, A. R., Selvaraj, C., Singh, S. K., Sheu, J. R., Thomas, P. A., and Geraldine, P. (2015) Structure-Based Virtual Screening and Biological Evaluation of a Calpain Inhibitor for Prevention of Selenite-Induced Cataractogenesis in an in Vitro System. J. Chem. Inf. Model. 55, 1686-1697. Fukunishi, Y. (2010) Structural ensemble in computational drug screening. Expert Opin. Drug Metab. Toxicol. 6, 835-849. Ellingson, S. R., Miao, Y., Baudry, J., and Smith, J. C. (2015) Multi-conformer ensemble docking to difficult protein targets. J. Phys. Chem. B 119, 1026-1034. Sun, X. Q., Chen, L., Li, Y. Z., Li, W. H., Liu, G. X., Tu, Y. Q., and Tang, Y. (2014) Structure-based ensemble-QSAR model: a novel approach to the study of the EGFR tyrosine kinase and its inhibitors. Acta Pharmacol. Sin. 35, 301-310. Amaro, R. E., Baron, R., and McCammon, J. A. (2008) An improved relaxed complex scheme for receptor flexibility in computer-aided drug design. J. Comput. Aided Mol. Des. 22, 693-705. Bashford, D., and Case, D. A. (2000) Generalized born models of macromolecular solvation effects. Annu. Rev. Phys. Chem. 51, 129-152. Du, J., Sun, H., Xi, L., Li, J., Yang, Y., Liu, H., and Yao, X. (2011) Molecular modeling study of checkpoint kinase 1 inhibitors by multiple docking strategies and prime/MMGBSA calculation. J. Comput. Chem. 32, 2800-2809. Harvey, C. B., and Williams, G. R. (2002) Mechanism of thyroid hormone action. Thyroid 12, 441-446. Flamant, F., Baxter, J. D., Forrest, D., Refetoff, S., Samuels, H., Scanlan, T. S., Vennstrom, B., and Samarut, J. (2006) International Union of Pharmacology. LIX. The pharmacology and classification of the nuclear receptor superfamily: thyroid hormone receptors. Pharmacol. Rev. 58, 705-711. Ibhazehiebo, K., Iwasaki, T., Kimura-Kuroda, J., Miyazaki, W., Shimokawa, N., and Koibuchi, N. (2011) Disruption of thyroid hormone receptor-mediated transcription and thyroid hormone-induced Purkinje cell dendrite arborization by polybrominated diphenyl ethers. Environ. Health Perspect. 119, 168-175. RCSB-PDB. (2015) RCSB Protein Data Bank (PDB). http://www.rcsb.org/pdb/. Jorgensen, W. L., Maxwell, D. S., and TiradoRives, J. (1996) Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225-11236. Schrödinger. (2015) Protein Preparation Wizard 2015-2, Schrodinger, LLC, New York, NY.

ACS Paragon Plus Environment

24

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

(47)

(48)

(49)

(50)

(51)

(52) (53) (54) (55) (56)

(57)

(58)

(59) (60) (61) (62)

(63)

(64)

Verdonk, M. L., Chessari, G., Cole, J. C., Hartshorn, M. J., Murray, C. W., Nissink, J. W., Taylor, R. D., and Taylor, R. (2005) Modeling water molecules in protein-ligand docking using GOLD. J. Med. Chem. 48, 6504-6515. Osterberg, F., Morris, G. M., Sanner, M. F., Olson, A. J., and Goodsell, D. S. (2002) Automated docking to multiple target structures: incorporation of protein mobility and structural water heterogeneity in AutoDock. Proteins 46, 34-40. Li, Y., Shen, J., Sun, X., Li, W., Liu, G., and Tang, Y. (2010) Accuracy assessment of protein-based docking programs against RNA targets. J. Chem. Inf. Model. 50, 11341146. Bortolato, A., Tehan, B. G., Bodnarchuk, M. S., Essex, J. W., and Mason, J. S. (2013) Water network perturbation in ligand binding: adenosine A(2A) antagonists as a case study. J. Chem. Inf. Model. 53, 1700-1713. Mysinger, M. M., Carchia, M., Irwin, J. J., and Shoichet, B. K. (2012) Directory of useful decoys, enhanced (DUD-E): better ligands and decoys for better benchmarking. J. Med. Chem. 55, 6582-6594. ChemAxon. (2014) JChem 6. 2. 0, 2014, ChemAxon (http://www.chemaxon.com). Schrödinger. (2014) Maestro, Schrödinger, LLC, New York, NY. Schrödinger. (2015) Glide, version 6.7, Schrödinger, LLC, New York, NY. Jain, A. N., and Nicholls, A. (2008) Recommendations for evaluation of computational methods. J. Comput. Aided Mol. Des. 22, 133-139. Wei, B. Q., Baase, W. A., Weaver, L. H., Matthews, B. W., and Shoichet, B. K. (2002) A model binding site for testing scoring functions in molecular docking. J. Mol. Biol. 322, 339-355. Li, J., Abel, R., Zhu, K., Cao, Y., Zhao, S., and Friesner, R. A. (2011) The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling. Proteins 79, 2794-2812. Iakovleva, I., Brannstrom, K., Nilsson, L., Gharibyan, A. L., Begum, A., Anan, I., Walfridsson, M., Sauer-Eriksson, A. E., and Olofsson, A. (2015) Enthalpic Forces Correlate with the Selectivity of Transthyretin-Stabilizing Ligands in Human Plasma. J. Med. Chem. 58, 6507-6515. Putcha, B. D., and Fernandez, E. J. (2009) Direct interdomain interactions can mediate allosterism in the thyroid receptor. J. Biol. Chem. 284, 22517-22524. Tellinghuisen, J. (2008) Isothermal titration calorimetry at very low c. Anal. Biochem. 373, 395-397. Turnbull, W. B., and Daranas, A. H. (2003) On the value of c: can low affinity systems be studied by isothermal titration calorimetry? J. Am. Chem. Soc. 125, 14859-14866. Borngraeber, S., Budny, M. J., Chiellini, G., Cunha-Lima, S. T., Togashi, M., Webb, P., Baxter, J. D., Scanlan, T. S., and Fletterick, R. J. (2003) Ligand selectivity by seeking hydrophobicity in thyroid hormone receptor. Proc. Natl. Acad. Sci. U. S. A. 100, 1535815363. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983) Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 79, 926-935. Brooks, B. R., Brooks, C. L., Mackerell, A. D., Nilsson, L., Petrella, R. J., Roux, B., Won, Y., Archontis, G., Bartels, C., Boresch, S., Caflisch, A., Caves, L., Cui, Q., Dinner, A. R., Feig, M., Fischer, S., Gao, J., Hodoscek, M., Im, W., Kuczera, K., Lazaridis, T.,

ACS Paragon Plus Environment

25

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(65)

(66) (67)

(68)

(69)

(70)

(71)

(72)

(73)

(74)

(75) (76)

(77)

Page 26 of 32

Ma, J., Ovchinnikov, V., Paci, E., Pastor, R. W., Post, C. B., Pu, J. Z., Schaefer, M., Tidor, B., Venable, R. M., Woodcock, H. L., Wu, X., Yang, W., York, D. M., and Karplus, M. (2009) CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 30, 1545-1614. Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R. D., Kale, L., and Schulten, K. (2005) Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 1781-1802. Petersen, H. G. (1995) Accuracy and Efficiency of the Particle Mesh Ewald Method. J. Chem. Phys. 103, 3668-3679. Best, R. B., Zhu, X., Shim, J., Lopes, P. E., Mittal, J., Feig, M., and Mackerell, A. D., Jr. (2012) Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone phi, psi and side-chain chi(1) and chi(2) dihedral angles. J. Chem. Theory Comput. 8, 3257-3273. Vanommeslaeghe, K., and MacKerell, A. D., Jr. (2012) Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing. J. Chem. Inf. Model. 52, 3144-3154. Warren, G. L., Andrews, C. W., Capelli, A. M., Clarke, B., LaLonde, J., Lambert, M. H., Lindvall, M., Nevins, N., Semus, S. F., Senger, S., Tedesco, G., Wall, I. D., Woolven, J. M., Peishoff, C. E., and Head, M. S. (2006) A critical assessment of docking programs and scoring functions. J. Med. Chem. 49, 5912-5931. Shelley, J. C., Cholleti, A., Frye, L. L., Greenwood, J. R., Timlin, M. R., and Uchimaya, M. (2007) Epik: a software program for pKa prediction and protonation state generation for drug-like molecules. J. Comput. Aided Mol. Des. 21, 681-691. Kitamura, S., Jinno, N., Ohta, S., Kuroki, H., and Fujimoto, N. (2002) Thyroid hormonal activity of the flame retardants tetrabromobisphenol A and tetrachlorobisphenol A. Biochem. Biophys. Res. Commun. 293, 554-559. Ren, X. M., Zhang, Y. F., Guo, L. H., Qin, Z. F., Lv, Q. Y., and Zhang, L. Y. (2015) Structure-activity relations in binding of perfluoroalkyl compounds to human thyroid hormone T3 receptor. Arch. Toxicol. 89, 233-242. Du, G., Hu, J., Huang, H., Qin, Y., Han, X., Wu, D., Song, L., Xia, Y., and Wang, X. (2013) Perfluorooctane sulfonate (PFOS) affects hormone receptor activity, steroidogenesis, and expression of endocrine-related genes in vitro and in vivo. Environ. Toxicol. Chem. 32, 353-360. Zhu, T., Cao, S., Su, P. C., Patel, R., Shah, D., Chokshi, H. B., Szukala, R., Johnson, M. E., and Hevener, K. E. (2013) Hit identification and optimization in virtual screening: practical recommendations based on a critical literature analysis. J. Med. Chem. 56, 6560-6572. Schneider, G. (2010) Virtual screening: an endless staircase? Nat. Rev. Drug Discov. 9, 273-276. Wagner, R. L., Huber, B. R., Shiau, A. K., Kelly, A., Cunha Lima, S. T., Scanlan, T. S., Apriletti, J. W., Baxter, J. D., West, B. L., and Fletterick, R. J. (2001) Hormone selectivity in thyroid hormone receptors. Mol. Endocrinol. 15, 398-410. Putcha, B. D., Wright, E., Brunzelle, J. S., and Fernandez, E. J. (2012) Structural basis for negative cooperativity within agonist-bound TR:RXR heterodimers. Proc. Natl. Acad. Sci. U. S. A. 109, 6084-6087.

ACS Paragon Plus Environment

26

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

(78)

(79)

(80)

(81)

(82) (83)

(84)

(85)

(86) (87)

(88)

(89)

(90)

(91)

Martinez, L., Nascimento, A. S., Nunes, F. M., Phillips, K., Aparicio, R., Dias, S. M., Figueira, A. C., Lin, J. H., Nguyen, P., Apriletti, J. W., Neves, F. A., Baxter, J. D., Webb, P., Skaf, M. S., and Polikarpov, I. (2009) Gaining ligand selectivity in thyroid hormone receptors via entropy. Proc. Natl. Acad. Sci. U. S. A. 106, 20717-20722. Cheek, A. O., Kow, K., Chen, J., and McLachlan, J. A. (1999) Potential mechanisms of thyroid disruption in humans: interaction of organochlorine compounds with thyroid receptor, transthyretin, and thyroid-binding globulin. Environ. Health Perspect. 107, 273278. Chevrier, J., Gunier, R. B., Bradman, A., Holland, N. T., Calafat, A. M., Eskenazi, B., and Harley, K. G. (2013) Maternal urinary bisphenol a during pregnancy and maternal and neonatal thyroid function in the CHAMACOS study. Environ. Health Perspect. 121, 138-144. Bunevicius, R., Kazanavicius, G., Zalinkevicius, R., and Prange, A. J., Jr. (1999) Effects of thyroxine as compared with thyroxine plus triiodothyronine in patients with hypothyroidism. N. Engl. J. Med. 340, 424-429. Burns, C. J., and Swaen, G. M. (2012) Review of 2,4-dichlorophenoxyacetic acid (2,4-D) biomonitoring and epidemiology. Crit. Rev. Toxicol. 42, 768-786. Asimakopoulos, A. G., Thomaidis, N. S., and Kannan, K. (2014) Widespread occurrence of bisphenol A diglycidyl ethers, p-hydroxybenzoic acid esters (parabens), benzophenone type-UV filters, triclosan, and triclocarban in human urine from Athens, Greece. Sci. Total Environ. 470-471, 1243-1249. Wang, L., Wu, Y., Zhang, W., and Kannan, K. (2012) Widespread occurrence and distribution of bisphenol A diglycidyl ether (BADGE) and its derivatives in human urine from the United States and China. Environ. Sci. Technol. 46, 12968-12976. Morgan, M. K., Wilson, N. K., and Chuang, J. C. (2014) Exposures of 129 preschool children to organochlorines, organophosphates, pyrethroids, and acid herbicides at their homes and daycares in North Carolina. Int. J. Environ. Res. Public Health 11, 3743-3764. Ansbaugh, N., Shannon, J., Mori, M., Farris, P. E., and Garzotto, M. (2013) Agent Orange as a risk factor for high-grade prostate cancer. Cancer 119, 2399-2404. Clemens, M. W., Kochuba, A. L., Carter, M. E., Han, K., Liu, J., and Evans, K. (2014) Association between Agent Orange exposure and nonmelanotic invasive skin cancer: a pilot study. Plast. Reconstr. Surg. 133, 432-437. Goldner, W. S., Sandler, D. P., Yu, F., Shostrom, V., Hoppin, J. A., Kamel, F., and LeVan, T. D. (2013) Hypothyroidism and pesticide use among male private pesticide applicators in the agricultural health study. J. Occup. Environ. Med. 55, 1171-1178. Saldana, T. M., Basso, O., Hoppin, J. A., Baird, D. D., Knott, C., Blair, A., Alavanja, M. C., and Sandler, D. P. (2007) Pesticide exposure and self-reported gestational diabetes mellitus in the Agricultural Health Study. Diabetes Care 30, 529-534. Lemaire, G., Mnif, W., Mauvais, P., Balaguer, P., and Rahmani, R. (2006) Activation of alpha- and beta-estrogen receptors by persistent pesticides in reporter cell lines. Life Sci. 79, 1160-1169. Yi, S. W., Hong, J. S., Ohrr, H., and Yi, J. J. (2014) Agent Orange exposure and disease prevalence in Korean Vietnam veterans: the Korean veterans health study. Environ. Res. 133, 56-65.

ACS Paragon Plus Environment

27

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(92)

(93)

(94)

(95)

(96) (97)

(98)

(99)

(100)

(101)

(102)

(103)

(104) (105)

Page 28 of 32

Ensminger, M. P., Budd, R., Kelley, K. C., and Goh, K. S. (2013) Pesticide occurrence and aquatic benchmark exceedances in urban surface waters and sediments in three urban areas of California, USA, 2008-2011. Environ. Monit. Assess. 185, 3697-3710. Coady, K. K., Kan, H. L., Schisler, M. R., Gollapudi, B. B., Neal, B., Williams, A., and LeBaron, M. J. (2014) Evaluation of potential endocrine activity of 2,4dichlorophenoxyacetic acid using in vitro assays. Toxicol. In Vitro 28, 1018-1025. Morgan, M. K., Sheldon, L. S., Thomas, K. W., Egeghy, P. P., Croghan, C. W., Jones, P. A., Chuang, J. C., and Wilson, N. K. (2008) Adult and children's exposure to 2,4-D from multiple sources and pathways. J. Expo. Sci. Environ. Epidemiol. 18, 486-494. Deziel, N. C., Colt, J. S., Kent, E. E., Gunier, R. B., Reynolds, P., Booth, B., Metayer, C., and Ward, M. H. (2015) Associations between self-reported pest treatments and pesticide concentrations in carpet dust. Environ. Health 14, 27-38. Landrigan, P. J., and Benbrook, C. (2015) GMOs, Herbicides, and Public Health. N. Engl. J. Med. 373, 693-695. Wang, L., Liao, C., Liu, F., Wu, Q., Guo, Y., Moon, H. B., Nakata, H., and Kannan, K. (2012) Occurrence and human exposure of p-hydroxybenzoic acid esters (parabens), bisphenol A diglycidyl ether (BADGE), and their hydrolysis products in indoor dust from the United States and three East Asian countries. Environ. Sci. Technol. 46, 11584-11593. Rannar, S., and Andersson, P. L. (2010) A novel approach using hierarchical clustering to select industrial chemicals for environmental impact assessment. J. Chem. Inf. Model. 50, 30-36. Liao, C., and Kannan, K. (2014) Widespread occurrence of benzophenone-type UV light filters in personal care products from China and the United States: an assessment of human exposure. Environ. Sci. Technol. 48, 4103-4109. Nakajima, D., Asada, S., Kageyama, S., Yamamoto, T., Kuramochi, H., Tanaka, N., Takeda, K., and Goto, S. (2006) Activity related to the carcinogenicity of plastic additives in the benzophenone group. J. UOEH. 28, 143-156. Wang, L., Asimakopoulos, A. G., Moon, H. B., Nakata, H., and Kannan, K. (2013) Benzotriazole, benzothiazole, and benzophenone compounds in indoor dust from the United States and East Asian countries. Environ. Sci. Technol. 47, 4752-4759. Kunisue, T., Chen, Z., Buck Louis, G. M., Sundaram, R., Hediger, M. L., Sun, L., and Kannan, K. (2012) Urinary concentrations of benzophenone-type UV filters in U.S. women and their association with endometriosis. Environ. Sci. Technol. 46, 4624-4632. Schmutzler, C., Bacinski, A., Gotthardt, I., Huhne, K., Ambrugger, P., Klammer, H., Schlecht, C., Hoang-Vu, C., Gruters, A., Wuttke, W., Jarry, H., and Kohrle, J. (2007) The ultraviolet filter benzophenone 2 interferes with the thyroid hormone axis in rats and is a potent in vitro inhibitor of human recombinant thyroid peroxidase. Endocrinology 148, 2835-2844. OECD. (2008) Report of the validation of the amphibian metamorphosis assay (Phase 3). https://www.oecd.org/chemicalsafety/testing/41620749.pdf (accessed May 6,2016). Molina-Molina, J. M., Escande, A., Pillon, A., Gomez, E., Pakdel, F., Cavailles, V., Olea, N., Ait-Aissa, S., and Balaguer, P. (2008) Profiling of benzophenone derivatives using fish and human estrogen receptor-specific in vitro bioassays. Toxicol. Appl. Pharmacol. 232, 384-395.

ACS Paragon Plus Environment

28

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

Table 1. Studied compounds with chemical abstract services (CAS) registry numbers, abbreviations, industrial application, prediction from the virtual screening protocol and binding affinities (Kd) determined using isothermal titration calorimetry. Compounds (abbreviations)

CAS No.

Applications

Prediction Kd (µM)a

triiodothyronine (T3)

689302-3

Positive control

-

1.6

2,4,5-trichlorophenoxyacetic acid (2,4,5-T)

93-76-5

Herbicide

Active

60

bisphenol A (3-chloro-2-hydroxypropyl) (2,3- 227947- Plasticizer dihydroxypropyl) ether (BADGE-HCl-H2O) 06-0

Active

87.5b

2,2',4,4'-tetrahydroxybenzophenone (BP2)

131-555

UV filter

Active

200

2,4-dichlorophenoxyacetic acid (2,4-D)

94-75-7

Herbicide

Active

463

1-((2,4113136- Herbicide dichlorophenyl)carbamoyl)cyclopropanecarboxylic 77-9 acid (Cyclanilide)

Active

ND

5-chloro-2-(2,4-dichlorophenoxy)phenol (Triclosan)c

338034-5

Antimicrobial Inactive agent (Negative control)

ND

a

ND refers to non-detected affinity; bThe binding affinity of BADGE-HCl-H2O was only an approximate value, probably due to poor solubility in the ITC buffer; cTriclosan was tested as a negative control. It was identified among the 131 initial hits but showed no interaction with the four hydrophilic residues (Figure 1).

ACS Paragon Plus Environment

29

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32

Figure 1. Scheme adopted in this study with details of the virtual screening (VS) model development process, the isothermal titration calorimetry (ITC) validation study and the molecular dynamics (MD) study of ligand-protein interactions. ssDock model: single structure docking model; esDock model: ensemble structure docking model. In brackets are the numbers of hits at different steps.

ACS Paragon Plus Environment

30

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

Figure 2. Molecular structures of active compounds and thermograms of their binding to THRβ1 determined by isothermal titration calorimetry. The abbreviations of compounds are defined in Table 1. For BADGE-HCl-H2O, the dissociation constant (Kd) was an approximate value likely due to its poor solubility in the ITC buffer.

ACS Paragon Plus Environment

31

Chemical Research in Toxicology

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 32

Figure 3. Comparison of docking and MD-averaged poses of potential THR disrupters: (B) 2,4,5T, (C) BP2, and (D) BADGE-HCl-H2O. In (A), the THRβ1 ligand binding domain is shown in cartoon, and the location of the ligand binding pocket is indicated with the two helices (H3 and H5), a binder (2,4,5-T, pink ball-and-stick) and two arginine residues (Arg282 and Arg316). In (B)-(D), MD-averaged poses (from the last 10 ns) are shown in stick with different colors, and VS proposed poses in gray stick, respectively. In (B), the electrostatic interactions between ligands and residues are shown in purple dash.

ACS Paragon Plus Environment

32