预测海岸带持久性有机污染物消除速率的方法

摘要

本发明涉及预测海岸带持久性有机污染物消除速率的方法，包括以下步骤：采用偏最小二乘方法构建QSAR模型，根据QSAR模型得到有机污染物的消除速率。本发明的模型应用域明确，具有良好的拟合效果、稳健性和预测能力。在获得化合物分子结构的基础上，仅通过计算表征结构特征的描述符，应用所建QSAR模型，即能快速高效预测POPs在贻贝体内的消除速率，其成本低廉，简便快捷，节省测试所需的大量人力费用和时间，为持久性有机污染物的生态风险评价和管理提供重要的数据支持，具有重要意义。

说明书

技术领域

本发明涉及一种通过建立定量构效关系模型(QSAR)预测持久性有机污染物 在贻贝(Elliptio complanata)体内消除速率的方法，属于生态风险评价的定量结构 与活性关系技术领域。

背景技术

多环芳烃(PAHs)、多溴联苯醚(PBDEs)和多氯联苯(PCBs)是典型的持久性有 机污染物(POPs)，具有持久性、生物积累性和三致效应。贻贝作为一类世界性分 布的沿岸底栖双壳类海洋动物，已被广地用于海洋污染的生物监测中。贻贝具 备了作为指示生物的4个基本条件：地理分布广、是沿海和河口海区的优势种、 生物体能累积污染物、能对多种环境污染物产生响应。此外，对贻贝的基础生 物学特点研究已较透彻，因而将贻贝作为海洋污染监测的指示生物已普遍被人 们接受。贻贝分布广泛，可生物富集POPs，目前已在大连湾、胶州湾、鸭儿湖 及珠三角地区的贻贝体内发现了PAHs、PBDEs和PCBs。消除速率(k_d)是一个重 要的动力学参数，可用于估计污染物被贻贝的摄取速率和污染物在贻贝体内的 生物富集因子。

实验测定是目前获得化合物k_d值的主要途径。但是由于现有化学品数量已 超过14万种，据欧盟于2007年6月开始全面实施的化学品管理新法规《化学 品注册、评估、授权和限制法规》(Registration,Evaluation,Authorization and  Restriction of Chemicals,简称REACH法规)”估算，每一种化学物质的基本检测 费用约需8.5万欧元(不含长期环境影响的评估费用)，每一新物质全面检测费用 约需57万欧元，这意味着如果对每种化学品都开展实验测定，将会带来巨额的 财政开支。进行全面的实验测试，也不符合化学品风险管理中的动物实验伦理 的“3R原则”，即：减少(Reduce)试验动物的个数、改进(Refine)动物实验的方法、 替代(Replace)实验动物的模型。目前，“3R原则”逐渐被欧美各国的研究机构所 采纳，并在美国、英国、法国乃至整个欧洲的一些立法和政府规章中得以明确。 因此，如果全部通过试验测试获得所有相关信息，不但存在测试费用昂贵和动 物伦理方面的问题，而且有的化学品的标样难以获得，在时间上是滞后的，难 以满足有机化学品生态风险评价和监管的数据需求。因此，通过定量结构活性 关系(QSAR)方法发展一种能够快速高效获取有机污染物生物消除速率的模型具 有重要的应用意义。

QSAR是指关联有机污染物的分子结构与其理化性质、环境行为和毒理学参 数(统称为活性)的定量预测模型，能够反映和揭示有机污染物的分子结构与其环 境行为和毒理效应之间的内在联系，比试验测试更快、更廉价，该方法已经在 有机污染物的环境光化学行为、光致毒性参数、环境分配参数和生态毒理学效 应等方面得以应用。因此，化合物的生物消除速率QSAR模型的建立，相对于 基于传统的实验方法来确定k_d数值，更有利于快捷高效的评价化学品的生物富 集能力及进行生态风险评估。

经济合作和发展组织(OECD)围绕化学品的安全性问题，开展了QSAR技术 的应用研究。2004年，为规范QSAR技术的应用，OECD提出了QSAR模型构 建和应用导则。该导则要求建立的QSAR应满足如下标准：①针对明确定义的 环境指标；②具有清晰和明确的数学算法；③定义了模型的应用域；④模型 具有适当的拟合度、稳定性和预测能力；⑤最好能够进行机理解释。2007年， OECD发布了关于确认和验证QSAR模型的指导文件。

