冯长君，堵锡华

徐州工程学院化学化工学院，江苏徐州221111

生物富集一般指生物体从周围环境中吸收不降解或难降解的化学物质，在体内积累，使这些物质在生物体内的浓度超过其生长环境中浓度的现象.生物富集程度常用生物富集因子(bioconcentration factors，BCF)表示，即生物体内污染物的平衡浓度(cS)与其生存环境中该污染物浓度(cH)的比值

目前，对于非离子性有机物在鱼体内生物富集实验数据的处理多采用“二室模型”［1］，即将实验测定的模拟池分为水室和鱼室，假定有机物在池内的迁移转化为一级动力学过程，可推得类似公式

其中，cw和cf分别为平衡时有机物在水体及鱼体内浓度(单位:mol/L).

污染物在不同生态系统中的含量，随食物链营养级的升高而增加，其富集系数在各营养级中均可达极其惊人的数值.处于食物链顶端的人类，便成为生物富集的最终受害者.如生物富集已引起包括水俣病和痛痛病在内的多起生态公害事件.

卤代苯(halogenated benzenes，HB)属于环境外来品，是化工、医药、制革和电子等行业广泛应用的化工原料、有机合成中间体或有机溶剂难降解而在环境中高度滞留.此类物质通过食物链进入人体，具有强烈的致癌、致畸和致突变作用，对人体造成严重危害.因此，其中的氯苯、邻二氯苯、对二氯苯和六氯苯已被列入美国及欧洲共同体环境保护部门所确定的129种优先控制污染物名单.通过研究生物富集作用，利于阐明物质在生态系统内的迁移和转化规律，评价和预测污染物进入生物体后可能造成的危害，在利用生物体对环境进行监测和净化及制定排放标准等方面，也具有重要的意义.

由于实测BCF成本高、周期长，在实际工作中通常需要采用估算方法来获取所需数据.已进行较多研究的估算法主要包括BCF与正辛醇-水分配系数(Kow)、水溶解度(SW)和吸附系数(Koc)的各种经验关系式及分子拓扑法等［2-3］.虽然 Kow、SW和Koc等关系式精度较高，但不足之处在于参数本身是实验测定值，不便应用.分子拓扑法建立BCF模型虽不依赖于实验，但所用的拓扑参数无明确物理意义，不便探讨影响生物富集的结构因素及富集机理.

物质定量构效关系(quantitative structure-activity relationship，QSAR)研究［4-6］不仅是数据建模，显示变量间的数学关系;更重要的是模型的有效性和实用性.有效性即模型质量，可用一系列统计回归指标来衡量;实用性是在有效性基础上，模型所蕴含的物理化学意义.因此，选用具有明确物理化学意义的描述符对于建立良好的QSAR模型至关重要.目前研究多卤联苯、多卤联苯醚(包含卤代烃和卤代苯)生物富集因子构效关系的较多［7-8］，单独探讨卤代苯BCF的QSAR模型尚未见报道.因此，本研究利用密度泛函理论(density functional theory，DFT)［9-11］对卤代苯分子进行量化计算，获得不需要实验测定，无实验误差且具有明确物理意义的量子化学参数(SL).然后基于最佳变量子集回归方法建立上述量子化学参数与其lg BCF［8］的最佳数学模型，用其探寻影响卤代苯生物富集行为的结构因素，揭示生物富集微观机理，提供理论诠释参考.

1 理论方法

1.1 卤代苯及其生物富集因子

从文献［8］中选取21种卤代苯，其名称和生物富集因子数据见表1.测试物种的鱼类有虹鳟鱼、孔雀鱼、黑头呆鱼和蓝鳃太阳鱼等.因不同卤代苯间BCF值相差比较大，且据热力学原理，采用其对数值lg BCF表示更为科学.lg BCF数值越大，表示该物质越易被富集.

1.2 量化参数计算

运用DFT的Becke三参数混合泛函(B3LYP)方法，在较高基组6-311G＊＊(d，p)水平上对卤代苯类化合物分子结构进行几何优化和理论计算.对所选化合物分子进行结构优化及频率计算，验证皆无虚频，表明其构型存在真实极小点;然后对以上优势构型进行单点能计算，从中提取一些相关的量化参数及热力学数据，均以SL表示.

1.3 变量筛选与模型验证

