CN105069315A - 基于金属形态和有效性的水生生物毒性预测方法 - Google Patents

Abstract

本发明涉及环境中有毒物质预测领域，具体为基于金属形态和有效性的水生生物毒性预测方法，根据金属的结构特征、外部环境条件与水生生物毒性效应的定量关系预测未知金属的“原位”毒性终点；是综合金属的生物有效性和特殊致毒机理建立的金属毒性耦合预测模型，并将其应用于预测未知金属毒性效应终点的一种方法。本发明基于金属的生物有效性和致毒机理，构建金属毒性预测模型，通过QSAR模型方法预测未知毒性值，快速、简单，依赖较少的试验测试数据完成多种毒性数据缺乏金属的“原位”毒性终点预测。

Description

基于金属形态和有效性的水生生物毒性预测方法

技术领域

本发明涉及环境中有毒物质预测领域，具体为基于金属形态和有效性的水生生物毒性预测方法。

背景技术

定量结构活性相关(QSAR)理论已取得重要进展，并在有机污染物的构效关系和毒性预测方面得到广泛应用。金属的定量构效关系研究始于上世纪七八十年代，但由于金属的形态和生物效应的复杂性，导致相关研究遇到了瓶颈，也是金属毒性预测领域的难点和挑战。例如，有效反映致毒机理的结构描述符难于获得，制约着QSAR在金属毒性预测中的应用。目前，获得金属毒性效应数据的唯一途径依然是通过传统毒理学试验测定，不仅耗费大量的人力、物力和财力，也必须面对生态伦理的挑战。进行金属定量构效关系和建模方法的研究对预测和评估生物效应具有广阔的应用前景。

近年来，随着金属形态分析技术的进步和对致毒机理的认识逐渐深入，为金属QSAR研究提供了克服瓶颈的新途径。基于致毒机理相同的科学假设，在相对理想的体系下，可以基于已有的QSAR模型经济、高效地实现一系列金属生物效应的纵向预测。目前的金属毒性预测技术基于金属的物理化学参数初步构建了一系列定量离子特征-毒性模型。然而，已报道的研究工作仅关注于金属自身的结构特点，而忽视了外界水化学条件对金属生物效应的影响，QSAR模型不能实现对不同环境条件下生物效应的横向预测。考虑到金属的形态和生物有效性对金属毒性的影响，已开发出生物配位体(BioticLigandmodel，BLM)模型。模型的主要思想是金属毒性以金属在生物敏感受体-鱼鳃上的累积过程模拟，包括金属形态和竞争性阳离子的防护效应。影响毒性的所有水化学特征包括在内，整合了金属离子与Ca²⁺、Na⁺、Mg²⁺、H⁺的竞争作用，与非生物配体DOC、氯化物、碳酸盐和硫化物结合。该模型比较全面地考虑了影响金属生物有效性的因素，和建立在硬度基础上的模型相比，生物配体模型具有更好的预测能力，是一个能够替代生物毒性试验预测金属毒性的有用工具。到目前为止，模型中最全面的数据是关于铜对虹鳟鱼的毒性。Playle等发现鳃膜上的钠离子通道是铜急性毒性效应的生物配位体，通过在宽泛水质条件范围内测定铜在鱼鳃表面的累积来校准BLM模型中铜的吸收。RebortSantore发展了BLM模型，能够预测五种金属对五种水生生物的毒性效应。马义兵等对环境中金属BLM模型研究进展进行了较为全面的综述，并应用BLM模型预测我国多个水体的重金属毒性。Balistrieri等借鉴了生物配位体模型的原理，进一步预测了金属混合物对鳟鱼的复合毒性，同时指出温度、pH、主要离子浓度和溶解有机碳共同影响复合毒性。

但是BLM模型的应用也存在一些限制，主要包括以下方面：多数物种的BLM模型还未开发出来，相关的平衡系数有待率定和评估；针对不同物种的环境条件范围存在差异，超出范围时BLM模型不再适用；痕量金属的吸收要受生物学的控制，急性毒性的预测优于慢性。自然界的重金属并不是单独存在的，模型没有考虑多种重金属的复合效应对这种金属毒性的影响；针对特定金属开发的BLM模型，不能直接通过参数校正预测其他金属的毒性。因此，在BLM模型中引入金属结构特征参数，将有效扩展模型的应用范围。在不同外部环境条件下预测的毒性效应更贴近实际环境背景，使水生生物的保护更有针对性和科学性。

