CN104915563A - 基于金属定量构效关系的淡水慢性基准预测方法 - Google Patents

Abstract

本发明涉及水质污染评价领域，具体为基于金属定量构效关系的淡水慢性基准预测方法，根据重金属离子的结构特征与水生生物慢性毒性效应的定量关系预测未知金属的毒性终点，结合不同物种的敏感度分布分析推导保护不同比例的水生生物的危险浓度；是综合重金属理化结构参数和不同水生生物的致毒机理建立QSAR金属毒性预测模型，并将其应用于预测未知基准连续浓度的一种方法。本发明基于生态学原理，系统筛选六门八科水生物种作为最小生物预测集，分别构建单参数的毒性预测模型，提高模型精度和预测能力。

Description

基于金属定量构效关系的淡水慢性基准预测方法

技术领域

本发明涉及水质污染评价领域，具体为基于金属定量构效关系的淡水慢性基准预测方法。

背景技术

重金属污染日益严重。过量的重金属进入自然环境，给生态环境和水生生物造成有害影响。基准连续浓度是污染物风险评估中广泛采用的水环境质量基准指标，有效表征污染物的慢性毒性效应和生态风险。目前，美国环保局颁布了10种金属的基准连续浓度，包括铜，锌，汞，镍，镉，铬，铅，铝，铁，砷。其他金属由于毒性数据的匮乏尚未颁布，制约了科学评价水质，应急环境时间处置，污染控制和风险管理等方面的工作开展。DeForest等基于生物配位体模型开发出锌的急慢性毒性预测模型，并尝试应用于美国的基准制定中。但是总体而言，对于金属的慢性毒性效应和相关水生态基准的研究尚处于起步阶段。通过标准化的生物毒性测试手段，是目前获得基准值的最主要途径。但是，由于重金属的种类繁多，结构和形态复杂，用于基准推导的大量毒性测试需要耗费人力，物力和财力，因此阻碍了重金属水质基准研究的发展。相比于短期暴露，由于标准测试方法和模式生物等相关研究的滞后，慢性毒性终点信息更难获得。特别是针对特定物种，由于自身尺寸和缺少驯养的标准化流程，不可能进行重复可控的实验室测试。

定量结构活性相关(QSAR)方法采用统计分析手段寻找目标污染物的结构与生物活性间的内在联系，作为毒理机制研究的有效手段已被广泛应用于各类毒性效应的预测评价中。QSAR方法不受实验条件和测试仪器的限制，采用各种计算化学和数据挖掘技术来研究和预测污染物的生物活性，因而在面对批量污染物和传统毒理学测试难以开展的情况时，QSAR在毒性预测与风险评价方面展现出独特魅力。发明人围绕重金属的急性毒性和水质基准预测开展了相关创新性研究工作，提出了一种基于金属定量构效关系的淡水急性基准预测方法，根据重金属离子的结构特征与水生生物急性毒性效应的定量关系预测未知金属的毒性终点，结合不同物种的敏感度分布分析推导保护不同比例的水生生物的危险浓度。以上方法都是基于水生物种的急性毒性预测模型，缺乏对慢性毒性预测和分析，模型的预测能力和应用范围非常有限。与急性预测模型相比，在建模方法上存在差异，表现在数据收集和筛选方法，敏感生物种类，模型参数选择，模型结构，模型评价方法等不同环节。目前，关于重金属慢性毒性预测模型的研究尚少，并且未检索到与之相关的专利公布。

发明内容

本发明的目的在于提供一种基于金属定量构效关系的淡水慢性基准预测方法，用以克服上述技术缺陷。

为实现上述目的，本发明提供的基于金属定量构效关系的淡水慢性基准预测方法，根据重金属离子的结构特征与水生生物慢性毒性效应的定量关系预测未知金属的毒性终点，结合不同物种的敏感度分布分析推导保护不同比例的水生生物的危险浓度；

该具体过程包括以下步骤：

步骤a，建模毒性数据采集，筛选，运算和汇总；

