范德玲，汪 貞，王 蕾，周林軍，古 文，劉濟(jì)寧，石利利

(生態(tài)環(huán)境部南京環(huán)境科學(xué)研究所，江蘇 南京 210042)

化學(xué)品在給人們生產(chǎn)生活帶來(lái)便利的同時(shí)，也給人體健康帶來(lái)極大威脅。大量的揮發(fā)性有機(jī)化學(xué)品可通過(guò)直接排放或從土壤和水相揮發(fā)進(jìn)入大氣層。大氣中有機(jī)化學(xué)品可通過(guò)物理過(guò)程移除，如干、濕沉降；也可通過(guò)化學(xué)過(guò)程降解，如直接光解或與大氣氧化劑(OH、NO3和臭氧)反應(yīng)。在大氣對(duì)流層中有機(jī)物與臭氧的反應(yīng)是其轉(zhuǎn)化的重要途徑。表征有機(jī)化學(xué)品與O3自由基反應(yīng)的速率常數(shù)(KO3，cm3·mol-1·s-1)是反映有機(jī)污染物在大氣中持久性能力的重要參數(shù)，是進(jìn)行有機(jī)污染物生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)的基礎(chǔ)指標(biāo)[1-2]。然而，目前KO3的實(shí)驗(yàn)數(shù)據(jù)較少，且實(shí)驗(yàn)耗時(shí),費(fèi)力,成本高，不能滿足有機(jī)化學(xué)品生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)的需求。因此，有必要開發(fā)出快速有效預(yù)測(cè)KO3的方法[3]。目前，定量結(jié)構(gòu)-活性關(guān)系(QSAR)模型是用于獲取KO3的重要方法。為指導(dǎo)各國(guó)構(gòu)建滿足化學(xué)品風(fēng)險(xiǎn)管理需求的QSAR模型，經(jīng)濟(jì)合作與發(fā)展組織于2007年發(fā)布了QSAR模型構(gòu)建與驗(yàn)證的導(dǎo)則(簡(jiǎn)稱導(dǎo)則)[4]，導(dǎo)則全面闡述了符合管理要求的QSAR模型應(yīng)滿足的標(biāo)準(zhǔn)。

目前，已有研究開發(fā)了關(guān)于臭氧反應(yīng)速率常數(shù)的QSAR模型[5-6]。2007年，REN等[7]報(bào)道了116種有機(jī)化合物基于KO3的QSAR模型，該模型采用DUPLEX分類算法劃分模型訓(xùn)練集和測(cè)試集，采用啟發(fā)式方法(Heuristic Algorigthm)篩選最優(yōu)描述符，并采用多元線性回歸、支持向量機(jī)和投影尋蹤回歸方法構(gòu)建預(yù)測(cè)模型。GRAMATICA等[8]基于遺傳算法篩選最優(yōu)描述符，采用多元線性回歸算法構(gòu)建了125種有機(jī)化合物基于KO3的QSAR模型，結(jié)果顯示留一法交叉驗(yàn)證系數(shù)(QLOO2)達(dá)到82%～88%，外部驗(yàn)證決定系數(shù)(QEXT2)達(dá)到90%，均方根誤差(RMSE)達(dá)到0.73。由美國(guó)國(guó)家環(huán)境保護(hù)局有毒物質(zhì)污染防治辦公室和Syracuse Research Corporation (SRC)公司共同開發(fā)的EPI(estimation programs interface)Suite軟件，采用基團(tuán)貢獻(xiàn)法構(gòu)建了112個(gè)烯烴和炔烴不飽和有機(jī)物基于KO3的QSAR模型，相關(guān)系數(shù)達(dá)到0.94，絕對(duì)平均殘差達(dá)到0.35[9]。但上述模型并不滿足導(dǎo)則要求，如缺少模型穩(wěn)健性和預(yù)測(cè)能力表征，或未定義模型應(yīng)用域[6-9]，不利于模型使用者評(píng)估需預(yù)測(cè)的有機(jī)化合物是否處于模型應(yīng)用域內(nèi)。因此，根據(jù)導(dǎo)則要求采用簡(jiǎn)單透明的遺傳算法-多元線性回歸(GA-MLR)算法構(gòu)建基于KO3的新QSAR模型，并對(duì)模型進(jìn)行擬合優(yōu)度、穩(wěn)健性、預(yù)測(cè)能力、應(yīng)用域表征和機(jī)制解釋。所構(gòu)建的模型對(duì)實(shí)現(xiàn)環(huán)境行為參數(shù)預(yù)測(cè)軟件化具有重要意義。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源與處理

