有機(jī)物對月牙藻的定量結(jié)構(gòu)-毒性關(guān)系模型研究

鄭霄晗，賈青竹，閆方友

(1.天津科技大學(xué)海洋與環(huán)境學(xué)院，天津 300457；2.天津科技大學(xué)化工與材料學(xué)院，天津 300457)

近幾十年來，工業(yè)、農(nóng)業(yè)以及個人護(hù)理領(lǐng)域使用了大量有機(jī)化學(xué)物質(zhì)[1-3].這些物質(zhì)經(jīng)過消費(fèi)、使用，通過地表徑流以及廢水排放進(jìn)入河流湖泊、淺層地下水、近海海域等水體環(huán)境中[4-5]，對生物多樣性和生態(tài)系統(tǒng)造成損害.

藻類是能量的初級生產(chǎn)者，水體中毒性物質(zhì)先被藻類吸收，然后通過食物鏈富集，最終危害人類和整個生態(tài)系統(tǒng)[6-7].月牙藻(Pseudokirchneriella subcapitata)作為單細(xì)胞藻類廣泛分布在淡水水域中，因其易培養(yǎng)、生命周期短、接觸表面積大、生理反應(yīng)快、對大量有機(jī)毒物敏感[8-10]，被經(jīng)濟(jì)合作與發(fā)展組織(OECD)推薦為用于生態(tài)毒理學(xué)生物檢測的物種之一.

定量構(gòu)效關(guān)系(QSAR)已成為當(dāng)前在化學(xué)和生物化學(xué)中具有重要理論和實(shí)踐意義的方法[11].當(dāng)前急劇增加的污染壓力要求進(jìn)行快速準(zhǔn)確的危害和風(fēng)險評估，避免有機(jī)物危及水體營養(yǎng)網(wǎng)絡(luò).已有研究成功利用 QSAR模型預(yù)測少量有機(jī)化合物對月牙藻的毒性，常用 10%抑制的有效濃度(EC10)和半數(shù)效應(yīng)濃度(EC50)表示毒性.Lee等[12]用苯甲酸及其衍生物對月牙藻毒性(EC50)建模，模型統(tǒng)計分析結(jié)果表明苯甲酸毒性模型 R2為 0.921，苯甲酸衍生物模型 R2為0.965.在丁腈化合物對月牙藻的毒性研究中，Huang等[13]利用正辛醇/水分配系數(shù)的對數(shù)(logKow)和最低非占據(jù)分子軌道能量作為描述符，建立基于溶解氧的響應(yīng)終點(diǎn)和藻類生長速率的 QSAR模型，兩個模型均表現(xiàn)出良好的穩(wěn)定性和預(yù)測能力(前者R2為0.92，Q2為 0.81；后者 R2為 0.92，Q2為 0.51).Khan 等[14]和 Yu[15]為了解決有機(jī)污染物樣本量小、訓(xùn)練集結(jié)構(gòu)相似的問題，開發(fā)了結(jié)構(gòu)復(fù)雜多樣的有機(jī)污染物模型.Khan等[14]采用偏最小二乘回歸技術(shù)建立最終的QSAR模型，并對月牙藻生態(tài)毒性模型的應(yīng)用性進(jìn)行驗(yàn)證；結(jié)果顯示模型 R2(pEC10)為 0.70，Q2(pEC10)為0.68 和 R2(pEC50)為 0.72，Q2(pEC50)為 0.70.Yu[15]利用支持向量機(jī)(SVM)算法對 pEC10和 pEC50構(gòu)建了兩個定量結(jié)構(gòu)-毒性關(guān)系(quantitative structuretoxicity relationship，QSTR)模型，在測試集中使用了更多的樣本(訓(xùn)練集樣本數(shù)為 167，測試集樣本數(shù)為167)，pEC10的訓(xùn)練集 R2為 0.76，測試集 R2為 0.75，pEC50的訓(xùn)練集R2為0.75，測試集R2為0.74，統(tǒng)計學(xué)參數(shù)顯示預(yù)測結(jié)果令人滿意.

