CN118800357B - 有机化学品的低营养级食物链生物放大预测模型构建方法 - Google Patents

Abstract

本发明公开了一种有机化学品的低营养级食物链生物放大预测模型，该模型为首次基于QSAR模型构建的PBDEs和HBCDs的低营养级食物链生物放大预测模型。所述模型通过样本采集及筛选、分子描述符计算、模型构建、模型验证等步骤获得。利用本方法构建出的有机化学品的低营养级食物链生物放大预测模型能准确预测出PBDEs和HBCDs类化学品有机污染物生物放大因子，提高了预测结果的准确率，节省了人力、物力和时间，简单、快速有效，并且严格按照OECD规定的QSAR模型使用规则，从分子描述符结构上解释影响生物放大因子的关键因素，对PBDEs和HBCDs等毒害化学物质的风险管控和环境安全具有重要意义。

Description

有机化学品的低营养级食物链生物放大预测模型构建方法

技术领域

本发明涉及有机化学品环境风险评估领域，具体涉及一种有机化学品的低营养级食物链生物放大预测模型构建方法。

背景技术

有毒有害化学物质具有亲脂性、难降解性，且在生物体内积累和富集。在生态系统中，随着营养级的增加呈现生物放大效应，对高等生物及人体产生毒性。因此研究水生生物体内有机污染物的含量特征为当地居民的健康和食品安全具有重要意义。

化合物通过食物链在高营养级生物和人体中的富集程度对评价该化合物对生态、环境毒性意义重大。多年来，科学家一直从事相关工作，通过长期的研究，他们发现并建立各种介质间化合物积累模型。评价有机污染物是否存在生物富集效应一般有两个标准。第一是化合物的K_OW(正辛醇-水分配系数)，一般地，当logK_OW>4-5时，该化合物可能具有生物富集效应，而logK_OW在5-7时，化合物具有最大的生物富集效应；第二是生物富集因子(BAF)。生物富集因子(BAF)可以表征一个化合物的相对生物可富集能力。

水生食物链中生物富集因子BAF的计算公式如下：

C_生物表示水生生物体中污染物的浓度，单位pg/kg lw；C_{水中溶解相}表示水体中溶解的污染物的浓度，单位pg/L。生物中的浓度为脂肪归一化浓度。

在BAF模型中，当某种污染物BAF值高于5000(或LogBAF>3.7)时，可认为该污染物在食物链中具有生物富集效应；当BAF值在2000～5000(或Lo gBAF>3.7)时，认为其有潜在的生物富集效应。

PBDEs和HBCDs作为斯德哥尔摩公约新增列的一类新型持久性有机污染物，具有典型的生物富集性。由于PBDEs和HBCDs具有高亲脂性和代谢能力差等特点，所以具有沿食物链放大的潜力。而最近的研究也证实，PBDEs和HBCDs可以在食物链中随营养级传递。但是仅通过实验方法获取化学品生物放大因子(BMF)成本高、费时费力，难以满足化学物质生态风险性评价的需要。目前，针对化学物质在食物链上的生物放大效应模型还空缺。所以，迫切需要发展科学快捷有效的食物链传递模型的理论计算方法。经济合作与发展组织(OECD)于2007年发布了QSAR模型构建与验证的导则，提出了QSAR模型应满足的标准：①具有明确定义的环境指标；②具有清晰和明确的数学算法；③定义了模型的应用域；④模型具有适当的拟合优度、稳健性和预测能力；⑤尽可能进行模型机理解释。

发明内容

本发明首次基于QSAR模型，构建PBDEs和HBCDs的低营养级食物链生物放大预测模型，如下所示：

PEC_{oral,predator} ＝ PEC_water*BCF_fish*BMF (1)

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v (2)

其中各参数含义如下：

本发明提供了上述低营养级食物链生物放大模型的构建方法，具体如下：

⑴样本采集及筛选

实验室数据获得包含了9个有机化合物的生物放大因子数据，这些化合物涵盖了PBDEs和HBCDs的有机物。为了建立有效的QSAR模型，首先把数据集分成训练集和测试集。为保证训练集化合物的代表性，本次工作所用的分组方法是Kennard&Stone方法，这种方法在一定程度上能够避免训练集样本分布不均匀，能够很好的将数据集划分为训练集和测试集，本数据集划分为6个训练集和3个验证集。

⑵分子描述符计算