烷烴、烯烴、芳香烴、含氧揮發(fā)性有機(jī)物和酚類等152種有機(jī)化學(xué)品的KO3數(shù)據(jù)來(lái)源于文獻(xiàn)[10]。選擇-lgKO3作為模型指標(biāo)。為避免樣本分布不均勻，采用KENNARD等[11]分組方法將數(shù)據(jù)集劃分為訓(xùn)練集和驗(yàn)證集，將結(jié)構(gòu)差異較大的樣本選入訓(xùn)練集，其他與之相近的樣本選入驗(yàn)證集，從而使代表性樣本全部進(jìn)入訓(xùn)練集。訓(xùn)練集有107種化學(xué)品，驗(yàn)證集有45種化學(xué)品。

1.2 分子結(jié)構(gòu)描述符的計(jì)算

分子結(jié)構(gòu)描述符是用于反映分子結(jié)構(gòu)信息的參數(shù)，根據(jù)分子結(jié)構(gòu)按照一定理論或規(guī)則計(jì)算得到。筆者構(gòu)建的模型采用的分子結(jié)構(gòu)描述符為Dragon描述符。采用Hyperchem 7.0軟件中MM+和AM1方法對(duì)選取的152種有機(jī)化學(xué)品結(jié)構(gòu)進(jìn)行優(yōu)化[12]，采用Dragon 5.4軟件計(jì)算優(yōu)化后結(jié)構(gòu)的描述符[13]，并對(duì)得到的1 664個(gè)描述符進(jìn)行初步篩選，去掉常數(shù)項(xiàng)、近似常數(shù)項(xiàng)和高度相關(guān)的分子結(jié)構(gòu)描述符，最終得到488個(gè)分子結(jié)構(gòu)描述符。

1.3 模型的建立

采用MobyDigs軟件中遺傳算法選擇與-lgKO3高度相關(guān)的描述符[14]。由遺傳算法變量篩選得到最優(yōu)描述符，并采用多元線性回歸(MLR)方法構(gòu)建預(yù)測(cè)模型。遺傳算法相關(guān)參數(shù)設(shè)置為種群大小(population size)為100，初始模型允許的最大變量數(shù)(maximum allowed variables)為7，變異均衡值(mutation trade-off，T)為0.5，交叉(crossover)和變異(mutation)概率均基于T值。當(dāng)增加變量數(shù)目對(duì)結(jié)果影響不大時(shí)，得到8個(gè)最優(yōu)描述符。

1.4 模型的表征與評(píng)價(jià)

根據(jù)導(dǎo)則要求，對(duì)構(gòu)建的QSAR模型進(jìn)行內(nèi)部驗(yàn)證(訓(xùn)練集的擬合優(yōu)度和穩(wěn)健性評(píng)估)和外部驗(yàn)證(驗(yàn)證集的預(yù)測(cè)能力評(píng)估)。采用實(shí)驗(yàn)值與預(yù)測(cè)值之間校正后的決定系數(shù)(Radj2)和均方根誤差(RMSE，ERMS)表征模型擬合優(yōu)度，采用留一法交叉驗(yàn)證系數(shù)(QLOO2)表征模型穩(wěn)定性，采用外部檢驗(yàn)參數(shù)(QEXT2)、驗(yàn)證集相關(guān)系數(shù)(REXT2)和驗(yàn)證集均方根誤差(ERMS，EXT)等外部驗(yàn)證決定系數(shù)表征模型預(yù)測(cè)能力，基于杠桿值(leverage，hi)的Williams圖定義模型應(yīng)用域[15]。外部驗(yàn)證決定系數(shù)計(jì)算公式為

(1)

(2)