步骤b，六门八科水生模式生物筛选；

步骤c，构建金属离子结构描述符数据集，通过各金属对应的结构参数为自变量进行线性相关性分析，通过相关系数排序，获得最佳结构描述符；

其中，构建金属离子结构描述符数据集，包括软指数σp、最大配合物稳定常数log-β_n、鲍林(Pauling)电负性X_m、共价指数X_m ²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR、Z/AR²。

以单物种的毒性终点为因变量，各金属对应的结构参数为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

$r=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2{(y_i-\stackrel{\‾}{y})}^2}}---\left(1\right)$

式中，分别表示各结构参数和毒性值的平均值，x_i和y_i分别表示第i种金属对应的结构参数和毒性值；相关系数r＞0.8为显著相关参数；

步骤d，构建毒性预测模型及稳健性检验；建立一元线性回归方程，对参数进行估计，采用F统计量对应的P值进行检验；

具体过程包括：

步骤d1，一元回归方程的构建与参数估计；

以步骤d中确定的最佳结构参数为自变量X，金属活性值为因变量Y，利用一元线性回归分析方法构建各模式生物的QICAR方程Y＝XB+E，具体为公式(2)：

$Y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_n\end{array}\right),X=\left(\begin{array}{cc}1& x_{11}\\ 1& x_{21}\\ .& .\\ .& .\\ .& .\\ 1& x_{n1}\end{array}\right),B=\left(\begin{array}{c}{\mathrm{\β}}_0\\ {\mathrm{\β}}_1\\ {\mathrm{\β}}_2\end{array}\right),E=\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_n\end{array}\right)---\left(2\right)$

n为观测值个数；

采用最小二乘法对方程中参数进行估计，X′为X的转置矩阵：

步骤d2，拟合优度检验和回归方程的显著性检验，采用F检验；

模型的拟合优度检验指标为：相关系数的平R²和自由度校正的相关系数标准偏差RMSE；

F检验的指标为单因子方差分析计算得到的F值和相关概率p；采用F统计量对应的P值进行检验；

步骤d3，判别标准：根据毒性数据获取途径，体外实验R²≥0.81，体内试验R²≥0.64；显著水平为α，当p＜α时，回归方程显著。

步骤d3按照下述公式计算，

$R^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}---\left(4\right)$

$\stackrel{\‾}{R^2}=1-\left({1-R}^2\right)\frac{n-1}{n-3}$

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

$F=\frac{[\mathrm{SS}\left(\mathrm{total}\right)-\mathrm{SS}\left(\mathrm{residual}\right)]/2}{\mathrm{SS}\left(\mathrm{residual}\right)/(n-3)}=\frac{[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2](n-3)}{2\×\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}---\left(6\right)$

式中，R²表示相关系数的平方，R²表示自由度校正的相关系数，RMSE表示标准偏差。

步骤e，QSAR模型的内部验证；

具体过程为：

步骤e1，在给定的建模样本中，抽取一个样本作为预测集，其余样本作为训练集建模，并计算该样本的预测误差；

步骤e2，记录每个方程中预测误差的平方加和，直到所有的样本都被预报了一次而且仅被预报一次；

步骤e3，计算交叉验证相关系数Q² _cv和交叉验证均方根误差，判别依据：Q² _cv＞0.6，R²-Q² _cv≤0.3。

步骤e3采用的计算公式为：

$Q_{\mathrm{CV}}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-{\stackrel{\‾}{y}}^{\mathrm{obs}})}^2}---\left(7\right)$

$\mathrm{RMSECV}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{n}}---\left(8\right)$

式中，y_i ^obs表示第i个化合物毒性的实测值，代表第i个化合物毒性的预测值，代表训练集毒性的平均值，n表示训练集中化合物数。

步骤f，模型适用范围计算；经过校验的模型，以杠杆值h为横坐标，以各数据点的标准残差为纵坐标，绘制Williams图；

杠杆值h_i的计算公式为：

h_i＝x_i ^T(X^TX)^-1x_i    (9)

