鹵代苯在魚體內(nèi)生物富集因子的理論模型

馮長君，堵錫華

徐州工程學院化學化工學院，江蘇徐州221111

生物富集一般指生物體從周圍環(huán)境中吸收不降解或難降解的化學物質(zhì)，在體內(nèi)積累，使這些物質(zhì)在生物體內(nèi)的濃度超過其生長環(huán)境中濃度的現(xiàn)象.生物富集程度常用生物富集因子(bioconcentration factors，BCF)表示，即生物體內(nèi)污染物的平衡濃度(cS)與其生存環(huán)境中該污染物濃度(cH)的比值

目前，對于非離子性有機物在魚體內(nèi)生物富集實驗數(shù)據(jù)的處理多采用“二室模型”［1］，即將實驗測定的模擬池分為水室和魚室，假定有機物在池內(nèi)的遷移轉(zhuǎn)化為一級動力學過程，可推得類似公式

其中，cw和cf分別為平衡時有機物在水體及魚體內(nèi)濃度(單位:mol/L).

污染物在不同生態(tài)系統(tǒng)中的含量，隨食物鏈營養(yǎng)級的升高而增加，其富集系數(shù)在各營養(yǎng)級中均可達極其驚人的數(shù)值.處于食物鏈頂端的人類，便成為生物富集的最終受害者.如生物富集已引起包括水俁病和痛痛病在內(nèi)的多起生態(tài)公害事件.

鹵代苯(halogenated benzenes，HB)屬于環(huán)境外來品，是化工、醫(yī)藥、制革和電子等行業(yè)廣泛應用的化工原料、有機合成中間體或有機溶劑難降解而在環(huán)境中高度滯留.此類物質(zhì)通過食物鏈進入人體，具有強烈的致癌、致畸和致突變作用，對人體造成嚴重危害.因此，其中的氯苯、鄰二氯苯、對二氯苯和六氯苯已被列入美國及歐洲共同體環(huán)境保護部門所確定的129種優(yōu)先控制污染物名單.通過研究生物富集作用，利于闡明物質(zhì)在生態(tài)系統(tǒng)內(nèi)的遷移和轉(zhuǎn)化規(guī)律，評價和預測污染物進入生物體后可能造成的危害，在利用生物體對環(huán)境進行監(jiān)測和凈化及制定排放標準等方面，也具有重要的意義.

由于實測BCF成本高、周期長，在實際工作中通常需要采用估算方法來獲取所需數(shù)據(jù).已進行較多研究的估算法主要包括BCF與正辛醇-水分配系數(shù)(Kow)、水溶解度(SW)和吸附系數(shù)(Koc)的各種經(jīng)驗關(guān)系式及分子拓撲法等［2-3］.雖然 Kow、SW和Koc等關(guān)系式精度較高，但不足之處在于參數(shù)本身是實驗測定值，不便應用.分子拓撲法建立BCF模型雖不依賴于實驗，但所用的拓撲參數(shù)無明確物理意義，不便探討影響生物富集的結(jié)構(gòu)因素及富集機理.

物質(zhì)定量構(gòu)效關(guān)系(quantitative structure-activity relationship，QSAR)研究［4-6］不僅是數(shù)據(jù)建模，顯示變量間的數(shù)學關(guān)系;更重要的是模型的有效性和實用性.有效性即模型質(zhì)量，可用一系列統(tǒng)計回歸指標來衡量;實用性是在有效性基礎(chǔ)上，模型所蘊含的物理化學意義.因此，選用具有明確物理化學意義的描述符對于建立良好的QSAR模型至關(guān)重要.目前研究多鹵聯(lián)苯、多鹵聯(lián)苯醚(包含鹵代烴和鹵代苯)生物富集因子構(gòu)效關(guān)系的較多［7-8］，單獨探討鹵代苯BCF的QSAR模型尚未見報道.因此，本研究利用密度泛函理論(density functional theory，DFT)［9-11］對鹵代苯分子進行量化計算，獲得不需要實驗測定，無實驗誤差且具有明確物理意義的量子化學參數(shù)(SL).然后基于最佳變量子集回歸方法建立上述量子化學參數(shù)與其lg BCF［8］的最佳數(shù)學模型，用其探尋影響鹵代苯生物富集行為的結(jié)構(gòu)因素，揭示生物富集微觀機理，提供理論詮釋參考.

