一种基于核受体共调节因子的甲状腺激素干扰物虚拟筛选及其干扰活性的定量计算方法

摘要

本发明公开了一种基于核受体共调节因子的甲状腺激素干扰物虚拟筛选及其干扰活性的定量计算方法，它基于核受体共调节因子调节机制，采用分子动力学模拟和结合自由能计算方法，通过12号螺旋的稳定位置、甲状腺激素的后踢与共激活因子、共抑制因子相互关系判断受测化合物的甲状腺激素干扰效应，区分积极的抗性、消极的抗性和拟性效应，通过计算配体‑受体结合自由能预测相关活性大小。与现有技术相比，本发明相比于传统的体外筛查方法，成本低廉、操作简单、效率更高；相比于其他虚拟筛选方法，此方法更全面地考虑了分子作用机制，可以区分拟性和抗性，定性和定量预测结果更可靠。

说明书

技术领域

本发明属于预测毒理学领域，具体涉及一种基于核受体共调节因子的甲状腺激素干扰物虚拟筛选及其干扰活性的定量计算方法。

背景技术

甲状腺激素是生物体内重要的激素，对生物体生长、分化、代谢和能量平衡等具有重要作用。甲状腺激素干扰物质可以影响胎儿大脑、中枢神经系统的发育，影响心脏的发育，造成各种疾病，其危害不容忽视。随着越来越多人工合成的或自然形成的化学物质被检测出具有甲状腺干扰活性，例如多溴联苯醚、多氯联苯、双酚A等，这些化学物质受到人们的广泛关注。为了筛查可疑的甲状腺激素干扰物质，人们发展了各种有效的体内和体外实验方法，其中主要为包括竞争结合、双杂交酵母和报告基因等体外试验方法。然而，一方面，采用这些方法费时费力，而且还相当昂贵；另一方面，环境中有成千上万的化学物质，很难逐一筛查。

面对这样的挑战，科学家们发展了基于计算机模拟的虚拟筛选方法，这些方法也逐渐得到人们的认可，如定量构效关系(Quantitative Structure-ActivityRelationship,QSAR)、分子对接和分子动力学(Molecular Dynamics,MD)模拟等。然而，定量构效关系需要大量活性数据用于建模，且忽视了小分子与受体间的相互作用；分子对接也有忽视受体柔性等缺点，难以准确反映化合物的活性状态；分子动力学模拟虽然能很好地反映分子作用机制，但通常的分子动力学模拟仅考虑配体与受体的作用，忽略了其他共调节因子，如共抑制因子和共激活因子的作用，容易导致假阳性/阴性结果的产生，且没有引入有针对性的预测参数达到定量预测的效果。

作为核受体家族的一员，甲状腺激素受体(Thyroid Hormone Receptor,TR)的配体结合区域(Ligand binding domain,LBD)在甲状腺激素受体作用过程起着至关重要的作用。激动剂如三碘甲腺原氨酸T3等通过与配体结合区域结合，诱导其产生构象变化，激活共激活因子结合位置，进而促进相关基因的转录，相关效应为拟甲状腺激素效应。而拮抗剂可以通过与配体结合区域结合，诱导其产生构象变化，使第12号螺旋阻挡共激活因子结合位置，从而产生抗甲状腺激素效应，称为消极的拮抗剂(Passive Antagonists)，相关抗性效应称为消极的抗性效应；也可以通过诱导TR捕获共抑制因子，从而抑制相关基因的转录，称为积极的拮抗剂(Active Antagonists)，相关抗性效应称为积极的抗性效应。在此，核受体介导作用过程中共调节因子(共激活因子、共抑制因子)对活性的产生起了至关重要的作用。因此，我们可以基于以上三种作用方式，通过模拟共调节因子与受体的相互作用，鉴别化合物的拟/抗甲状腺激素效应。