总体而言，以上方法只考虑金属离子结构特征与毒性的定量关系，忽视了水化学特征对毒性的影响。对于金属的生物有效性和原位金属毒性的定量预测，缺乏系统的研究和可靠的预测方法。

发明内容

本发明的目的在于提供一种基于金属形态和有效性的水生生物毒性预测方法，用以解决忽视了水化学特征对毒性的影响而导致预测结果不准确的问题。

为实现上述目的，本发明提出了一种基于金属形态和有效性的水生生物毒性预测方法，包括以下步骤：

步骤a，建立QSAR耦合模型采集、筛选、运算和汇总毒性数据；

步骤b，构建金属结构描述符数据集和对应的环境条件指标，以金属的结构参数和环境指标的耦合参数为自变量进行线性相关性分析和主成分分析，获得最优的金属结构描述符集合与环境指标耦合参数；

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

步骤d，QSAR耦合模型的内部验证；

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

步骤f，对未知金属和特定环境条件的水生生物毒性进行快速筛选与预测。

作为优选，所述步骤b中，构建金属结构描述符数据集，包括金属离子软指数σp、单位电荷软指数σP/Z、原子序数AN、离子半径r、IP：O_N态离子的离子势、IP(N+1)：O_N+1态离子的离子势、IP(N+1)和IP的差值△IP、原子半径R、原子量AW、鲍林电负性X_m、共价指数X_m ²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、实测电负性x、极化率z/rx、离子价Z、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR和Z/AR²；

其中，环境条件指标包括温度、酸度pH、硬度和盐度。

作为优选，所述步骤b具体包括以下过程：

步骤b1，以毒性终点为因变量，各金属对应的结构参数和环境指标为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

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

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

当相关系数r≥0.8判定为结构参数与毒性显著相关；

步骤b2，在显著相关的前提下，通过主成分分析，获得最优的金属结构描述符集合，具体公式：

F＝a_1i*Z_X1+a_2i*Z_X2+……+a_pi*Z_Xp(2)

其中a_1i，a_2i，……，a_pi(i＝1，……，m)为X的协方差阵Σ的特征值所对应的特征向量，m为主成分数，p为观测指标数，Z_X1，Z_X2，……，Z_Xp是原始变量经过标准化处理的值；

A＝(a_ji)_p×m＝(a₁，a₂，…，a_m)(3)

Ra_i＝λ_ia_i(4)

其中，A为特征向量矩阵(j＝1，……，p)，R为相关系数矩阵，λ_i、a_i是相应的特征值和单位特征向量，λ₁≥λ₂≥…≥λ_p≥0。

作为优选，所述的步骤c包括以下过程：

步骤c1，多元回归方程的构建与耦合参数估计；

步骤b中确定的最优的金属结构描述符集合和环境指标为自变量矩阵X，金属水生毒性值为因变量矩阵Y，利用多元线性回归分析方法构建QSAR耦合模型方程Y＝XB+E，见公式(5)，其中：

$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(5\right)$

n为观测值个数，Y为被解释变量的观测值向量，X为解释变量的观测值矩阵，B为总体回归参数向量，E随机误差项向量；

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

$B=\left(\begin{array}{c}{\mathrm{\β}}_0\\ {\mathrm{\β}}_1\\ \mathrm{...}\\ {\mathrm{\β}}_m\end{array}\right)={\left(X^{\′}X\right)}^{-1}{X}^{\′}Y---\left(6\right)$

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

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

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

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

作为优选，所述步骤c3按照下述公式计算，

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

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

$RMSE=\sqrt{\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}{n-3}}---\left(8\right)$

$F=\frac{\&lsqb;SS\left(total\right)-SS\left(residual\right)\&rsqb;/2}{SS\left(residual\right)/(n-3)}=\frac{\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2\&rsqb;(n-3)}{2\&times;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}---\left(9\right)$

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

作为优选，所述步骤d的具体过程包括：

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

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

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

