CN104820873A - 一种基于金属定量构效关系的淡水急性基准预测方法 - Google Patents

Abstract

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

Description

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

技术领域

本发明涉及淡水水质基准模型领域，尤其涉及一种基于金属定量构效关系的淡水急性基准预测方法。

背景技术

金属污染是本世纪最具挑战性的环境问题之一。过量的重金属进入自然环境，破坏生物多样性，给生态环境和人体健康造成有害影响。为金属制定科学的环境基准是环境保护和风险评估的基础。国际上最先开展基准研究的国家是美国，我国现有基准体系主要照搬或借鉴国外成果，缺乏科学性。在最新的基准文件中，15种金属被列入优控污染物和非优控污染物名录，但只有10种金属具有基准参考值。大部分金属的水质基准值缺失，主要原因是生物毒性数据不足，其次受到环境要素的影响，只有铜，镍等金属基准的研究较为深入。目前，通过标准化的生物毒性测试手段，是目前获得基准值的唯一途径。但是，由于重金属的种类繁多，结构和形态复杂，用于基准推导的大量毒性测试需要耗费人力，物力和财力，并且复杂生物体系中金属形态难于准确测定，因此阻碍了重金属水质基准研究的发展。虽然也有研究者采用计算手段对各种毒性终点进行预测，但真正用于毒性和水质基准预测的还未见报道。发展不依赖于试验测定的基准预测方法，更符合我国的国情，节省大量的人力、物力和财力。

定量结构活性相关(QSAR)方法采用统计分析手段寻找目标污染物的结构与生物活性间的内在联系，作为毒理机制研究的有效手段已被广泛应用于各类毒性效应的预测评价中。QSAR方法不受实验条件和测试仪器的限制，采用各种计算化学和数据挖掘技术来研究和预测污染物的生物活性，因而在面对批量污染物和多受试物种时具有尤为明显的优势，QSAR正在毒性预测与风险评价方面展现出独特魅力。众所周知，离子形态是金属最为活跃的形态，溶解态金属的生物活性与自由离子浓度密切相关。在理想体系下，研究者尝试开展金属离子的定量构效关系研究，提出定量离子特征-活性相关预测其生物活性的方法。Newman等利用海洋发光细菌(V.fischeri)的毒性实验数据建立了QSAR方程并预测了金属毒性。结果发现，第一水解常数|logK_OH|与金属离子对生物体的毒性效应有很强的作用关系。Bogaerts等在评价原生动物(T.pyriformis)的毒性效应与金属离子理化特征之间的作用关系时也指出金属离子软指数σp为毒性预测方程最佳建模参数。

以上方法都是基于单一物种的单参数预测模型，缺乏对生态系统中多物种的系统性毒性预测和分析，模型的预测能力和应用范围非常有限。

鉴于上述缺陷，本发明创作者经过长时间的研究和实践终于获得了本创作。

发明内容

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

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

该具体过程为：

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

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

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

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

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

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

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

进一步地，在上述步骤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为显著相关参数。

进一步地，在上述步骤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²。

进一步地，上述步骤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}{ccc}1& x_{11}& x_{12}\\ 1& x_{21}& x_{22}\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& x_{n1}& x_{n2}\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检验的指标为多因子方差分析(Multi-ANOVA)计算得到的F值和相关概率p(Significance F)；采用F统计量对应的P值进行检验；

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

进一步地，上述步骤d3按照下述公式计算，

$\begin{array}{c}{R}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2}\\ \stackrel{\&OverBar;}{R^2}=1-(1-R^2)\frac{n-1}{n-3}\end{array}---\left(4\right)$

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(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{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2](n-3)}{2\&times;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2}---\left(6\right)$

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

进一步地，上述步骤e的具体过程为：

步骤e1，在给定的建模样本中，选取大部分样本进行建模型，留小部分样本用建立的模型进行预测，并计算这小部分样本的预测误差；

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

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

进一步地，上述步骤e3采用的计算公式为：

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

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

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

进一步地，在上述步骤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}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\\ .& .\\ .& .\\ .& .\\ x_{n1}& x_{n2}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

进一步地，在上述步骤f中，临界值h*的计算公式为：

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

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

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

进一步地，上述步骤g的具体过程为：

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