式中，x_i代表第i个金属的结构参数组成的列向量；对于单参数模型， $x_i=\left(\begin{array}{c}{x}_{i1}\\ x_{i2}\end{array}\right),$ $X=\left(\begin{array}{c}{x}_{11}\\ x_{21}\\ .\\ .\\ .\\ x_{n1}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

临界值h*的计算公式为：

$h^{*}=\frac{3(p+1)}{n}---\left(10\right)$

式中，p代表模型中变量数，单参数模型中p＝1，n代表模型训练集化合物的数量，根据步骤a-e校验过后各QSAR方程中训练集金属个数决定；

在Williams图中h＜h*的坐标空间为模型的适用范围。

步骤g，采用获得的毒性预测值和物种敏感度分析对未知金属的毒性和基准连续浓度值进行快速筛选与预测。

具体过程为：

步骤g1，按照上述步骤a-f所述的方法，依次获得优选六门八科水生生物的单参数QSAR预测方程；

步骤g2，对于步骤b中筛选出的八个物种，重复步骤c-f的过程，获得八个典型模式生物的QSAR方程；搜集并整理待预测金属以上八个方程中出现的所有结构描述符的值，依次代入每个方程计算待预测金属对各物种的慢性毒性终点；

步骤g3，每种金属对应的各物种金属毒性数据由低到高排序后，以累积百分率为纵坐标构建物种敏感度分布图；

步骤g4，采用非线性Sigmoidal-Logistic拟合方程对曲线进行拟合，根据拟合方程计算累积百分率为0.05，0.1和0.2时对应的危险浓度HC₅，HC₁₀和HC₂₀。

本发明提供的基于金属定量构效关系的淡水慢性基准预测方法，基于生态学原理，系统筛选六门八科水生物种作为最小生物预测集，分别构建单参数的慢性毒性预测模型，提高模型精度和预测能力。通过QSAR模型方法预测多种金属的毒性值，快速、简单，依赖较少的试验测试数据完成多种毒性数据缺乏的金属的基准预测。弥补现有技术只对急性毒性终点基于多参数QSAR模型进行预测、预测结果不够准确的问题。

附图说明

图1为本发明的流程示意图；

图2为实施例的Williams图；

图3为实施例的汞的慢性毒性预测值的物种敏感度分布曲线。

具体实施方式

以下结合附图，对本发明上述的和另外的技术特征和优点作更详细的说明。

本发明的原理是根据重金属离子的结构特征与水生生物慢性毒性效应的定量关系预测未知金属的毒性终点，结合不同物种的敏感度分布分析推导保护5％、10％和20％水生生物的危险浓度。是综合重金属理化结构参数和不同水生生物的致毒机理建立QSAR金属毒性预测模型，并将其应用于预测未知基准参考值的一种方法。

如图1所示，其为本发明基于金属定量构效关系的淡水慢性基准预测方法的流程示意图，该具体过程为：

步骤a，建模毒性数据采集，筛选，运算和汇总；

步骤a1，数据采集过程：

步骤a2，数据筛选过程；数据筛选满足的条件为：

1)每个物种包含至少6种金属的毒性数据；

2)毒性终点数据类型包括致死率，生长率和繁殖率，表现为最大无效应剂量(NOAEL)，最低有效应剂量(LOAEL)或毒物最大可接受浓度(MATC)；

3)毒性测试必须在一定范围的环境条件下以标准的操作流程进行；

4)生物测试暴露时间7～14天。

步骤a3，数据运算过程；在本发明实施例中的运算方法为：

以金属自由离子浓度为数据的衡量指标，如单位为质量浓度除以分子量统一转化为摩尔浓度，即mol/L。LOAEL和MATC统一转换成NOAEL，换算公式为NOAEL＝2*LOAEL；NOAEL＝1.414*MATC。

步骤a4，数据汇总过程：

最终得到的数据集包括金属化合物分子式，毒性效应类型，终点指标，受试生物类型。最终的物种毒性终点为多个实验条件下NOAEL的几何平均值。