本方法首先在ChemDraw软件中构建出9个有机化合物的分子结构，然后导入HyperChem程序对分子进行优化。优化分为两步：首先是MM+分子力场方法进行初步的能量优化，然后使用半经验量子力学AM1方法对结构进行更加准确的构型优化。优化后的结构导入到DRAGON5.4软件中计算1664个不同类型的理论分子描述符。建模前对这些描述符进行了预处理，即将常数项、接近常数的项和具有高度相关(相关系数大于0.96的两个分子描述符中与目标值相关系数较小的)的分子描述符删除。最终剩余1169个描述符用于后面的变量选择过程。

⑶模型构建

采用遗传算法来选择与生物富集具有高度相关的描述符集，这个过程在MobyDigs中实现。经过遗传算法变量选择后，用多元线性回归(MLR)方法建立线性QSAR模型，即GA-MLR模型。模型评价函数选择留一法交互检验(leave-one-out cross validation)，即当增加一个描述符后模型的性能没有明显变化时(增加一个描述符Q2增加小于0.02)，即达到最佳描述符个数。本方法中，最佳描述符个数为7。建模中的相关参数设置为：种群大小(population size)为100，初始模型允许的做大变量数(maximum allowed variables)为7，变异均衡值(mutation trade-off，T)为0.5，交叉(crossover)和变异(mutation)概率均基于T参数。

(4)模型验证

经过遗传算法变量选择后，用多元线性回归方法建立线性QSAR模型，即MLR模型。线性MLR方程如下：

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

n_tr＝6Q² _LOO＝0.9981R² _fitting＝0.9996R² _adj＝0.9994RMSE_tr＝0.0665R² _boot＝0.7175

n_ext＝3R² _ext＝0.836，Q² _ext＝0.8662R² _adj＝0.9994RMSE_ext＝0.0289

其中，GGI3表示3阶拓扑电荷指数，Mor21v表示3D-MoRSE-加权原子范德华体积，GGI3与生物放大因子呈正相关性，Mor21v与生物放大因子呈负相关性，n_tr和n_ext分别为训练集和验证集化合物数量；R² _adj为经自由度校正的决定系数；RMSE_tr为训练集均方根误差；Q² _LOO为去一法交叉验证系数；R² _boot为Bootstrapping方法验证系数；R² _ext为验证集实验值和预测值之间的决定系数，Q² _ext为外部验证决定系数，RMSE_ext为验证集均方根误差,SE为标准误差。R² _fitting是训练集实验值和预测值之间的决定系数。训练集和验证集RMSE分别为0.0665和0.0289，模型预测效果较好。

表1PBDEs和HBCDs BMF实验值和预测值

物质	CAS	训练集/验证集	实验值	预测值
					BDE100	189084-64-8	训练集	4.07	4.07
BDE153	68631-49-2	验证集	1.14	1.81
					BDE154	207122-15-4	训练集	3.05	3.05
BDE183	207122-16-5	验证集	1.5	1.33
					BDE28	41318-75-6	训练集	3.6	3.65
BDE47	5436-43-1	训练集	3.46	3.42
					BDE99	60348-60-9	验证集	3.85	4.53
BPE209	145538-74-5	训练集	6.04	6.05
					HBCD	3194-55-6	训练集	7.26	7.25

本方法建立的有机化学品的低营养级食物链生物放大预测模型的优点在于：由实验手段测定水生生物体中污染物的浓度PBDEs和HBCDs含量水平、水体中溶解的污染物的浓度，再计算得到化学品在食物链上的生物放大时间长、成本高。利用本方法构建出的有机化学品的低营养级食物链生物放大预测模型能准确预测出PBDEs和HBCDs类化学品有机污染物生物放大因子，提高了预测结果的准确率，节省了人力、物力和时间，简单、快速有效，并且严格按照OECD规定的QSAR模型使用规则，从分子描述符结构上解释影响生物放大因子的关键因素，对PBDEs和HBCDs等毒害化学物质的风险管控和环境安全具有重要意义。

附图说明

图1BMF预测模型拟合图。

图2BMF预测模型预测模型应用域表征图。

具体实施方式

以下结合具体实施例对本发明作进一步说明。

实施例1

有机化学品的低营养级食物链生物放大预测模型构建具体步骤如下：

(1)数据收集，设置训练集和验证集样本化合物

实验室获得9种有机化学品生物放大因子BMF。训练集共选取6个样本化合物，验证集共选取3个样本化合物。

(2)计算描述符

采用Hyperchem 7.0软件中的MM+分子力学对化合物结构进行预优化，半经验AM1方法对化合物结构进行优化，基于优化的结构，使用Dragon 5.4软件计算描述符，并对计算的1664个描述符进行初步筛选，

(3)模型构建

采用MobyDigs软件中遗传算法(GA)进行变量选择，基于筛选的变量，采用多元线性回归(MLR)方法建立预测模型，即GA-MLR模型：

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