步骤g1，搜集并整理待预测金属在八个方程中出现的所有结构描述符的值，代入方程计算待预测金属对各物种的急性毒性终点；

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

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

与现有技术相比本发明的有益效果为：1、现有技术只对单一物种的毒性终点进行预测，模型预测不够准确，预测误差在两个数量级左右。本发明基于生态学原理，系统筛选五门八科水生物种作为最小生物预测集，分别构建多参数的毒性预测模型，提高模型精度和预测能力。

2.QSAR模型与SSD分析结合预测基准最大浓度(CMCs)

现有技术通过实验测试手段获得毒性终点值，再进行物种敏感度分析进而推导基准值。本专利通过QSAR模型方法预测多种金属的毒性值，快速、简单，依赖较少的试验测试数据完成多种毒性数据缺乏的金属的基准预测。

附图说明

图1为本发明的基于金属定量构效关系的淡水急性基准预测方法的流程示意图；

图2a为本发明模型适用范围评价的Williams图一；

图2b为本发明模型适用范围评价的Williams图一；

图3为本发明的汞毒性预测值的物种敏感度分布曲线。

具体实施方式

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

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

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

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

步骤a1，数据采集过程；

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

1)每个物种的急性毒性数据须来自同一试验来源，同一研究组和相同试验条件；

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

3)毒性终点数据类型包括致死率，生长率和繁殖率，表现为EC₅₀或LC₅₀；

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

5)生物测试暴露时间48～96小时。

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

以金属自由离子浓度为数据的衡量指标，如单位为质量浓度除以分子量统一转化为摩尔浓度，即mol/L。

步骤a4，数据汇总过程：

最终得到的数据集包括金属化合物分子式，受试生物类型，毒性效应类型，终点指标，试验条件，暴露时间，数据来源。

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

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

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

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

步骤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{\&Sigma;}}}(x_i-\stackrel{\&OverBar;}{x})(y_i-\stackrel{\&OverBar;}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-\stackrel{\&OverBar;}{x})}^2{(y_i-\stackrel{\&OverBar;}{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}{ccc}1& x_{11}& x_{12}\\ 1& x_{21}& x_{22}\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& x_{n1}& x_{n2}\end{array}\right),B=\left(\begin{array}{c}{\mathrm{\&beta;}}_0\\ {\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\end{array}\right),E=\left(\begin{array}{c}{\mathrm{\&epsiv;}}_1\\ {\mathrm{\&epsiv;}}_2\\ .\\ .\\ .\\ {\mathrm{\&epsiv;}}_n\end{array}\right)---\left(2\right)$

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

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

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

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

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

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

$\begin{array}{c}{R}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2}\\ \stackrel{\&OverBar;}{R^2}=1-(1-R^2)\frac{n-1}{n-3}\end{array}---\left(4\right)$

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(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{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2](n-3)}{2\&times;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2}---\left(6\right)$

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

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

步骤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{\&Sigma;}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i^{\mathrm{obs}}-{\stackrel{\&OverBar;}{y}}^{\mathrm{obs}})}^2}---\left(7\right)$

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

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

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

步骤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}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\\ .& .\\ .& .\\ .& .\\ x_{n1}& x_{n2}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

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

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

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

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

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

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

步骤g1，搜集并整理待预测金属在八个方程中出现的所有结构描述符的值，代入方程计算待预测金属对各物种的急性毒性终点。

步骤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(x-x_c)}}---\left(11\right)$

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

以下通过实施例结合附图对本发明进一步说明。

实施例1：

采用本发明步骤a所述的方法，对大型蚤的急性毒性数据进行汇总，如表1所示

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

实施例2：

采用本发明步骤b所述的方法，优选五门八科水生生物信息，如表2所示

表2 优选急慢性模式生物

实施例3：

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

按照步骤a-d所述的方法，分别构建八科模式生物的毒性预测方程，如表3所示。计算汞的最优结构参数σp＝0.065，log-β_n＝21.7，X_m ²r＝4.08，AN/ΔIP＝9.62，Z/r＝1.96，|logK_OH|＝3.4，ΔE₀＝0.91。依次代入方程获得各物种的毒性预测值。

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

实施例4：