(3)

(4)

Williams圖是標(biāo)準(zhǔn)殘差(δ)和hi值定義的模型應(yīng)用域，其計(jì)算公式為

(5)

hi=xiT(XTX)-1xi。

(6)

式(5)～(6)中，xi為第i種化合物分子結(jié)構(gòu)描述符的行向量；X為n×m的矩陣，構(gòu)成訓(xùn)練集化合物的描述符空間。

當(dāng)訓(xùn)練集中化合物hi值大于警戒值(h*)時(shí)，說(shuō)明在數(shù)據(jù)集中該物質(zhì)的子結(jié)構(gòu)出現(xiàn)較少，會(huì)對(duì)模型預(yù)測(cè)結(jié)果有顯著影響。h*值計(jì)算公式為

h*=3(m+1)/n。

(7)

2 結(jié)果與討論

2.1 QSAR模型的構(gòu)建

模型描述符意義、回歸系數(shù)、回歸系數(shù)偏差和標(biāo)準(zhǔn)回歸系數(shù)見(jiàn)表1。構(gòu)建的GA-MLR回歸方程為Y=17.898-0.371X1+0.334X2+0.215X3+0.193X4+0.426X5-0.453X6-0.260X7-0.308X8，n訓(xùn)練集=107，Radj，訓(xùn)練集2=0.784，QLOO2=0.744，ERMS，訓(xùn)練集=1.127，P<0.000 1，nEXT=45，REXT2=0.664，QEXT2=0.761，ERMS，EXT=1.039。

GOLBRAIKH等[20]研究認(rèn)為，QSAR模型可接受標(biāo)準(zhǔn)為Q2>0.50和R2>0.60。由圖1可知，筆者構(gòu)建的模型擬合優(yōu)度和穩(wěn)健性較好，預(yù)測(cè)能力也較好。

表1 臭氧自由基反應(yīng)速率常數(shù)模型描述符物理化學(xué)意義和相應(yīng)系數(shù)

Table 1 Physical chemistry meaning and corresponding coefficients of the descriptors used in the -lgKO3MLR model

變量描述符定義回歸系數(shù)回歸系數(shù)偏差標(biāo)準(zhǔn)回歸系數(shù)常數(shù)項(xiàng)17.8980.480X1PW3Randic形狀指數(shù)-7.9701.471-0.371X2HOMA分子芳香性指數(shù)2.4600.4020.334X3RDF035u徑向分布函數(shù)描述符0.0830.0260.215X4G1s原子電拓?fù)浼訖?quán)指數(shù)1.4550.3800.193X5HATS2e加權(quán)原子Sanderson電負(fù)性信息2.1150.2800.426X6Nr=Cs脂肪族化合物分子中仲碳原子個(gè)數(shù)-0.9240.107-0.453X7Nr=Ct脂肪族化合物分子中叔碳原子個(gè)數(shù)-1.0950.247-0.260X8H-050與雜原子相連的氫原子個(gè)數(shù)-1.4070.228-0.308

圖1 臭氧自由基反應(yīng)速率常數(shù)MLR模型的預(yù)測(cè)值與實(shí)驗(yàn)值的擬合圖

2.2 應(yīng)用域表征

利用杠桿方法制作Williams圖分析和評(píng)價(jià)模型應(yīng)用范圍，可以圖形方式量化模型應(yīng)用范圍。模型對(duì)應(yīng)用域內(nèi)物質(zhì)預(yù)測(cè)性能較好，而對(duì)應(yīng)用域外物質(zhì)預(yù)測(cè)性能差。采用Williams圖表征的QSAR模型應(yīng)用域見(jiàn)圖2。

虛線為警戒值(h*=0.252)。

由圖2可知，數(shù)據(jù)集152種化合物中只有肼的h值≥h*(h*=0.252)，位于應(yīng)用域范圍外，為X離群點(diǎn)。所有化合物標(biāo)準(zhǔn)殘差在-3～3范圍內(nèi)，即無(wú)Y離群點(diǎn)。因此，構(gòu)建的QSAR模型可用于預(yù)測(cè)應(yīng)用域內(nèi)其他化合物-lgKO3值。