通过DFT方法获得的卤代苯类化合物分子的量化描述符SL之间一定存在或多或少的信息重叠.为遵循奥卡姆剃刀原则，尽可能建立简单的QSAR模型，须对其进行筛选，以建立具有良好稳定性和较高预测能力的QSAR模型.目前对变量压缩与筛选的方法较多，如逐步回归、最佳变量子集回归和遗传算法等.本研究基于最佳变量子集回归方法对变量SL进行压缩建立QSAR模型.模型质量采用留一法交叉验证(leave-one-out cross validation，LOOCV)和留多法交叉验证(leave-multiple-out cross validation，LMOCV)的相关系数(Rcv-o2和Rcv-m2)作为模型质量的判断标准.为确立最终的QSAR模型，引入Akaike信息判据AC(Akaike's information criterion)和Kubinyi函数 KT(Kubinyi function)［12-13］，其计算公式为

其中，RSS为方差和;f为化合物数;b为变量数.AC值越小，KT值越大，所建的模型越稳定，预测能力越高.

2 结果讨论

2.1 lg BCF的QSAR模型

采用DFT方法，在B3LYP/6-311G＊＊(d，p)基组上对21种卤代苯分子进行结构优化及频率计算，从中提取一些用于定量描述卤代苯-受体分子间作用力的量化及热力学参数:①能量参数有最高占有轨道能(the energy of the highest occupied molecular orbital，EHOMO)、次最高占有轨道能(the energy ofthe nexthighestoccupied molecular orbital，ENHOMO)、最低空轨道能(the energy of lowest unoccupied molecular orbital，ELUMO)、次最低空轨道能(the energy of the next lowest unoccupied molecular orbital，ENLUMO)及两者轨道能级间隙 ΔE1(ΔE1=ELUMO-EHOMO)和ΔE2(ΔE2=ENLUMO-ENHOMO)，单位均为Hartree;②电子性质包括分子偶极距μ(单位:Debye)、分子的可极化率α(单位:a.u.)、苯环上各碳原子的净电荷(指原子核带的正电荷与聚集在该原子上的电子所带的负电荷之和)Q1～Q6及其碳原子净电荷之和Qsum;③热力学性质有分子的零点振动能(Ev)、热力学能(E)、热容(Cv)、熵(Sm)、电子运动空间(VE)及分子配分函数Q(取Q的常用对数为lg Q)等.

基于最佳变量子集回归方法对21个卤代苯的生物富集因子与SL进行变量压缩建立QSAR模型，结果见表2，其中，R2、Radj2、Rcv-o2、S、F、AC和KT分别为判定系数(亦称削减误差比例)、校正判定系数(以消除自变量个数及样本容量对判定系数的影响)、逐一剔除法的交叉验证系数、估计标准误差、Fisher统计值、Akaike信息判据和Kubinyi函数.表2为建立高稳定性、高预测能力的QSAR模型提供指导.

表2 lg BCF与SL的最佳变量子集回归结果Table2 The results of SLand lg BCFwith Leaps-and-Bounds regression

由表2可见，随着模型变量增多，对应R2和Radj2持续增大，S持续减小，Rcv-o2、F和KT均呈先增后减趋势，而AC变化趋势则完全反之.但其出现转折点的位置不一，对应最大Rcv-o2的模型为

从式(5)可见，描述符ENLUMO、lg Q和Q2的回归系数均是其标准偏差的2倍以上，而VE的回归系数不足其标准偏差的2倍，说明VE的引入对模型的稳定是不利的.在置信度为95%时，标准t值(tα/2)为 2.080.模型 lg BCF中 ENLUMO、lg Q、Q2和VE的t值依次为 -11.356、 -7.711、2.495 和1.681.可见，VE的绝对值小于标准t值，进一步证明该方程是不稳健的.扣除VE后的三元数学模型

从式(6)可见，描述符ENLUMO、lg Q和Q2的回归系数分别是其标准偏差的10.6、8.0和2.6倍，这说明它们的引入对模型的稳定性是有利的;结合AC与KT均在此处发生转折，可认为该模型具有良好的稳定性与预测能力.对比式(5)和(6)，可以发现其常数项和描述符ENLUMO、lg Q和Q2的回归系数变化均很小，表明VE因素对lg BCF的贡献很小，可以舍弃.因此采用3个变量建模既保证良好的稳定性与预测能力，又符合奥卡姆剃刀原则，即尽可能选择简单的模型.

2.2 模型预测能力和稳定性的再验证

在置信度为95%时，标准t值(tα/2)为2.0796.模型 lg BCF中 ENLUMO、lg Q和Q2的t值依次为-10.674、 -7.931 和2.694，它们的绝对值均大于tα/2，证明方程(6)是稳健的.按此模型给出的估算值见表1，与相应实验值吻合.若要建立一个具有良好预测能力的QSAR模型，其R2≥0.9［14］.模型(6)的R2=0.922 ＞0.9，且 Rcv-o2=0.866 ＞ 0.5［15］，可认为该模型具有良好的稳健性与预测能力.

将模型(6)的样本随机分为2个样本:第4、8、12、16和18个化合物为检验集，余下为训练集.对于训练集的模型为此模型的R2、和S与模型(6)很接近.按照模型(7)给出训练集的估算值和检验集的预测值，均与相应实验值很接近，如图1.

图1 卤代苯化合物生物富集因子(lg BCF)的实验值与计算值的关系Fig.1 Relationship between experimental and calculated bioconcentration factors(lg BCF)of halogenated benzenes

为进一步验证所建QSAR模型的预测能力，把21个卤代苯分成训练集(16个化合物)和测试集(5个化合物)，共形成20 349个训练集与测试集(模型(7)只是其中之一).通过Matlab编程计算这些训练集、测试集的相关系数及逐五剔除法(leave-five-out，LFO;即每次剔除5个卤代苯分子，用余下16个分子建模)的交叉验证相关系数(Rcv-f2).对于20 349个测试集的Rcv-f2=0.855，与=0.866很接近;相应20 349个训练集的R2的平均值为0.925，与模型(6)的0.922非常接近.

为定量反映模型的稳健性，作者曾定义相关系数(R或R2)的相对标准偏差(RS)予以衡量

其中，Ri2为模型i的相关系数;Ra2为相关系数的平均值;w为自由度数，等于相关系数的个数减1.显然，RS越小，该模型的稳健性越高.根据相关回归分析中的显著性检验，通常确定显著性水平α=0.05.为此，建议RS＜0.05，说明Ri2变异性很小，呈现明显的集中趋势，反映所建模型具有良好的稳健性.20 349个训练集的Ri2的RS=0.022 ＜0.05，表明模型(6)是稳健的，不存在异常值.

2.3 生物富集机理分析

根据进入模型(6)中3个变量ENLUMO、lg Q和Q2的标准化回归系数(standardized regression coefficients，SR)依次为 -2.531、 -1.904 和0.203;从SR的绝对值可见，对lg BCF影响顺序是ENLUMO＞lg Q ＞ Q2.lg BCF分别与ENLUMO、lg Q、Q2回归，它们的判定系数 R2依次为 0.624、0.372 和0.172，此与SR结果一致.据此可知:

1)根据分子轨道理论，前线轨道(HOMO和LUMO)及其附近的分子轨道对分子的化学性质及生物活性影响最大.其中，LUMO及其附近的未占据轨道(LUMO+x)等具有优先接受电子的作用，能级越低，越易接受电子.模型(6)中ENLUMO越小，越易接受鱼体中生物大分子的电子，发生亲电反应而被富集，lg BCF越大.即ENLUMO与lg BCF呈负相关，故ENLUMO前系数小于0.