Lee 等[12]、Huang 等[13]、Khan 等[14]和 Yu[15]將實(shí)驗(yàn)值參數(shù) logKow直接或間接作為描述符，成功預(yù)測月牙藻生態(tài)毒性值，并證明 logKow有助于提高這些模型的質(zhì)量，然而使用計算的 logKow[16]數(shù)據(jù)，不能完全避免誤差放大的可能性.

近年來，本課題組根據(jù)分子結(jié)構(gòu)定義一系列范數(shù)描述符，并據(jù)此建立定量構(gòu)效關(guān)系模型、定量結(jié)構(gòu)性質(zhì)關(guān)系模型，預(yù)測有機(jī)物的物理化學(xué)性質(zhì)和有機(jī)物及離子液體的毒性.先前的工作已經(jīng)成功預(yù)測離子液體的生態(tài)毒性[17-18]、殺蟲劑對虹鱒魚的急性毒性[19]、有機(jī)化合物對斑馬魚胚胎的急性毒性[20]、農(nóng)藥對大型水蚤的毒性[21]等.這說明本課題組提出的范數(shù)描述符能夠準(zhǔn)確描述有機(jī)物結(jié)構(gòu)與毒性之間的關(guān)系.

本工作的目的是建立廣義的、應(yīng)用于多種有機(jī)物的月牙藻pEC10和pEC50毒性模型，并通過留一交叉驗(yàn)證(LOO-CV)、外部驗(yàn)證、Y隨機(jī)驗(yàn)證和應(yīng)用性域分析，對模型的魯棒性和預(yù)測能力進(jìn)行評價.

1 數(shù)據(jù)集與研究方法

1.1 數(shù)據(jù)集

月牙藻生態(tài)毒性數(shù)據(jù)(EC10和 EC50)來自 Kusk等[22]的生長抑制實(shí)驗(yàn).使用 298種有機(jī)物建立月牙藻毒性定量構(gòu)效關(guān)系模型，在數(shù)據(jù)收集過程中利用Chemicalbook數(shù)據(jù)庫和 NIST數(shù)據(jù)庫，根據(jù) CAS號對涉及的物質(zhì)結(jié)構(gòu)進(jìn)行數(shù)據(jù)核對，確保建立模型依據(jù)的分子結(jié)構(gòu)準(zhǔn)確無誤.其中，訓(xùn)練集和測試集約按1∶1的比例隨機(jī)劃分，分別包含 148和 150種物質(zhì).按照 QSAR分析的慣例，將實(shí)驗(yàn)生態(tài)毒性(EC10和 EC50)進(jìn)行單位轉(zhuǎn)換，即將以“mg/L”為單位的數(shù)值轉(zhuǎn)換為以“mmol/L”為單位的數(shù)值，并對數(shù)值取負(fù)對數(shù)，標(biāo)記為pEC10或pEC50，298種有機(jī)物的生態(tài)毒性數(shù)據(jù)及物質(zhì)數(shù)字識別號碼列在附表 S1(本文所有附屬文件均已上傳至天津科技大學(xué)學(xué)報網(wǎng)站，網(wǎng)址為http://xuebao.tust.edu.cn)中.

1.2 分子結(jié)構(gòu)優(yōu)化

分子幾何結(jié)構(gòu)優(yōu)化工作在軟件Gaussian 16中完成，用 B3LYP雜化泛函計算了密度泛函理論(DFT)能級上的能量.由于電荷密度、軌道能級對分子性質(zhì)有較大影響，因此使用彌散和極化函數(shù)增廣分裂價6-311+G(d，p)基組以及自然布居分析(NPA)等對電荷影響進(jìn)行分析.

1.3 原子分布矩陣和范數(shù)描述符

根據(jù)分子中的原子間連接關(guān)系，得到如式(1)—式(5)所示的步長矩陣，包括步長矩陣S、相鄰步長矩陣 SA、相間步長矩陣 SB、相跳步長矩陣 SC、相鄰-相間-相跳步長矩陣 SABC.由各個原子在分子中的空間位置推導(dǎo)得出如式(6)—式(11)所示的距離矩陣，包括距離矩陣 D、相鄰距離矩陣 DA、相間距離矩陣DB、相跳距離矩陣 DC、相鄰-相間距離矩陣 DAB、相鄰-相間-相跳距離矩陣DABC.

式(1)—式(11)中：sij為原子 i和 j之間的步長，dij為原子i和j之間的歐幾里得空間距離.

同時，為了進(jìn)一步區(qū)分每個原子還引入了原子性質(zhì)，如電子層數(shù)、最外層電子數(shù)、電離能、NPA里德堡電子數(shù)和NPA價電子數(shù)等，性質(zhì)矩陣列于表1.

表1 性質(zhì)矩陣Tab.1 Property matrices

步長矩陣或距離矩陣與性質(zhì)矩陣通過如式(12)—式(19)所示的8種方式進(jìn)行組合，建立原子分布矩陣(M).對原子分布矩陣進(jìn)行范數(shù)計算，本工作使用如式(20)—式(23)所示的 4種范數(shù)進(jìn)行數(shù)據(jù)挖掘.

1.4 模型驗(yàn)證

針對本工作建立的有機(jī)物對月牙藻定量結(jié)構(gòu)-毒性關(guān)系模型，選擇驗(yàn)證指標(biāo)如決定系數(shù)(R2)、基于留一交叉驗(yàn)證的、平均絕對誤差(eMA)等評價模型的穩(wěn)定性、魯棒性、適應(yīng)度和預(yù)測性.采用 Y隨機(jī)驗(yàn)證檢驗(yàn)?zāi)Ｐ偷臋C(jī)會相關(guān)性.模型的應(yīng)用域可用于評價模型的應(yīng)用范圍.