1 理論方法

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

從文獻［8］中選取21種鹵代苯，其名稱和生物富集因子數(shù)據(jù)見表1.測試物種的魚類有虹鱒魚、孔雀魚、黑頭呆魚和藍鰓太陽魚等.因不同鹵代苯間BCF值相差比較大，且據(jù)熱力學原理，采用其對數(shù)值lg BCF表示更為科學.lg BCF數(shù)值越大，表示該物質(zhì)越易被富集.

1.2 量化參數(shù)計算

運用DFT的Becke三參數(shù)混合泛函(B3LYP)方法，在較高基組6-311G＊＊(d，p)水平上對鹵代苯類化合物分子結(jié)構(gòu)進行幾何優(yōu)化和理論計算.對所選化合物分子進行結(jié)構(gòu)優(yōu)化及頻率計算，驗證皆無虛頻，表明其構(gòu)型存在真實極小點;然后對以上優(yōu)勢構(gòu)型進行單點能計算，從中提取一些相關(guān)的量化參數(shù)及熱力學數(shù)據(jù)，均以SL表示.

1.3 變量篩選與模型驗證

通過DFT方法獲得的鹵代苯類化合物分子的量化描述符SL之間一定存在或多或少的信息重疊.為遵循奧卡姆剃刀原則，盡可能建立簡單的QSAR模型，須對其進行篩選，以建立具有良好穩(wěn)定性和較高預測能力的QSAR模型.目前對變量壓縮與篩選的方法較多，如逐步回歸、最佳變量子集回歸和遺傳算法等.本研究基于最佳變量子集回歸方法對變量SL進行壓縮建立QSAR模型.模型質(zhì)量采用留一法交叉驗證(leave-one-out cross validation，LOOCV)和留多法交叉驗證(leave-multiple-out cross validation，LMOCV)的相關(guān)系數(shù)(Rcv-o2和Rcv-m2)作為模型質(zhì)量的判斷標準.為確立最終的QSAR模型，引入Akaike信息判據(jù)AC(Akaike's information criterion)和Kubinyi函數(shù) KT(Kubinyi function)［12-13］，其計算公式為

其中，RSS為方差和;f為化合物數(shù);b為變量數(shù).AC值越小，KT值越大，所建的模型越穩(wěn)定，預測能力越高.

2 結(jié)果討論

2.1 lg BCF的QSAR模型

采用DFT方法，在B3LYP/6-311G＊＊(d，p)基組上對21種鹵代苯分子進行結(jié)構(gòu)優(yōu)化及頻率計算，從中提取一些用于定量描述鹵代苯-受體分子間作用力的量化及熱力學參數(shù):①能量參數(shù)有最高占有軌道能(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;②電子性質(zhì)包括分子偶極距μ(單位:Debye)、分子的可極化率α(單位:a.u.)、苯環(huán)上各碳原子的凈電荷(指原子核帶的正電荷與聚集在該原子上的電子所帶的負電荷之和)Q1～Q6及其碳原子凈電荷之和Qsum;③熱力學性質(zhì)有分子的零點振動能(Ev)、熱力學能(E)、熱容(Cv)、熵(Sm)、電子運動空間(VE)及分子配分函數(shù)Q(取Q的常用對數(shù)為lg Q)等.

基于最佳變量子集回歸方法對21個鹵代苯的生物富集因子與SL進行變量壓縮建立QSAR模型，結(jié)果見表2，其中，R2、Radj2、Rcv-o2、S、F、AC和KT分別為判定系數(shù)(亦稱削減誤差比例)、校正判定系數(shù)(以消除自變量個數(shù)及樣本容量對判定系數(shù)的影響)、逐一剔除法的交叉驗證系數(shù)、估計標準誤差、Fisher統(tǒng)計值、Akaike信息判據(jù)和Kubinyi函數(shù).表2為建立高穩(wěn)定性、高預測能力的QSAR模型提供指導.

表2 lg BCF與SL的最佳變量子集回歸結(jié)果Table2 The results of SLand lg BCFwith Leaps-and-Bounds regression

由表2可見，隨著模型變量增多，對應R2和Radj2持續(xù)增大，S持續(xù)減小，Rcv-o2、F和KT均呈先增后減趨勢，而AC變化趨勢則完全反之.但其出現(xiàn)轉(zhuǎn)折點的位置不一，對應最大Rcv-o2的模型為

從式(5)可見，描述符ENLUMO、lg Q和Q2的回歸系數(shù)均是其標準偏差的2倍以上，而VE的回歸系數(shù)不足其標準偏差的2倍，說明VE的引入對模型的穩(wěn)定是不利的.在置信度為95%時，標準t值(tα/2)為 2.080.模型 lg BCF中 ENLUMO、lg Q、Q2和VE的t值依次為 -11.356、 -7.711、2.495 和1.681.可見，VE的絕對值小于標準t值，進一步證明該方程是不穩(wěn)健的.扣除VE后的三元數(shù)學模型

從式(6)可見，描述符ENLUMO、lg Q和Q2的回歸系數(shù)分別是其標準偏差的10.6、8.0和2.6倍，這說明它們的引入對模型的穩(wěn)定性是有利的;結(jié)合AC與KT均在此處發(fā)生轉(zhuǎn)折，可認為該模型具有良好的穩(wěn)定性與預測能力.對比式(5)和(6)，可以發(fā)現(xiàn)其常數(shù)項和描述符ENLUMO、lg Q和Q2的回歸系數(shù)變化均很小，表明VE因素對lg BCF的貢獻很小，可以舍棄.因此采用3個變量建模既保證良好的穩(wěn)定性與預測能力，又符合奧卡姆剃刀原則，即盡可能選擇簡單的模型.

2.2 模型預測能力和穩(wěn)定性的再驗證

在置信度為95%時，標準t值(tα/2)為2.0796.模型 lg BCF中 ENLUMO、lg Q和Q2的t值依次為-10.674、 -7.931 和2.694，它們的絕對值均大于tα/2，證明方程(6)是穩(wěn)健的.按此模型給出的估算值見表1，與相應實驗值吻合.若要建立一個具有良好預測能力的QSAR模型，其R2≥0.9［14］.模型(6)的R2=0.922 ＞0.9，且 Rcv-o2=0.866 ＞ 0.5［15］，可認為該模型具有良好的穩(wěn)健性與預測能力.

將模型(6)的樣本隨機分為2個樣本:第4、8、12、16和18個化合物為檢驗集，余下為訓練集.對于訓練集的模型為此模型的R2、和S與模型(6)很接近.按照模型(7)給出訓練集的估算值和檢驗集的預測值，均與相應實驗值很接近，如圖1.

圖1 鹵代苯化合物生物富集因子(lg BCF)的實驗值與計算值的關(guān)系Fig.1 Relationship between experimental and calculated bioconcentration factors(lg BCF)of halogenated benzenes

為進一步驗證所建QSAR模型的預測能力，把21個鹵代苯分成訓練集(16個化合物)和測試集(5個化合物)，共形成20 349個訓練集與測試集(模型(7)只是其中之一).通過Matlab編程計算這些訓練集、測試集的相關(guān)系數(shù)及逐五剔除法(leave-five-out，LFO;即每次剔除5個鹵代苯分子，用余下16個分子建模)的交叉驗證相關(guān)系數(shù)(Rcv-f2).對于20 349個測試集的Rcv-f2=0.855，與=0.866很接近;相應20 349個訓練集的R2的平均值為0.925，與模型(6)的0.922非常接近.

為定量反映模型的穩(wěn)健性，作者曾定義相關(guān)系數(shù)(R或R2)的相對標準偏差(RS)予以衡量

其中，Ri2為模型i的相關(guān)系數(shù);Ra2為相關(guān)系數(shù)的平均值;w為自由度數(shù)，等于相關(guān)系數(shù)的個數(shù)減1.顯然，RS越小，該模型的穩(wěn)健性越高.根據(jù)相關(guān)回歸分析中的顯著性檢驗，通常確定顯著性水平α=0.05.為此，建議RS＜0.05，說明Ri2變異性很小，呈現(xiàn)明顯的集中趨勢，反映所建模型具有良好的穩(wěn)健性.20 349個訓練集的Ri2的RS=0.022 ＜0.05，表明模型(6)是穩(wěn)健的，不存在異常值.

2.3 生物富集機理分析

根據(jù)進入模型(6)中3個變量ENLUMO、lg Q和Q2的標準化回歸系數(shù)(standardized regression coefficients，SR)依次為 -2.531、 -1.904 和0.203;從SR的絕對值可見，對lg BCF影響順序是ENLUMO＞lg Q ＞ Q2.lg BCF分別與ENLUMO、lg Q、Q2回歸，它們的判定系數(shù) R2依次為 0.624、0.372 和0.172，此與SR結(jié)果一致.據(jù)此可知:

1)根據(jù)分子軌道理論，前線軌道(HOMO和LUMO)及其附近的分子軌道對分子的化學性質(zhì)及生物活性影響最大.其中，LUMO及其附近的未占據(jù)軌道(LUMO+x)等具有優(yōu)先接受電子的作用，能級越低，越易接受電子.模型(6)中ENLUMO越小，越易接受魚體中生物大分子的電子，發(fā)生親電反應而被富集，lg BCF越大.即ENLUMO與lg BCF呈負相關(guān)，故ENLUMO前系數(shù)小于0.

2)配分函數(shù)lg Q反映分子大小和對稱性，分子越小、對稱性越大，lg Q越小.鹵代苯分子越大，越難擴散透過魚鱗表面的磷脂膜，即穿過魚鱗膜進入血液后溶解在肝臟內(nèi)的生物富集量較小.可見，lg Q與lg BCF亦呈負相關(guān).生物富集機理中可假定包含有機物分子鑲嵌在生物受體提供的“空穴”內(nèi)的過程，顯然分子的對稱性與空穴的“契合”程度越高，該有機物的lg BCF值越大.對于同分異構(gòu)體的lg BCF值不同，便是此“假定”的佐證.例如1，2，3-三氯苯、1，2，4-三氯苯和1，3，5-三氯苯的分子對稱性不同，與“空穴”的吻合程度不同，相應 lg BCF亦不同，依次為2.90、2.95 和3.26.