2)配分函数lg Q反映分子大小和对称性，分子越小、对称性越大，lg Q越小.卤代苯分子越大，越难扩散透过鱼鳞表面的磷脂膜，即穿过鱼鳞膜进入血液后溶解在肝脏内的生物富集量较小.可见，lg Q与lg BCF亦呈负相关.生物富集机理中可假定包含有机物分子镶嵌在生物受体提供的“空穴”内的过程，显然分子的对称性与空穴的“契合”程度越高，该有机物的lg BCF值越大.对于同分异构体的lg BCF值不同，便是此“假定”的佐证.例如1，2，3-三氯苯、1，2，4-三氯苯和1，3，5-三氯苯的分子对称性不同，与“空穴”的吻合程度不同，相应 lg BCF亦不同，依次为2.90、2.95 和3.26.

3)形成细胞膜的磷脂属于两亲性分子，亲水性头部(磷酸基)带负电荷指向细胞膜外，疏水性尾部(脂肪酸)带正电荷指向细胞膜内.Q2是苯环第2个碳原子上的净电荷量，其值越高即正电荷量越大，与带负电荷的亲水性磷酸基的静电引力越强，对鱼鳞膜的亲和力越强，越易进入鱼体内而生物富集量越大.Q2值越低即带负电荷量越大，与带负电荷的亲水性磷酸基的静电斥力越大，且与鱼鳞膜外围水分子形成氢键能力越强，越难穿过鱼鳞膜进入鱼体，相应lg BCF越小.即Q2与lg BCF成正相关，此与模型(6)中Q2前系数大于零是一致的.由此可见，C2是卤代苯分子中的活性部位，因其受取代基影响最为敏感.

