有机磷类化合物大鼠急性毒性QSAR模型构建与毒性机制研究

郑子廷，闫赛红,#，查金苗,*

1. 中国科学院生态环境研究中心，环境水质学国家重点实验室，北京 100085 2. 中国科学院大学，北京 100049

有机磷类化合物(organophosphorus compounds, OPs)包括有机磷阻燃剂(organophosphorus flame retardants, OPFRs)、有机磷杀虫剂(organophosphorus pesticides, OPPs)和其他含磷的有机化合物。工业生产上，OFPRs作为溴化阻燃剂的低毒替代品被应用于家电、建筑和电子领域；农业上，OPPs则被广泛用于杀虫和害虫控制[1-2]。此外在医学领域，一些抗病毒的药物中也含有OPs[3]。因此，OPs容易通过磨损、挥发和溶解等方式进入各种环境介质中。目前科研人员已在诸如河流、地下水、海水、空气和室内等环境介质中检测到了OPs[4-6]。甚至在东亚地区的人类母乳样本中检测到10余种OPs，最高浓度为70 ng·g-1[7]。越来越广泛存在的OPs可能造成的危害和健康风险引起了研究者们的关注。

OPs具有神经毒性[8-9]、内分泌毒性[10]和致癌性[11]。研究表明，接触OPPs(如毒死蜱)会引起人类的神经系统发育缺陷[12]。OPFRs如磷酸三甲苯(TMPP)、磷酸三苯酯(TPHP)和磷酸三(2,3-二溴丙基)酯(TDBPP)则显示出对核受体(nuclear receptors, NRs)的亲和力和内分泌的干扰作用[13-14]。上述研究及更多OPs的毒性毒理研究往往都是通过传统的体内(invivo)或者体外(invitro)实验进行的，这些方法需要耗费大量的时间与金钱。此外，实验手段无法完整覆盖一整类化合物的毒性数据，并且测试在实验室间存在误差[15]。因此，迫切需要一种有效的替代方法，对化合物毒性数据进行评估。

随着计算科学的发展，计算毒理学已经成为一种具有成本效益且快速的实验替代方法[16]。计算毒理学整合各种来源的数据，开发基于数学、统计科学和机器学习的模型，对化学品进行高效和高通量的定性或定量预测[17]。定量结构活性关系(quantitative structure-activity relationship, QSAR)模型便是计算毒理学领域被广泛利用的模型之一，其原理是建立化合物的化学结构与各种毒理学终点效应和生物活性之间的函数关系[18]。QSAR模型被研究者们广泛应用于药物发现、毒理学和生物学等领域[19]。除了建立化学结构与毒性数据之间的关系模型外，QSAR模型还能通过分子描述符的解释从一定程度进行毒理机制的阐述。Alves等[20]使用QSAR方法为COVID-19病毒建立了抑制性药物预测模型，并应用该模型从小分子数据库中高通量地筛选了11个候选药物，其中3个被美国转化科学推动中心(National Center for Advancing Translational Sciences, NCATS)实验证实为有效的病毒抑制剂。Jeon等[21]则收集大量实验数据，为鱼胚胎急性毒性建立了一个通用的基线毒性QSAR模型，并使用该模型筛选毒性化合物。大量的研究表明QSAR模型用于化合物的毒性预测行之有效。

分子对接是研究化合物与目标蛋白结合部位行为的工具，能够通过计算手段探究化合物与特定靶点蛋白的结合特异性，从而从一定程度上反映化合物的毒性。Zou等[22]使用分子对接得到的配体-受体相互作用能量(binding energy)来反映化学品的毒性作用。Li等[23]通过对接羟基多溴联苯醚(HO-PBDEs)与人类甲状腺激素受体，从一定程度上揭示了HO-PBDEs的毒性机制。