3)形成細胞膜的磷脂屬于兩親性分子，親水性頭部(磷酸基)帶負電荷指向細胞膜外，疏水性尾部(脂肪酸)帶正電荷指向細胞膜內(nèi).Q2是苯環(huán)第2個碳原子上的凈電荷量，其值越高即正電荷量越大，與帶負電荷的親水性磷酸基的靜電引力越強，對魚鱗膜的親和力越強，越易進入魚體內(nèi)而生物富集量越大.Q2值越低即帶負電荷量越大，與帶負電荷的親水性磷酸基的靜電斥力越大，且與魚鱗膜外圍水分子形成氫鍵能力越強，越難穿過魚鱗膜進入魚體，相應lg BCF越小.即Q2與lg BCF成正相關(guān)，此與模型(6)中Q2前系數(shù)大于零是一致的.由此可見，C2是鹵代苯分子中的活性部位，因其受取代基影響最為敏感.

綜上研究，鹵代苯分子在魚體內(nèi)生物富集的微觀歷程可表述為:溶解于水體的鹵代苯，首先要穿過魚鱗膜，即從水室進入魚室，顯然其分子體積及所帶負電荷越大，不利于透過此生物膜;穿過魚鱗膜隨血液進入肝臟的鹵代苯分子，其分子對稱性與生物受體提供的“空穴”的“契合”程度越高，以及ENLUMO越低，越易發(fā)生親電反應而使lg BCF越大.

