QSAR模型預(yù)測(cè)石化廢水中芳香族物質(zhì)對(duì)厭氧菌群的綜合毒性

吳少奇，王黎,*，夏正海，孫義，張春雨，魯逸飛

1. 武漢科技大學(xué)資源與環(huán)境工程學(xué)院，環(huán)境污染綠色控制與修復(fù)技術(shù)研究中心，武漢 430081 2. 中韓(武漢)石油化工有限公司發(fā)展技術(shù)部，武漢 430082

石化、焦化和鋼鐵制造等行業(yè)產(chǎn)生的工業(yè)廢水成分復(fù)雜，且普遍存在著芳香化合物[1]。芳香化合物是一類具有持久性、生物累積性和毒性的有機(jī)污染物[2]，對(duì)自然環(huán)境和人體健康有潛在威脅。含芳香化合物的有機(jī)廢水處理方式主要包括物理化學(xué)法(如超聲波處理、光化學(xué)氧化和電化學(xué)氧化等)[3]和生物法(如好氧降解和厭氧消化等)。其中，厭氧消化很好地克服了物理化學(xué)法和好氧降解所面臨的高成本、污染物降解不充分以及二次污染等問(wèn)題[4]。但厭氧微生物對(duì)環(huán)境較為敏感，高濃度的芳香化合物廢水會(huì)對(duì)其活性產(chǎn)生明顯的抑制作用，從而影響出水水質(zhì)。因此，確定芳香族化合物對(duì)厭氧微生物的抑制閾值很有現(xiàn)實(shí)意義。

通過(guò)傳統(tǒng)的逐一毒性試驗(yàn)獲取有機(jī)物毒性數(shù)據(jù)，不僅工作量大而且費(fèi)時(shí)費(fèi)力，已無(wú)法滿足當(dāng)下芳香化合物廢水處理的實(shí)際需求。定量構(gòu)效關(guān)系(quantitative structure activity relationship, QSAR)模型是一種通過(guò)統(tǒng)計(jì)學(xué)手段建立物質(zhì)結(jié)構(gòu)參數(shù)與有機(jī)體活性定量關(guān)系的分析方法，可實(shí)現(xiàn)對(duì)有機(jī)物性質(zhì)的有效預(yù)測(cè)[5]。同時(shí)，芳香化合物的毒性大小在很大程度上受其結(jié)構(gòu)參數(shù)(如物理化學(xué)參數(shù)和量子化學(xué)參數(shù))的影響[6]。因此，通過(guò)建立芳香化合物結(jié)構(gòu)與其毒性關(guān)系的QSAR模型，可探究影響芳香化合物毒性大小的主要結(jié)構(gòu)參數(shù)，從而合理預(yù)測(cè)芳香化合物的毒性。

為預(yù)測(cè)芳香化合物對(duì)厭氧微生物的毒性作用，筆者進(jìn)行了厭氧微生物毒理試驗(yàn)，以19種芳香化合物的48 h半數(shù)效應(yīng)濃度(48 h-EC50)值為樣本，基于20個(gè)理化和量子化學(xué)結(jié)構(gòu)描述符，利用多元逐步線性回歸(multiple stepwise linear regression, MSLR)法建立了芳香化合物毒性與結(jié)構(gòu)參數(shù)相關(guān)聯(lián)的QSAR模型。

1 材料與方法(Materials and methods)

1.1 實(shí)驗(yàn)材料

1.1.1 模擬廢水

采用模擬廢水，其成分根據(jù)石化廢水處理某工段進(jìn)水水質(zhì)，化學(xué)需氧量(COD)為4 000 mg·L-1，用葡萄糖、NH4Cl和KH2PO4按m(COD)∶m(N)∶m(P)=200∶5∶l的比例配制而成，同時(shí)加入微量元素Zn2+、Cu2+、Mn2+、Co2+和Ni2+。

1.1.2 接種污泥

試驗(yàn)所用的接種污泥來(lái)源于武漢某污水處理廠的厭氧污泥，接種污泥的揮發(fā)性固體(VS)、總固體(TS)和pH值分別為2.29%、4.32%和7.1。采集的厭氧污泥先在1.1.1所述的模擬廢水中馴化20 d，之后用于厭氧微生物毒理試驗(yàn)。

1.1.3 芳香化合物

