SCREENING METHODS FOR THYROID HORMONE DISRUPTORS BASED ON CO-REGULATOR INVOLVED SIMULATIONS
20170285007 · 2017-10-05
Inventors
- Hongxia Yu (Nanjing, CN)
- Wei Shi (Nanjing, CN)
- Qinchang Chen (Nanjing, CN)
- Xiaoxiang Wang (Nanjing, CN)
Cpc classification
G16B15/00
PHYSICS
G01N2333/723
PHYSICS
G16B5/00
PHYSICS
G16B15/30
PHYSICS
International classification
Abstract
The present patent relates to a method for qualitative identification and quantitative prediction of thyroid hormone disrupting chemicals base on the interaction between thyroid hormone receptor and co-regulators (coactivator and corepressor).The method identifies chemicals as passive antagonists, active antagonists and agonists by means of co-regulator involved molecular dynamics simulations, and predicts the relative disrupting potencies by use of binding free energy, therefore, may be used for screening of thyroid hormone disruptors among environmental pollutants. Upon more comprehensive consideration of the functioning mechanism of thyroid hormone receptor, the present invention is able to sufficiently identify thyroid hormone disruptors as agonists and antagonists, and gives more accurate prediction of the disrupting potency. Further, since nuclear receptors, just as thyroid hormone receptor, are strongly associated with co-regulators, the method may be expanded to the screening of nuclear receptor mediated endocrine disruptors.
Claims
1. A screening method for thyroid hormone disruptors based on co-regulator involved simulations, comprising: 1) apo conformations of thyroid hormone (TH) receptor (TR) ligand binding domain (LBD) are constructed by used of the reported protein crystal as templates, followed by quality evaluation and structural optimization; ligand structures are built and optimized; the ligands are docked into the TRs to form ligand-TR complexes; Molecular dynamics (MD) simulations are performed for the ligand-TR complexes; 2) equilibrated conformations are extracted from the MD trajectories; the relative ligands are identified as passive antagonists for those helix-12s (H12s)are induced to block the coactivator binding site; 3) coactivator and corepressor are docked to the equilibrated conformations extracted in step 2) whose co-regulator binding sites are not blocked and thus be exposed to co-regulators; whether the equilibrated conformations are more attractive to coactivator or corepressoris estimated according to the docking scores; the relative ligands are identified as agonists if the equilibrated conformations are more selectively bind to coactivators; the relative ligands are identified as activeantagonists if the equilibrated conformations are more selectively bind to corepressors; 4) for those identified as agonists in step 3), MD simulations are performed again for the obtained ligand-TR-coactivator complexes; for those identified as active antagonists in step 3), MD simulations are performed again for the obtained ligand-TR-corepressor complexes; 5) bioassays are performed for the representative TH disrupting chemicals to obtain the relative disrupting potencies;binding free energies with TR are calculated for the passive antagonists, active antagonists and agonists in step 2) and4); the predicted TH disrupting potencies are calculated by the prediction models achieved by regression analysis between binding free energy and relative TH disrupting potencies.
2. The method of claim 1, wherein the construction of the apo conformations in step 1) is based on homology modeling: the loop between H11 and H12 is constructed according to the template coded 1A52 in the protein database; the rest part is constructed based on the reported crystal structures of TRs in the protein database.
3. The method of claim 1, wherein the optimization of ligand structures in step 1)is first performedwith MM2 (Molecular Mechanics, Allinger Force Field version 2), followed by Powell gradient algorithm with the Tripos force field.
4. The method of claim 1, wherein the MD simulations of ligand-TR complexes in step 1) are performed following the protocol: CHARMM force fields are applied to the receptors and ligands, respectively; each complex is immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm; Na+ and Cl− ions are added to neutralize the system; all MD simulations are performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions; electrostatic interactions are calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm; all simulations are performed for 20-22 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps.
5. The method of claim 1, wherein the MD simulations of ligand-TR-coactivator/corepressor complexes in step 4) are performed following the protocol: CHARMM force fields are applied to the receptors with coactivator/corepressor and ligands, respectively; each complex is immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm; Na+ and Cl− ions are added to neutralize the system; all MD simulations are performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions; electrostatic interactions are calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm; all simulations are performed for 20-22 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps.
6. The method of claim 1, wherein the blocking of coactivator binding site in step 2) is evaluated according to the distance between H12 and key residues: for TRα, H12 is considered blocking coactivator binding site if the distances with V284 and K306 are less than 3 Å, and the distances with K288 and 1302 are less than 5 Å; for TRβ, H12 is considered blocking coactivator binding site if the distances with V230 and K252 are less than 3 Å, and the distances with K234 and 1248 are less than 5 Å.
7. The method of claim 1, wherein calculation of binding free energy in step 5) is based on MM-PBSA (molecular mechanics Poisson-Boltzmann surface area) method.
8. The method of claim 1, wherein the binding free energies between ligands and TR used for predictions of thyroid disrupting potencies in step 5) are following the rules: for prediction of passive antagonistic potencies, binding free energies are calculated by use of MD trajectories of passive antagonist-TR complexes in step 1); for prediction of active antagonistic and agonistic potencies, binding free energies are calculated by use of MD trajectories of active antagonist-TR-corepressor and agonist-TR-coactivator complexes, respectively, in step 4).
Description
BRIEF DESCRIPTION OF THE DRAWINGS
[0009]
[0010]
[0011]
[0012]
[0013]
[0014]
[0015]
DETAILED DESCRIPTION
Technical solutions
[0016] The present screening method for thyroid hormone disruptors based on co-regulator involved simulations can be understood by reference to the following description and the examples included therein.
[0017] The protocol of the present screening method for thyroid hormone disruptors based on co-regulator involved simulations is depicted in
[0023] In an embodiment of the step 1), the reported crystal structures are those can be searched and downloaded from the protein database (http://www.rcsb.org/pdb/home/home.do), which includes all published structures that are approval.
[0024] In an embodiment of the step 1), the construction of the apo conformations is based on homology modeling, the loop between H11 and H12 is constructed according to the template coded 1A52 in the protein database, the rest part is constructed based on the reported crystal structures (PDB codes: 4LNW and 1NQ0) of TRs in the protein database.
[0025] In an embodiment of the step 1), the quality of the constructed apo conformations are evaluated by Ramachandran plot.
[0026] In an embodiment of the step 1),the optimization of ligand structures is first performed with MM2 (Molecular Mechanics, Allinger Force Field version 2), followed by Powell gradient algorithm with the Tripos force field.
[0027] In an embodiment of the step 1), the MD simulations of ligand-TR complexes are performed following the protocol: CHARMM force fields are applied to the receptors and ligands, respectively; each complex is immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm; Na+ and Cl− ions are added to neutralize the system; all MD simulations are performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions; electrostatic interactions are calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm; all simulations are performed for 20-22 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps.
[0028] In an embodiment of the step 4), the MD simulations of ligand-TR-coactivator/corepressor complexes are performed following the protocol: CHARMM force fields are applied to the receptors with coactivator/corepressor and ligands, respectively; each complex is immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm; Na+ and Cl− ions are added to neutralize the system; all MD simulations are performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions; electrostatic interactions are calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm; all simulations are performed for 20-22 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps.
[0029] In an embodiment of the step 2), equilibrated conformations are obtained by calculation of root-mean-square deviation (RMSD) of H12, reference to backbone of the apo structures. Conformation extracted when the RMSD of H12 is equilibrated is considered equilibrated conformation.
[0030] In an embodiment of the step 2), equilibrated conformations can also be obtained by extracting conformation when H12 is in the average position of the fluctuation.
[0031] In an embodiment of the step 2), the blocking of coactivator binding site is evaluated according to the distance between H12 and key residues: for TRα, H12 is considered blocking coactivator binding site if the distances with V284 and K306 are less than 3 Å, and the distances with K288 and 1302 are less than 5 Å; for TRβ, H12 is considered blocking coactivator binding site if the distances with V230 and K252 are less than 3 Å, and the distances with K234 and 1248 are less than 5 Å.
[0032] In an embodiment of the step 3), docking of co-regulators to the equilibrated conformations are performed as follows:
[0033] Reference structures of receptors containing corepressor (PDB code: 2OVM) or coactivator (PDB code: 2B1V) are aligned to the extracted conformations, and the positions of co-regulators are set as the reference positions of the input co-regulators. Co-regulators were then docked to the extracted TRs using shape-based 3D fast Fourier transform (FFT) docking methods. The receptor and ligand range angles are all set to 15 degrees to make sure the co-regulators not rotate far away from the reference positions. The docking result with least Edock(docking score) out of 500 results is chosen as the final docking result.
[0034] In an embodiment of the step 3), conformations withEdockof coactivator greater than corepressor are considered to selectively bind to corepressor; conformations with E.sub.dock of corepressor greater than coactivator are considered to selectively bind to coactivator.
[0035] In an embodiment of the step 5), binding free energy is calculated based on snapshots of every 100 ps extracted from the MD trajectories using MM-PBSA (molecular mechanics Poisson-Boltzmann surface area) method.
[0036] In an embodiment of the step 5), the binding free energies between ligands and TR used for predictions of thyroid disrupting potencies are following the rules: for prediction of passive antagonistic potencies, binding free energies are calculated by use of MD trajectories of passive antagonist-TR complexes in step 1); for prediction of active antagonistic and agonistic potencies, binding free energies are calculated by use of MD trajectories of active antagonist-TR-corepressor and agonist-TR-coactivator complexes, respectively, in step 4).
[0037] In an embodiment of the step 5), the prediction model by use of GH3 cell proliferation assay, with 11 HO-PBDEs as training set and 2 as validation set, is built and given in Equation. 1
-logRIC.sub.20=0.4115−0.0244ΔG.sub.pas/act (1)
where RIC.sub.20 (mol/L) is 20% inhibition of proliferation of GH3 cells induced by 0.5 nM T3 (triiodothyronine); ΔG.sub.pas/act is the combined binding free energy combining binding free energies of TRα and TRβ with passive and active antagonists.
Beneficial Effect
[0038] In compare with previous methods, the present invention has the following significant results.
[0039] (1) The present invention considers the function of coactivator and corepressor, which are essential for the functioning of TR, in the toxicity simulations. The involvement of co-regulators gives more comprehensive consideration of the functioning mechanism of TR.
[0040] (2) The MD simulations perfectly exhibited the interaction among ligand, TR and co-regulators, which are able to effect the functioning the others.
[0041] (3) Based on the mechanism comprehension, thyroid hormone disruptors are identified and classified as passive antagonists, active antagonists and agonists, which leads to more accurate potency prediction.
[0042] (4) Agonistic and antagonistic potencies are quantitatively predicted by prediction model using binding free energy.
[0043] (5) The present invention gives more simple and sufficient prediction with less cost than in vivo and in vitro screening.
EXAMPLES
[0044] The following examples further illustrate the method of identifying TH disrupting chemicals and building of quantitative prediction model. It will be understood, however, that the examples are for better comprehension of the present invention, which should not limit what is claimed in the claim section. Many variations and modifications of the methods can be made while remaining within the scope and spirit of the present invention. For example, human and rat TRα and TRβ are used as the nuclear receptor, HO-PBDEs are used as the potential thyroid hormone disruptors in the following examples. Other nuclear receptors and other endocrine disruptors can be used for construction of virtual screen models base on the protocol.
Example 1
[0045] Apo structures of human TRα- and TRβ-LBD were built by homology modeling. Previously reported human TRα- (PDB code: 4LNX) and TRβ-LBD (PDB code: 1NQ0) were used as templates of the main bodies of TRα and TRβ, respectively, and apo structure of estrogen receptor (PDB code: 1A52) was used as the template of the loops between H11 and H12 of both receptors. Qualities of the constructed apo structures were further evaluated with Ramachandran plot generated by the Structure Analysis and Verification Server (SAVES, http://services.mbi.ucla.edu/SAVES/).
[0046] Eight HO-PBDEs with T3 are selected as ligands for this example. Ligand structures are first optimized with MM2 (Molecular Mechanics, Allinger Force Field version 2), and then optimized using Powell gradient algorithm with the Tripos force field. Then the Surflex-Dock program interfaced with SYBYL 7.3 was used to dock the optimized ligands into the docking cavities of TR-LBDs. Docking scores of the ligands were calculated and the ligand-TR complexes were used for MD simulations.
[0047] MD simulations were performed by use of GROMACS software package. The complexes were immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm. Na+ and Cl− ions were added to neutralize the system. All MD simulations were performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions. Electrostatic interactions were calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm. All simulations were performed for 20 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps. After the MD simulations, snapshots of every 100ps were extracted for calculation of binding free energy by use of g_mmpbsa program.
[0048] Binding affinities of ligands binding with TR were detected in the competitive binding assays, by use of Fluorescein-labeled T3 as probe. The results (
[0049] Lower binding free energy and higher docking score reveal greater binding affinity. However, the results of docking score and binding free energy were dramatically different (
Example 2
[0050] Apo structures of rat TRα- and TRβ-LBD were constructed by homology modeling. Because no rat TR crystal structures was reported, human TRα- (PDB code: 4LNX) and TRβ-LBD (PDB code: 1NQ0) were used as templates of the main bodies of TRα and TRβ, respectively, and apo structure of estrogen receptor (PDB code: 1A52) was used as the template of the loops between H11 and H12 of both receptors. Qualities of the constructed apo structures were further evaluated with Ramachandran plot generated by the Structure Analysis and Verification Server (SAVES, http://services.mbi.ucla.edu/SAVES/).
[0051] Sixteen HO-PBDEs with T3 are selected as ligands for this example. Ligand structures are first optimized with MM2 (Molecular Mechanics, Allinger Force Field version 2), and then optimized using Powell gradient algorithm with the Tripos force field. Then the Surflex-Dock program interfaced with SYBYL 7.3 was used to dock the optimized ligands into the docking cavities of TR-LBDs. Docking scores of the ligands were calculated and the ligand-TR complexes were used for MD simulations.
[0052] MD simulations were performed by use of GROMACS software package. The complexes were immersed with TIP3P water molecules in a box, keeping the minimum distance between the complex and the boundary larger than 1.4 nm. Na+ and Cl− ions were added to neutralize the system. All MD simulations were performed in the NPT (constant pressure and temperature) ensemble in 1 atm, 300 K, with periodic boundary conditions. Electrostatic interactions were calculated with the particle mesh Ewald method, with van der Waals interactions cutoff of 1.0 nm. All simulations were performed for 20 ns using a 2 fs time step, and snapshots for analysis are saved every 2 ps.
[0053] The rat pituitary tumor cell line GH3 was cultured for determination of the relative TH disrupting potencies. Thirteen out of sixteen HO-PBDEs were determined to be thyroid hormone antagonists, with the relative anti-TH potencies (″logRIC20) ranging from 6.51 to 8.42 (Tab. 1). Due to cytotoxicity, the other 3 HO-PBDEs were not tested for anti-TH potencies.
[0054] Snapshots of equilibrated conformations were extracted to investigate the reposition of H12 and identify chemicals as passive antagonists. For T3-bound TRβ, AF-2 was exposed to the co-regulators, so that coactivators were able to bind it (
[0055] For those co-regulator binding site was exposed to co-regulators, coactivator and corepressor were docked to the co-regulator binding site by use of Hex 8.0.0 software. Reference structures of receptors containing corepressor (PDB code: 20VM) or coactivator (PDB code: 2B1V) are aligned to the equilibrated conformations, and the positions of co-regulators are set as the reference positions of the input co-regulators. Co-regulators were then docked to the extracted TRs using shape-based 3D fast Fourier transform (FFT) docking methods. The receptor and ligand range angles are all set to 15 degrees to make sure the co-regulators not rotate far away from the reference positions. The docking result with least Edock out of 500 results is chosen as the final docking result.
[0056] Docking of co-regulators to TR was performed to identify chemicals as active antagonists and agonists. Most of the HO-PBDE-bound TRs bound to corepressor with Edock value less than that of the coactivator except 6-HO-BDE-47-bound TRβ (
[0057] MD simulations were performed again for the active antagonist-TR-corepressor complexes. The procedures of simulations were the same as above.Ligand-receptor binding free energies of ligand-TR complexes (ΔG.sub.lig-rTRα and ΔG.sub.lig-rTRβ) and active antagonist-TR-corepressor complexes (ΔG.sub.lig-rTRα/cor and ΔG.sub.lig-rTRβ/cor) were then calculated (Tab. 1). Among the thirteen HO-PBDEs that were identified as (active or passive) antagonists by MD simulations and were detected antagonists in the GH3 cell proliferation assay, two (2′-HO-BDE-66 and 4-HO-BDE-90) were selected as validation set and the other eleven were considered training set. Because the GH3 cells contain both TRα and TRβ, the binding free energies were combined before regression analysis. In this example, ΔG.sub.sum,cor-rTR represented the sum of ΔG.sub.lig-rTRα and ΔG.sub.lig-rTRβ, while ΔG.sub.pas/act was the combination of ΔG.sub.lig-rTRα or ΔG.sub.lig-rTRβ of passive antagonist and ΔG.sub.lig-rTRα/coror ΔG.sub.lig-rTRβ/cor of active antagonist.
TABLE-US-00001 TABLE 1 Binding free energy and the relative anti-TH potencies TRα (kJ/mol) TRβ (kJ/mol) Combined (kJ/mol) Ligand −log RIC.sub.20 ΔG.sub.lig-rTRα ΔG.sub.lig-rTRα/cor ΔG.sub.lig-rTRβ ΔG.sub.lig-rTRβ/cor ΔG.sub.sum, cor-rTR ΔG.sub.pas/act 4-HO-BDE-188 8.42 −175.2 −170.1 −169.7 −169.4 −344.9 −339.5 3-HO-BDE-100 7.83 −139.7 −150.7 −137.5 −140.6 −277.2 −291.3 4′-HO-BDE-101 8.11 −150.7 −150.6 −150.6 −146.5 −301.3 −297.1 6-HO-BDE-157 8.14 −149.5 −150.6 −156.8 −154.3 −306.3 −304.9 6-HO-BDE-90 7.76 −143.3 −145.6 −151.7 −152.9 −295.1 −298.5 4′-HO-BDE-49 7.30 −141.5 −142.1 −141.4 −141.6 −282.9 −283.7 2′-OH-BDE-66 7.09 −142.5 −139.0 −136.4 −133.1 −278.8 −272.1 6′-HO-BDE-17 6.78 −128.1 −136.4 −135.0 −133.7 −263.1 −270.1 2′-OH-BDE-68 7.15 −136.2 −132.1 −142.5 −147.3 −278.7 −279.5 4′-OH-BDE-17 6.51 −126.9 −126.2 −131.4 −131.0 −258.3 −257.2 6-HO-BDE-137 8.03 −161.3 −156.3 −162.2 — −323.5 −318.6 2-HO-BDE-123 7.44 −152.5 −149.8 −149.5 — −302.0 −299.3 4-HO-BDE-90 7.59 −147.8 −151.8 −146.2 — −293.9 −297.9 6-HO-BDE-47 — −145.3 −140.8 −135.4 — −280.8 — 6-HO-BDE-85 — −138.3 — −150.8 −151.3 −289.1 −289.6 6′-HO-BDE-99 — −139.1 −137.7 −144.3 −148.5 −283.4 −286.2
[0058] As shown in
−logRIC.sub.20=0.4115−0.0244ΔG.sub.pas/act (2)
[0059] The determination coefficient R.sup.2=0.826, and the external explained variance Q.sup.2.sub.ext=0.926. The results indicated that the prediction model had good predictive ability. This prediction model can be utilized to predict the relative anti-TH potencies of potential TR antagonists, especially among HO-PBDEs, after qualitative identification following the procedure above.