有機化學品與臭氧反應速率常數的定量預測模型研究

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

(生態環境部南京環境科學研究所，江蘇 南京 210042)

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

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

1 材料與方法

1.1 數據來源與處理

烷烴、烯烴、芳香烴、含氧揮發性有機物和酚類等152種有機化學品的KO3數據來源于文獻[10]。選擇-lgKO3作為模型指標。為避免樣本分布不均勻，采用KENNARD等[11]分組方法將數據集劃分為訓練集和驗證集，將結構差異較大的樣本選入訓練集，其他與之相近的樣本選入驗證集，從而使代表性樣本全部進入訓練集。訓練集有107種化學品，驗證集有45種化學品。

1.2 分子結構描述符的計算

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

1.3 模型的建立

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

1.4 模型的表征與評價

根據導則要求，對構建的QSAR模型進行內部驗證(訓練集的擬合優度和穩健性評估)和外部驗證(驗證集的預測能力評估)。采用實驗值與預測值之間校正后的決定系數(Radj2)和均方根誤差(RMSE，ERMS)表征模型擬合優度，采用留一法交叉驗證系數(QLOO2)表征模型穩定性，采用外部檢驗參數(QEXT2)、驗證集相關系數(REXT2)和驗證集均方根誤差(ERMS，EXT)等外部驗證決定系數表征模型預測能力，基于杠桿值(leverage，hi)的Williams圖定義模型應用域[15]。外部驗證決定系數計算公式為

(1)

(2)

(3)

(4)

Williams圖是標準殘差(δ)和hi值定義的模型應用域，其計算公式為

(5)

hi=xiT(XTX)-1xi。

(6)

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

當訓練集中化合物hi值大于警戒值(h*)時，說明在數據集中該物質的子結構出現較少，會對模型預測結果有顯著影響。h*值計算公式為

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

(7)

2 結果與討論

2.1 QSAR模型的構建

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

GOLBRAIKH等[20]研究認為，QSAR模型可接受標準為Q2>0.50和R2>0.60。由圖1可知，筆者構建的模型擬合優度和穩健性較好，預測能力也較好。

表1 臭氧自由基反應速率常數模型描述符物理化學意義和相應系數

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

變量描述符定義回歸系數回歸系數偏差標準回歸系數常數項17.8980.480X1PW3Randic形狀指數-7.9701.471-0.371X2HOMA分子芳香性指數2.4600.4020.334X3RDF035u徑向分布函數描述符0.0830.0260.215X4G1s原子電拓撲加權指數1.4550.3800.193X5HATS2e加權原子Sanderson電負性信息2.1150.2800.426X6Nr=Cs脂肪族化合物分子中仲碳原子個數-0.9240.107-0.453X7Nr=Ct脂肪族化合物分子中叔碳原子個數-1.0950.247-0.260X8H-050與雜原子相連的氫原子個數-1.4070.228-0.308

圖1 臭氧自由基反應速率常數MLR模型的預測值與實驗值的擬合圖

2.2 應用域表征

利用杠桿方法制作Williams圖分析和評價模型應用范圍，可以圖形方式量化模型應用范圍。模型對應用域內物質預測性能較好，而對應用域外物質預測性能差。采用Williams圖表征的QSAR模型應用域見圖2。

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

由圖2可知，數據集152種化合物中只有肼的h值≥h*(h*=0.252)，位于應用域范圍外，為X離群點。所有化合物標準殘差在-3～3范圍內，即無Y離群點。因此，構建的QSAR模型可用于預測應用域內其他化合物-lgKO3值。

2.3 有機化學品與臭氧反應機制解釋

通過解釋線性化合物臭氧反應速率的QSAR模型中所選描述符的物理化學意義，可以獲得決定化合物臭氧反應速率的結構信息。描述符的相對重要程度由模型中每個描述符的標準回歸系數決定。標準回歸系數絕對值的大小表示對應描述符對臭氧反應速率影響程度的強弱，正負號表示對應描述符與臭氧反應速率相關性的正負。在模型的8個描述符中，nR=Cs和nR=Ct均為官能團數目描述符，且標準回歸系數均為負值(表1)，這表明nR=Cs和nR=Ct與臭氧的反應速率常數呈正相關。HATS2e為GETAWAY(geometry，topology and atom-weights assembly)類描述符，與-lgKO3呈較大負相關。HOMA為幾何描述符，PW3為拓撲描述符，可通過計算分子中每個原子的path數目與walk數目的比值，再將這些比值求和后除以分子中的原子數目得到。由于path/walk獨立于分子大小，所以PW3可以較好地表征分子形狀。RDF035u為徑向分布函數描述符，表示在一個半徑為R的球形體內發現特定類似原子的概率。G1s為WHIM描述符，在模型中表征分子靜電拓撲狀態。H-050為以原子為中心的碎片描述符，表征與雜原子相連的氫原子個數。

2.4 模型驗證和應用案例

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

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

表2 20種有機化學品-lgKO3實驗和預測數據比較

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

CAS編號實驗值EPI Suite軟件該模型預測值殘差預測值殘差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 結論

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