采用本发明步骤e所述的方法对模型进行内部验证。以大型蚤的急性毒性预测方程log-EC₅₀＝(-0.272±18.674)σp+(-0.360±0.136)log-β_n+(6.604±4.093)为例，对模型进行抽一法内部验证，相关拟合参数见表4。根据步骤e中的公式(7)和(8)，计算Q² _cv＝0.63，RMSECV＝1.139，R²-Q² _cv＝0.239。满足模型稳健性判别依据Q² _cv>0.6，R²-Q² _cv≤0.3，该模型通过内部验证。

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

实施例5：

采用本发明步骤f所述的方法计算模型适用范围，绘制Williams图。以鲤鱼的急性毒性预测方程log LC₅₀＝(33.439±6.256)σp+(0.412±0.137)Z/r+(-3.159±0.559)为例，训练集各金属的结构参数和毒性终点为如表5所示。临界值h*＝3*(2+1)/6＝1.5。

表5 鲤鱼的急性毒性预测方程适用范围的计算

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

实施例5：

根据本发明步骤g所述，获得金属汞的QSAR-SSDs曲线拟合方程：

$y=\frac{(0.9273\&PlusMinus;0.0593)}{1+e^{-(3.3526\&PlusMinus;0.5625)[x-(-0.7808\&PlusMinus;0.0647)]}}$

评价拟合优度的各参数分别为：Adj.r²＝0.9594，RSS＝0.019，F＝231.176，P＝1.18×10^-5。根据SSDs曲线(如图3)，当y等于0.05，0.10和0.20时，对应的logHC₅，logHC₁₀和logHC₂₀的值为-1.6352，-1.4022，-1.1658。美国环保局1985年发布的水质基准指南中，基于实验室测定推导出汞的危害浓度为-1.8560，预测误差为0.119。

上述详细说明是针对本发明其中之一可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本发明技术方案的范围内。

Claims (10)

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

该具体过程为：

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

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

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

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

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

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

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

2.根据权利要求1所述的基于金属定量构效关系的淡水急性基准预测方法，其特征在于，在上述步骤c中，以单物种的毒性终点为因变量，各金属对应的结构参数为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

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

式中，分别表示各结构参数和毒性值的平均值，x_i和y_i分别表示第i种金属对应的结构参数和毒性值；

相关系数r>0.8为显著相关参数。

3.根据权利要求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²。

4.根据权利要求2所述的基于金属定量构效关系的淡水急性基准预测方法，其特征在于，上述步骤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}{ccc}1& x_{11}& x_{12}\\ 1& x_{21}& x_{22}\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& x_{n1}& x_{n2}\end{array}\right),B=\left(\begin{array}{c}{\mathrm{\&beta;}}_0\\ {\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\end{array}\right),E=\left(\begin{array}{c}{\mathrm{\&epsiv;}}_1\\ {\mathrm{\&epsiv;}}_2\\ .\\ .\\ .\\ {\mathrm{\&epsiv;}}_n\end{array}\right)---\left(2\right)$

n为观测值个数；

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

$\hat{B}=\left(\begin{array}{c}\widehat{{\mathrm{\&beta;}}_0}\\ \widehat{{\mathrm{\&beta;}}_1}\\ ...\\ \widehat{{\mathrm{\&beta;}}_m}\end{array}\right)={\left(X^{\&prime;}X\right)}^{-1}{X}^{\&prime;}Y---\left(3\right)$

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

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

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

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

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

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

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

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(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{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2](n-3)}{2\&times;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-\hat{y})}^2}---\left(6\right)$

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

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

步骤e1，在给定的建模样本中，选取大部分样本进行建模型，留小部分样本用建立的模型进行预测，并计算这小部分样本的预测误差；

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

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

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

$Q_{\mathrm{CV}}^2=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i^{\mathrm{obs}}-y_i^{\mathrm{predcv}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_o^{\mathrm{obs}}-{\stackrel{\&OverBar;}{y}}^{\mathrm{obs}})}^2}---\left(7\right)$

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

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

8.根据权利要求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}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\\ .& .\\ .& .\\ .& .\\ x_{n1}& x_{n2}\end{array}\right);$ X^T表示矩阵X的转置矩阵，(X^TX)^-1表示对X^TX矩阵的逆矩阵。

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

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

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

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

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

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

步骤g1，搜集并整理待预测金属在八个方程中出现的所有结构描述符的值，代入方程计算待预测金属对各物种的急性毒性终点；

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

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