作为优选，所述步骤d3采用的计算公式为：

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

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

式中，y_i ^obs表示第i个及金属毒性的实测值，代表第i个金属毒性的预测值，代表训练集毒性的平均值。

作为优选，所述步骤e中，杠杆值h_i的计算公式为：

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

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

作为优选，所述步骤e中，当杠杆值h_i超过了某个预警值后，这个化合物的结构与用于标定的化合物存在显著差异，表现在化合物位于最优预测空间之外，以这个预警值作为计算模型预测范围的临界值h*，计算公式为：

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

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

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

作为优选，其特征在于，所述步骤f中具体过程为，按照上述步骤a-步骤e所述的方法获得基于金属生物有效性的QSAR耦合模型预测方程，搜集并整理待预测金属的所有结构描述符和对应的环境条件指标，即可代入方程计算出待预测金属的“原位”毒性终点。

本发明提供的基于金属形态和生物有效性的毒性预测方法，基于最优离子结构参数和致毒机理，合理引入影响金属形态的水化学特征，建立定量构效关系-生物有效性耦合模型。通过QSAR模型方法预测未知“原位”毒性值，快速、简单，依赖较少的试验测试数据完成多种毒性数据缺乏金属的毒性终点预测。

附图说明

图1为本发明的预测方法的流程示意图；

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

图3为本发明模型适用范围评价的Williams图之二

图4为本发明模型适用范围评价的Williams图之三。

具体实施方式

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

本发明的原理是根据金属的结构特征、外部环境条件与水生毒性建立定量QSAR耦合模型并进一步预测未知金属的“原位”毒性终点。是综合金属的生物有效性和致毒机理建立金属“原位”毒性预测模型，并将其应用于预测未知金属的毒性终点的一种方法。

如图1所示，其为本发明一种基于金属形态和有效性的水生生物毒性预测方法的流程示意图，该具体过程为：

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

步骤a1，数据采集过程；

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

1)急慢性毒性数据自近10年文献报道和美国环保局ECOTOX数据库；

2)每个毒性数据必须包括科学，完整的试验条件参数，优选温度，酸度，硬度和盐度等参数，不同水化学参数的范围是温度10～30℃之间，酸度5.5～8之间，硬度(CaCO₃mg/L)20～5000mg/L之间。

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

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

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

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

步骤a4，数据汇总过程：

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

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

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

以杂色鳉鱼(Cyprinodonvariegatus)的急性毒性终点为例，进行数据筛选，运算和汇总，结果如表1所示。

表1杂色鳉鱼急性毒性数据筛选，运算和汇总范例

步骤b，构建金属结构描述符数据集和对应的环境条件指标，以结构参数和环境指标为自变量进行线性相关性分析和主成分分析，获得最优的结构描述符与环境指标耦合参数；

其中，构建金属结构描述符数据集，分别为金属离子软指数σp、单位电荷软指数σP/Z、原子序数AN、离子半径r、IP：O_N态离子的离子势、IP(N+1)：O_N+1态离子的离子势、IP(N+1)和IP的差值△IP、原子半径R、原子量AW、鲍林(Pauling)电负性X_m、共价指数X_m ²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、实测电负性x、极化率z/rx、离子价Z、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR、Z/AR²。

环境条件指标包括温度，酸度(pH)，硬度(CaCO₃，mg/L)和盐度。

步骤b1，以毒性终点为因变量，各金属对应的结构参数和环境指标为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

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

式中，x_i和y_i分别表示第i种金属对应的结构参数和实测毒性值，分别表示各结构参数和实测毒性值的平均值。相关系数r>0.8为显著相关参数。采用皮尔逊相关可以简便，客观的度量两个因子之间的关联程度。该步骤中通过相关系数r，筛选出与毒性显著相关的结构参数，避免了伪相关参数引入模型。

按照步骤b1所述的方法，分别计算各结构参数的皮尔逊相关系数r，如表2所示，由于表格较大，拆分为表2(上)和表2(下)，但其实是一个表格。

表2各参数之间的皮尔逊相关系数(上)

表2各参数之间的皮尔逊相关系数(下)

步骤b2，在显著相关的前提下，通过主成分分析，获得最优的结构描述符组合。具体公式：

