CN103425872A - 通过定量构效关系模型预测大气中有机物与羟基反应速率常数的方法 - Google Patents

Abstract

本发明公开了一种预测大气中有机物与羟基反应速率常数(k_OH)的模型方法。在获得化合物分子结构的基础上，通过多元线性回归方法，构建了298K和温度依附性的k_OH模型，可快捷、高效地预测有机化合物在298K或其它温度下k_OH值；还依照经济合作与发展组织(OECD)关于QSAR模型构建和验证导则，对模型拟合优度、稳健性和预测能力进行了表征；模型应用域明确，涵盖不同种类化合物近900种；预测模型简明、易于程序化，能为有机化学品的环境持久性评价提供重要数据支持。

Description

通过定量构效关系模型预测大气中有机物与羟基反应速率常数的方法

技术领域

本发明涉及一种通过建立定量构效关系模型(QSAR)预测大气中有机物与羟基反应速率常数的方法，属于生态风险评价测试策略领域。

背景技术

与大气中羟基自由基(·OH)反应，是有机化合物在大气中去除的重要途径之一。有机化合物与羟基反应速率常数(k_OH)是表征有机化合物与·OH反应强度与能力的参数，其值越大说明有机化合物与·OH反应能力越强，化合物越容易通过与·OH反应而得以去除。因此，k_OH常被用来评估化学品环境持久性，是化学品生态风险评估的一个基础数据。

k_OH可通过激光闪光光解法，放电流技术等实验直接测定，也可通过相对速率方法间接测定。尽管近几年开展了大量的实验工作用于测定有机化合物的k_OH值，但是由于实验方法具有耗时、昂贵，依赖于实验仪器和设备等缺点，目前仅900余种化学品具有可获取的气相k_OH数据。截止2013年7月份，美国化学文摘社(Chemical Abstracts Service,CAS)登记数据库包含了7千万多种有机和无机化学物质，其中大约14万种有机化学品已投入商业化使用。这表明具有k_OH数据的化学品仅占需要进行环境持久性评估的化学品的一小部分。大部分化学品缺少k_OH数据，将不利于对其进行环境持久性评估，从而不能满足对有机化学品风险评价和管理工作的数据需求。因此，通过发展非实验技术获取化学品环境迁移转化行为参数如k_OH数据是社会发展的需要。

基于分子结构性质来预测化合物的理化性质及环境行为或毒理学参数(统称为活性)的方法，即定量构效关系(QSAR)方法是其中一种重要的预测环境迁移转化行为参数技术。它是一种比实验测试技术更快、更廉价，不依赖于实验标准品的方法。QSAR反映和揭示有机污染物的分子结构与其环境行为和毒理效应之间的内在联系，具有弥补测试数据的缺失、降低测试费用、评估数据的不确定性的优点，并在有机化学品的生态风险评价领域得到了越来越广泛的重视和应用。因此，化合物与羟基反应速率常数QSAR模型的建立，相对于基于传统的实验方法来确定k_OH数值，更有利于快捷高效地评价化学品的环境持久性及进行生态风险评估。

QSAR方法，作为一种可信的技术工具，在欧盟《关于化学品注册、评估、许可和限制制度》里也有体现，制度明确规定QSAR方法可以为化学品的注册提供信息支持。2004年OECD提出的QSAR模型构建和使用准则规定，具有如下5个标准的QSAR模型可以应用于化学品的风险评价与管理：(1)具有明确定义的环境指标；(2)具有明确的算法；(3)定义了模型的应用域；(4)模型具有适当的拟合度、稳定性和预测能力；(5)最好能够进行机理解释。