2.3 有機(jī)化學(xué)品與臭氧反應(yīng)機(jī)制解釋

通過(guò)解釋線性化合物臭氧反應(yīng)速率的QSAR模型中所選描述符的物理化學(xué)意義，可以獲得決定化合物臭氧反應(yīng)速率的結(jié)構(gòu)信息。描述符的相對(duì)重要程度由模型中每個(gè)描述符的標(biāo)準(zhǔn)回歸系數(shù)決定。標(biāo)準(zhǔn)回歸系數(shù)絕對(duì)值的大小表示對(duì)應(yīng)描述符對(duì)臭氧反應(yīng)速率影響程度的強(qiáng)弱，正負(fù)號(hào)表示對(duì)應(yīng)描述符與臭氧反應(yīng)速率相關(guān)性的正負(fù)。在模型的8個(gè)描述符中，nR=Cs和nR=Ct均為官能團(tuán)數(shù)目描述符，且標(biāo)準(zhǔn)回歸系數(shù)均為負(fù)值(表1)，這表明nR=Cs和nR=Ct與臭氧的反應(yīng)速率常數(shù)呈正相關(guān)。HATS2e為GETAWAY(geometry，topology and atom-weights assembly)類描述符，與-lgKO3呈較大負(fù)相關(guān)。HOMA為幾何描述符，PW3為拓?fù)涿枋龇?，可通過(guò)計(jì)算分子中每個(gè)原子的path數(shù)目與walk數(shù)目的比值，再將這些比值求和后除以分子中的原子數(shù)目得到。由于path/walk獨(dú)立于分子大小，所以PW3可以較好地表征分子形狀。RDF035u為徑向分布函數(shù)描述符，表示在一個(gè)半徑為R的球形體內(nèi)發(fā)現(xiàn)特定類似原子的概率。G1s為WHIM描述符，在模型中表征分子靜電拓?fù)錉顟B(tài)。H-050為以原子為中心的碎片描述符，表征與雜原子相連的氫原子個(gè)數(shù)。

2.4 模型驗(yàn)證和應(yīng)用案例

近年來(lái)，計(jì)算毒理學(xué)技術(shù)在歐美、日本和OECD得到大力發(fā)展。美國(guó)國(guó)家環(huán)境保護(hù)局研發(fā)了化學(xué)品理化性質(zhì)/環(huán)境行為指標(biāo)參數(shù)與預(yù)測(cè)模型軟件EPI Suite，其中的AOPWIN模塊采用基團(tuán)貢獻(xiàn)法預(yù)測(cè)有機(jī)化學(xué)品臭氧自由基反應(yīng)常數(shù)。OECD允許使用QSAR方法彌補(bǔ)數(shù)據(jù)缺失，并于2008年發(fā)布了第1版QSAR Toolbox工具包。其中的臭氧反應(yīng)速率數(shù)據(jù)主要來(lái)源于EPI Suite軟件。與發(fā)達(dá)國(guó)家相比，我國(guó)在計(jì)算毒理學(xué)技術(shù)研發(fā)和應(yīng)用方面具有較大差距。近年來(lái)我國(guó)已經(jīng)啟動(dòng)化學(xué)品環(huán)境安全信息預(yù)測(cè)技術(shù)研究，在一定程度上填補(bǔ)了我國(guó)化學(xué)品固有屬性預(yù)測(cè)技術(shù)的空白。其中，生態(tài)環(huán)境部南京環(huán)境科學(xué)研究所基于簡(jiǎn)化分子線性輸入規(guī)范(SMILES)解析碎片拆分技術(shù)，開發(fā)了具有我國(guó)自主知識(shí)產(chǎn)權(quán)的化學(xué)品定量構(gòu)效預(yù)測(cè)軟件[17]。