在TR介导的甲状腺激素干扰效应产生过程中，只有配体与TR结合才能产生干扰效应，配体与TR的结合效力是衡量配体甲状腺激素干扰活性大小的重要因子。已有研究采用结合自由能对化合物的干扰效应进行虚拟筛选，如Slynko等用结合自由能筛选蛋白激酶相关激酶1抑制剂，得到很好的预测效果(Slynko I,Scharfe M,Rumpf T,et al.Virtualscreening of PRK1inhibitors:ensemble docking,rescoring using binding freeenergy calculation and QSAR model development[J].Journal of chemicalinformation and modeling,2014,54(1):138-150.)，因此可以通过计算配体与TR结合自由能预测其甲状腺激素活性大小。自由能包括分子力学势能、熵和溶剂化自由能，而结合自由能是整体自由能与配体、受体结合自由能之和的差。其中溶剂化自由能，包括极性溶剂化自由能和非极性溶剂化自由能，是结合自由能计算过程中最难计算也是最耗时的。g_mmpbsa软件集合了gromacs和apbs等软件，采用分子力学/泊松玻尔兹曼表面积(MolecularMechanics/Poisson Boltzmann Surface Area，MM/PBSA)方法，极大简化了结合自由能的计算过程。因此，我们采用g_mmpbsa软件计算结合自由能，达到预测甲状腺激素干扰效应的目的。同时，通过考虑消极的抗性、积极的抗性和拟性效应等甲状腺激素干扰效应作用机制，使预测结果更准确。

于红霞等(于红霞,史薇,王小享.基于分子动力学模拟的核受体介导内分泌干扰物质的虚拟筛选方法:CN,CN103324861A[P].2013.)和张爱茜等(张爱茜,蔺远,彭素芬,刘磊,高常安,韩朔睽.一种有机物雌激素受体激动和拮抗作用的识别方法:CN,CN101381894A[P].2009.)都曾通过核受体H12的稳定情况判断化合物是否具有内分泌干扰效应，然而，这两种方法中都没有引入对核受体作用至关重要的共调节因子，且不能有效预测干扰活性当量。文献检索结果表明，在本发明完成之前，还未发现基于共调节因子作用过程产生消极的抗性、积极的抗性和拟性效应鉴别甲状腺激素干扰物质，或考虑以上作用机制的条件下，用结合自由能预测干扰活性的报道。

发明内容

本发明要解决的技术问题是提供一种基于核受体共调节因子的甲状腺激素干扰物虚拟筛选方法，以解决现有技术存在的筛选和预测不准确等问题。

本发明还要解决的技术问题是提供上述筛选出的甲状腺激素干扰物虚拟的干扰活性的定量计算方法。

为解决上述技术问题，本发明采用的技术方案如下：

一种基于核受体共调节因子的甲状腺激素干扰物虚拟筛选及其干扰活性的定量计算方法，它包括以下步骤：

(1)以已报道的受体晶体结构为模板构建甲状腺激素受体(Thyroid hormoneReceptor,TR)初始结构，并对初始结构进行验证和优化得到甲状腺激素受体；构建并优化配体分子，将化合物配体分子对接到受体的对接口袋中，得到配体-受体复合体；最后将所得配体-受体复合体进行分子动力学模拟(Molecular Dynamics,MD)；

(2)从分子动力学模拟得到的分子运动轨迹中提取平衡状态下的分子构象，即平衡构象；对于促使第十二号螺旋阻挡共激活因子结合位置的配体，判定为消极的甲状腺激素受体拮抗剂；

(3)分别将共激活因子和共抑制因子对接到步骤(2)中提取到的平衡构象中，并通过对接打分判断平衡构象更趋向于与共激活因子结合还是与共抑制因子结合；若平衡构象趋向于与共激活因子结合，则将配体鉴定为甲状腺激素受体激活剂；若平衡构象趋向于与共抑制因子结合，则将配体鉴定为积极的甲状腺激素受体拮抗剂；

(4)对于步骤(3)中鉴定得到的甲状腺激素受体激活剂，将其与共激活因子对接得到的配体-受体-共激活因子复合体进行分子动力学模拟；对于步骤(3)中鉴定得到的积极的甲状腺激素受体拮抗剂，将其与共抑制因子对接得到的配体-受体-共抑制因子复合体进行分子动力学模拟；

(5)选取已有报道的具有甲状腺激素干扰活性的代表性化合物进行活性测试，对其与甲状腺激素受体之间的结合自由能和其活性大小进行回归建模；计算步骤(2)～(4)中筛选出的消极的甲状腺激素受体拮抗剂、甲状腺激素受体激活剂和积极的甲状腺激素受体拮抗剂与受体之间的结合自由能，并代入回归建模所得公式进行计算，即可定量预测筛选出的甲状腺激素干扰物的活性。