2.4 熱力學分析

根據(jù)熱力學原理，任何物理過程、化學反應還是生物過程的推動力都源于焓作用和熵作用.對于生物反應中的焓作用包含靜電作用、立體作用(如分子大小和對稱性等);熵作用有疏水作用、轉(zhuǎn)動熵、平動熵以及構(gòu)象熵等.靜電作用又分為離子-離子作用、離子-偶極作用、偶極-偶極作用、氫鍵以及電荷轉(zhuǎn)移作用等，其中蘊含軌道反應及電子得失.由機理討論可見，模型(6)和模型(7)涉及過程的焓作用，表面上似乎沒有熵作用.

鹵代苯分子(HB)在魚體內(nèi)生物受體(biological receptors，BRs)的生物富集機理可表示為

該過程的第1步是鹵代苯在水體中的溶解過程，此是物理過程;第2步是鹵代苯在水體中與在魚體中達平衡的富集過程，其平衡常數(shù)(Kbio)表達式為

其中，［BRs］為純生物受體濃度，取值為1 mol/L.定溫下，此生物反應的吉布斯自由能變(ΔGbio)

其中，ΔHbio為富集過程焓變;ΔSbio為富集過程熵變;R為氣體常數(shù);T為反應溫度.比較模型(6)、(7)與式(13)可見，其ENLUMO、lg Q和Q2主要對應ΔHbio;而常數(shù)項則呈現(xiàn)ΔSbio作用，即疏水性作用及分子平動、轉(zhuǎn)動和構(gòu)象的熵變.由于常數(shù)項是負數(shù)，即ΔSbio＜0，此是熵減過程.即各自獨立運動的水合鹵代苯分子、生物受體大分子相互結(jié)合成HB-BRs而損失的熵(ΔSbio'＜0)，大于體系中增加無序水分子而獲得的熵(ΔSw＞0):ΔSbio=ΔSbio'+ΔSw＜0.實際上，此ΔSw在鹵代苯溶解過程中是被抵消的.因此，整個生物富集過程是熵減過程:ΔSbio=ΔSbio'＜0.

綜上所述，鹵代苯類化合物在魚體內(nèi)的生物富集過程中是焓變起到主要及正向作用，熵變(主要是疏水性作用)則起次要及負向作用，這也表明其生物富集是有序性提高的熵減過程.

結(jié) 語

基于DFT方法，在B3LYP/6-311G＊＊(d，p)基組上對21種鹵代苯分子進行結(jié)構(gòu)優(yōu)化及頻率計算，從中提取了量化及熱力學參數(shù)SL.基于最佳變量子集回歸方法，建立了它們在魚體內(nèi)lg BCF的QSAR模型，經(jīng)R2、Rcv-o2、Rcv-f2、F、AC和FT等質(zhì)量指標檢驗，顯示良好的相關(guān)性、穩(wěn)健性及預測能力.

根據(jù)進入模型中的量化參數(shù)分析，鹵代苯分子在魚體內(nèi)生物富集的微觀機理可能是:鹵代苯分子體積及所帶負電荷越大，不利于穿過魚鱗膜;進入肝臟的鹵代苯分子與生物受體內(nèi)“空穴”的“契合”程度越高，以及ENLUMO越低，越易發(fā)生親電反應而使lg BCF越大.

經(jīng)熱力學分析，鹵代苯類化合物在魚體內(nèi)的生物富集過程中的焓變起到主要及正向作用，熵變(主要是疏水性作用)則起次要及負向作用，這也表明其生物富集是有序性提高的熵減過程.

/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)侯 玲，胡長敏，趙曉明，等.江河中有機污染物在鯉魚體內(nèi)富集與釋放的模擬研究 ［J］.環(huán)境科學，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激酶抑制劑的三維定量構(gòu)效關(guān)系 ［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)堵錫華.多溴代二苯并呋喃/噻吩熱力學性質(zhì)的定量構(gòu)效關(guān)系 ［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)曹紅翠，吳啟勛，王愛麗，等.多氯聯(lián)苯醚類(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)劉亞偉，牛芳芳，曾鵬舉，等.取代并五苯電子結(jié)構(gòu)與光譜的密度泛函研究［J］.深圳大學學報理工版，2010，27(3):267-272.

［10］Song Jie，Wang Zongde，F(xiàn)indlater 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-三唑類化合物抑菌活性的定量構(gòu)效關(guān)系和結(jié)構(gòu)修飾的理論研究 ［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，F(xiàn)erná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.