目前，已有很多研究者应用QSAR方法成功地建立了有机化合物与·OH反应的k_OH的预测模型。如文献“Journal of Computational Chemistry.2011,32:2386-2396”基于DRAGON描述符和一些免费在线描述符应用多元线性回归(MLR)方法建立了涵盖460种均相有机化合物的k_OH预测模型,具有很好的预测能力和稳健性；文献“Atmospheric Environment.2009,43:1131-1135”运用偏最小二乘(PLS)方法，建立了298K温度下包含22个分子结构描述符的QSAR模型，数据集扩展到722种有机化合物，模型的相关系数达到了0.878。文献“Atmospheric Environment.2005,39:2189-2200”则采用333个描述符，建立了涵盖733种有机化合物的羟基反应速率常数QSAR预测模型，相关系数达到了0.90以上(r²=0.906)。虽然模型具有良好的预测性能，但是由于变量数目过多，导致使用起来过于麻烦。文献“Journal of Atmospheric Chemistry.2012,69:201-213”采用遗传算法与MLR相结合的方法最终筛选出3个量子化学描述符，成功建立了161种烷烃化合物的支持向量机(SVM)模型，其相关系数达到0.919。文献“Environmetrics.1994,5:197-208”通过PLS方法建立了58种卤代脂肪烃化合物的k_OH预测模型。但是上述模型只适用于单一结构类型的化合物(烷烃化合物，卤代脂肪烃化合物)，即模型的应用域小。文献“Journal ofChemical Information and Computer Sciences.2004,44:1794-1802”采用基因算法与MLR相结合的方法，选取4个描述符，成功建立了460种不同种类化合物的羟基反应速率常数QSAR模型(r²=0.828)。上述研究建立的羟基反应速率常数QSAR模型，虽然应用域较广并且预测性能良好，但是只能预测室温(即25⁰C)下的k_OH。反应速率常数具有较强的温度依附性，可以用阿伦尼乌斯方程即k=Ae^--E/RT来表示，有时也会用三参数的公式即k=ATⁿe^-E/RT来表示，对于大气对流层而言，温度分布介于200-300K，获得涵盖这一温度范围的速率常数对于评价处于对流层不同高度的有机化合物的归趋具有重要意义。因此本发明除了发展一个298K温度下的QSAR模型，拓展原有模型的应用域外，还建立一个应用域覆盖多种不同种类化合物的温度依附性羟基反应速率常数QSAR模型，即将温度T及1/T作为变量加入模型。同时，应OECD的要求，建立模型后，对模型进行应用域的表征及机理解释。

发明内容

本发明提供了一种简便、快捷、高效预测有机化学品与羟基反应速率常数的方法，该方法可以根据化合物分子结构预测其k_OH数值，进而可以对其环境持久性及大气归趋进行预测和评价，为化学品风险评价和管理提供必要的基础数据。

本发明的技术方案如下：

首先，收集得到化合物数据集涵盖温度T为206K～1364K的1543个k_OH数值，其中在温度为298K的数据872个；然后，对收集得到的1543个k_OH数值进行建立模型，包括温度依附性的k_OH QSAR模型和298K的k_OH QSAR模型；将上述收集得到的1543个和872个化合物数据集分别按照4:1的比例拆分为训练集和验证集，其中，温度依附性的k_OH QSAR模型训练集化合物个数为n=1234，验证集化合物个数n_ext=309；298K的k_OH QSAR模型训练集化合物个数为n=696，验证集化合物个数n_ext=176。

上述两种模型选取2种来源的描述符：

(1)量子化学描述符，选取了6个量子化学描述符；

(2)DRAGON描述符，基于优化后的构型得到了2000多种分子结构描述符；建模过程中，将温度T及1/T作为温度依附性的k_OH QSAR模型描述符；分别对训练集中的k_OH数据和上述所有描述符进行MLR回归分析；

两个模型的线性关系式如下：

模型1:温度依附性的k_OH QSAR模型

logk_OH=-8.613-0.02100×X%+14.38×E_HOMO-0.6430×Mor29u+0.5870×NdsCH+0.5870×GATS1e+0.5770×X3A-0.2450×SdsCH-167.0×(1/T)+1.103×BIC1+0.1170×RDF015m-1.044×SpMin8_Bh(p)+0.2390×nR=Cp-0.1980×NssssC-0.5080×F02[F-Br]

