有机化学品与臭氧反应速率常数的定量预测模型研究

范德玲，汪 贞，王 蕾，周林军，古 文，刘济宁，石利利

(生态环境部南京环境科学研究所，江苏 南京 210042)

化学品在给人们生产生活带来便利的同时，也给人体健康带来极大威胁。大量的挥发性有机化学品可通过直接排放或从土壤和水相挥发进入大气层。大气中有机化学品可通过物理过程移除，如干、湿沉降；也可通过化学过程降解，如直接光解或与大气氧化剂(OH、NO3和臭氧)反应。在大气对流层中有机物与臭氧的反应是其转化的重要途径。表征有机化学品与O3自由基反应的速率常数(KO3，cm3·mol-1·s-1)是反映有机污染物在大气中持久性能力的重要参数，是进行有机污染物生态风险评价的基础指标[1-2]。然而，目前KO3的实验数据较少，且实验耗时,费力,成本高，不能满足有机化学品生态风险评价的需求。因此，有必要开发出快速有效预测KO3的方法[3]。目前，定量结构-活性关系(QSAR)模型是用于获取KO3的重要方法。为指导各国构建满足化学品风险管理需求的QSAR模型，经济合作与发展组织于2007年发布了QSAR模型构建与验证的导则(简称导则)[4]，导则全面阐述了符合管理要求的QSAR模型应满足的标准。

目前，已有研究开发了关于臭氧反应速率常数的QSAR模型[5-6]。2007年，REN等[7]报道了116种有机化合物基于KO3的QSAR模型，该模型采用DUPLEX分类算法划分模型训练集和测试集，采用启发式方法(Heuristic Algorigthm)筛选最优描述符，并采用多元线性回归、支持向量机和投影寻踪回归方法构建预测模型。GRAMATICA等[8]基于遗传算法筛选最优描述符，采用多元线性回归算法构建了125种有机化合物基于KO3的QSAR模型，结果显示留一法交叉验证系数(QLOO2)达到82%～88%，外部验证决定系数(QEXT2)达到90%，均方根误差(RMSE)达到0.73。由美国国家环境保护局有毒物质污染防治办公室和Syracuse Research Corporation (SRC)公司共同开发的EPI(estimation programs interface)Suite软件，采用基团贡献法构建了112个烯烃和炔烃不饱和有机物基于KO3的QSAR模型，相关系数达到0.94，绝对平均残差达到0.35[9]。但上述模型并不满足导则要求，如缺少模型稳健性和预测能力表征，或未定义模型应用域[6-9]，不利于模型使用者评估需预测的有机化合物是否处于模型应用域内。因此，根据导则要求采用简单透明的遗传算法-多元线性回归(GA-MLR)算法构建基于KO3的新QSAR模型，并对模型进行拟合优度、稳健性、预测能力、应用域表征和机制解释。所构建的模型对实现环境行为参数预测软件化具有重要意义。

1 材料与方法

1.1 数据来源与处理

烷烃、烯烃、芳香烃、含氧挥发性有机物和酚类等152种有机化学品的KO3数据来源于文献[10]。选择-lgKO3作为模型指标。为避免样本分布不均匀，采用KENNARD等[11]分组方法将数据集划分为训练集和验证集，将结构差异较大的样本选入训练集，其他与之相近的样本选入验证集，从而使代表性样本全部进入训练集。训练集有107种化学品，验证集有45种化学品。

1.2 分子结构描述符的计算

分子结构描述符是用于反映分子结构信息的参数，根据分子结构按照一定理论或规则计算得到。笔者构建的模型采用的分子结构描述符为Dragon描述符。采用Hyperchem 7.0软件中MM+和AM1方法对选取的152种有机化学品结构进行优化[12]，采用Dragon 5.4软件计算优化后结构的描述符[13]，并对得到的1 664个描述符进行初步筛选，去掉常数项、近似常数项和高度相关的分子结构描述符，最终得到488个分子结构描述符。

1.3 模型的建立

采用MobyDigs软件中遗传算法选择与-lgKO3高度相关的描述符[14]。由遗传算法变量筛选得到最优描述符，并采用多元线性回归(MLR)方法构建预测模型。遗传算法相关参数设置为种群大小(population size)为100，初始模型允许的最大变量数(maximum allowed variables)为7，变异均衡值(mutation trade-off，T)为0.5，交叉(crossover)和变异(mutation)概率均基于T值。当增加变量数目对结果影响不大时，得到8个最优描述符。