將筆者研究數(shù)據(jù)集之外的20種有機(jī)化學(xué)品-lgKO3實(shí)驗(yàn)數(shù)據(jù)[18-19]與該模型和EPI Suite軟件中AOPWIN模塊的預(yù)測(cè)結(jié)果進(jìn)行比較發(fā)現(xiàn)，20種有機(jī)化學(xué)品-lgKO3實(shí)驗(yàn)值與筆者模型預(yù)測(cè)值的決定系數(shù)(R2)達(dá)到0.794，與EPI Suite預(yù)測(cè)值的R2為0.695(表2)。其中該模型15種化學(xué)品預(yù)測(cè)結(jié)果優(yōu)于EPI Suite，EPI Suite軟件5種化學(xué)品預(yù)測(cè)結(jié)果較好。SMILES是化學(xué)物質(zhì)1維結(jié)構(gòu)的線性表達(dá)，而2維和3維結(jié)構(gòu)描述符可更全面地表達(dá)化學(xué)物質(zhì)立體結(jié)構(gòu)的空間形態(tài)。由于EPI Suite軟件基于SMILES碼碎片拆分，選取的結(jié)構(gòu)碎片也許不能完全表達(dá)分子結(jié)構(gòu)信息，同時(shí)也未給出模型應(yīng)用域，因此筆者構(gòu)建的模型彌補(bǔ)了EPI Suite軟件的不足[20-21]。

表2 20種有機(jī)化學(xué)品-lgKO3實(shí)驗(yàn)和預(yù)測(cè)數(shù)據(jù)比較

Table 2 Comparison of predicted results with experimental results for -lgKO3of 20 organic chemicals

CAS編號(hào)實(shí)驗(yàn)值EPI Suite軟件該模型預(yù)測(cè)值殘差預(yù)測(cè)值殘差000096-33-317.978 17.756 0.22218.078 -0.100 000116-14-319.036 19.348 -0.312 18.755 0.281 000123-73-918.045 17.739 0.30418.174 -0.129000497-23-418.657 16.943 1.714 17.784 0.873001630-77-918.585 18.853 -0.26818.4740.111001630-78-017.677 18.552 -0.875 18.474 -0.797 017559-81-817.744 17.195 0.549 17.724 0.020 018409-46-617.001 17.410 -0.409 17.018 -0.017 000074-86-219.318 19.522 -0.204 19.124 0.194 000075-01-418.619 18.602 0.017 18.616 0.003 000075-38-718.721 18.552 0.16918.341 0.380000078-94-417.321 17.325 -0.00417.289 0.032000108-05-417.494 17.756 -0.262 18.396 -0.902 000140-88-517.244 17.756 -0.512 17.213 0.031000359-11-518.853 18.950 -0.09718.837 0.016 000463-51-419.154 19.057 0.09719.181 -0.027006728-26-317.698 17.739 -0.04117.687 0.011 000109-92-215.812 17.057 -1.24516.022 -0.210000087-44-513.935 15.354 -1.41915.652-1.717 000360-89-420.167 18.552 1.61518.417 1.750

3 結(jié)論

該研究建立了包括烷烴、烯烴、芳香烴、含氧揮發(fā)性有機(jī)物和酚類152種有機(jī)化合物與臭氧反應(yīng)速率常數(shù)預(yù)測(cè)模型。根據(jù)經(jīng)濟(jì)合作與發(fā)展組織關(guān)于QSAR模型構(gòu)建與驗(yàn)證的導(dǎo)則要求，構(gòu)建的有機(jī)化學(xué)品與臭氧反應(yīng)速率預(yù)測(cè)模型擬合能力、穩(wěn)健性和預(yù)測(cè)能力均較好，Williams圖定義模型應(yīng)用域(AD)結(jié)果也表明該模型應(yīng)用域較廣。模型機(jī)理研究結(jié)果表明分子芳香性、電負(fù)性和仲碳原子數(shù)目是影響有機(jī)化學(xué)品與臭氧自由基反應(yīng)速率(KO3)的關(guān)鍵因素。綜上所述，構(gòu)建的有機(jī)化合物與臭氧自由基反應(yīng)速率常數(shù)QSAR模型可以用于預(yù)測(cè)應(yīng)用域內(nèi)難以測(cè)定或未知有機(jī)化合物與臭氧自由基反應(yīng)速率常數(shù)，評(píng)估其持久性，進(jìn)而對(duì)有機(jī)污染物進(jìn)行生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)。