现有的一些针对大鼠急性口服毒性的QSAR模型存在没有使用测试集[24]、选用的分子描述符数量不规范[25]等问题。在本研究中，完善地建立了基于逐步算法的线性回归(SW-MLR)和基于遗传算法的线性回归(GA-MLR)QSAR模型用于预测大鼠的急性口服毒性。研究收集的数据集中含有9种没有毒性实验数据的OPs，最终开发完成的QSAR模型将用于预测这9种OPs的急性毒性。分子对接也将作为一个辅助工具来评估这9种OPs潜在的神经毒性。

1 材料与方法(Materials and methods)

1.1 数据集的收集和划分

本研究首先收集了OPs对大鼠急性口服毒性的数据。检索并从文献中提取了62种环境中常见的OPs清单[26-27]。通过查询EPA CompTox数据库(https://comptox.epa.gov/dashboard)，收集了同一实验条件下获得的53种OPs急性毒性数据LD50(导致一组实验动物中刚好死亡一半的化合物浓度)。获取的OPs的急性毒性数值使用log函数进行规范化处理。整个数据集按照4∶1的比例随机划分为训练集(n=43)和测试集(n=10)，没有毒性实验数据的9种OPs作为最后的预测集(n=9)。在之后的计算实验中，训练集仅用于建立模型，测试集则用于评估模型的预测能力，这种划分方式有利于验证QSAR模型的拟合和泛化能力。

1.2 分子描述符的计算

数据集中所有OPs的结构文件(.sdf)来自PubChem(https://pubchem.ncbi.nlm.nih.gov/)数据库。并使用PaDEL开源软件[28]计算所有化合物的分子结构描述符，共得到1 444个分子结构描述符，包括(1)一维结构描述符：如原子和化学键的数量、分子量等；(2)二维拓扑描述符：如拓扑指数、拓扑电荷和自相关系数等；(3)电荷描述符、原子电状态描述符等。本研究中使用了PaDEL计算的一维和二维描述符，因为这些类型的描述符计算速度快，可重复性好，并且不需要分子结构能量最小化过程，该过程中不同实验室间可能会在计算上产生误差。计算结束后，我们对得到的分子描述符数据进行了预处理，剔除了方差为0的单列分子描述符数据，相关系数>0.99的多组描述符最终仅保留一组。经过预处理后，最终保留了833组分子描述符用于后续分析。

1.3 模型的构建

分别使用逐步算法(SW)和遗传算法(GA)进行分子描述符的选择。2种算法选择的分子描述符变量通过基于普通最小二乘法(OLS)的多元线性回归(MLR)方法关联毒性终点并建立QSAR模型，之后进行模型的内部和外部评估。研究中使用SPSS 21.0软件建立了SW-MLR模型，使用DTC-QSAR开源软件[29]建立了GA-MLR模型。所有参数的统计分析过程在EXCEL软件中完成。本研究详细的QSAR模型构建及OPs的LD50预测工作流程如图1所示。

图1 定量结构活性关系(QSAR)模型构建及有机磷类化合物(OPs)的急性毒性半致死量(LD50)预测的工作流程图注：GA-MLR表示基于遗传算法的多元线性回归模型；SW-MLR表示基于逐步算法的多元线性回归模型；R2表示决定系数；MAE表示平均绝对误差。Fig. 1 Workflow diagram for quantitative structure-activity relationship (QSAR) models construction and acute toxicity semi-lethal dose (LD50) prediction of organophosphorus compounds (OPs)Note: GA-MLR represents the multiple linear regression model based on genetic algorithm; SW-MLR represents the multiple linear regression model based on stepwise algorithm; R2 represents the coefficient of determination; MAE represents mean absolute error.

1.4 模型评估

1.5 模型应用域

使用Williams图(又称杠杆方法)[30]对建立的QSAR模型进行应用域分析。基于帽子矩阵(一类投影矩阵)得出的值(h)与标准化残差(δ)的比较。h大于警告杠杆值(h*)的化合物被认为是结构异常值，标准化残差>2.5的化合物则被认为是计算异常值。Williams图的公式定义如下:

(1)

表1 本研究中用于模型评价的公式及参数Table 1 Equations and parameters used for model evaluation in this study

(2)

(3)

1.6 分子对接

我们先前的实验研究证实了稀有鮈鲫的丁酰胆碱酯酶(BChE)能够被一些OPs显著抑制，从而引起神经毒性[31]。这项研究中，借助分子对接工具探究预测集中9种OPs对哺乳动物类的BChE的结合与抑制能力。从蛋白质数据库中检索到BChE的人类蛋白质结构(PDB ID: 5K5E)作为对接蛋白。通过AutoDockTools4程序进行预处理，去除水分、添加氢并分配电荷。对接网格的大小设置为52 Å×52 Å×52 Å，并使用AutoDock Vina程序[32]进行分子对接。对接前，从BChE蛋白中提取了共结晶化合物并重新对接至原蛋白进行对接程序的验证。最终使用PLIP工具[33]分析了对接后化合物与蛋白质的结合模式。

2 结果(Results)

2.1 预测模型

逐步算法的参数设置为进入F值：<0.010，逐出F值：>0.050。建立的SW-MLR模型如下：

logLD50=-0.154(±1.5150)+0.008(±0.0018)ATS3s-2.783(±0.5247)GATS3e+2.054(±0.4083)MATS8i-88.853(±13.1459)JGI7-2.871(±0.6230)MATS3i+1.257(±0.2907)maxssO-(±0.0553)0.180topoDiameter

遗传算法的参数设置为最大迭代次数：100，方程参数数量：7，初始方程数：100，每代方程数：30，适应函数：MAE。建立的GA-MLR模型如下；

logLD50=5.608(±0.7354)+0.0067(±0.0049)AATSC7m-0.012(±0.0088)ATSC3i+0.0133(±0.0087)VR3_Dzp-3.1038(±0.9454)AVP_4-85.9765(±14.6601)JGI7+0.263(±0.1021)nsssCH-1.2293(±0.3834)GATS4p

对于每个模型还执行了2 000次Y-随机测试。随机测试结果显示SW-MLR和GA-MLR的最大伪R2值分别为0.488和0.504，最大伪Q2值分别为0.299和0.234。所有的伪R2和伪Q2都很低，这表明原始模型的结果并不是偶然关联的。SW-MLR和GA-MLR模型的拟合关系结果如图2所示。

图2 SW-MLR(a)和GA-MLR(b)模型中OPs实验值与预测值拟合关系Fig. 2 Fitting relationship between experimental and predicted values of OPs in SW-MLR (a) and GA-MLR (b) models

2.2 分子描述符解释

在SW-MLR方程中，分子描述符ATS3s、MATS8i和maxssO与毒性终点logLD50正相关。ATS3s、MATS8i分别是基于原子的内在状态和电离电位来表征化合物拓扑结构的分子描述符[34]。maxssO则代表分子结构中—O—基团的最大电性拓扑状态，前人的研究中也选用maxssO来描述一些无毒的抗疟化合物[35]，这与我们研究中maxssO描述符在SW-MLR方程中的正系数一致。分子描述符GATS3e、JGI7、MATS3i和topoDiameter则与毒性终点logLD50负相关。GATS3e和JGI7描述符与化合物中高电负性的原子相关[36-37]。MATS3i和topoDiameter是特定拓扑结构的分子描述符。

在GA-MLR方程中，分子描述符AATSC7m、nsssCH和VR3_Dzp与毒性终点logLD50正相关。AATSC7m和nsssCH分别对应OPs分子结构中R基上支链的长度[38]、支链或取代基烷基化的程度[39]。VR3_Dzp则是比较抽象的二维描述符，代表极化率加权的Barysz矩阵的特征向量[40]。分子描述符ATSC3i、AVP_4、JGI7和GATS4p则与毒性终点logLD50负相关。这4个描述符与分子的疏水性、极性和原子的电负性相关[36,41]。对于2个模型分子描述符的解释从一定程度上揭示了OPs的化学结构与毒性相关的内在因素，对于识别OPs的毒性起到关键作用。

2.3 模型应用域

建立的SW-MLR模型和GA-MLR模型应用域表征Willimas图如图3所示。SW-MLR模型中，所有OPs类化合物都很好地落入了应用域范围内。GA-MLR模型中，仅有一种化合物(dimethoate, CAS: 60-51-5,h=0.570>h*=0.558)为离群点。结论表明本研究的2种模型能够较好地预测OPs类化合物对大鼠的急性毒性数据。表2中给出了使用2种QSAR模型预测的没有实验数据的9种OPs的急性毒性结果。

图3 SW-MLR(a)和GA-MLR(b)的Williams应用域图Fig. 3 Williams plotting for SW-MLR (a) and GA-MLR (b)

2.4 模型比较

2.5 分子对接

AutoDock Vina重新对接BChE(5K5E)的原配体的结合能结果为-11.9 kcal·mol-1，均方根偏差(RMSD)为0.89 Å，说明对接程序适用于进一步的对接与研究。使用AutoDock Vina软件对接预测集(9种OPs)与BChE(5K5E)的最低对接能量结果如表2所示。PLIP工具对OPs与BChE的结合模式可视化如图4所示。

AutoDock Vina对接的结果显示其中8个OPs的对接结合能低于阈值[43]-6 kcal·mol-1，表明了潜在的结合能力以及进一步可能造成的神经毒性。此外，取代基或支链烷基化程度较高的化合物(如磷酸三(3-氯丙基)酯)对接结合能也较高，这表明更低的结合能力与更低的神经毒性，这一结果与GA-MLR模型中分子描述符nsssCH的含义一致。因此我们认为OPs的R基及取代基的烷基化能够降低毒性。由图4可知，BChE与OPs结合的关键氨基酸残基有Gly116和Ser198(氢键)；His428、Trp332和Phe329(盐桥)；Trp231、Phe329、Trp82、Tyr332和His438(π-堆叠)。这些关键氨基酸残基与前人的研究较为一致[44-45]。

3 讨论(Discussion)

研究采用了SW和GA算法，基于PaDEL分子描述符，对OPs化合物的大鼠口服急性毒性建立了多元线性回归QSAR模型，SW-MLR模型具有更好的拟合能力，GA-MLR模型有着更好的预测能力。对模型分子描述符的分析显示了OPs的R基的高烷基化和支链的延长能够降低毒性。此外，通过模型预测了数据集中没有实验数据的9种OPs的毒性，并以分子对接方法检测其与BChE的潜在结合能力。所构建模型能够方便有效地预测OPs类化合物的毒性，并为新化学品的合成与环境保护政策的制定提供帮助。

表2 SW-MLR和GA-MLR对9种OPs的LD50预测值、没有实验值的9种OPs与丁酰胆碱酯酶(BChE)对接的能量与关键残基Table 2 Predicted LD50 values by SW-MLR and GA-MLR for nine OPs, docking energy and key residue results by docking to butyrylcholinesterase (BChE)

图4 9种OPs与BChE(5K5E)分子对接结果的可视化注：A～I对应表2中化合物序号；疏水性相互作用为灰色虚线，平行π-堆叠为亮绿色虚线，垂直π-堆叠为深绿色虚线，氢键为蓝色实线，盐桥为黄色虚线。Fig. 4 Ligand interaction diagram of the main interactions of nine OPs in the active site of BChE (5K5E)Note: A～I represents the compounds in Table 2; hydrophobic interactions are represented by gray dashes; π-stackings are represented by bright green (parallel) and dark green (perpendicular) dashes; hydrogen bonds are represented by blue solid lines; salt bridges are represented by yellow dashes.

通讯作者简介：查金苗(1975—)，男，博士，研究员，主要研究方向为水生模型生物体系的构建与发展、水环境生物毒性测试方法、环境内分泌干扰物的筛选技术研究、环境污染物对水生生物分子毒理机制和水生态系统完整性评估方法等。

共同通讯作者简介：闫赛红(1988—)，女，博士，助理研究员，主要研究方向为水生态毒理学。