(1)

模型2:T=298K的k_OH QSAR模型

log(k_OH)=-6.511+15.85×E_HOMO-0.03800×AMW+0.1300×NdsCH+0.1630×Mor14i+0.3170×nR=Cp+0.7790×nP+0.3930×nRCHO-0.01900×X%-0.4550×SpMaxA_AEA(dm)+0.5890×(C-020)-0.05600×nCbH+0.1410×CATS2D_03_DL

(2)

其中，X%表示卤素原子在分子中所占的百分比，E_HOMO表示最高占据分子轨道能量，Mor29u表示未加权的3D分子结构描述符，NdsCH表示分子中具有=CH-结构的数目,GATS1e表示Sanderson电负性加权的lag1的Geary自相关指数，X3A表示3阶平均连接指数，SdsCH表示=CH-结构电性拓扑态的总和，1/T表示温度的倒数，BIC1表示1阶临近对称键信息内容指标，RDF015m表示0.15质量加权的径向分布函数，SpMin8_Bh(p)表示极化率加权的burden矩阵的最小特征值，nR=Cp表示末端sp2杂化的主碳数目，NssssC表示分子中含>C<结构的个数，F02[F-Br]表示在拓扑距离为2处F-Br出现的频率，AMW表示平均分子质量，Mor14i表示离子化势加权的3D分子结构描述符，nP表示分子中磷原子个数，nRCHO表示分子中脂肪族醛的个数，SpMaxA_AEA(dm)表示偶极矩加权的来自于扩增边缘临界处的标准化的主要特征值，C-020表示分子中=CX2结构信息，nCbH表示非取代苯环上sp2杂化碳的数目，CATS2D_03_DL表示位于lag03处的CATS2D亲脂性供体。

所得模型1训练集化合物个数为n=1234，每个描述符的变量膨胀因子(VIF)均小于10，自变量与因变量组成的矩阵M_YX以及自变量矩阵M_X的K相关指数K_XY和K_XX满足K_XX(0.374)<K_XY(0.405)，都表明模型不存在多重相关性；模型的拟合能力由R²和均方根误差(RMSE)表征，R²=0.873，RMSE=0.369，表明该模型具有良好的拟合能力；模型的稳健性由内部验证的交叉验证系数(Q² _LOO)和Bootstrapping方法所得Q² _BOOT评价，Q² _LOO=0.871，Q² _BOOT=0.797，R²和Q²之差远小于0.3，可认为该模型不存在过拟合现象，具有良好的稳健性；在模型的外部验证过程中，验证集数据个数n_ext=309，外部预测相关系数的R² _ext=0.838,Q² _ext=0.835，RMSE_ext=0.452，表面该模型具有良好的外部预测能力。

所得模型2训练集化合物个数为n=696，每个描述符的变量膨胀因子(VIF)均小于10，自变量与因变量组成的矩阵M_YX以及自变量矩阵M_X的K相关指数K_XY和K_XX满足K_XX(0.312)<K_XY(0.358)，都表明模型不存在多重相关性；模型的拟合能力由R²和均方根误差(RMSE)表征，R²=0.883，RMSE=0.419，表明该模型具有良好的拟合能力；模型的稳健性由内部验证的交叉验证系数(Q² _LOO)和Bootstrapping方法所得Q² _BOOT评价，Q² _LOO=0.879，Q² _BOOT=0.795，R²和Q²之差远小于0.3，可认为该模型不存在过拟合现象，具有良好的稳健性；在模型的外部验证过程中，验证集数据个数n_ext=176，外部预测相关系数的R² _ext==0.858,Q² _ext=0.851，RMSE_ext=0.489，表面该模型具有良好的外部预测能力；该模型较前人的研究成果比，涵盖化合物的数目更多，种类也从原有的烷烃，烯烃，炔烃，芳香类化合物，醛，酮，醚，酸，酯，卤代化合物，硝基化合物，含硫化合物和有机硅化合物拓展到有机磷酸酯，二甲基萘，有机硒和有机汞，拓展了模型的应用域；且模型采用了较少数量的描述符，在保证预测结果准确可靠的情况下，使得使用起来更加方便快捷。