1.4 模型的表征与评价

根据导则要求，对构建的QSAR模型进行内部验证(训练集的拟合优度和稳健性评估)和外部验证(验证集的预测能力评估)。采用实验值与预测值之间校正后的决定系数(Radj2)和均方根误差(RMSE，ERMS)表征模型拟合优度，采用留一法交叉验证系数(QLOO2)表征模型稳定性，采用外部检验参数(QEXT2)、验证集相关系数(REXT2)和验证集均方根误差(ERMS，EXT)等外部验证决定系数表征模型预测能力，基于杠杆值(leverage，hi)的Williams图定义模型应用域[15]。外部验证决定系数计算公式为

(1)

(2)

(3)

(4)

Williams图是标准残差(δ)和hi值定义的模型应用域，其计算公式为

(5)

hi=xiT(XTX)-1xi。

(6)

式(5)～(6)中，xi为第i种化合物分子结构描述符的行向量；X为n×m的矩阵，构成训练集化合物的描述符空间。

当训练集中化合物hi值大于警戒值(h*)时，说明在数据集中该物质的子结构出现较少，会对模型预测结果有显著影响。h*值计算公式为

h*=3(m+1)/n。

(7)

2 结果与讨论

2.1 QSAR模型的构建

模型描述符意义、回归系数、回归系数偏差和标准回归系数见表1。构建的GA-MLR回归方程为Y=17.898-0.371X1+0.334X2+0.215X3+0.193X4+0.426X5-0.453X6-0.260X7-0.308X8，n训练集=107，Radj，训练集2=0.784，QLOO2=0.744，ERMS，训练集=1.127，P<0.000 1，nEXT=45，REXT2=0.664，QEXT2=0.761，ERMS，EXT=1.039。

GOLBRAIKH等[20]研究认为，QSAR模型可接受标准为Q2>0.50和R2>0.60。由图1可知，笔者构建的模型拟合优度和稳健性较好，预测能力也较好。

表1 臭氧自由基反应速率常数模型描述符物理化学意义和相应系数

Table 1 Physical chemistry meaning and corresponding coefficients of the descriptors used in the -lgKO3MLR model

变量描述符定义回归系数回归系数偏差标准回归系数常数项17.8980.480X1PW3Randic形状指数-7.9701.471-0.371X2HOMA分子芳香性指数2.4600.4020.334X3RDF035u径向分布函数描述符0.0830.0260.215X4G1s原子电拓扑加权指数1.4550.3800.193X5HATS2e加权原子Sanderson电负性信息2.1150.2800.426X6Nr=Cs脂肪族化合物分子中仲碳原子个数-0.9240.107-0.453X7Nr=Ct脂肪族化合物分子中叔碳原子个数-1.0950.247-0.260X8H-050与杂原子相连的氢原子个数-1.4070.228-0.308

图1 臭氧自由基反应速率常数MLR模型的预测值与实验值的拟合图

2.2 应用域表征

利用杠杆方法制作Williams图分析和评价模型应用范围，可以图形方式量化模型应用范围。模型对应用域内物质预测性能较好，而对应用域外物质预测性能差。采用Williams图表征的QSAR模型应用域见图2。

虚线为警戒值(h*=0.252)。

由图2可知，数据集152种化合物中只有肼的h值≥h*(h*=0.252)，位于应用域范围外，为X离群点。所有化合物标准残差在-3～3范围内，即无Y离群点。因此，构建的QSAR模型可用于预测应用域内其他化合物-lgKO3值。

2.3 有机化学品与臭氧反应机制解释

通过解释线性化合物臭氧反应速率的QSAR模型中所选描述符的物理化学意义，可以获得决定化合物臭氧反应速率的结构信息。描述符的相对重要程度由模型中每个描述符的标准回归系数决定。标准回归系数绝对值的大小表示对应描述符对臭氧反应速率影响程度的强弱，正负号表示对应描述符与臭氧反应速率相关性的正负。在模型的8个描述符中，nR=Cs和nR=Ct均为官能团数目描述符，且标准回归系数均为负值(表1)，这表明nR=Cs和nR=Ct与臭氧的反应速率常数呈正相关。HATS2e为GETAWAY(geometry，topology and atom-weights assembly)类描述符，与-lgKO3呈较大负相关。HOMA为几何描述符，PW3为拓扑描述符，可通过计算分子中每个原子的path数目与walk数目的比值，再将这些比值求和后除以分子中的原子数目得到。由于path/walk独立于分子大小，所以PW3可以较好地表征分子形状。RDF035u为径向分布函数描述符，表示在一个半径为R的球形体内发现特定类似原子的概率。G1s为WHIM描述符，在模型中表征分子静电拓扑状态。H-050为以原子为中心的碎片描述符，表征与杂原子相连的氢原子个数。