详细的毒性数据获取过程如下：

建模的慢性毒性数据优先采集自美国环保局ECOTOX毒性数据库(http://cfpub.epa.gov/ecotox/)。如果毒性数据不足，以近30年SCI科学引文索引查询的有效数据(ISI Web of Knowledge)作为补充。通过数据库和文献检索引擎，输入金属名称、待测物种名称和慢性毒性等关键词，导出满足条件的毒性数据集。在满足步骤a2条件的前提下，筛选出合格的毒性数据。以金属自由离子浓度为数据的衡量指标，如果原始数据以离子化合物质量为毒性终点指标。需除以分子量统一转化为微摩尔浓度，即μmol/L。LOAEL和MATC统一转换成NOAEL，换算公式为NOAEL＝2*LOAEL；NOAEL＝1.414*MATC。最终的物种毒性终点为多个实验条件下NOAEL的几何平均值。在数据汇编过程中，记录金属原子或分子式，原子或分子量，毒性效应类型，终点指标，受试生物类型等信息，整理成Excel表格作为建模依据。

对大型蚤的慢性毒性数据进行汇总，如表1所示

表1.慢性毒性数据筛选，运算和汇总范例

步骤b，六门八科水生模式生物筛选；

慢性模式生物以美国环保局推荐的推导水质基准的三门八科生物为基础，筛选出淡水中对重金属敏感的六门八科模式生物，包括浮游甲壳类节肢动物2种，脊索动物2种，软体动物，轮虫，绿藻和浮萍各1种。对于每一类模式生物，对应的毒性数据需严格按照数据采集和筛选的要求，依次汇总各物种慢性毒性数据。如果满足要求的物种数超过最少物种数要求，选择受试金属种类丰富的进行建模。例如，通过数据收集，浮游甲壳类生物中满足条件的物种类型有5种，按照受试金属元素的数量进行排序，选取前两种作为模式生物。在进行模式生物筛选后，确定八种生物的科学命名，所属门和科。

优选六门八科水生生物信息，如表2所示

表2优选重金属慢性效应评价的模式生物

物种	门	科
			摇蚊	节足动物门	摇蚊科
大型蚤	节足动物门	溞科
			鲤鱼	脊索动物门	鲤科
斑马鱼	脊索动物门	鲤科
			小球藻	绿藻门门	小球藻科
轮虫	轮形动物门	臂尾轮虫科
			贻贝	软体动物门	贻贝科
浮萍	被子植物门	天南星科

步骤c，构建金属离子结构描述符数据集；

构建金属离子结构描述符集合，包括软指数σp、最大配合物稳定常数log-β_n、鲍林(Pauling)电负性X_m、共价指数X_m ²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR、Z/AR²。

步骤c1，以单物种的毒性终点为因变量，各金属对应的结构参数为自变量进行线性相关性分析，根据下述公式(1)计算皮尔逊相关系数r；

$r=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2{(y_i-\stackrel{\‾}{y})}^2}}---\left(1\right)$

式中，x_i和y_i分别表示第i种金属对应的结构参数和实测毒性值，分别表示各结构参数和实测毒性值的平均值。相关系数r＞0.8为显著相关参数。采用皮尔逊相关可以简便，客观的度量两个因子之间的关联程度。

步骤c2，在显著相关的前提下，通过相关系数排序，获得最优的结构描述符。该步骤中通过相关系数r，筛选出与毒性显著相关的结构参数，避免了伪相关参数引入模型。

步骤d，构建毒性预测模型及稳健性检验；

步骤d1，一元回归方程的构建与参数估计；

以上述步骤d中确定的两最佳结构参数为自变量X，金属活性值为因变量Y，利用一元线性回归分析方法构建各模式生物的QICAR方程Y＝XB+E，请参阅下述公式(2)，其中：

