预测水相中有机化合物与羟基自由基反应速率常数的方法

技术领域

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

背景技术

羟基自由基(·OH)是水环境中一类常见的活性物种，能够通过水中溶解性物 质(如溶解性有机质)的光化学转化过程或人为的高级氧化过程(如Fenton体系) 产生。·OH的标准氧化还原电势为1.9V，是一种重要的强氧化剂，能够与有机 污染物发生氧化降解反应。该反应的二级反应速率常数(k_OH)是表征有机污染物 与·OH反应强度与能力的参数，也是评估有机污染物在环境中的持久性和进行 生态风险评价的一个重要指标。

化合物的水相k_OH可通过实验方法获得。但是实验的方法存在耗时、昂贵、 依赖于设备的缺点。截止到2014年11月，美国化学文摘社(Chemical Abstracts  Service,CAS)登记的化学品已超过9千万种(http://www.cas.org/)，其中投入商业 化使用的有机化学品超过14万种。如此数量庞大的有机化学品，无法通过实验 方法对其k_OH值进行一一测定，从而限制了对其环境持久性的评估。目前仅有 1000余种化合物具有可获取的水相k_OH值，因此，迫切需要发展非实验技术以 便高效快捷地获取k_OH数值，以满足有机化学品生态风险评价和管理的需求。

定量结构-活性关系(QSAR)是基于分子结构性质来预测化合物的理化性质、 环境行为或毒理学参数(统称为活性)的方法。它具有弥补测试数据的缺失、降低 测试费用的优点，在有机化学品的生态风险评价领域得到了越来越广泛的重视 和应用。通过构建QSAR模型，可以实现高效、快捷地获取有机化学品的k_OH数据，为化学品的环境持久性和生态风险评估提供支撑。2004年OECD提出了 QSAR模型构建和使用的准则，指出符合以下5个要求的QSAR模型可以应用 于化学品的风险评价与管理：(1)具有明确定义的环境指标；(2)具有明确的算 法；(3)定义了模型的应用域；(4)模型具有适当的拟合度、稳定性和预测能力； (5)最好能够进行机理解释。

目前，国内外用于预测k_OH的QSAR方法多集中在气相方面，水相k_OH的研 究较少。k_OH的气相和水相模型虽有一些相似和关联之处，但也有所区别，如水 相中氢键的作用、水分子的极性以及溶剂化效应等因素会显著影响有机物与·OH 的反应机制、反应自由能变等，从而改变其反应速率，同一化合物在水相和气 相中的k_OH差别非常显著。因此，对气相k_OH具有较好预测效果的QSAR模型并 不适用于水相k_OH的预测，非常有必要建立针对水相k_OH的预测模型。关于水相 k_OH的模型构建，前人开展了少许研究。目前已构建的水相k_OH预测模型虽然有 其自身的特点，也存在一些不足之处。这些不足主要体现在以下几方面：第一， 大部分模型涵盖化合物种类和数量较少，应用域较窄，如文献“QSAR and  Combinatorial Science.2009,28(11-12):1309-1316”中基于多元线性回归(MLR)方 法建立的包括酚类、烷烃和醇类的水相QSAR预测模型，涵盖化合物仅55种； 文献“Atmospheric Environment.2008,42(33):7611-7622”中基于SAR方法建立的 用于烷烃、醇、有机酸等脂肪族小分子的水相预测模型，仅包含72种化合物。 在大部分模型中，很多化合物未被涵盖，如含氮、磷、硫等原子的化合物，这 些化合物多为生产和生活中重要的化学品(如农药、抗生素等)，因此非常有必要 在模型应用域中涵盖这些类别的化合物。第二，有些模型采用实验测定的参数 作为描述符，限制了模型的推广应用，如文献“Physical Chemistry Chemical  Physics.2004,6(16):4118-4126”,“International Journal of Chemical Kinetics.2008, 40(4),174-188.”中利用实验测定的键离解能(BDE)来预测化合物的水相k_OH，其 推广应用很大程度上受到BDE实验值缺乏的限制。第三，模型不够透明，如文 献“Atmospheric Environment.2003,37(2):269-276”中基于分子结构碎片和人工 神经网络(ANN)建立的水相k_OH预测模型，但由于神经网络模型属于“黑箱模型”， 其无法将k_OH与模型参数的关系用明确的数学表达式呈现出来，也难以对模型进 行机理解释，因此不宜被接受和推广使用。第四，一些模型，如文献“Atmospheric  Environment.2005,39(40):7667-7688”构建的模型虽然具有较好的拟合性能，但 模型验证和表征并不全面，如欠缺内部和外部验证、模型应用域表征等，不符 合OECD对QSAR模型构建和使用准则的要求。此外，还有的模型根据化合物 结构和类别进行分类预测，如文献“Environmental Science & Technology.2009, 43(16):6220-6227”中应用基团贡献法基于反应官能团的特征对水相k_OH进行分 类预测，文献“环境化学.1999,18(3):232-237”“哈尔滨工业大学学报2002,34(4): 521-528”等中根据化合物结构和官能团将其划分为烷烃、醇、酚、有机酸等类分 别进行预测。分类模型中一般部分化合物的预测模型较优，但也存在预测效果 不理想(如拟合R值较低)的类别，且模型形式复杂，不便于应用。基于上述 研究现状，迫切需要发展应用域涵盖化合物数量较多且结构种类丰富、描述符 易获得、具有明确算法且透明度高、易于机理解释和便于应用推广的QSAR综 合模型。另外，根据OECD对QSAR模型构建和使用导则的要求，对模型进行 全面的模型性能评估、应用域的表征及机理解释，为有机化学品环境持久性评 估及生态风险评价提供基础数据。

发明内容

本发明提供了一种简便、快捷、高效预测有机化学品水相羟基反应速率常 数的方法，该方法可以根据化合物分子结构预测其k_OH数值，进而可以评估其环 境持久性，为化学品风险评价和管理提供必要的基础数据。在建模过程中参照 OECD对QSAR模型构建和使用导则，不仅进行了内、外部验证考察模型的预 测能力和稳健性，而且对模型应用域进行了表征。

通过查阅大量文献，搜集了水相k_OH实验数据，建立了一个涵盖526种有机 化合物的数据集，化合物种类不仅包括常见的烃类、芳香类、醇类、酸类化合 物，而且涵盖了以往研究中很少研究的含氮、磷、硫等杂原子的化合物。

本发明的技术方案如下：

预测水相中有机化合物与羟基自由基反应速率常数的方法，步骤如下：

首先收集526种有机化合物的水相k_OH数值，将其按照4:1随机拆分为训练 集和验证集；训练集中的421种化合物用于构建模型，验证集中的105种化合 物用于评估模型的外部预测能力；对上述526种有机化合物的分子结构进行优 化，得到相应有机化合物的稳定构型并选取9个量子化学描述符；通过上述稳 定构型得到2418个Dragon描述符；采用MLR回归分析方法筛选分子描述符和 构建QSAR模型；

筛选出的最优QSAR模型如下：

logk_OH＝11.566+6.233×E_HOMO-0.074×HATS2s-0.183×Mor23u+0.238 ×GATS1e-0.099×N-075+0.107×nR＝Cp-0.230×nRCONH2-0.070×C-001 +0.080×MLOGP+0.113×nS-0.265×nBR+0.651×q_H⁺+0.119× Eig03_EA(dm)；

其中，E_HOMO表示最高占据分子轨道能量，HATS2s是与内蕴状态相关的 GETAWAY描述符，Mor23u表示未加权的3D分子结构描述符，GATS1e表示 Sanderson电负性加权的lag 1的Geary自相关指数，N-075表示苯环上的N原子 或与O、N、S、卤素等电负性原子形成离域键的N原子碎片数，nR＝Cp表示末 端sp2杂化的主碳数目，nRCONH2表示分子中含RCONH₂结构的数目，C-001 表示分子中-CH₃/CH₄结构信息，MLOGP表示Moriguchi辛醇-水分配系数，nS 表示分子中含S原子的数目，nBR表示分子中含Br原子的数目，q_H⁺表示H原 子最正净电荷，Eig03_EA(dm)表示偶极矩加权的本征值的边界邻接指数。

所述有机化合物为烷烃类化合物、烯烃类化合物、炔烃类化合物、芳香烃 类化合物、醇类化合物、醛类化合物、酮类化合物、醚类化合物、酸类化合物、 脂类化合物、卤代类化合物、含氮化合物、含硫化合物或含磷化合物。

所得模型中每个描述符的方差膨胀因子(VIF)均小于10，且自变量与因变量 组成的矩阵M_YX以及自变量矩阵M_X的K相关指数满足K_XX(0.209)<K_XY(0.249)， 表明模型不存在多重相关性。模型的拟合能力由R²_adj和均方根误差(RMSE)表征， R²_adj＝0.805，RMSE＝0.165，表明该模型具有良好的拟合能力；模型的稳健性由内 部验证的交叉验证系数(Q²_LOO)和Bootstrapping方法所得Q²_BOOT评价， Q²_LOO＝0.797，Q²_BOOT＝0.791，R²和Q²之差远小于0.3，可认为该模型不存在过拟 合现象，具有良好的稳健性；在模型的外部验证过程中，外部预测相关系数的 R²_ext＝0.802,Q²_ext＝0.801,RMSE_ext＝0.232，表明该模型具有良好的外部预测能力。 采用Williams图表征模型的应用域。结果表明所建的模型能够有效地用于烷烃 类化合物、烯烃类化合物、炔烃类化合物、芳香烃类化合物、醇类化合物、醛 类化合物、酮类化合物、醚类化合物、酸类化合物、脂类化合物、卤代类化合 物、含氮化合物、含硫化合物、含磷化合物等的水相k_OH预测。

本发明的有益效果是：所建模型可以用于预测多种类有机化合物的水相 k_OH。该方法简便快捷、成本低廉。水相k_OH预测方法符合OECD规定的QSAR 模型发展和使用导则，因此，使用该发明专利的k_OH预测结果，可以为化学品监 管提供数据支持，对化学品的生态风险性评价具有重要意义。

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

1.模型数据集涵盖烃类、醇类、酮类、酚类、酸类、芳香类等多种结构的有机 化合物，尤其包含了其他模型中很少研究的含氮、磷、硫元素的化合物，应 用域广，是目前涵盖化合物数目最大、种类最丰富的k_OH预测模型；

2.建模过程中采用OECD对QSAR模型构建和使用导则推荐的透明算法 ——MLR算法，所建模型包含2个量子化学描述符和11个Dragon描述符， 模型算法透明，机理易于解释，便于应用推广；

3.依照OECD关于QSAR模型的构建和使用导则构建和评估模型，所建模型具 有良好的拟合能力、稳健性和预测能力，可以用于化学品的风险评价与管理。

附图说明

图1为训练集log k_OH的实测值与预测值的拟合图，训练集化合物为421种。

图2为验证集log k_OH的实测值与预测值的拟合图，验证集化合物为105种。

图3为模型应用域的Williams图。

具体实施方式

以下结合附图和技术方案进一步说明本发明的具体实施方式。

本发明的实施例，对构建的化合物分子结构，先进行能量最小化，之后再 进行能量优化，然后基于优化结构，提取量子化学描述符并计算Dragon描述符。

实施例1

正庚醇，实验测定的水相logk_OH值为9.87，计算的量子化学和Dragon描述符为： [E_HOMO]＝-0.394,[HATS2s]＝0.818,[Mor23u]＝-1.293,[GATS1e]＝1.049,[N-075]＝0, [nR＝Cp]＝0,[nRCONH2]＝0,[C-001]＝1,[MLOGP]＝1.940,[nS]＝0,[nBR]＝0, [q_H⁺]＝0.326,[Eig03_EA(dm)]＝0。

由MATLAB计算的h＝0.0171<h^*＝0.0998,所以该化合物在应用域内，由模型计算 如下：

logk_OH＝6.233[E_HOMO]-0.074[HATS2s]-0.183[Mor23u]+0.238[GATS1e]-0.099[N-07 5]+0.107[nR＝Cp]-0.230[nRCONH2]-0.070[C-001]+0.080[MLOGP]+0.113[nS]-0.26 5[nBR]+0.651[q_H⁺]+0.119[Eig03_EA(dm)]+11.566

＝6.233×(-0.394)-0.074×0.818-0.183×(-1.293)+0.238×1.049-0.07+0.080×1.940+0.65 1×0.326+11.566

＝9.83

实施例2

苯甲腈，芳香族含氮化合物，实验测定的水相logk_OH值为9.64，计算的量子化 学和Dragon描述符为：[E_HOMO]＝-0.373,[HATS2s]＝0.585,[Mor23u]＝-0.500, [GATS1e]＝0.476,[N-075]＝0,[nR＝Cp]＝0,[nRCONH2]＝0,[C-001]＝0, [MLOGP]＝1.769,[nS]＝0,[nBR]＝0,[q_H⁺]＝0.167,[Eig03_EA(dm)]＝0。

由MATLAB计算的h＝0.0171<h^*＝0.0998,所以该化合物在应用域内，由模型计算 如下：

logk_OH＝6.233[E_HOMO]-0.074[HATS2s]-0.183[Mor23u]+0.238[GATS1e]-0.099[N-07 5]+0.107[nR＝Cp]-0.230[nRCONH2]-0.070[C-001]+0.080[MLOGP]+0.113[nS]-0.26 5[nBR]+0.651[q_H⁺]+0.119[Eig03_EA(dm)]+11.566

＝6.233×(-0.373)-0.074×0.585-0.183×(-0.500)+0.238×0.476+0.080×1.769+0.651×0.1 67+11.566

＝9.65

实施例3

二异丙基亚砜，含S化合物，实验测定的水相logk_OH值为9.83，计算的量子化 学和Dragon描述符为：[E_HOMO]＝-0.314,[HATS2s]＝0.864,[Mor23u]＝-0.296, [GATS1e]＝0.400,[N-075]＝0,[nR＝Cp]＝0,[nRCONH2]＝0,[C-001]＝4, [MLOGP]＝1.587,[nS]＝1,[nBR]＝0,[q_H⁺]＝0.171,[Eig03_EA(dm)]＝0。

由MATLAB计算的h＝0.0277<h^*＝0.0998,所以该化合物在应用域内，由模型计算 如下：

logk_OH＝6.233[E_HOMO]-0.074[HATS2s]-0.183[Mor23u]+0.238[GATS1e]-0.099[N-07 5]+0.107[nR＝Cp]-0.230[nRCONH2]-0.070[C-001]+0.080[MLOGP]+0.113[nS]-0.26 5[nBR]+0.651[q_H⁺]+0.119[Eig03_EA(dm)]+11.566

＝6.233×(-0.314)-0.074×0.864-0.183×(-0.296)+0.238×0.400+0.070×4+0.080×1.587+ 0.113+0.651×0.171+11.566

＝9.77

实施例4

葡萄糖-1-磷酸，含P化合物，实验测定的水相logk_OH值为9.15，计算的量子化 学和Dragon描述符为：[E_HOMO]＝-0.391,[HATS2s]＝3.643,[Mor23u]＝0.039, [GATS1e]＝1.310,[N-075]＝0,[nR＝Cp]＝0,[nRCONH2]＝0,[C-001]＝0, [MLOGP]＝-2.438,[nS]＝0,[nBR]＝0,[q_H⁺]＝0.389,[Eig03_EA(dm)]＝0。

由MATLAB计算的h＝0.0192<h^*＝0.0998,所以该化合物在应用域内，由模型计算 如下：

logk_OH＝6.233[E_HOMO]-0.074[HATS2s]-0.183[Mor23u]+0.238[GATS1e]-0.099[N-07 5]+0.107[nR＝Cp]-0.230[nRCONH2]-0.070[C-001]+0.080[MLOGP]+0.113[nS]-0.26 5[nBR]+0.651[q_H⁺]+0.119[Eig03_EA(dm)]+11.566

＝6.233×(-0.391)-0.074×3.643-0.183×0.039+0.238×1.310+0.080×(-2.438)+0.651×0.3 89+11.566

＝9.22

实施例5

5-溴吲哚，含溴化合物，实验测定的水相logk_OH值为10.20，计算的量子化学和 Dragon描述符为：[E_HOMO]＝-0.329,[HATS2s]＝0.644,[Mor23u]＝-0.574, [GATS1e]＝0.851,[N-075]＝0,[nR＝Cp]＝0,[nRCONH2]＝0,[C-001]＝0, [MLOGP]＝2.573,[nS]＝0,[nBR]＝1,[q_H⁺]＝0.294,[Eig03_EA(dm)]＝0。

由MATLAB计算的h＝0.1038>h^*＝0.0998,所以该化合物不在应用域内，由模型计 算如下：

logk_OH＝6.233[E_HOMO]-0.074[HATS2s]-0.183[Mor23u]+0.238[GATS1e]-0.099[N-07 5]+0.107[nR＝Cp]-0.230[nRCONH2]-0.070[C-001]+0.080[MLOGP]+0.113[nS]-0.26 5[nBR]+0.651[q_H⁺]+0.119[Eig03_EA(dm)]+11.566

＝6.233×(-0.329)-0.074×0.644-0.183×(-0.574)+0.238×0.851+0.080×2.573-0.265+0.6 51×0.294+11.566

＝9.91

模型预测值9.91与实验值10.20较为接近，说明模型能够较好预测应用域外的 化合物，具有一定的外推能力。