2.4 模型验证和应用案例

近年来，计算毒理学技术在欧美、日本和OECD得到大力发展。美国国家环境保护局研发了化学品理化性质/环境行为指标参数与预测模型软件EPI Suite，其中的AOPWIN模块采用基团贡献法预测有机化学品臭氧自由基反应常数。OECD允许使用QSAR方法弥补数据缺失，并于2008年发布了第1版QSAR Toolbox工具包。其中的臭氧反应速率数据主要来源于EPI Suite软件。与发达国家相比，我国在计算毒理学技术研发和应用方面具有较大差距。近年来我国已经启动化学品环境安全信息预测技术研究，在一定程度上填补了我国化学品固有属性预测技术的空白。其中，生态环境部南京环境科学研究所基于简化分子线性输入规范(SMILES)解析碎片拆分技术，开发了具有我国自主知识产权的化学品定量构效预测软件[17]。

将笔者研究数据集之外的20种有机化学品-lgKO3实验数据[18-19]与该模型和EPI Suite软件中AOPWIN模块的预测结果进行比较发现，20种有机化学品-lgKO3实验值与笔者模型预测值的决定系数(R2)达到0.794，与EPI Suite预测值的R2为0.695(表2)。其中该模型15种化学品预测结果优于EPI Suite，EPI Suite软件5种化学品预测结果较好。SMILES是化学物质1维结构的线性表达，而2维和3维结构描述符可更全面地表达化学物质立体结构的空间形态。由于EPI Suite软件基于SMILES码碎片拆分，选取的结构碎片也许不能完全表达分子结构信息，同时也未给出模型应用域，因此笔者构建的模型弥补了EPI Suite软件的不足[20-21]。

表2 20种有机化学品-lgKO3实验和预测数据比较

Table 2 Comparison of predicted results with experimental results for -lgKO3of 20 organic chemicals

CAS编号实验值EPI Suite软件该模型预测值残差预测值残差000096-33-317.978 17.756 0.22218.078 -0.100 000116-14-319.036 19.348 -0.312 18.755 0.281 000123-73-918.045 17.739 0.30418.174 -0.129000497-23-418.657 16.943 1.714 17.784 0.873001630-77-918.585 18.853 -0.26818.4740.111001630-78-017.677 18.552 -0.875 18.474 -0.797 017559-81-817.744 17.195 0.549 17.724 0.020 018409-46-617.001 17.410 -0.409 17.018 -0.017 000074-86-219.318 19.522 -0.204 19.124 0.194 000075-01-418.619 18.602 0.017 18.616 0.003 000075-38-718.721 18.552 0.16918.341 0.380000078-94-417.321 17.325 -0.00417.289 0.032000108-05-417.494 17.756 -0.262 18.396 -0.902 000140-88-517.244 17.756 -0.512 17.213 0.031000359-11-518.853 18.950 -0.09718.837 0.016 000463-51-419.154 19.057 0.09719.181 -0.027006728-26-317.698 17.739 -0.04117.687 0.011 000109-92-215.812 17.057 -1.24516.022 -0.210000087-44-513.935 15.354 -1.41915.652-1.717 000360-89-420.167 18.552 1.61518.417 1.750

3 结论

该研究建立了包括烷烃、烯烃、芳香烃、含氧挥发性有机物和酚类152种有机化合物与臭氧反应速率常数预测模型。根据经济合作与发展组织关于QSAR模型构建与验证的导则要求，构建的有机化学品与臭氧反应速率预测模型拟合能力、稳健性和预测能力均较好，Williams图定义模型应用域(AD)结果也表明该模型应用域较广。模型机理研究结果表明分子芳香性、电负性和仲碳原子数目是影响有机化学品与臭氧自由基反应速率(KO3)的关键因素。综上所述，构建的有机化合物与臭氧自由基反应速率常数QSAR模型可以用于预测应用域内难以测定或未知有机化合物与臭氧自由基反应速率常数，评估其持久性，进而对有机污染物进行生态风险评价。