$Y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_n\end{array}\right),X=\left(\begin{array}{cc}1& x_{11}\\ 1& x_{21}\\ .& .\\ .& .\\ .& .\\ 1& x_{n1}\end{array}\right),B=\left(\begin{array}{c}{\mathrm{\β}}_0\\ {\mathrm{\β}}_1\\ {\mathrm{\β}}_2\end{array}\right),E=\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_n\end{array}\right)---\left(2\right)$

n为观测值个数；B代表未知参数，是方程中需要通过最小二乘法进行估计的；E代表随机误差项，反映了除x₁，x₂对y的线性关系之外的随机因素对y的影响。与一元线性回归相比，方程(2)采用一元线性回归建立了两种不同结构参数与金属毒性值的关系，完整、准确地表达预测对象与相关因素的关系。

采用最小二乘法对方程中参数进行估计，X′为X的转置矩阵：

最小二乘回归是从误差拟合角度对回归模型进行参数估计，是一种标准的一元建模工具，尤其适用于预测分析。

步骤d2，拟合优度检验和回归方程的显著性检验(F检验)；

模型的拟合优度检验指标为：相关系数的平方(R²)和自由度校正的相关系数标准偏差(RMSE)。F检验的指标为单因子方差分析(ANOVA)计算得到的F值和相关概率p(Significance F)。通常采用F统计量对应的P值进行检验。

步骤d3，判别标准：根据毒性数据获取途径，体外实验R²≥0.81，体内试验R²≥0.64。显著水平为α，当p＜α时，回归方程显著。

$R^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}---\left(4\right)$

$\stackrel{\‾}{R^2}=1-\left({1-R}^2\right)\frac{n-1}{n-3}$

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

$F=\frac{[\mathrm{SS}\left(\mathrm{total}\right)-\mathrm{SS}\left(\mathrm{residual}\right)]/2}{\mathrm{SS}\left(\mathrm{residual}\right)/(n-3)}=\frac{[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2](n-3)}{2\×\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}---\left(6\right)$

式中，y_i表示第i种金属实测的毒性值，y表示第i种金属预测的毒性值，表示各毒性值的平均值，n为训练集中金属的个数。

方程(4)、(5)的相关系数和标准偏差可以度量回归直线的拟合优度；方程(6)是检验因变量与自变量的线性关系是否显著的通用方法。

预测金属汞对八科模式生物的毒性值，结合SSDs曲线预测基准参考阈值。

按照步骤a-d所述的方法，分别构建八科模式生物的毒性预测方程，如表3所示。计算汞的最优结构参数σp＝0.065。依次代入方程获得各物种的毒性预测值。

表3八科模式生物的QSAR毒性预测方程

步骤e，QSAR模型的内部验证；

每个物种的QSAR模型还应采用抽一法进行验证，方法的核心思想是随机从训练集中抽出一个数据，用其他的毒性数据和步骤c获得的最佳结构描述符建立一元回归模型，根据抽出数据的预测值与实验值的比较，来校验所建立的网络模型。为了减少交叉验证结果的可变性，对一个样本数据集进行多次不同的划分，得到不同的互补子集，进行多次交叉验证。本步骤中，取多次验证的平均值作为验证结果。

此内部验证方法的优势在于用几乎所有的样本来训练模型，最接近样本，这样评估所得的结果比较可靠；实验没有随机因素，整个过程是可重复的。

具体步骤如下：

步骤e1，在给定的建模样本中，抽取一个样本作为预测集，其余样本作为训练集建模，并计算该样本的预测误差；

步骤e2，记录每个方程中预测误差的平方加和，直到所有的样本都被预报了一次而且仅被预报一次；

步骤e3，计算交叉验证相关系数Q² _cv和交叉验证均方根误差RMSECV，计算公式如下所述；判别依据：Q² _cv＞0.6，R²-Q² _cv≤0.3；

$Q_{\mathrm{CV}}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-{\stackrel{\‾}{y}}^{\mathrm{obs}})}^2}---\left(7\right)$

$\mathrm{RMSECV}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{n}}---\left(8\right)$

式中，y_i ^obs表示第i个化合物毒性的实测值，代表第i个化合物毒性的预测值，代表训练集毒性的平均值，n表示训练集中化合物数。