目前，已有许多研究者应用QSAR方法建立了有机污染物生物消除速率k_d的预测模型。如文献“QSAR&Combinatorial Science.2011,28:537–541”基于量 子化学描述符和偏最小二乘(PLS)方法，建立了PAHs在贻贝(Elliptio complanata) 体内k_d的QSAR模型。模型具有较好的稳健性，但仅适用于单一结构类型化合 物(多环芳烃化合物)，即模型应用域小。文献“Archives of Environment  Contamination and Toxicology.2004,47:74–83”基于logK_OW值构建了PCBs的k_d值模型，但模型过于简单，没有考虑作用机理和性能评价。文献“Chemosphere. 2007,69:362–370”建立了8个PBDEs和5个PCBs在贻贝体内k_d值与K_OW值的 QSAR模型，但训练集化合物的数目太少，也没有对模型的机理解释。文献 “Science in China,Series B Chemistry.2009,52:1281–1286”基于PLS和理论分子 描述符建立了k_d的QSAR模型，但没有对模型的应用域进行表征。因此有必要 建立一个涵盖多种类化合物，并且模型结构简单、预测规则透明、易于理解和 实际应用的QSAR模型，同时按照OECD导则对模型进行应用域表征和机理解 释。

发明内容

针对上述技术不足，本发明的目的是发展一种简洁、快捷、高效预测持久 性有机污染物在贻贝体内消除速率的方法。

本发明解决其技术问题所采用的技术方案是：预测海岸带持久性有机污染 物消除速率的方法，包括以下步骤：通过采用偏最小二乘方法构建的QSAR模 型得到有机污染物的消除速率。

所述QSAR模型为

${\log k}_d=1.34\times {10}^{-1}-1.07\times {10}^3\omega -4.86\times {10}^{-3}\mathrm{polar}-7.65\times {\overline{V_s}}^{-}-4.31\times {10}^{-2}\tau$

其中，k_d为有机污染物的消除速率，ω表示分子亲电性指数，polar表示平 均分子极化率；表示分子表面上负静电势的平均值；τ表示分子表面静电势的 平衡参数。

所述分子亲电性指数通过以下公式得到

$\omega =\frac{{\mu }^2}{2\eta },$化学势$\mu =\frac{E_{\mathrm{LUMO}}+E_{\mathrm{HOMO}}}{2},$化学硬度$\eta =\frac{E_{\mathrm{LUMO}}-E_{\mathrm{HOMO}}}{2},$E_HOMO是有机污 染物的最高占据轨道能；E_LUMO是有机污染物的最低未占据轨道能。

所述分子表面上负静电势的平均值通过以下公式得到

${\overline{V_S}}^{-}=\frac{1}{\beta }\underset{j=1}{\overset{\beta }{\Sigma }}{V_S}^{-}\left(r_j\right)$

其中，β为分子表面负的静电势计算的点数，V_S^-(r_j)为第j个点的分子表面 负的静电势。

所述分子表面静电势的平衡参数通过以下公式得到

$\tau =\frac{{\sigma }_{+}^2{\sigma }_{-}^2}{{\left({\sigma }_{\mathrm{tot}}^2\right)}^2},$${\sigma }_{\mathrm{tot}}^2={\sigma }_{+}^2+{\sigma }_{-}^2=\frac{1}{\alpha }\underset{i=1}{\overset{\alpha }{\Sigma }}{[V^{+}\left(r_i\right)-{\overline{V_s}}^{+}]}^2+\frac{1}{\beta }\underset{j=1}{\overset{\beta }{\Sigma }}{[V^{-}\left(r_j\right)-{\overline{V_s}}^{-}]}^2$

其中，和分别为分子Bader面静电势为正值和负值的区域静电势分布的 方差；α和β分别为分子表面正的静电势和负的静电势计算的点数；V⁺(r_i)和V^-(r_j) 分别表示分子表面上正的静电势和负的静电势；和分别为分子表面正静电 势和负静电势的平均值。

所述有机污染物包括多环芳烃、多溴联苯醚、多氯联苯。

本发明具有以下有益效果及优点：

1.本发明采用偏最小二乘回归算法，基于机理分析选取描述符构建了预测 模型，建立的QSAR模型具有良好的拟合优度、稳健性和预测能力。模型简洁， 透明度强，便于理解和实际应用。

2.本发明的QSAR模型，明确了应用域范围，模型涵盖多种结构的有机化 合物，可以为持久性有机污染物的风险评价和管理工作提供基础数据。

3.采用本发明方法可以预测不同类型持久性有机污染物在贻贝体内的消除 速率，该方法成本低廉、简便而快速，能够大量节省实验测试所需的人力、费 用和时间。