用于厭氧微生物毒理試驗(yàn)的芳香化合物(分析純)共19種，可分為3類：苯類、苯酚類以及苯胺類。苯類包括苯、鄰二氯苯、氯苯和對(duì)二甲苯；苯酚類包括苯酚、鄰甲基苯酚、間甲基苯酚、對(duì)甲基苯酚、鄰硝基苯酚、對(duì)硝基苯酚、間硝基苯酚、3，5-二甲基苯酚、2，4-二氯苯酚、2,6-二氯苯酚、鄰(對(duì))苯二酚和五氯苯酚；苯胺類包括苯胺和2,6-二乙基苯胺。

1.2 厭氧微生物毒理試驗(yàn)

在250 mL的消化瓶中加入200 mL模擬廢水和接種污泥所構(gòu)成的混合液，使初始的揮發(fā)性懸浮固體(VSS)濃度為15 g·L-1，同時(shí)調(diào)節(jié)pH為7.1±0.1。在消化瓶中分別添加1.1.3中所述的芳香化合物，使瓶中芳香化合物的濃度分別為100、500、1 000、2 000和4 000 mg·L-1。用氮?dú)庠谝好嫔喜看祾? min，然后用膠塞密封瓶口，保持瓶?jī)?nèi)處于厭氧狀態(tài)。將消化瓶置于(35±1) ℃的恒溫震蕩器(ZQTY-70N，上海知楚，中國(guó))中厭氧消化48 h，測(cè)定其最終累積產(chǎn)氣量。每組實(shí)驗(yàn)設(shè)置2個(gè)平行樣。

1.3 分析方法

1.3.1 測(cè)定方法

TS、VS、VSS以及pH按照美國(guó)公共衛(wèi)生協(xié)會(huì)(APHA)的標(biāo)準(zhǔn)方法進(jìn)行測(cè)定[7]；累積產(chǎn)氣量采用排水法進(jìn)行測(cè)定[8]。

1.3.2 數(shù)據(jù)分析方法

通過(guò)式(1)計(jì)算48 h芳香化合物對(duì)厭氧微生物的抑制率，然后將得到的48 h各組不同濃度的抑制率代入Logistic方程進(jìn)行模擬，獲得48 h-EC50值[9]。抑制率的計(jì)算式如下[10]。

(1)

式中：I為芳香化合物對(duì)厭氧微生物活性的抑制率(%)；GA為實(shí)驗(yàn)組t時(shí)刻厭氧微生物累積產(chǎn)氣量(mL)；GB為對(duì)照組t時(shí)刻厭氧微生物累積產(chǎn)氣量(mL)。

1.4 QSAR模型構(gòu)建

1.4.1 描述符

為了鑒定影響芳香化合物對(duì)厭氧微生物毒性的結(jié)構(gòu)參數(shù)，選擇了20個(gè)具有代表性的描述符，涉及理化參數(shù)包括正辛醇/水分配系數(shù)(logKow)、折射率(n)、水溶解度(logS)、生成焓(heat of form, HOF)、摩爾折射率(molar refractivity, MR)[11]，同時(shí)涉及量子化學(xué)參數(shù)，包括最高占據(jù)軌道能(EHOMO)、最低空軌道能(ELUMO)、次最高占據(jù)軌道能(ENHOMO)、次最低空軌道能(ENLUMO)、分子偶極矩(dipole moment,μ)、分子總能量(total energy, TE)、分子摩爾體積(molecular volume, MV)、極化率(polarizability,α)、自旋極化率(spin polarization, SP)、零點(diǎn)(zero point energy, ZPE)、超極化率(hyper polarizability,β)、熱力學(xué)能(thermodynamic energy, ThE)、熱容量(heat capacity,HC)、熵(entropy,En)、自洽場(chǎng)能(self-consistent field energy, SCFE)[12]。理化參數(shù)通過(guò)Chemoffice查詢獲得，量子化學(xué)參數(shù)通過(guò)Gaussian 09W軟件的密度泛函理論雜化泛函(DFT-B3LYP)方法，在6-311G(d, p)水平下對(duì)芳香化合物結(jié)構(gòu)優(yōu)化所得[13]。

1.4.2 QSAR模型構(gòu)建與分析

在上述3類芳香化合物中隨機(jī)各抽取一種組成預(yù)測(cè)集，用于檢測(cè)模型的預(yù)測(cè)能力，其余的芳香化合物作為訓(xùn)練集。為了減少數(shù)據(jù)的跨度，將EC50值轉(zhuǎn)換為10為底的對(duì)數(shù)值-lgEC50。采用MSLR建立-lgEC50與結(jié)構(gòu)參數(shù)的芳香化合物毒性預(yù)測(cè)模型。

(2)

(3)