F＝a_1i*Z_X1+a_2i*Z_X2+……+a_pi*Z_Xp(2)

其中a_1i，a_2i，……，a_pi(i＝1，……，m)为X的协方差阵Σ的特征值所对应的特征向量，m为主成分数，p为观测指标数，Z_X1，Z_X2，……，Z_Xp是原始变量经过标准化处理的值；

A＝(a_ji)_p×m＝(a₁，a₂，…，a_m)(3)

Ra_i＝λ_ia_i(4)

其中，A为特征向量矩阵(j＝1，……，p)，R为相关系数矩阵，λ_i、a_i是相应的特征值和单位特征向量，λ₁≥λ₂≥…≥λ_p≥0。

主成分分析是对于原先提出的所有变量，将关系紧密的变量删去多余，建立尽可能少的新变量，使得这些新变量是两两不相关的，而且这些新变量尽可能保持原有的信息。

按照步骤b2所述的方法，分别计算各结构参数的主成分贡献，如表3所示。

表3各结构参数的主成分分析

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

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

步骤b中确定的优选结构参数和环境指标为自变量矩阵X，金属水生毒性值为因变量矩阵Y，利用多元线性回归分析方法构建QSAR方程Y＝XB+E，见

公式(5)，其中： $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(5\right)$

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

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

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

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

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

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

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

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

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

$RMSE=\sqrt{\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}{n-3}}---\left(8\right)$

式中，R²表示相关系数的平方，表示自由度校正的相关系数，RMSE表示标准偏差。方程(7)、(8)的相关系数和标准偏差可以度量回归直线的拟合优度；

$F=\frac{\&lsqb;SS\left(total\right)-SS\left(residual\right)\&rsqb;/2}{SS\left(residual\right)/(n-3)}=\frac{\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2\&rsqb;(n-3)}{2\&times;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}---\left(9\right)$

式中，y_i表示第i种金属实测的毒性值，表示第i种金属预测的毒性值，表示各毒性值的平均值，n为训练集中金属的个数。方程(9)是检验因变量与多个自变量的线性关系是否显著的通用方法。

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

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

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

具体步骤如下：

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

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

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

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

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

式中，y_i ^obs表示第i个及金属毒性的实测值，代表第i个金属毒性的预测值，代表训练集毒性的平均值。

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

采用本发明步骤d所述的方法对模型进行内部验证。以预测方程log-C＝(5.4539±2.8414)+(-0.2109±0.0709)T+(-0.0308±0.0078)S+(13.7692±10.3215)σp为例，对模型进行抽一法内部验证，相关拟合参数见表4。根据步骤d3中的公式(7)和(8)，计算Q² _cv＝0.8773，RMSECV＝0.2935，R²-Q² _cv＝0.946-0.8773＝0.0687。满足模型稳健性判别依据Q² _cv>0.6，R²-Q² _cv≤0.3，该模型通过内部验证。

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

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

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

杠杆值h_i的计算公式为：

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

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

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

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

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

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

训练集各金属化合物的结构参数和毒性终点为如表5所示。临界值h*＝3*(3+1)/10＝1.2。

表5金属化合物“原位”毒性预测模型适用范围的计算

以各金属的温度，盐度和软指数σp的杠杆值为横坐标，预测残差为纵坐标绘制Williams图，如图2、图3和图4所示。图中三条虚线内部的空间为模型的适用范围，计算结果显示训练集的5种金属化合物和外部实验条件在模型的适用范围之内。

步骤f，对未知金属和特定环境条件的水生生物毒性进行快速筛选与预测。

具体过程为，按照上述步骤a-步骤e所述的方法获得QSAR预测方程，搜集并整理待预测金属的所有结构描述符和目标环境条件，代入方程计算待预测的毒性终点。

Claims (10)

1.一种基于金属形态和有效性的水生生物毒性预测方法，其特征在于，包括以下步骤：

步骤a，建立QSAR耦合模型采集、筛选、运算和汇总毒性数据；

步骤b，构建金属结构描述符数据集和对应的环境条件指标，以金属的结构参数和环境指标的耦合参数为自变量进行线性相关性分析和主成分分析，获得最优的金属结构描述符集合与环境指标耦合参数；

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