4.本发明涉及的生物消除速率预测方法的建立和验证严格依照OECD规定 的QSAR模型发展和使用导则，因此使用本发明专利的生物降解性预测结果， 可以为有机化学品生态风险评价和管理提供重要的数据支持，具有重要的理论 和现实意义。

附图说明

图1为本发明的方法流程图；

图2为logk_d的实测值与预测值的拟合图；

图3为模型应用域表征图。

具体实施方式

下面结合实施例对本发明做进一步的详细说明。

该方法可以直接根据化合物分子结构预测其消除速率，进而对目标化合物 的环境持久性进行预测和评价，为化学品风险评价和管理提供必要的基础数据。 如图1所示，包括如下步骤：

(1)数据搜集和划分：搜集文献中实验测定的28个多环芳烃(PAHs)、8个 多溴联苯醚(PBDEs)和28个多氯联苯(PCBs)用于预测在贻贝(Elliptio complanata) 体内的消除速率k_d值，将上述数据随机按照4:1划分为训练集和验证集；

(2)QSAR模型构建：根据作用机理分析，计算13种量子化学描述符，采 用偏最小二乘方法构建logk_d和训练集化合物分子描述符的定量构效关系QSAR 模型，表达式为：

${\log k}_d=1.34\times {10}^{-1}-1.07\times {10}^3\omega -4.86\times {10}^{-3}\mathrm{polar}-7.65\times {\overline{V_s}}^{-}-4.31\times {10}^{-2}\tau$

其中ω表示分子亲电性指数；polar表示平均分子极化率；表示分子表面 上负静电势的平均值；τ表示分子表面静电势的平衡参数。

(3)QSAR模型的验证与应用域表征：QSAR模型的验证结果用外部预测 相关系数的平方Q²_EXT和均方根误差RMSE表示，QSAR模型的化合物应用域采 用Williams图的方法进行表征。

本发明共搜集了28个多环芳烃(PAHs)、8个多溴联苯醚(PBDEs)和28个多 氯联苯(PCBs)在贻贝(Elliptio complanata)体内的消除速率kd值。

本发明采用的技术方案包括如下步骤：

(1)为保证用于建模的数据准确性，对文献收集的实验测定值进行评估和 分析。尽量选取采用相同实验方法的实验数据。化合物包括多环芳烃、多溴联 苯醚和多氯联苯。以4:1的比例随机划分训练集和验证集。

(2)机理剖析及分子结构描述符的选取。据推测，化合物在贻贝体内的消 除速率可能与其平衡分配和相互作用有关。这个过程包括：化合物在生物相中 的分配能力、偶极-偶极相互作用、氢键和静电相互作用。

基于此，选取13个理论量子化学描述符用于QSAR建模：电子密度体积(V)、 平均分子极化率(polar)用于表征空穴相互作用；分子最高占据轨道能(E_HOMO)、 分子最低未占据轨道能(E_LUMO)、分子中氢原子的最正净电荷(qH⁺)、电负性指数 (ω)、分子表面上的最正和最负静电势(V_s,max,V_s,min)、分子表面上正静电势和负静 电势的平均值分子表面静电势的分散度(П)、静电势的平衡参数(τ)被 选用来表征氢键和静电相互作用。

(3)分子描述符的计算。V定义为电子密度为0.001电子/立方波尔的空间 内的电子密度体积。量子化学描述符均采用Gaussian09程序软件计算得到。具 体过程是：首先在Chemoffice中生成初始的分子结构，并利用其中的PM3方法 进行初步优化，将得到的图形文件转化成Gaussian输入文件。然后运用Gaussian 09程序包中的密度泛函理论(DFT)中的B3LYP方法，对化合物分子在6-31G(d,p) 基组水平上进行结构优化，获得其稳定的分子结构。接着对优化好的几何构型 进行频率分析，以确保体系无虚频。在上述几何优化和DFT的计算中都使用自 洽场(SCRF)和积分方程形式极化连续介质模型(IEFPCM)来考虑溶剂(水)效应。 SCRF模型中，将溶质分子放置于真空空腔内，而溶剂被视为连续的无结构但具 有一定介电常数ε(水的ε=78.4)的介质。ω由下面的公式计算得到：

$\omega =\frac{{\mu }^2}{2\eta }---\left(1\right)$

$\mu =\frac{E_{\mathrm{LUMO}}+E_{\mathrm{HOMO}}}{2}---\left(2\right)$