步骤(1)中，所述的已报道的受体晶体结构是指能从蛋白质晶体结构库(http://www.rcsb.org/pdb/home/home.do)中搜索并下载的晶体结构，所有研究解析得到并被认可的蛋白质晶体结构都会被收录到该库中。

步骤(1)中，以受体晶体结构为模板构建甲状腺激素受体初始结构的方法为：采用同源建模的方法，第十一号螺旋到第十二号螺旋之间的柔性链以蛋白数据库编号为1A52的受体为模板，其余部分以已解析的甲状腺激素受体晶体结构为模板，构建初始构象。

步骤(1)中，验证方法为：用拉氏图检验结构的合理性。

步骤(1)中，优化方法为：采用Chemi3D软件，将所述的配体分子结构先用分子力学阿林格力场2(Molecular Mechanics,Allinger Force Field version 2，MM2)进行初步优化，再采用鲍威尔(Powell)梯度算法和特里波(Tripos)力场进行再次优化。

步骤(1)中，将配体和受体分别赋予CHARMM(Chemistry at HARvardMacromolecular Mechanics)力场后，将配体-受体复合体浸入TIP3P模型水中，加入钠离子或氯离子平衡体系电荷后，经能量收敛和升温，使体系维持在300K和1个标准大气压下；在蛋白受约束的情况下对配体-受体复合体或配体-受体-共激活因子复合体进行第一次分子动力学模拟优化，再撤销对蛋白的约束，进行第二次分子动力学模拟优化；最后进行20～22纳秒的分子动力学模拟。

步骤(4)中，将配体和受体分别赋予CHARMM(Chemistry at HARvardMacromolecular Mechanics)力场后，将配体-受体-共抑制因子复合体或配体-受体-共激活因子复合体浸入TIP3P模型水中，加入钠离子或氯离子平衡体系电荷后，经能量收敛和升温，使体系维持在300K和1个标准大气压下；在蛋白受约束的情况下对配体-受体-共抑制因子复合体或配体-受体-共激活因子复合体进行第一次分子动力学模拟优化，再撤销对蛋白的约束，进行第二次分子动力学模拟优化；最后进行20～22纳秒的分子动力学模拟。

步骤(2)中，所述的平衡构象由如下方法得到：计算第十二号螺旋相对于初始结构的均方根偏差，将均方根偏差处于稳定状态时的构象作为甲状腺激素受体的平衡构象。

步骤(2)中，所述的平衡构象还可由如下方法得到：提取第十二号螺旋摆动的平均位置为平衡构象。

步骤(2)中，第十二号螺旋是否阻挡共激活因子结合位置的情况通过第十二号螺旋与关键氨基酸的最短距离进行判断；对于甲状腺激素受体α亚型，若第十二号螺旋与V284和K306的距离均小于且与K288和I302的距离均小于则认为第十二号螺旋阻挡共激活因子结合位置；对于甲状腺激素受体β亚型，若第十二号螺旋与V230和K252的距离均小于且与K234和I248的距离均小于则认为第十二号螺旋阻挡共激活因子结合位置。

步骤(3)中，将共激活因子和共抑制因子对接到平衡构象中的方法为：

将参考受体叠合到平衡构象上，将参考受体中共激活因子(或共抑制因子)的位置作为对接时的位置；采用Hex 8.0.0软件基于蛋白质形状进行3D对接；对接时，将受体和共激活因子(或共抑制因子)的偏转范围均设定为15°，最终从得到的500个结果中选取对接打分值E_dock最低的情况作为对接结果。

其中，共抑制因子的参考受体为PDB编号为2OVM的受体，共激活因子的参考受体为PDB编号为4ZO1的受体。

其中，若平衡构象与共抑制因子的对接打分值低于平衡构象与共激活因子的对接打分值之差，则平衡构象趋向于与共抑制因子结合；反之，平衡构象趋向于与共激活因子结合。

步骤(5)中，用于计算消极的甲状腺激素受体拮抗剂的活性的结合自由能为步骤(1)中分子动力学模拟时，消极的甲状腺激素受体拮抗剂与受体之间的结合自由能；用于计算甲状腺激素受体激活剂和积极的甲状腺激素受体拮抗剂活性的结合自由能为步骤(4)中分子动力学模拟时，甲状腺激素受体激活剂和积极的甲状腺激素受体拮抗剂与受体之间的结合自由能。

步骤(5)中，回归建模所得公式如下：

-logRIC₂₀＝0.502-0.024ΔG_pas/act；

其中，ΔG_pas/act为区分消极的拮抗剂和积极的拮抗剂后所得甲状腺激素干扰物与甲状腺激素受体之间的结合自由能的组合，-logRIC₂₀为甲状腺激素干扰物的活性。

有益效果

本发明基于核受体共调节因子调节机制，采用分子动力学模拟和结合自由能计算方法，通过12号螺旋的稳定位置、甲状腺激素受体与共激活因子、共抑制因子相互关系判断受测化合物的甲状腺激素干扰效应，区分积极的抗性、消极的抗性和拟性效应，通过计算配体-受体结合自由能预测相关活性大小。

本发明首次在分子动力学模拟中引入共抑制因子、共激活因子，利用不同作用机制鉴别甲状腺激素干扰物，并用结合自由能量预测相关甲状腺激素干扰活性的大小。

与现有技术相比，本发明具有如下优势：

(1)在分子动力学模拟过程中引入共抑制因子、共激活因子等作用机制，较为全面地考虑甲状腺激素受体作用过程；

(2)利用消极的抗性、积极的抗性和拟性效应鉴别环境中的甲状腺激素干扰物，并在区分不同干扰机制后，用结合自由能预测干扰活性大小；

(3)相比于传统的体外筛查方法，此方法成本低廉、操作简单、效率更高；相比于其他虚拟筛选方法，此方法更全面地考虑了分子作用机制，可以区分拟性和抗性，定性和定量预测结果更可靠。

附图说明

图1为本发明虚拟筛选甲状腺激素干扰物并预测干扰效应的流程图；

图2为实施例1中结合自由能和对接打分结果及其对结合效力预测效果的比较；

图3A为实施例2中，甲状腺激素受体在三碘甲腺原氨酸T3诱导下的结构变化示意图；

图3B为实施例2中，甲状腺激素受体在3-OH-BDE-100诱导下的结构变化示意图；

图4A为实施例2中TRα与共激活因子和共抑制因子的对接打分结果；

图4B为实施例2中TRβ与共激活因子和共抑制因子的对接打分结果；

图5A为实施例2中未考虑分子作用机制下结合自由能与活性的线性回归方程；

图5B为实施例2中考虑分子作用机制后结合自由能与活性的线性回归方程。

具体实施方式

根据下述实施例，可以更好地理解本发明。然而，本领域的技术人员容易理解，实施例所描述的内容仅用于说明本发明，而不应当也不会限制权利要求书中所详细描述的本发明。

下面将结合本发明实施例和附图，对本发明技术方案清楚、完整地描述，显然，所描述的实施例仅是本发明的一个具体的实施例，而不是全部。基于本发明中的实施例，本研究领域其他普通技术人员在没有做出创造性劳动的前提下所获得的其他实施例，都处于本发明保护的范围。

以下实施例所采用的受体为人类和大鼠的TRα、TRβ，配体为羟基化多溴联苯醚(OH-PBDE)。

实施例1

采用同源建模的方法构建人的TRα、TRβ初始结构，主体部分分别以已有文献通过实验得到的人的TRα和TRβ(PDB编号分别为4LNX和1NQ0)为模板，H11与H12之间柔链以PDB编号为1A52的受体为模板，构建所得结构通过拉氏图检验。挑选8个OH-PBDE分子和三碘甲腺原氨酸(T3)为配体，采用SYBYL 7.3中的Surflex-Dock模块将配体对接于受体中形成配体-受体复合体，并计算对接打分值。将复合体用gromacs软件包进行20ns的MD模拟。MD模拟过程为：将配体和受体分别赋予CHARMM力场后，将复合体浸入TIP3P模型水中，加入钠离子或氯离子平衡体系电荷后，经能量收敛和升温，使体系维持在300K和1个标准大气压下；先后在蛋白受约束和蛋白不受约束的情况下进行分子动力学模拟优化，最后进行MD模拟。

在20ns的MD模拟过程中，每隔100ps抽取一个构象，利用g_mmpbsa程序，采用MM/PBSA的方法计算配体与受体的结合自由能。OH-PBDE配体的结合自由能和对接打分值如图2所示，结合自由能越低、对接打分越高说明结合越强，然而对接打分和结合自由能预测的结果有很大差异。

以荧光标记的三碘甲腺原氨酸(T3)为探针进行竞争结合实验，结果TRα和TRβ分别有3个和2个OH-PBDE检测出结合效力，且结合效力远低于T3，而其他OH-PBDE都未检测到效力(图2)。将模拟结果和实验结果比较发现结合自由能的预测效果明显优于对接打分结果，这为我们采用结合自由能预测甲状腺激素干扰活性奠定基础。

实施例2

采用同源建模的方法构建大鼠的TRα、TRβ初始结构，由于目前学者通过实验获得的TR受体只有人类的TR，因此主体部分分别以人的TRα和TRβ(PDB编号分别为4LNX和1NQ0)为模板，H11与H12之间柔链以PDB编号为1A52的受体为模板，构建所得结构通过拉氏图检验。挑选16个OH-PBDE分子和T3为配体，将配体对接于受体中形成配体-受体复合体，将复合体用gromacs软件包进行至少20ns的MD模拟。MD模拟过程为：将配体和受体分别赋予CHARMM力场后，将复合体浸入TIP3P模型水中，加入钠离子或氯离子平衡体系电荷后，经能量收敛和升温，使体系维持在300K和1个标准大气压下；先后在蛋白受约束和蛋白不受约束的情况下进行分子动力学模拟优化，最后进行MD模拟。

计算H12相对于初始结构的RMSD，对于RMSD趋于稳定的情况，将RMSD稳定状态称为TR的平衡状态，并提取代表性构象为平衡构象；对于RMSD不稳定的情况，通过动态变化轨迹的观察，提取H12摆动的平均位置为平衡构象。对提取的平衡构象进行观察，结果发现部分OH-PBDE可能为消极的拮抗剂，如3-OH-BDE-100等能通过诱导H12部分占据共激活因子结合位点AF-2，阻挡共激活因子的结合，从而产生抗甲状腺激素效应(如图3所示)。因此，此类HO-PBDE被鉴别为消极的甲状腺激素受体拮抗剂。

对于暴露共激活/抑制因子结合位置的情况，引入PDB编号为20VM受体中的共抑制因子和PDB编号为4ZO1受体中的共激活因子，将共抑制因子和共激活因子分别对接到受体的共调节因子结合表面。比较共抑制因子和共激活因子的对接结果发现大部分HO-PBDE，如4-HO-BDE-188，作用下，TR趋向于与共抑制因子结合，而非共激活因子，从而抑制转录，导致抗甲状腺激素效应的发生(图4)。因此，此类HO-PBDE被鉴别为积极的甲状腺激素受体拮抗剂。将上述情况得到的积极的拮抗剂-受体-共抑制因子复合体再次进行20ns的MD模拟。MD模拟过程同上。

采用大鼠垂体瘤细胞GH3细胞系的细胞增殖实验对甲状腺激素干扰效应进行定性和定量检测，结果所有16个OH-PBDE都不具有拟甲状腺激素效应，而除了3个细胞毒性太强外，其他13个OH-PBDE均检测出抗甲状腺激素效应，验证了通过模拟鉴别的结果。通过计算得到抗甲状腺激素活性大小，用-logRIC₂₀表示。

分别计算各个配体第一次MD模拟总的结合自由能ΔG_lig-rTRα和ΔG_lig-rTRβ，和加入共抑制因子后的结合自由能ΔG_{lig-rTRα/cor}和ΔG_{lig-rTRβ/cor}(如表1所示)。并将α、β两个受体的结合自由能进行组合，得到不区分消极的拮抗剂和积极的拮抗剂的结合自由能组合ΔG_sum,cor-rTR，和鉴别出消极的拮抗剂和积极的拮抗剂后的结合自由能组合ΔG_pas/act。用origin软件分别与-logRIC₂₀进行回归分析建模，结果ΔG_sum,cor-rTR与-logRIC₂₀的相关系数为0.754，而ΔG_pas/act与-logRIC₂₀的相关系数达到0.834(如图5所示)。说明ΔG_pas/act的预测效果更佳，也说明考虑分子作用机制有助于更好预测活性。

表1.结合自由能和相关的抗甲状腺激素活性

ΔG_pas/act与-logRIC₂₀回归分析建模的结果如下：

-logRIC₂₀＝0.502-0.024ΔG_pas/act

结合上述化合物甲状腺激素干扰效应的鉴别方法，该模型可用于定量预测化合物的甲状腺激素干扰活性。