2 結(jié)果和討論

2.1 模型建立

采用多元線性回歸分析法，提出有機(jī)物對月牙藻生態(tài)毒性的模型.為了避免模型過擬合，對建模所用描述符進(jìn)行優(yōu)化篩選，最后各選取了 18個范數(shù)描述符，建立月牙藻的毒性 pEC10和 pEC50模型如式(24)、式(25)所示，模型中系數(shù) K和范數(shù)描述符I表達(dá)式見附表S2、S3.模型的統(tǒng)計學(xué)參數(shù)見表2.

表2 pEC10和pEC50模型的統(tǒng)計學(xué)參數(shù)Tab.2 Statistical parameters of pEC10 and pEC50 models

2.2 內(nèi)部驗(yàn)證

有機(jī)物對月牙藻生態(tài)毒性數(shù)據(jù) pEC10和 pEC50模型的實(shí)驗(yàn)值和計算值的散點(diǎn)圖如圖1所示.由圖1可知：各點(diǎn)沿對角線分布，表明 pEC10和 pEC50模型的實(shí)驗(yàn)值和計算值基本一致，說明模型擬合良好.本模型與LOO-CV法模型的誤差分布如圖2所示.

圖1 pEC10和pEC50模型的實(shí)驗(yàn)值與計算值散點(diǎn)圖Fig.1 Scatter plot of experimental and calculated values of pEC10 and pEC50 models

圖2 本模型與LOO-CV法模型的誤差分布Fig.2 Error distribution of this model and the LOO-CV method model

訓(xùn)練集 pEC10和 pEC50的分別為0.7552和0.7461，表明本模型與 LOO-CV法的誤差分布是一致的；pEC10和 pEC50的較高，表明該模型穩(wěn)定、魯棒性強(qiáng).

2.3 外部驗(yàn)證

pEC10和pEC50模型中訓(xùn)練集與測試集的散點(diǎn)圖如圖3所示，圖中計算值與實(shí)驗(yàn)值距對角線越近，說明數(shù)據(jù)點(diǎn)擬合結(jié)果越好.pEC10模型的和分別為0.8068和0.8115，pEC50模型的和分別為0.8106和0.7954.較高的決定系數(shù)表明該模型具有較好的預(yù)測能力，能夠?qū)Σ煌Y(jié)構(gòu)的有機(jī)化合物進(jìn)行毒性預(yù)測.