模型的適用范圍由Williams圖(即標(biāo)準(zhǔn)化殘差與杠桿值)來(lái)表征[15]。計(jì)算標(biāo)準(zhǔn)化殘差(δ)、杠桿值(h)以及警告杠桿(h*)的公式如下。

(4)

(5)

h*=3(p+1)/n

(6)

式中：p為模型中包含的描述符數(shù)量；xi為化合物i的描述符行向量；X為訓(xùn)練集的描述符矩陣。

標(biāo)準(zhǔn)化殘差的絕對(duì)值>3的化合物被視為異常值。若預(yù)測(cè)的芳香化合物殘差在此范圍內(nèi)，則說(shuō)明模型的預(yù)測(cè)是可靠的。采用Origin 2017和IBM SPSS Statistics 25進(jìn)行上述統(tǒng)計(jì)分析。

2 結(jié)果與討論(Results and discussion)

2.1 芳香化合物毒性比較

芳香化合物對(duì)厭氧微生物的毒性大小可由48 h-EC50值表示，上述芳香化合物的48 h-EC50值范圍為486～2 224 mg·L-1(圖1)。EC50值反映了厭氧微生物所能忍受的毒性大小，EC50值越低，則芳香化合物的毒性越高。因此可知，五氯苯酚對(duì)厭氧微生物的毒性作用最強(qiáng)，對(duì)苯二酚則最弱。

圖1 芳香化合物的48 h半數(shù)效應(yīng)濃度(48 h-EC50)值Fig. 1 48 h median effect concentration (48 h-EC50) value of aromatic compounds

從毒理試驗(yàn)結(jié)果來(lái)看，苯類物質(zhì)對(duì)厭氧微生物的毒性高于苯酚類物質(zhì)，這與前人所得結(jié)論相近[16]。一方面，苯類物質(zhì)比苯酚類物質(zhì)具有更強(qiáng)的疏水性，導(dǎo)致前者能更輕易地與微生物細(xì)胞相結(jié)合，對(duì)其產(chǎn)生毒性作用。另一方面，有機(jī)化合物的生物降解過(guò)程可以間接影響其毒性的作用效果。苯酚類物質(zhì)的結(jié)構(gòu)比苯類物質(zhì)更為復(fù)雜，導(dǎo)致其難以被厭氧微生物降解。因此，最先與微生物作用的是苯類物質(zhì)，其毒性作用在反應(yīng)初期就顯現(xiàn)出來(lái)，苯酚類物質(zhì)則與之相反。

2.2 芳香化合物毒性預(yù)測(cè)的QSAR模型

2.2.1 描述符的篩選和QSAR模型的構(gòu)建

為了避免忽略具有低相關(guān)性但高重要性的描述符，將1.4.1中的20個(gè)描述符與訓(xùn)練集中16種芳香化合物進(jìn)行MSLR分析，回歸結(jié)果如表1所示。

表1 多元逐步線性回歸的芳香化合物定量構(gòu)效關(guān)系(QSAR)模型Table 1 The quantitative structure activity relationship (QSAR) model of aromatic compounds for multiple stepwise linear regression

-lgEC50=1.065+0.196logKow+0.154MR+0.028μ(7)

n=16，R2=0.928，F(xiàn)=47.581，P=0.000

2.2.2 模型檢驗(yàn)

對(duì)于內(nèi)部驗(yàn)證，式(1)的R2為0.928，逐步回歸的顯著性P=0.000<0.05，F(xiàn)=47.581>臨界值F0.05(3,14)=3.34，說(shuō)明該模型與訓(xùn)練集擬合優(yōu)度較好，-lgEC50與logKow、MR和μ所構(gòu)成的回歸方程可靠度很高。從表2可知，logKow、MR和μ三者之間的相關(guān)系數(shù)分別為0.745、0.116和0.188，均<0.8[17]，且模型的VIF值介于l～10之間，這說(shuō)明3個(gè)自變量之間沒(méi)有潛在的共線性，故從統(tǒng)計(jì)學(xué)角度所構(gòu)模型可以接受。模型的R2(>0.6)，F(xiàn)(>15)，P(<0.05)和VIF(<10)值均滿足最小統(tǒng)計(jì)穩(wěn)健性標(biāo)準(zhǔn)[18-19]。

2.2.3 外部預(yù)測(cè)能力

圖2 -lgEC50的實(shí)驗(yàn)值與計(jì)算值相關(guān)性圖Fig. 2 Correlation between experimental and calculated values of -lgEC50