$\eta =\frac{E_{\mathrm{LUMO}}-E_{\mathrm{HOMO}}}{2}---\left(3\right)$

这里，μ是化学势；η是化学硬度；E_HOMO是化合物的最高占据轨道能，E_LUMO是化合物的最低未占据轨道能。

分子表面静电势的计算在Gaussian09优化好的结构上进行。静电势计算采用 格点法，立方格精度设定为Cube=fine。这样对每个分子而言，可获得近100³个 点的静电势。然后，对这些点的静电势进行统计计算便获得所需参数。分别计 算了V_s,max、V_s,min、П和τ。具体参数定义如下：

分子表面正负静电势的平均值

${\overline{V_S}}^{+}=\frac{1}{\alpha }\underset{i=1}{\overset{\alpha }{\Sigma }}{V_S}^{+}\left(r_i\right)---\left(4\right)$

${\overline{V_S}}^{-}=\frac{1}{\beta }\underset{j=1}{\overset{\beta }{\Sigma }}{V_S}^{-}\left(r_j\right)---\left(5\right)$

平均分散度(П)

$\Pi =\frac{1}{\alpha +\beta }\underset{i=1}{\overset{\alpha +\beta }{\Sigma }}|V\left(r_i\right)-\overline{V_s}|---\left(6\right)$

$\overline{V_s}=\frac{\alpha {\overline{V_s}}^{+}+\beta {\overline{V_s}}^{-}}{\alpha +\beta }---\left(7\right)$

平衡常数(τ)

$\tau =\frac{{\sigma }_{+}^2{\sigma }_{-}^2}{{\left({\sigma }_{\mathrm{tot}}^2\right)}^2}---\left(8\right)$

${\sigma }_{\mathrm{tot}}^2={\sigma }_{+}^2+{\sigma }_{-}^2=\frac{1}{\alpha }\underset{i=1}{\overset{\alpha }{\Sigma }}{[V^{+}\left(r_i\right)-{\overline{V_s}}^{+}]}^2+\frac{1}{\beta }\underset{j=1}{\overset{\beta }{\Sigma }}{[V^{-}\left(r_j\right)-{\overline{V_s}}^{-}]}^2---\left(9\right)$

其中，s表示分子表面；α和β分别为分子表面正的静电势和负的静电势计算 的点数；V⁺(r_i)和V^-(r_j)分别表示分子表面上正的静电势和负的静电势；和分 别为分子表面正静电势和负静电势的平均值；表示分子表面的平均静电势；和分别为分子Bader面静电势为“正值”和“负值”的区域静电势分布的方差。 V_S⁺(r_i)、V_S^-(r_j)分别为第i个点的分子表面正的静电势、第j个点的分子表面负的静 电势，V(r_i)为第i个点的分子表面静电势。

(4)QSAR模型的建立。运用偏最小二乘(PLS)和步骤(3)中计算得到的分 子结构描述符，获得如下模型：

${\log k}_d=1.34\times {10}^{-1}-1.07\times {10}^3\omega -4.86\times {10}^{-3}\mathrm{polar}-7.65\times {\overline{V_s}}^{-}-4.31\times {10}^{-2}\tau$

其中ω表示分子亲电性指数；polar表示平均分子极化率；表示分子表面 上负静电势的平均值；τ表示分子表面静电势的平衡参数。

在模型中，logk_d表示为4个描述符变量的函数，训练集数据n=51。相关 系数平方(R²)为0.952，标准偏差(SE)=0.119，显著性水平(p)<0.0001，说明模 型具有良好的拟合优度。累积交叉验证系数(Q²_CUM)=0.947，表明模型具有很好 的稳健性(图2)。

(5)建立的QSAR模型的验证和应用域。

QSAR模型的预测能力需要通过外部验证来检验。验证集n=13个数据。外 部验证的结果可以通过外部预测相关系数的平方(Q²_EXT)以及外部验证结果的均 方根误差RMSE来表示。这两个参数定义为：

${Q_{\mathrm{EXT}}}^2=1-\frac{\underset{i=1}{\overset{n}{\Sigma }}{(y_i-{\hat{y}}_i)}^2}{\underset{i=1}{\overset{n}{\Sigma }}{(y_i-{\bar{y}}_{\mathrm{EXT}})}^2}---\left(10\right)$

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Sigma }}{(y_i-{\hat{y}}_i)}^2}{n}}---\left(11\right)$

其中，n代表化合物的个数，y_i和分别表示第i个化合物的实测值和预测 值；为化合物活性实测值的平均值；表示外部验证集预测值的平均值。