综上研究，卤代苯分子在鱼体内生物富集的微观历程可表述为:溶解于水体的卤代苯，首先要穿过鱼鳞膜，即从水室进入鱼室，显然其分子体积及所带负电荷越大，不利于透过此生物膜;穿过鱼鳞膜随血液进入肝脏的卤代苯分子，其分子对称性与生物受体提供的“空穴”的“契合”程度越高，以及ENLUMO越低，越易发生亲电反应而使lg BCF越大.

2.4 热力学分析

根据热力学原理，任何物理过程、化学反应还是生物过程的推动力都源于焓作用和熵作用.对于生物反应中的焓作用包含静电作用、立体作用(如分子大小和对称性等);熵作用有疏水作用、转动熵、平动熵以及构象熵等.静电作用又分为离子-离子作用、离子-偶极作用、偶极-偶极作用、氢键以及电荷转移作用等，其中蕴含轨道反应及电子得失.由机理讨论可见，模型(6)和模型(7)涉及过程的焓作用，表面上似乎没有熵作用.

卤代苯分子(HB)在鱼体内生物受体(biological receptors，BRs)的生物富集机理可表示为

该过程的第1步是卤代苯在水体中的溶解过程，此是物理过程;第2步是卤代苯在水体中与在鱼体中达平衡的富集过程，其平衡常数(Kbio)表达式为

其中，［BRs］为纯生物受体浓度，取值为1 mol/L.定温下，此生物反应的吉布斯自由能变(ΔGbio)