方程(7)、(8)是抽一法内部验证的指示参数，可有效降低模型对训练集数据的过拟合，测定训练集中有无特定金属对模型稳健性的影响。

对模型进行内部验证。以大型蚤的慢性毒性预测方程log-NOEC＝(37.073±3.983)σp+(-4.312±0.408)为例，对模型进行抽一法内部验证，相关拟合参数见表4。根据步骤e中的公式(7)和(8)，计算Q² _cv＝0.880，RMSECV＝0.3482，R²-Q² _cv＝0.055。满足模型稳健性判别依据Q² _cv＞0.6，R²-Q² _cv≤0.3，该模型通过内部验证。

表4模型内部验证抽一法相关参数

金属	σp	观测值	截距	斜率	预测值
						汞	0.104	-0.203	-4.3225	36.7462	-0.501
银	0.131	0.229	-4.6276	40.9784	0.741

铜	0.065	-1.847	-4.4074	37.8879	-1.945
						锌	0.126	0.333	-4.3313	37.3171	0.371
铅	0.081	-1.511	-4.1604	35.9176	-1.251
						铁	0.115	0.283	-4.2165	35.5108	-0.133
镍	0.074	-1.668	-4.2062	36.2077	-1.527

步骤f，模型适用范围计算；

经过校验的模型，采用杠杆值法计算模型的适用范围，以Williams图直观表示。此方法可保证模型在预测过程中具有最佳的可靠性。

杠杆值h_i的计算公式为：

h_i＝x_i ^T(X^TX)^-1x_i    (9)

式中，x_i代表第i个金属的结构参数组成的列向量；对于单参数模型， $x_i=\left(\begin{array}{c}{x}_{i1}\\ x_{i2}\end{array}\right),$ $X=\left(\begin{array}{c}{x}_{11}\\ x_{21}\\ .\\ .\\ .\\ x_{n1}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

临界值h*的计算公式为：

$h^{*}=\frac{3(p+1)}{n}---\left(10\right)$

式中，p代表模型中变量数，单参数模型中p＝1，n代表模型训练集化合物的数量，根据步骤a-e校验过后各QSAR方程中训练集金属个数决定。

以杠杆值h为横坐标，以各数据点的标准残差为纵坐标，绘制Williams图。在图中h＜h*的坐标空间为模型的适用范围。

以斑马鱼的慢性毒性预测方程log-NOEC＝(41.156±4.432)σp+(-4.468±0.476)为例，训练集各金属的结构参数和毒性终点为如表5所示。临界值h*＝3*(1+1)/7＝0.86。

表5斑马鱼的慢性毒性预测方程适用范围的计算

种类	σp	h(σp)	观测值	预测值	残差
						铜	0.104	0.143	0.037	-0.188	-0.225
六价铬	0.157	0.599	2.125	1.800	-0.325

汞	0.065	0.372	-1.989	-1.793	0.196
						镍	0.126	0.225	0.083	0.718	0.635
镉	0.081	0.220	-1.067	-1.134	-0.067
						锌	0.115	0.165	0.585	0.265	-0.32
银	0.074	0.276	-1.334	-1.422	-0.088

以各金属最优结构参数的杠杆值为横坐标，预测残差为纵坐标绘制Williams图，如图2所示。图中三条虚线内部的空间为模型的适用范围，计算结果显示训练集的7种金属在模型的预测范围之内。

步骤g，采用获得的毒性预测值和物种敏感度分析对未知金属的毒性和基准预测值进行快速筛选与预测。

步骤g1，按照步骤a-f所述的方法，依次获得优选六门八科水生生物的单参数QSAR预测方程。

步骤g2，对于步骤b中筛选出的八个物种，重复步骤c-f的过程，获得八个典型模式生物的QSAR方程。搜集并整理待预测金属以上八个方程中出现的所有结构描述符的值，依次代入每个方程计算待预测金属对各物种的慢性毒性终点；

步骤g3，每种金属对应的各物种金属毒性数据由低到高，即最敏感物种到最不敏感物种排序后，以累积百分率为纵坐标(P＝(R-0.5)/N，R物种序号，N物种数)构建物种敏感度分布图。

步骤g4，采用非线性Sigmoidal-logistic拟合方程对曲线进行拟合(公式)，根据拟合方程计算累积百分率为0.05，0.1和0.2时对应的危险浓度HC₅，HC₁₀和HC₂₀。

曲线拟合优度的判别指标包括 F和P。计算方法见方程(4)-(6)。

$y=\frac{a}{1+e^{-k\left({x-x}_c\right)}}---\left(11\right)$

式中，a代表拟合曲线的振幅，x_c代表中心值，k代表曲线斜率。大量研究证实，非线性Sigmoidal-Logistic拟合模型对物种敏感度曲线的拟合效果最佳。故本发明采用此方法作为推导金属基准连续浓度值的方法。

具体的，获得金属汞的QSAR-SSDs曲线拟合方程：

$y=\frac{(2.338\±1.974)}{{1+e}^{-(2.075\±0.588)[x-(-0.834\±0.709)]}}$

评价拟合优度的各参数分别为：Adj.r²＝0.980，RSS＝0.009，F＝446.960，P＝2.07×10^-6。根据SSDs曲线，如图3。当y等于0.05，0.10和0.20时，对应的logHC₅，logHC₁₀和logHC₂₀的值为-2.678，-2.326，-1.965。美国环保局1985年发布的水质基准指南中，基于实验室测定推导出汞的危害浓度为-2.417，预测误差为0.261。

Claims (9)

1.基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，根据重金属离子的结构特征与水生生物慢性毒性效应的定量关系预测未知金属的毒性终点，结合不同物种的敏感度分布分析推导保护不同比例的水生生物的危险浓度；

该具体过程包括以下步骤：

步骤a，建模毒性数据采集，筛选，运算和汇总；

步骤b，六门八科水生模式生物筛选；

步骤c，构建金属离子结构描述符数据集，通过各金属对应的结构参数为自变量进行线性相关性分析，通过相关系数排序，获得最佳结构描述符；

以单物种的毒性终点为因变量，各金属对应的结构参数为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

$r=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2{(y_i-\stackrel{\‾}{y})}^2}}---\left(1\right)$