结果表明模型1和2可以有效预测烷烃类化合物、烯烃类化合物、炔烃类化合物、芳香烃类化合物、醇类化合物、醛类化合物、酮类化合物、醚类化合物、酸类化合物、脂类化合物、卤代类化合物、含氮化合物、含硫化合物、含磷化合物、有机硅化合物、有机汞化合物、有机铅化合物、有机硒化合物等化合物的k_OH数值；以上两个模型的应用域表征均采用欧几里德距离方法。某一化合物特征向量(描述符参数)的欧几里德距离是其到描述符空间中心点的特征向量的距离；公式如下：

$X_k=\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_{\mathrm{jk}}---\left(3\right)$

其中X_k是中心点特征向量的第k个参数(描述符)，X_jk为化合物j的第k个描述符值，n为训练集化合物总数。则某一化合物i在M维空间(M为所采用分子结构描述符个数)的特征向量的欧几里德距离d_i可表示为：

$d_i=\sqrt{{\mathrm{\&Sigma;}}_{k=1}^M{(X_{\mathrm{ik}}-X_k)}^2}---\left(4\right)$

式中X_ik为化合物i的第k个描述符值。对于模型1和模型2，训练集化合物特征向量到中心点特征向量的欧几里德距离范围分别为0.1870-1.599和0.1530-1.593，因此特征向量欧几里德距离不大于1.599和1.593的化合物分别适用于模型1和模型2。

本发明的有益效果是：模型1可以用于有效预测不同温度下，多种类有机化合物与·OH反应的k_OH数值；模型2可以用于有效预测温度298K下有机化合物与·OH反应的k_OH数值。该方法不仅成本低廉、简便而快速，而且节省大量的人力、物力和财力；该发明涉及的k_OH预测方法的建立和验证严格依照OECD规定的QSAR模型发展和使用的导则，因此使用该发明专利的k_OH预测结果，可以为化学品监管工作提供重要的基础数据，并对生态风险评价具有重要的指导意义。

本发明提供的方法具有如下特点：

1、模型的应用域涵盖脂肪族，芳香族类，有机硅、铅、硒类等多种结构的有机化合物，可用于预测多种不同种类化合物在不同温度下的k_OH数值，且温度覆盖范围广泛，可为有机化学品环境持久性评估及生态风险评价提供基础数据。

2、建模过程中采用OECD推荐的透明算法-MLR算法，对于温度依附性模型1，筛选出14个描述符，模型2筛选出12个描述符，模型简洁，便于分析理解和应用。

3、按照OECD关于QSAR模型构建和使用的导则，建立的模型具有良好的拟合能力，稳健性和预测能力。

附图说明

图1为温度依附性k_OH预测模型的训练集和验证集log k_OH的实测值与预测值的拟合图。

图2为298K温度下k_OH预测模型的训练集和验证集log k_OH的实测值与预测值的拟合图。

图3为建立的温度依附性k_OH预测模型应用域表征图。

图4为建立的298K温度下k_OH预测模型应用域表征图。

图中：○表示训练集化合物，△表示验证集化合物。

具体实施方式

实施例1