本发明采用基于leverage的Williams图评价了建立的QSAR模型的应用域， 如图3所示。Williams图是综合应用分子结构描述符的影响(leverage)值(以h_i表 示，i代表不同的化合物)和经过标准化的交叉验证残差来表征QSAR模型的描 述符域，这种方法已经被用于QSAR模型应用域的评价。

训练集化合物的H值(帽子矩阵)可以通过下面的公式求得：

H=X(X^TX)^–1X^T   (12)

式中，X为n×k的矩阵，n为训练集化合物的个数，k为预测变量(分子结构描 述符)的个数，X矩阵构成了训练集化合物的描述符空间。

每个化合物的leverage(h_i)值（即杠杆值）是H所对应的对角线的值，可由公 式(13)计算得到：

h_i=x_i^T(X^TX)^–1x_i   (13)

式中，x_i为第i个化合物分子结构描述符的行向量。

警戒值(h^*)定义为：

h^*=3(k+1)/n   (14)

其中，k为描述符的个数。

对于本发明的QSAR模型来说，在应用域内验证集化合物k_d预测的结果为： n=13，Q²_EXT=0.892，RMSE=0.160。表面该QSAR模型具有良好的预测能力(图 2)。

实施例1

2,6-二甲基萘：采用Williams图法计算得到其杠杆值为0.2349<h^*(预警值) =0.294，标准残差(SE)=-0.0288>-3，说明此化合物在QSAR模型应用域内。 基于作用机理，采用Gaussian09软件，按照发明中叙述方法计算得到模型中的 四个描述符值。

2,6-二甲基萘在贻贝体内的污染消除速率logk_d的实测值为：-0.577。基于 QSAR模型预测步骤如下：

logk_d=0.134–1070×(0.000674)–0.00486×(174.4643)–76.5×(-0.01197)–0.0431 ×(0.018006)=-0.527。预测值与实测值十分相符。

实施例2

二苯骈呋喃：采用Williams图法计算得到其杠杆值为0.1335<h^*(预警值)= 0.294，标准残差(SE)=0.2374<3，说明此化合物在QSAR模型应用域内。基于 作用机理，采用Gaussian09软件，按照发明中叙述方法计算得到模型中的四个 描述符值。

二苯骈呋喃在贻贝体内的污染消除速率logk_d的实测值为：-0.635。基于 QSAR模型预测步骤如下：

logk_d=0.134–1070×(0.000815)–0.00486×(169.273)–76.5×(-0.012547)–0.0431× (0.208105)=-0.527。预测值与实测值十分相符。

实施例3

1-甲基菲：采用Williams图法计算得到其杠杆值为0.0414<h^*(预警值)= 0.294，标准残差(SE)=0.3421<3，说明此化合物在QSAR模型应用域内。基于 作用机理，采用Gaussian09软件，按照发明中叙述方法计算得到模型中的四个 描述符值。

1-甲基菲在贻贝体内的污染消除速率logk_d的实测值为：-0.858。基于QSAR 模型预测步骤如下：

logk_d=0.134–1070×(0.000699)–0.00486×(224.2747)–76.5×(-0.0117)–0.0431× (0.120829)=-0.860。预测值与实测值十分相符。

实施例4

BDE-47：采用Williams图法计算得到其杠杆值为0.0320<h^*(预警值)= 0.294，标准残差(SE)=-0.3495>-3，说明此化合物在QSAR模型应用域内。基 于作用机理，采用Gaussian09软件，按照本发明中叙述方法计算得到模型中的 四个描述符值。

BDE-47在贻贝体内的污染消除速率logk_d的实测值为：-1.721。基于QSAR 模型预测步骤如下：

logk_d=0.134–1070×(0.000954)–0.00486×(264.171)–76.5×(-0.00788)–0.0431 ×(0.187566)=-1.648。预测值与实测值十分相符。

实施例5

PCB-97：采用Williams图法计算得到其杠杆值为0.0577<h^*(预警值)= 0.294，标准残差(SE)=-0.1619>-3，说明此化合物在QSAR模型应用域内。基 于作用机理，采用Gaussian09软件，按照本发明中叙述方法计算得到模型中的 四个描述符值。

PCB-97在贻贝体内的污染消除速率logk_d的实测值为：-1.721。基于QSAR 模型预测步骤如下：

logk_d=0.134–1070×(0.001234)–0.00486×(216.0167)–76.5×(-0.00904)–0.0431 ×(0.190156)=-1.626。预测值与实测值十分相符。