式中，分别表示各结构参数和毒性值的平均值，x_i和y_i分别表示第i种金属对应的结构参数和毒性值；相关系数r＞0.8为显著相关参数；

步骤d，构建毒性预测模型及稳健性检验；建立一元线性回归方程，对参数进行估计，采用F统计量对应的P值进行检验；

步骤e，QSAR模型的内部验证；

步骤f，模型适用范围计算；经过校验的模型，以杠杆值h为横坐标，以各数据点的标准残差为纵坐标，绘制Williams图；

步骤g，采用获得的毒性预测值和物种敏感度分析对未知金属的毒性和基准连续浓度值进行快速筛选与预测。

2.根据权利要求1所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述的步骤c中，构建金属离子结构描述符数据集，包括软指数σp、最大配合物稳定常数log-β_n、鲍林(Pauling)电负性X_m、共价指数Xm²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR、Z/AR²。

3.根据权利要求1所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述的步骤d的具体过程包括：

步骤d1，一元回归方程的构建与参数估计；

以步骤d中确定的最佳结构参数为自变量X，金属活性值为因变量Y，利用一元线性回归分析方法构建各模式生物的QICAR方程Y＝XB+E，具体为公式(2)：

$Y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_n\end{array}\right),X=\left(\begin{array}{cc}1& x_{11}\\ 1& x_{21}\\ .& .\\ .& .\\ .& .\\ 1& x_{n1}\end{array}\right),B=\left(\begin{array}{c}{\mathrm{\β}}_0\\ {\mathrm{\β}}_1\\ {\mathrm{\β}}_2\end{array}\right),E=\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_n\end{array}\right)---\left(2\right)$