步骤d，QSAR耦合模型的内部验证；

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

步骤f，对未知金属和特定环境条件的水生生物毒性进行快速筛选与预测。

2.根据权利要求1所述的基于金属形态和生物有效性的毒性预测方法，其特征在于，所述步骤b中，构建金属结构描述符数据集，包括金属离子软指数σp、单位电荷软指数σP/Z、原子序数AN、离子半径r、IP：O_N态离子的离子势、IP(N+1)：O_N+1态离子的离子势、IP(N+1)和IP的差值△IP、原子半径R、原子量AW、鲍林电负性X_m、共价指数X_m ²r、原子电离势AN/ΔIP、第一水解常数|logK_OH|、电化学势ΔE₀、原子大小AR/AW、实测电负性x、极化率z/rx、离子价Z、极化力参数Z/r、Z/r²、Z²/r、似极化力参数Z/AR和Z/AR²；

其中，环境条件指标包括温度、酸度pH、硬度和盐度。

3.根据权利要求1或2所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤b具体包括以下过程：

步骤b1，以毒性终点为因变量，各金属对应的结构参数和环境指标为自变量进行线性相关性分析，根据下述公式(1)计算相关系数r；

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

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

当相关系数r≥0.8判定为结构参数与毒性显著相关；

步骤b2，在显著相关的前提下，通过主成分分析，获得最优的金属结构描述符集合，具体公式：

F＝a_1i*Z_X1+a_2i*Z_X2+……+a_pi*Z_Xp(2)

其中a_1i，a_2i，……，a_pi(i＝1，……，m)为X的协方差阵Σ的特征值所对应的特征向量，m为主成分数，p为观测指标数，Z_X1，Z_X2，……，Z_Xp是原始变量经过标准化处理的值；

A＝(a_ji)_p×m＝(a₁，a₂，…，a_m)(3)

Ra_i＝λ_ia_i(4)

其中，A为特征向量矩阵(j＝1，……，p)，R为相关系数矩阵，λ_i、a_i是相应的特征值和单位特征向量，λ₁≥λ₂≥…≥λ_p≥0。

4.根据权利要求1所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述的步骤c包括以下过程：

步骤c1，多元回归方程的构建与耦合参数估计；

步骤b中确定的最优的金属结构描述符集合和环境指标为自变量矩阵X，金属水生毒性值为因变量矩阵Y，利用多元线性回归分析方法构建QSAR耦合模型方程Y＝XB+E，见公式(5)，其中：

$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(5\right)$

n为观测值个数，Y为被解释变量的观测值向量，X为解释变量的观测值矩阵，B为总体回归参数向量，E随机误差项向量；

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

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

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

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

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

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

5.根据权利要求4所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤c3按照下述公式计算，

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

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

$RMSE=\sqrt{\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}{n-3}}---\left(8\right)$

$F=\frac{\&lsqb;SS\left(total\right)-SS\left(residual\right)\&rsqb;/2}{SS\left(residual\right)/(n-3)}=\frac{\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-\stackrel{\&OverBar;}{y})}^2-\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2\&rsqb;(n-3)}{2\&times;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-y)}^2}---\left(9\right)$

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

6.根据权利要求1所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤d的具体过程包括：

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

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

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

7.根据权利要求6所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤d3采用的计算公式为：

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

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

式中，y_i ^obs表示第i个及金属毒性的实测值，代表第i个金属毒性的预测值，代表训练集毒性的平均值。

8.根据权利要求1所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤e中，杠杆值h_i的计算公式为：

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

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

9.根据权利要求8所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤e中，当杠杆值h_i超过了某个预警值后，这个化合物的结构与用于标定的化合物存在显著差异，表现在化合物位于最优预测空间之外，以这个预警值作为计算模型预测范围的临界值h*，计算公式为：

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

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

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

10.根据权利要求1所述的基于金属形态和有效性的水生生物毒性预测方法，其特征在于，所述步骤f中具体过程为，按照上述步骤a-步骤e所述的方法获得基于金属生物有效性的QSAR耦合模型预测方程，搜集并整理待预测金属的所有结构描述符和对应的环境条件指标，即可代入方程计算出待预测金属的“原位”毒性终点。