其中，GGI3表示3阶拓扑电荷指数，Mor21v表示3D-MoRSE-加权原子范德华体积，GGI3与生物放大因子呈正相关性，Mor21v与生物放大因子呈负相关性，训练集和验证集RMSE分别为0.0665和0.0289。

(4)模型验证

根据OECD关于QSAR模型的导则，需要对构建的模型进行内部验证(拟合优度和稳健性评估)和外部验证(预测能力评估)。采用校正后的实验值与拟合值之间的相关系数平方(R² _adj)、均方根误差(RMSE)来表征模型的拟合优度：

其中，n代表化合物的个数，m为预测变量的个数，y_i和分别表示第i个化合物活性指标的实验值和预测值；为化合物活性指标实验值的平均值。

采用去一法交叉验证系数(Q² _LOO)和Bootstrapping法(Q² _BOOT)表征模型的稳定性：

其中，表示训练集化合物活性指标实验值的平均值。Bootstrapping法采用去1/5交叉验证，重复5000次。

采用外部验证相关系数(Q² _EXT),R² _EXT,RMSE_EXT表征模型预测能力：

其中n_EXT代表验证集化合物个数，表示验证集化合物活性指标实验值和预测值的平均值。得到模型的表征与评价参数：

n_tr＝6Q² _LOO＝0.9981R² _fitting＝0.9996R² _adj＝0.9994RMSE_tr＝0.0665R² _boot＝0.7175

n_ext＝3R² _ext＝0.836，Q² _ext＝0.8662R² _adj＝0.9994RMSE_ext＝0.0289

其中，n_tr和n_ext分别是训练集和验证集化合物数，p为显著性水平。R² _adj为经自由度校正的决定系数；RMSE为均方根误差；Q² _LOO为去一法交叉验证系数；Q² _BOOT为Bootstrapping方法验证系数；R² _ext为实验值和预测值相关系数，Q² _ext为外部验证决定系数，RMSE_ext为验证集均方根误差。

结果表明，模型具有较好的预测能力和稳健性。

实施例2

本实施例对上述预测模型进行应用域表征。

Williams图是由标准残差(δ)和杠杆值(以h_i表示，i代表不同的化合物)定义的一种模型应用域。δ采用下式计算：

训练集化合物的杠杆值(leverage，h_i)可以通过下面的公式求得：

h_i ＝ x_i ^T (X^TX)^–1 x_i (8)

式中，x_i为第i个化合物分子结构描述符的行向量。警戒值(h^*)定义为：

h^* ＝ 3(k + 1)/n (9)

其中，k为描述符的个数，n为训练集个数。

模型应用域表征结果如图1、图2所示。图1中h^*＝3(k+1)/n＝3(2+1)/6＝1.5。Williams图纵坐标用实验值和预测值的标准残差来表征实验值的离散程度，当化合物的标准残差δ的绝对值大于3.0时，被视为离群点。横坐标代表训练集中化合物h_i值，h_i大于警戒值(h*＝1.5)时，说明在训练集中该物质的子结构出现较少，会对模型预测结果有显著影响。

由图可见，所有化合物的杠杆值h在警戒杠杆值h^*内，表明这个化合物的结构与训练集化合物的结构有一定的相似性，标准残差均在(-3，+3)范围内，说明本模型适用于BDE153(CAS:68631-49-2),BDE 183(CAS:207122-16-5)和BDE-99(CAS:60348-60-9)的标准残差均落在(-3，+3)范围内，表明本模型适用于这三种物质的预测，能够很好地被预测，见图2。

实施例3

利用实施例1构建的模型对9种持久性有机污染物食物链放大因子进行预测，结果见表2。模型的R² _adj＝0.999，表明模型具有较强的拟合能力。Q² _LOO＝0.9981，Q² _BOOT＝0.7175，表明模型较为稳健。R² _ext＝0.9994，Q² _ext＝0.8662，Golbraikh等人研究认为，QSAR模型可接受的标准是Q²>0.50和R²>0.60。结果表明，模型具有较好的预测能力，能够成功应用到训练集以外的化合物中。图1为持久性有机化学品食物链放大因子BMF的预测值与实验值的拟合图，从图1中可见大部分物质的预测值和实验值拟合较好，BDE153(CAS:68631-49-2),BDE183(CAS:207122-16-5)和BDE-99(CAS:60348-60-9)能够被较好地预测。

实施例4

