张文灏,陈景文,徐童,王雅
工业生态与环境工程教育部重点实验室,大连理工大学环境学院,大连 116024
外源性化合物可通过鱼的摄食、呼吸和直接接触等途径进入到体内,经过吸收和代谢竞争,在体内蓄积[1-2]。蓄积的化学品可对鱼体产生毒害效应[3-5],因此有必要评价外源化合物的生物蓄积效应[6-8]。经济合作与发展组织(OECD)在2012年发布的“鱼体内生物蓄积:水和食物暴露”导则[9]指出,除生物富集因子(BCF)、生物放大因子(BAF)外,化合物在鱼体内的半减期(t1/2, d)和消除速率常数(KM, 1/d)也可用于评价化学品的生物蓄积效应。为方便统计和比较,通常取logt1/2进行记录和计算[10-12]。目前具有logt1/2实测值的化学品仅有几百种,实验测定t1/2的速度慢、成本高,难以满足化学品生态风险评价的需求,需要发展替代实验的模型预测方法。
定量构效关系(QSAR)模型可用于预测化合物的logt1/2值[13]。Arnot等[14]基于632种化合物在鱼体的实测logt1/2值,采用正辛醇/水分配系数(logKow)、分子量(Mw) 2种描述符,以及57个分子碎片构建了QSAR模型。在此基础上,Brown等[15]构建了包含34个分子碎片及logKow和Mw2种描述符的QSAR模型。Papa等[6]运用同一数据库,建立了包含9个2D分子描述符的QSAR模型。然而,这些模型的训练集和验证集中,均不包含近年来引起广泛关注的药物和个人护理用化学品(PPCPs)类物质。
近年来,PPCPs类污染物在各处水体和水生生物中被检出[16-18],具有潜在的生态风险[19-20]。如果QSAR模型能够预测包含抗抑郁药、降压药、麻醉剂、抗过敏药、抗病毒药和抗生素在内的PPCPs在鱼体内的t1/2,则模型有助于评价PPCPs类化合物的危害性和风险。本研究在以往数据库的基础上,整理搜集了包括PPCPs在内的653个化合物在鱼体内的logt1/2实测值,采用多元线性回归(MLR)[21-23]和支持向量机(SVM)[24-27]2种方法分别建立logt1/2的预测模型,并对模型进行验证、应用域表征和机理解释。
本研究从文献[28-46]和数据库(EPI Suite Package)中,共搜集653个化合物在鱼体内的logt1/2数据,涉及鱼种类包含鲤鱼、虹鳟、斑马鱼、罗非鱼、鲑鱼、鲈鱼、太阳鱼和青鳉鱼等十几种鱼类,它们体长、体重、身体构造、生长习性均不相同,实验温度、pH等实验条件也不尽相同。KM和t1/2关系如式(1)所示:
t1/2= ln2/KM
(1)
为了弱化鱼种和实验条件对模型的影响,本研究采用Arnot等[14]的方法,通过式(1)得到新收集数据的消除速率常数(KM,X,单位为d-1)值,再采用式(2)对数据进行规范化处理,即:
KM, N=KM, X(WN/WX)-0.25exp[0.01(TN-TX)]
(2)
这里(WN/WX)为规范化体重(0.01 kg)与实际体重之比;(TN-TX)为规范化温度(15 ℃)与实际温度的差。据式(2)得到规范化处理的消除速率常数(KM, N),最后根据式(1)重新得到t1/2进行建模。经过规范化处理后,得到含有653种化合物的新数据集,包括多环芳烃、多氯联苯、多溴联苯醚、有机磷农药和药物等典型化合物以及其他烷烃、环烷烃、烯烃、醇、醚、酸、酯、酮、卤代化合物、芳香族化合物、含硫、氮、磷化合物。以2∶1的比例将数据集随机拆分成训练集(n= 436)和验证集(n= 217)。
将得到的化合物3D结构使用Gaussian 09软件包[47]中的B3LYP/6-311+G(d, p)方法进行结构优化,其中I原子采用Lanl2DZ赝势基组[48]。基于优化后的分子结构使用DRAGON (6.0)软件[49]计算分子结构描述符,去掉常数项、近常数项以及数据不完整的描述符,得到2291个描述符。
采用杠杆值(h)和标准残差(δ)做Williams图进行模型应用域表征。h和δ计算公式如下:
(3)
h*=3(k+1)/n
(4)
(5)
MLR法建立的logt1/2的QSAR模型为:
logt1/2= 0.96 + 0.064MLOGP2 -0.09Mor04u+ 0.037RDF045p-0.045CATS2D_07_ll-2.82(R1e+) + 0.206Mor16m-0.139NaaaC-1.64SpMaxA_B(s)-0.11ATS7s+ 1.36B06[N-P] -0.728nP04
模型评价参数表明,2个模型均具有良好的预测能力和稳健性。对于MLR模型的外部验证也表明模型具有良好的外部预测能力。各描述符的VIF值均小于10,表明模型不存在多重共线性。描述符含义及VIF值见表1。MLR和SVM预测模型实验值和预测值拟合图如图1所示。
表1 logt1/2的QSAR模型中描述符意义、VIF值及t值Table 1 Explanation of descriptors, VIF and t values in the QSAR model of logt1/2
注:VIF值表示方差膨胀因子,t值表示t检验值。
Note:VIFstands for variance inflation factor;tvalue stands fort-test value.
图1 多元线性回归(MLR)和支持向量机(SVM)模型中logt1/2实测值和预测值拟合关系Fig. 1 Plots of the experimental versus predicted logt1/2 values by multiple linear regression (MLR) and support vector machine (SVM) models
2种建模方法表征应用域的Willimas图如图2所示,MLR模型中,邻苯二甲酸二异壬酯和多菌灵2种化合物为离群点,来自验证集。SVM模型中有6个离群点,包括训练集化合物三聚氰胺、异丙隆、六氯丁二烯、五氯苯甲醚、2-乙基己基乙烯醚和验证集化合物多菌灵。
2种建模方法中离群点化合物共7个,其中2个醚类化合物,数据集中醚类化合物共有134个,说明本研究模型可以预测大部分含-C-O-C-结构的化合物。除此之外,数据集中包含了10种邻苯二甲酸酯类化合物,只有一种未被准确预测,说明对大多数邻苯二甲酸酯类化合物具有较好的预测效果。多菌灵可以与无机酸反应生成盐,本研究中多菌灵实验数据来自鲈鱼,其为有胃鱼,可以分泌盐酸与多菌灵反应生成盐。同样,三聚氰胺在生物体内容易水解生成三聚氰酸等化合物,因此,参与体内循环的化合物并非本体化合物,进而导致其预测结果不准确。
表1给出了模型中涉及的11个描述符意义、VIF值及t值。从表1中数据可以看出,MLOGP2的t值明显大于其他描述符,说明MLOGP2是该模型中十分重要的描述符,这与前人的研究结果一致[14]。MLOGP2和CATS2D_07_ll2种描述符与化合物的疏水性(亲脂性)相关,前者与logt1/2正相关,后者为负相关。Mor04u和Mor16m为3D-MoRSE描述符[49],前者直接表征分子结构,与logt1/2负相关,后者基于质量表征分子结构,与logt1/2正相关。SpMaxA_B(s)是与分子原子连接有关的拓扑描述符,与logt1/2负相关。NaaaC表示::C:结构的数量[6],指的是苯环对接处C原子个数,与logt1/2负相关。R1e+与分子尺寸和电负性相关,而ATS7s也与分子的尺寸有关。有研究表明,分子的体积大小对其在生物体内的吸收分布具有显著影响[50]。RDF045p与logt1/2正相关,经过统计分析,硅氧烷、环烷烃以及含有2个及以上苯环的长链化合物的RDF045p值较大。B06[N-P]表示在拓扑距离6时,是否存在N-P结构,存在值为1,不存在为0,与logt1/2正相关。nP04表示分子中磷酸盐或者硫代磷酸盐基团的个数,与logt1/2正相关。
图2 MLR和SVM模型的Williams图Fig. 2 Williams plots of the MLR and SVM models
在所有描述符中,共有75个化合物NaaaC值不为零,均来自蒽、菲、萘、芴和噻吩等多环芳烃类化合物及其类似物,多氯联苯类化合物以及多溴联苯醚。这些化合物往往有2个以上的苯环相连,该描述符的加入,更准确地表达了这一结构特征对t1/2的影响。
比较前人研究[6,14-15]和本研究2个模型的统计学参数,结果汇总于表2,显然非线性模型预测效果优于线性模型。
与前人研究模型相比,本研究模型数据集增加了药物类化合物,使模型预测范围更广泛。其次,在Arnot等[14]和Brown等[15]的研究中,存在14个预测效果不好的化合物,在Papa等[6]的研究中部分化合物的预测结果得到优化。这些化合物包含在本研究模型的训练集和验证集中,表3给出了14种化合物的logt1/2实验值以及在此前各个模型中的预测最优值,可以看出,除了2,3,4,5-四氯硝基苯和六氯苯,其余化合物的预测准确度都有不同程度的提升,且SVM模型明显优于其他几种线性模型。另外,本研究的数据集中共包含28个氟化物,此前的研究[6,14-15]对这些化合物很难准确预测,而在本研究中,MLR法建立的模型中14种氟化物|δ|<1,预测较为准确,SVM法中21种氟化物也得到较准确的预测。
表2 不同预测模型对比Table 2 Comparison of different models
注:M代表描述符个数,N代表模型包含数据个数;R2表示校正后决定系数,RMSE表示均方根误差;GA表示遗传算法。
Note:Mrepresents the number of descriptors;Nrepresents the number of data in the model;R2represents the adjusted determination coefficient;RMSErepresents the root mean squared error; GA represents genetic algorithm.
表3 各模型对14种化合物的logt1/2预测值Table 3 Prediction of logt1/2 for 14 chemicals from different models
综上所述,本研究运用MLR和SVM这2种方法,基于Dragon分子描述符,构建了鱼体logt1/2值的QSAR预测模型,SVM模型具有更良好的预测能力和稳健性。结果表明,以下几种类型化合物不容易被生物体代谢,鱼体内生物半减期较长:分子体积大、亲脂性高的化合物,硅氧烷、环烷烃、含有2个及以上不直接相连苯环的复杂长链化合物以及含氮、磷的化合物。所构建的模型可以预测多环芳烃、多氯联苯、多溴联苯、多溴联苯醚、农药、药物以及其他烷烃、环烷烃、烯烃、醇、醚、酸、酯、酮、卤代化合物、芳香族化合物、含硫、氮、磷化合物等的鱼体内生物半减期。