给定一个化合物甲基丙烯酸丁酯，要预测其在287K，294K,298K，303K和313K温度下的logk_OH数值。首先需要根据甲基丙烯酸丁酯的结构信息，使用Gaussian09软件对其进行结构优化后，即可计算出E_HOMO值为-0.2752，通过Draogon6.0软件计算出X%，Mor29u，NdsCH，GATS1e，X3A，SdsCH，BIC1，RDF015m，SpMin8_Bh(p)，nR=Cp，NssssC和F02[F-Br]的值分别为0，-0.84，0，0.708，0.267，0，0.565，5.99，0，1，0和0。然后根据公式(4)计算得到特征向量的欧几里德距离为0.5316(<1.599)，在模型应用域范围内，故可以利用模型1进行预测。将T=287K，294K，298K，303K和313K，及上述描述符的值带入公式(1)，即可得出logk_OH在T=287K，294K，298K，303K和313K的预测值为-10.48，-10.47，-10.46，-10.45和-10.43。与其在T=287K，294K，298K，303K和313K下的logk_OH实验数据-10.12，-10.16，-10.18，-10.20和-10.25进行拟合，得出二者相关性系数R²=0.986，预测值与实验数据十分相符。

实施例2

给定一个实验数据较多的化合物2,3-二甲基丁烷，预测其在298K，843K，869K，872K，887K，893K，937K，938K，951K，977K，1009K，1038K，1039K，1043K，1082K，1107K，1144K，1156K，1161K，1183K，1206K，1231K，1247K，1249K和1292K温度下的logk_OH数值并与实验值进行对比。根据2,3-二甲基丁烷的结构信息，使用Gaussian09软件对其进行结构优化后，即可计算出E_HOMO值为-0.30974，通过Draogon6.0软件计算出X%，Mor29u，NdsCH，GATS1e，X3A，SdsCH，BIC1，RDF015m，SpMin8_Bh(p)，nR=Cp，NssssC和F02[F-Br]的值分别为0，-0.868，0，1.667，0.333，0，0.272，4.225，0，0，0，0。然后根据公式(4)计算得到特征向量的欧几里德距离为0.4967(<1.599)，在模型应用域范围内，故可以利用模型1进行预测。将T=298K，843K，869K，872K，887K，893K，937K，938K，951K，977K，1009K，1038K，1039K，1043K，1082K，1107K，1144K，1156K，1161K，1183K，1206K，1231K，1247K，1249K，1292K及上述描述符值带入公式(1),即可得出logk_OH在298K，843K，869K，872K，887K，893K，937K，938K，951K，977K，1009K，1038K，1039K，1043K，1082K，1107K，1144K，1156K，1161K，1183K，1206K，1231K，1247K，1249K，1292K温度下的logk_OH预测值为-11.11，-10.74，-10.74，-10.74，-10.73，-10.73，-10.72，-10.72，-10.72，-10.72，-10.71，-10.71，-10.71，-10.71，-10.7，-10.7，-10.69，-10.69，-10.69，-10.69，-10.68，-10.68，-10.68，-10.68，-10.67。其在T=298K，843K，869K，872K，887K，893K，937K，938K，951K，977K，1009K，1038K，1039K，1043K，1082K，1107K，1144K，1156K，1161K，1183K，1206K，1231K，1247K，1249K，1292K下的logk_OH实验数据-11.21，-10.67，-10.68，-10.66，-10.71，-10.67，-10.68，-10.61，-10.65，-10.64，-10.61，-10.62，-10.62，-10.56，-10.57，-10.53，-10.56，-10.51，-10.51，-10.5，-10.46，-10.47，-10.43，-10.4，-10.4进行拟合，得出二者相关性系数R²=0.8516，预测值与实验数据十分相符。

实施例3

给定一个化合物对二甲苯，预测其在T=298K下的logk_OH并与实验值进行对比。根据2,3-二甲基丁烷的结构信息，使用Gaussian09软件对其进行结构优化后，即可计算出E_HOMO值为-0.23683，通过Draogon6.0软件计算出AMW，NdsCH，Mor14i，nR=Cp，nP，nRCHO，X%，SpMaxA_AEA(dm)，C-020，nCbH，CATS2D_03_D分别为5.899，0，0.533，0，0，0，0，0.34，0，4，0。然后根据公式(4)计算得到特征向量的欧几里德距离为0.3876(<1.593)，在模型应用域范围内，故可以利用模型2进行预测。将上述描述符的值带入公式(2)，即可得出logk_OH的预测值为-10.78，。与其在T=298K下的logk_OH实验数据-10.86对比，差值为0.08。预测值与实验值十分相近，说明其预测性良好。