n为观测值个数；

采用最小二乘法对方程中参数进行估计，X′为X的转置矩阵：

步骤d2，拟合优度检验和回归方程的显著性检验，采用F检验；

模型的拟合优度检验指标为：相关系数的平R²和自由度校正的相关系数标准偏差RMSE；

F检验的指标为单因子方差分析计算得到的F值和相关概率p；采用F统计量对应的P值进行检验；

步骤d3，判别标准：根据毒性数据获取途径，体外实验R²≥0.81，体内试验R²≥0.64；显著水平为α，当p＜α时，回归方程显著。

4.根据权利要求3所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，上述步骤d3按照下述公式计算，

$R^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}---\left(4\right)$

$\stackrel{\‾}{R^2}=1-(1-R^2)\frac{n-1}{n-3}$

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

$F=\frac{[\mathrm{SS}\left(\mathrm{total}\right)-\mathrm{SS}\left(\mathrm{residual}\right)]/2}{\mathrm{SS}\left(\mathrm{residual}\right)/(n-3)}=\frac{[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2](n-3)}{2\×\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\hat{y})}^2}---\left(6\right)$

式中，R²表示相关系数的平方，R²表示自由度校正的相关系数，RMSE表示标准偏差。

5.根据权利要求1所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述步骤e的具体过程为：

步骤e1，在给定的建模样本中，抽取一个样本作为预测集，其余样本作为训练集建模，并计算该样本的预测误差；

步骤e2，记录每个方程中预测误差的平方加和，直到所有的样本都被预报了一次而且仅被预报一次；

步骤e3，计算交叉验证相关系数Q² _cv和交叉验证均方根误差，判别依据：Q² _cv＞0.6，R²-Q² _cv≤0.3。

6.根据权利要求5所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，上述步骤e3采用的计算公式为：

$Q_{\mathrm{CV}}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-{\stackrel{\‾}{y}}^{\mathrm{obs}})}^2}---\left(7\right)$

$\mathrm{RMSECV}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{n}}---\left(8\right)$

式中，表示第i个化合物毒性的实测值，代表第i个化合物毒性的预测值，代表训练集毒性的平均值，n表示训练集中化合物数。

7.根据权利要求1所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述的步骤f中，杠杆值h_i的计算公式为：

h_i＝x_i ^T(X^TX)^-1x_i   (9)

式中，x_i代表第i个金属的结构参数组成的列向量；对于单参数模型， $x_i=\left(\begin{array}{c}{x}_{i1}\\ x_{i2}\end{array}\right),$ $X=\left(\begin{array}{c}{x}_{11}\\ x_{21}\\ .\\ .\\ .\\ x_{n1}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

8.根据权利要求7所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述的步骤f中，临界值h*的计算公式为：

$h^{*}=\frac{3(p+1)}{n}---\left(10\right)$

式中，p代表模型中变量数，单参数模型中p＝1，n代表模型训练集化合物的数量，根据步骤a-e校验过后各QSAR方程中训练集金属个数决定；

在Williams图中h＜h*的坐标空间为模型的适用范围。

9.根据权利要求1所述的基于金属定量构效关系的淡水慢性基准预测方法，其特征在于，所述的步骤g的具体过程为：

步骤g1，按照上述步骤a-f所述的方法，依次获得优选六门八科水生生物的单参数QSAR预测方程；

步骤g2，对于步骤b中筛选出的八个物种，重复步骤c-f的过程，获得八个典型模式生物的QSAR方程；搜集并整理待预测金属以上八个方程中出现的所有结构描述符的值，依次代入每个方程计算待预测金属对各物种的慢性毒性终点；

步骤g3，每种金属对应的各物种金属毒性数据由低到高排序后，以累积百分率为纵坐标构建物种敏感度分布图；

步骤g4，采用非线性Sigmoidal-Logistic拟合方程对曲线进行拟合，根据拟合方程计算累积百分率为0.05，0.1和0.2时对应的危险浓度HC₅，HC₁₀和HC₂₀。