其中，ΔHbio为富集过程焓变;ΔSbio为富集过程熵变;R为气体常数;T为反应温度.比较模型(6)、(7)与式(13)可见，其ENLUMO、lg Q和Q2主要对应ΔHbio;而常数项则呈现ΔSbio作用，即疏水性作用及分子平动、转动和构象的熵变.由于常数项是负数，即ΔSbio＜0，此是熵减过程.即各自独立运动的水合卤代苯分子、生物受体大分子相互结合成HB-BRs而损失的熵(ΔSbio'＜0)，大于体系中增加无序水分子而获得的熵(ΔSw＞0):ΔSbio=ΔSbio'+ΔSw＜0.实际上，此ΔSw在卤代苯溶解过程中是被抵消的.因此，整个生物富集过程是熵减过程:ΔSbio=ΔSbio'＜0.

综上所述，卤代苯类化合物在鱼体内的生物富集过程中是焓变起到主要及正向作用，熵变(主要是疏水性作用)则起次要及负向作用，这也表明其生物富集是有序性提高的熵减过程.

结 语

基于DFT方法，在B3LYP/6-311G＊＊(d，p)基组上对21种卤代苯分子进行结构优化及频率计算，从中提取了量化及热力学参数SL.基于最佳变量子集回归方法，建立了它们在鱼体内lg BCF的QSAR模型，经R2、Rcv-o2、Rcv-f2、F、AC和FT等质量指标检验，显示良好的相关性、稳健性及预测能力.

根据进入模型中的量化参数分析，卤代苯分子在鱼体内生物富集的微观机理可能是:卤代苯分子体积及所带负电荷越大，不利于穿过鱼鳞膜;进入肝脏的卤代苯分子与生物受体内“空穴”的“契合”程度越高，以及ENLUMO越低，越易发生亲电反应而使lg BCF越大.

经热力学分析，卤代苯类化合物在鱼体内的生物富集过程中的焓变起到主要及正向作用，熵变(主要是疏水性作用)则起次要及负向作用，这也表明其生物富集是有序性提高的熵减过程.

/References:

［1］Hou Ling，Hu Changmin，Zhao Xiaoming，et al.Simulating studies on bioaccumulation and elimination of the organic compounds of river in the cyprinus carpio［J］.Environmental Science，1997，18(2):34-36.(in Chinese)侯 玲，胡长敏，赵晓明，等.江河中有机污染物在鲤鱼体内富集与释放的模拟研究 ［J］.环境科学，1997，18(2):34-36.

［2］Toropova Alla P，Toropov Andrey A，Benfenati Emilio，et al.CORAL:quantitative models for estimating bioconcentration factor of organic compounds［J］.Chemometrics and Intelligent Laboratory Systems，2012，118(3):70-73.

［3］Andrea Gissi，Orazio Nicolotti，Angelo Carotti，et al.Integration of QSAR models for bioconcentration suitable for REACH ［J］.Science of the Total Environment，2013，456/457:325-332.

［4］Luo Haibin，Chen Guowen，Shao Yongxian，et al.3D QSAR studies on pyrimidine analogues as cyclin-dependent kinase 1 inhibitors［J］.Journal of Shenzhen University Science and Engineering，2012，29(5):438-443.(in Chinese)罗海彬，陈国文，邵咏贤，等.嘧啶类CDK1激酶抑制剂的三维定量构效关系 ［J］.深圳大学学报理工版，2012，29(5):438-443.

［5］Feng Changjun.Theoretical studies on the action strength of DOM of phenyl-isopropylamine dopes［J］.Journal of Shenzhen University Science and Engineering，2012，29(1):67-72.(in Chinese)冯长君.苯异丙基胺类兴奋剂兴奋强度理论研究［J］.深圳大学学报理工版，2012，29(1):67-72.

［6］Du Xihua.QSPR study on thermodynamic properties of polybrominated dibenzofurans and polybrominated dibenzothiophenes［J］.Journal of Chemical Industry and Engineering，2010，61(12):3059-3066.(in Chinese)堵锡华.多溴代二苯并呋喃/噻吩热力学性质的定量构效关系 ［J］.化工学报，2010，61(12):3059-3066.

［7］Cao Hongcui，Wu Qixun，Wang Aili，et al. QSPR models for polychlorinated diphenyl ethers［J］.Journal of Henan Polytechnic University Natural Science，2012，31(6):749-752.(in Chinese)曹红翠，吴启勋，王爱丽，等.多氯联苯醚类(PCDEs)化合物的QSPR研究［J］.河南理工大学学报自然科学版，2012，31(6):749-752.

［8］Lu X X，Tao S，Hu H Y，et al.Estimation of bioconcen tration factors of nonionic organic compounds in fish by molecular connectivity indices and polarity correction factors［J］.Chemosphere，2000，41(10):1675-1688.

［9］Liu Yawei，Niu Fangfang，Zeng Pengju，et al.DFT study on electronic structures and spectra properties of substituent pentacene connected with electron donating and withdrawing groups［J］.Journal of Shenzhen University Science and Engineering，2010，27(3):267-272.(in Chinese)刘亚伟，牛芳芳，曾鹏举，等.取代并五苯电子结构与光谱的密度泛函研究［J］.深圳大学学报理工版，2010，27(3):267-272.

［10］Song Jie，Wang Zongde，Findlater A，et al.Terpenoid mosquito repellents:a combined DFT and QSAR study［J］.Bioorganic＆ Medicinal Chemistry Letters，2013，23(5):1245-1248.

［11］Feng Changjun.Theoretical studies on quantitative structure-activity relationship and structural modification for 3-substituted sulfur-5-(2-hydroxyphenyl)-4H-1，2，4-triazole compounds ［J］.Acta Chimica Sinica，2012，70(4):512-518.(in Chinese)冯长君.3-取代硫基-5-(2-羟基苯基)-4H-1，2，4-三唑类化合物抑菌活性的定量构效关系和结构修饰的理论研究 ［J］.化学学报，2012，70(4):512-518.

［12］Urra L S，Gonza'lez M P，Teijeira M.2D-autocorrelation descriptors for predicting cytotoxicity of naphtha-quinone ester derivatives against oral human epidermoid carcinoma［J］.Bioorganic＆ Medicinal Chemistry，2007，15(10):3565-3571.

［13］Urra L S，Gonza'lez M P，Teijeira M.QSAR studies about cytotoxicity of benzophenazines with dual inhibition toward both topoisomerases I and II:3D-MoRSE descriptors and statistical considerations about variable selection ［J］.Bioorganic＆ Medicinal Chemistry，2006，14(21):7347-7358.

［14］Yao S W，Lopes V H C，Fernández F，et al.Synthesis and QSAR study of the anticancer activity of some novel indane carbocyclic nucleosides［J］.Bioorganic ＆Medicinal Chemistry，2003，11(23):4999-5006.

［15］Tropsha A，Gramatica P，Gombar V K.The importance of being earnest:Validation is the absolute essential for successful application and interpretation of QSPR models［J］.QSAR ＆ Combinatorial Science，2003，22(1):69-77.