Claims (2)

1.一种通过定量构效关系模型预测大气中有机物与羟基反应速率常数的方法，其特征在于，

首先，收集得到化合物数据集涵盖温度T为206K～1364K的1543个k_OH数值，其中在温度为298K的数据872个；

然后，对收集得到的1543个k_OH数值进行建模，包括温度依附性的k_OH QSAR模型和298K的k_OH QSAR模型；将上述收集得到的1543个和872个化合物数据集分别按照4:1的比例拆分为训练集和验证集，其中，温度依附性的k_OH QSAR模型训练集化合物个数为n=1234，验证集化合物个数n_ext=309；298K的k_OH QSAR模型训练集化合物个数为n=696，验证集化合物个数n_ext=176；

上述两种模型选取2种来源的描述符：

(1)量子化学描述符，选取了6个量子化学描述符；

(2)DRAGON描述符，基于优化后的构型得到了2000多种分子结构描述符；建模过程中，将温度T及1/T作为温度依附性的k_OH QSAR模型描述符；分别对训练集中的k_OH数据和上述所有描述符进行MLR回归分析；

两个模型的线性关系式如下：

模型1:温度依附性的k_OH QSAR模型

logk_OH=-8.613-0.02100×X%+14.38×E_HOMO-0.6430×Mor29u+0.5870×NdsCH+0.5870×GATS1e+0.5770×X3A-0.2450×SdsCH-167.0×(1/T)+1.103×BIC1+0.1170×RDF015m-1.044×SpMin8_Bh(p)+0.2390×nR=Cp-0.1980×NssssC-0.5080×F02[F-Br]

(1)

模型2:T=298K的k_OH QSAR模型

log(k_OH)=-6.511+15.85×E_HOMO-0.03800×AMW+0.1300×NdsCH+0.1630×Mor14i+0.3170×nR=Cp+0.7790×nP+0.3930×nRCHO-0.01900×X%-0.4550×SpMaxA_AEA(dm)+0.5890×(C-020)-0.05600×nCbH+0.1410×CATS2D_03_DL

(2)

其中，X%表示卤素原子在分子中所占的百分比，E_HOMO表示最高占据分子轨道能量，Mor29u表示未加权的3D分子结构描述符，NdsCH表示分子中具有=CH-结构的数目,GATS1e表示Sanderson电负性加权的lag1的Geary自相关指数，X3A表示3阶平均连接指数，SdsCH表示=CH-结构电性拓扑态的总和，1/T表示温度的倒数，BIC1表示1阶临近对称键信息内容指标，RDF015m表示0.15质量加权的径向分布函数，SpMin8_Bh(p)表示极化率加权的burden矩阵的最小特征值，nR=Cp表示末端sp2杂化的主碳数目，NssssC表示分子中含>C<结构的个数，F02[F-Br]表示在拓扑距离为2处F-Br出现的频率AMW表示平均分子质量，Mor14i表示离子化势加权的3D分子结构描述符，nP表示分子中磷原子个数，nRCHO表示分子中脂肪族醛的个数，SpMaxA_AEA(dm)表示偶极矩加权的来自于扩增边缘临界处的标准化的主要特征值，C-020表示分子中=CX2结构信息，nCbH表示非取代苯环上sp2杂化碳的数目，CATS2D_03_DL表示位于lag03处的CATS2D亲脂性供体。

2.根据权利要求1所述的方法，其特征在于，所述的化合物包括烷烃类化合物、烯烃类化合物、炔烃类化合物、芳香烃类化合物、醇类化合物、醛类化合物、酮类化合物、醚类化合物、酸类化合物、脂类化合物、卤代类化合物、含氮化合物、含硫化合物、含磷化合物、有机硅化合物、有机汞化合物、有机铅化合物、有机硒化合物。

Images