圖3 pEC10和pEC50模型的訓(xùn)練集與測試集的散點(diǎn)圖Fig.3 Scatter plot of the training and test sets of the pEC10 and pEC50 models

2.4 Y隨機(jī)驗(yàn)證

為進(jìn)一步檢驗(yàn)?zāi)Ｐ筒淮嬖谂既恍院瓦^度擬合的可能，對 pEC10和 pEC50模型進(jìn)行 Y隨機(jī)驗(yàn)證.在本工作中，Y隨機(jī)驗(yàn)證重復(fù) 10000次，結(jié)果如圖4所示.pEC10模型的＝0.0067，＝0.0150；pEC50模型的為 0.0070，為 0.0154，說明 pEC10和pEC50模型是穩(wěn)定的，不存在偶然相關(guān).

圖4 10000次Y隨機(jī)驗(yàn)證結(jié)果Fig.4 Results of the 10000 times Y-randomization test

2.5 應(yīng)用域分析

根據(jù) OECD的原則，QSAR模型應(yīng)具有明確的應(yīng)用范圍，利用所建立的 QSTR模型，可以預(yù)測內(nèi)部化合物的應(yīng)用域，通常采用杠桿法(Williams圖)實(shí)現(xiàn)可視化[23-25]，結(jié)果如圖5所示.

圖5 pEC10和pEC50模型的應(yīng)用域Fig.5 Applicability domain of pEC10 and pEC50 model

由圖5可知：月牙藻的pEC10和pEC50兩個模型中，訓(xùn)練集和測試集中的大部分化學(xué)品都被標(biāo)準(zhǔn)殘差(-3，3)和臨界帽子值(h*＝0.3851)包圍.結(jié)果表明，兩種模型均具有廣泛的應(yīng)用性.對于 h大于 h*的化合物，如點(diǎn) 10(氯己定)和點(diǎn) 237(腺苷二磷酸)，通常被認(rèn)為是好的影響點(diǎn)，這些點(diǎn)會使模型更加穩(wěn)定.對于交叉驗(yàn)證的標(biāo)準(zhǔn)殘差大于 3個標(biāo)準(zhǔn)偏差單位的化合物，如點(diǎn) 132(環(huán)嗪酮)則認(rèn)為是屬于響應(yīng)異常的范疇.總之，模型具有廣泛的應(yīng)用性.

2.6 本文模型與文獻(xiàn)模型的對比

將本文模型的統(tǒng)計學(xué)數(shù)據(jù)與其他文獻(xiàn)模型的進(jìn)行比較，結(jié)果見表3.Lee等[12]和 Huang等[13]的模型有令人滿意的預(yù)測結(jié)果和較高的 R2，但模型中涉及的化學(xué)物質(zhì)較少，種類也不豐富.Khan等[14]和 Yu[15]的模型結(jié)果中，使用計算的 logKow數(shù)據(jù)(MLOGP)時，不能完全避免誤差放大的可能性.在有機(jī)物數(shù)量豐富的情況下，本工作模型有著更高的 R2和 Q2，表明此模型的預(yù)測值更準(zhǔn)確，模型更穩(wěn)定.

表3 本工作與文獻(xiàn)研究的比較Tab.3 Comparison of this research with references

3 結(jié) 論

基于本課題組提出的分子結(jié)構(gòu)的范數(shù)描述符建立 QSTR模型，預(yù)測有機(jī)物對月牙藻的生態(tài)毒性pEC10和pEC50.結(jié)果表明：本工作模型對月牙藻生態(tài)毒性(pEC10和 pEC50)的預(yù)測精度較高；統(tǒng)計結(jié)果和模型驗(yàn)證結(jié)果顯示了本工作模型的穩(wěn)定性和廣泛應(yīng)用性；基于原子分布矩陣構(gòu)建的范數(shù)描述符能夠準(zhǔn)確描述有機(jī)物的分子結(jié)構(gòu)，據(jù)此建立的模型對于有機(jī)物的生態(tài)風(fēng)險評價具有重要意義.