利用实施例1构建的模型，预测BDE 203(SMILES：BrC1＝C(OC2＝CC(Br)＝C(Br)C(Br)＝C2Br)C(Br)＝CC(Br)＝C1Br)的生物放大因子BMF。首先根据化学物质分子结构，使用Dragon软件计算出2种描述符GGI3和Mor21v，分别为1.937，-0.06，Hat为0.435，在模型应用域范围内。

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

BMF＝-5.04472+0.8374*(1.937)-35.46426*(-0.06)＝3.25

则BDE 203(BMF)预测值为3.25，与试验测定结果(1.73)接近。

实施例5

利用实施例1构建的模型，预测BDE 196(SMILES：BrC1＝CC＝C(OC2＝C(Br)C(Br)＝C(Br)C(Br)＝C2Br)C(Br)＝C1Br)的生物放大因子BMF。首先根据化学物质分子结构，使用Dragon软件计算出2种描述符GGI3和Mor21v，分别为2.375，-0.016，Hat为0.651，在模型应用域范围内。

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

BMF＝-5.04472+0.8374*(2.375)-35.46426*(-0.016)＝2.05

则BDE 196(BMF)预测值为2.05，与试验测定结果(1.43)接近。

Claims (1)

1.一种有机化学品的低营养级食物链生物放大预测模型的构建方法，其特征在于，所述低营养级食物链生物放大预测模型如下：

PEC_{oral,predator} ＝ PEC_water*BCF_fish*BMF (1)

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v (2)

其中，PEC_{oral,predator}指低营养级捕食者体内浓度，mg·kg_wetfish ^-1；PEC_water指水中预测浓度，mg/L；BCF_fish指鱼的生物富集因子，L·kg_wetfish ^-1；BMF指生物放大因子；GGI3指3阶拓扑电荷指数；Mor21v指3D-MoRSE-加权原子范德华体积；

所述有机化学品的低营养级食物链生物放大预测模型的构建方法，包括以下步骤：

(1)样本采集及筛选；

(2)分子描述符计算；

(3)模型构建；

(4)模型验证；

步骤(1)具体为：实验室数据获得包含了9个有机化合物的生物放大因子数据，这些化合物涵盖了PBDEs和HBCDs的有机物，首先把数据集分成训练集和测试集，所用的分组方法是Kennard&Stone方法，划分为6个训练集和3个验证集；

步骤(2)具体为：首先在ChemDraw软件中构建出9个有机化合物的分子结构，然后导入HyperChem程序对分子进行优化；优化分为两步：首先是MM+分子力场方法进行初步的能量优化，然后使用半经验量子力学AM1方法对结构进行更加准确的构型优化，优化后的结构导入到DRAGON5.4软件中计算1664个不同类型的理论分子描述符；建模前对这些描述符进行预处理，即将常数项、接近常数的项和具有高度相关的分子描述符删除，最终剩余1169个描述符用于后面的变量选择过程；

步骤(3)具体为：采用遗传算法来选择与生物富集具有高度相关的描述符集，这个过程在MobyDigs中实现；经过遗传算法变量选择后，用多元线性回归MLR方法建立线性QSAR模型，模型评价函数选择留一法交互检验，即当增加一个描述符后模型的性能没有明显变化时，即达到最佳描述符个数；本方法中，最佳描述符个数为7；建模中的相关参数设置为：种群大小population size为100，初始模型允许的最大变量数maximum allowed variables为7，变异均衡值mutation trade-off，T为0.5，交叉crossover和变异mutation概率均基于T参数；

步骤(4)具体为：经过遗传算法变量选择后，用多元线性回归方法建立线性QSAR模型，即MLR模型，线性MLR方程如下：

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

n_tr＝6Q² _LOO＝0.9981R² _fitting＝0.9996R² _adj＝0.9994RMSE_tr＝0.0665R² _boot＝0.7175

n_ext＝3R² _ext＝0.836，Q² _ext＝0.8662R² _adj＝0.9994RMSE_ext＝0.0289

其中，GGI3表示3阶拓扑电荷指数，Mor21v表示3D-MoRSE-加权原子范德华体积，GGI3与生物放大因子呈正相关性，Mor21v与生物放大因子呈负相关性，训练集和验证集RMSE分别为0.0665和0.0289；n_tr和n_ext分别为训练集和验证集化合物数量；R² _adj为经自由度校正的决定系数；RMSE_tr为训练集均方根误差；Q² _LOO为去一法交叉验证系数；R² _boot为Bootstrapping方法验证系数；R² _ext为验证集实验值和预测值之间的决定系数，Q² _ext为外部验证决定系数，RMSE_ext为验证集均方根误差,SE为标准误差，R² _fitting是训练集实验值和预测值之间的决定系数。