此外，利用Williams圖分析了QSAR模型對(duì)芳香化合物毒性預(yù)測(cè)的適用范圍(圖3)。水平虛線表示標(biāo)準(zhǔn)化的剩余異常值(±3)，垂直虛線表示警示杠桿值(h*=0.750)。由圖3可知，訓(xùn)練集和測(cè)試集中的所有數(shù)據(jù)點(diǎn)都屬于適用范圍，表明QSAR模型對(duì)于通過(guò)-lgEC50與結(jié)構(gòu)參數(shù)的線性關(guān)系預(yù)測(cè)芳香化合物對(duì)厭氧微生物的毒性是可靠的。

圖3 QSAR模型的Williams圖Fig. 3 Williams plot of the QSAR model

表2 回歸方程中各自變量系數(shù)相關(guān)矩陣及其方差膨脹因子(VIF)Table 2 Correlation matrix of each variable coefficient and variance inflation factor (VIF) in the regression equation

表3 模型參數(shù)以及模型驗(yàn)證結(jié)果Table 3 Model parameters and model verification results

2.3 芳香化合物的毒性與描述符的作用機(jī)理

QSAR模型表示為式(7)，表明芳香化合物對(duì)厭氧微生物的毒性影響與logKow、MR和μ有關(guān)，且-lgEC50均與三者成正比。從表2可知，logKow、MR和μ對(duì)此模型的貢獻(xiàn)程度分別為51.30%、27.63%和21.07%，說(shuō)明logKow可能是影響芳香化合物對(duì)厭氧微生物毒性大小的主要因素。

由圖4(a)可知，logKow與-lgEC50值相關(guān)性顯著(r2=0.8105)。芳香化合物的毒性隨logKow增大而增大。logKow反映了芳香化合物的溶脂性大小，而微生物的細(xì)胞膜是由磷脂雙分子層所構(gòu)成。芳香化合物的logKow越大，其溶脂能力越強(qiáng)，可能越容易通過(guò)細(xì)胞膜，從而表現(xiàn)出的毒性作用也就越大。

圖4 模型參數(shù)與-lgEC50的相關(guān)性圖注：(a) logKow與-lgEC50的相關(guān)性；(b) MR與-lgEC50的相關(guān)性；(c) μ與-lgEC50的相關(guān)性。Fig. 4 Correlation between model parameters and -lgEC50Note: (a) the correlation between logKow and -lgEC50; (b) the correlation between MR and -lgEC50; (c) the correlation between μ and -lgEC50.

MR是由于在光的照射下分子中電子云相對(duì)于分子骨架相對(duì)運(yùn)動(dòng)的結(jié)果，可作為分子中電子極化率的度量[21]。由圖4(b)可知，MR與-lgEC50值相關(guān)性顯著(r2=0.7128)，芳香化合物對(duì)厭氧生物的毒性隨MR增大而增大。

μ體現(xiàn)了芳香化合物影響微生物細(xì)胞極性的能力。由圖4(c)可知，μ與-lgEC50值相關(guān)性不顯著(r2=0.0695)。Su等[19]探究了四胺基大環(huán)配體(TAML)/H2O2)氧化芳香族有機(jī)污染物的定量構(gòu)效關(guān)系，發(fā)現(xiàn)所構(gòu)多元線性模型的自變量#H:C與因變量lnkobs的相關(guān)性并不顯著。因此，應(yīng)將μ與其他因素結(jié)合起來(lái)考慮。

芳香化合物的濃度會(huì)顯著影響厭氧消化過(guò)程中厭氧微生物活性，從而引起出水水質(zhì)的波動(dòng)。本研究通過(guò)厭氧微生物毒理試驗(yàn)，以芳香化合物48 h-EC50值為樣本，基于理化和量子化學(xué)結(jié)構(gòu)描述符，采用MSLR建立了QSAR模型，得出如下結(jié)論。

(1)QSAR模型為：-lgEC50=1.065+0.196logKow+0.154MR+0.028μ，R2=0.928，此模型可用于芳香化合物的毒性大小預(yù)測(cè)(h*=0.750)，有助于提高石化類工業(yè)廢水的厭氧生物處理效率。

(2)芳香化合物對(duì)厭氧微生物的毒性影響與logKow、MR和μ有關(guān)，三者對(duì)模型的貢獻(xiàn)程度分別為51.30%、27.63%和21.07%，logKow可能是影響芳香化合物對(duì)厭氧微生物毒性大小的主要因素。

(3)-lgEC50值與logKow和MR相關(guān)性顯著，r2分別為0.8105和0.7128；-lgEC50值與μ相關(guān)性不顯著，需要與其他因素綜合起來(lái)考慮。