本發(fā)明涉及一種通過建立定量構效關系模型(qsar)預測水相中有機物與硫酸自由基反應速率常數的方法，屬于生態(tài)風險評價測試策略技術領域。

背景技術：

高級氧化技術(advancedoxidationprocesses,aops)是水污染控制中常用的工藝，在去除水體含有的微量污染物方面具有良好的發(fā)展前景。其原理是通過激發(fā)產生羥基自由基、硫酸自由基等活性物種，進而攻擊水體中的各類污染物，將其降解為co2、水和礦物鹽類，從而實現水體中毒害物質的脫除。根據自由基激活方式的不同，高級氧化技術分為：光化學氧化、電化學氧化、光催化氧化、濕式催化氧化、聲化學氧化、濕式氧化及fenton氧化等。

傳統(tǒng)的高級氧化技術大多基于羥基自由基作為活性物種來完成降解反應，但這一體系對反應條件的要求(較低ph下進行反應)比較苛刻，故具有一定的局限性?；诹蛩嶙杂苫?so4^·-)為活性物種的高級氧化技術，具有較寬的ph適用性和較強的氧化能力，越來越受到研究者的關注。有機污染物與so4^·-的水相反應速率常數是表征污染物與so4^·-反應能力與強度的物理量，具有越大數值的化合物具有越強的so4^·-反應活性，越容易被so4^·-降解。因此，獲取污染物的數值，有助于污水處理中選擇合適的工藝參數(如設計反應時間，藥物的投加量等)，也有助于評價該污染物在環(huán)境中的持久性與生態(tài)風險。

化合物的數值的主要實驗測定方法是競爭動力學法，該方法雖然可以獲得比較準確的結果，但成本較高、耗費人力物力。目前有報道的有數值的化合物僅約200種，遠不能滿足so4^·-氧化技術工程參數的需求，也無法滿足化學品風險管理的需要。截止2017年5月，美國化學文摘社(chemicalabstractsservice,cas)登記的化學品已超過1.29億種，已經得到市場化的化學品超過14萬種，面對數量如此龐大的化學品，單靠實驗測定顯然無法滿足對的需要。因此，非常有必要發(fā)展可靠的實驗替代技術用于獲取有機物的數值。

定量結構-活性關系(quantitativestructureactivityrelationship,qsar)是指有機污染物的分子結構與其理化性質、環(huán)境行為或毒理學效應參數之間的定量關系，通過建立qsar模型對有機物的各項性質進行預測，具有降低測試費用、縮短測試時間、彌補實驗數據缺失、評估試驗數據不確定性等諸多優(yōu)勢，因此qsar技術在化學品風險管理中得到了廣泛應用。為此2007年世界經濟合作與發(fā)展組織(oecd)提出了qsar模型發(fā)展的使用準則：(1)具有明確定義的環(huán)境指標；(2)具有明確的算法；(3)定義了模型的應用域；(4)模型具有適當的擬合度、穩(wěn)定性和預測能力；(5)最好能夠進行機理解釋。該準則的提出為qsar模型的發(fā)展和使用指明了方向。

目前有報道的針對有機物與硫酸自由基水相反應速率常數的qsar模型共有三例。文獻“environ.sci.technol,2015,49：13394-13402”選取分子的氧碳比及l(fā)umo與homo的能量差作為分子描述符，預測了85種有機物分子(含有的基團包括碳碳雙鍵、碳碳三鍵、醇羥基、酚羥基、醛基、羰基、羧基、酯基、氨基、醚鍵、氯原子、溴原子、二硫鍵等)的數值，相關系數r²達到0.86以上，但所含的化合物數目較少，模型應用域偏窄。文獻“rscadv.2016,110：108448-108457”利用padel描述符基于單一決策樹和梯度提升決策樹兩種機器學習方法建立了預測同樣85種有機物分子數值的qsar模型，該文獻報道的模型參數雖然較好，但缺少明確的計算公式，算法不透明，很大程度上限制了該模型的應用。文獻“waterres.2017,116：106-115”基于基團貢獻方法建立了兩個預測113種有機物分子數值的qsar模型，該模型僅需分子結構式即可確定各描述符的值并計算出反應速率常數，但缺點在于模型的外部驗證能力不理想(外部驗證系數r²ext僅有0.42-0.62，即預測能力不佳)。

綜上，基于硫酸自由基的反應在高級氧化技術中具有良好的發(fā)展前景，但目前數值缺失還十分嚴重，限制了這一技術的發(fā)展?jié)摿?，而對該參數的實驗測定耗時費力，難以大批量進行，現有的qsar模型覆蓋化合物數量較少，缺少同時具備性能良好、使用方便等優(yōu)勢的模型。因此有必要基于現有實驗數據的基礎上，發(fā)展覆蓋化合物范圍更廣的、便于應用的、可準確預測的qsar模型，以滿足化學品風險管理的需求。

技術實現要素：

本發(fā)明提供了一種簡便、準確、高效預測有機化學品與硫酸自由基水相反應速率常數的方法，該方法可以根據有機物分子結構預測其從而得以評估其環(huán)境持久性及水環(huán)境行為和歸趨，進而為化學品風險評價與管理提供支持，并為基于硫酸自由基的高級氧化技術發(fā)展提供必要的數據支撐。

本發(fā)明的技術方案：

一種預測水相中有機物與硫酸自由基水相反應速率常數的定量結構活性關系模型，步驟如下：

收集得到197種有機化合物的實驗測定值，根據目標化合物名稱和cas號，確定目標化合物對應的結構，接下來對各分子結構進行mm2預優(yōu)化，在此基礎上，以pm6eps＝78.6charge＝0ef1scfgnorm＝0.0100mullikpolardebugshift＝80為關鍵詞，對三維mol文件進行pm6結構優(yōu)化，生成最終用于計算描述符的mol文件和out文件；

所述的有機化合物中含有碳碳雙鍵、碳碳三鍵、醇羥基、酚羥基、羰基、醛基、羧基、酯基、酰胺基、硝基、氨基、氰基、醚鍵、二硫鍵、氟原子、氯原子、溴原子、碘原子和/或砷原子。

有機化合物的結構通過量子化學描述符和dragon描述符進行描述，量子化學描述符從優(yōu)化后的結構中提取，dragon描述符根據優(yōu)化后的結構，計算得到4885個dragon描述符，對符合以下原則的dragon描述符進行舍棄：(1)標準偏差值低于0.0001的描述符；(2)常數或近常數描述符；(3)相關性不小于0.99的描述符；(4)存在缺失值的描述符；最終得到1241個dragon描述符。

將197種有機化合物隨機分成5組，按照4:1的比例分為訓練集和驗證集；用去多法(leave-many-out)將每組各做一次驗證集，同時其余四組為訓練集，使用逐步多元線性回歸(mlr)方法對訓練集數據進行擬合建模，用決定系數r²adj和均方根誤差rmse表征模型擬合度，按照擬合度的高低進行排序；然后進行內部驗證和外部驗證，用去一法交叉驗證系數q²loo和bootstrapping法交叉驗證系數q²boot表征模型的穩(wěn)健性；用r²ext、q²ext及rmseext表征模型預測能力；綜合模型的擬合度、穩(wěn)健性、預測能力等多項指標選出最合適的模型，所得模型如下：

logk＝2.383×(avs_b(p))+0.614×(ehomo)-1.043×(gats3m)-0.566×(saaac)+0.716×(narno2)-0.312×(mpc05)+0.482×(ic1)-1.625×(spmax_ea(dm))+8.317ntrain＝158,m＝8,r²adj＝0.819,q²loo＝0.805,q²boot＝0.787,rmsetrain＝0.486,next＝39,r²ext＝0.786,q²ext＝0.784，rmseext＝0.530

式中，k指有機物與硫酸自由基水相反應速率常數，單位為m^-1·s^-1。ntrain為訓練集包含的化合物個數，m為模型選用的描述符個數，r²adj為調整決定系數，q²loo為去一法交叉驗證系數，q²boot為bootstrapping法交叉驗證系數，rmsetrain為訓練集的均方根誤差，next為驗證集包含的化合物個數，r²ext為外部驗證的決定系數，q²ext為外部驗證系數，rmseext為外部驗證的均方根誤差。

模型建立采用了8個dragon描述符，分別為avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)。其中，avs_b(p)是極化率加權的載荷矩陣的平均頂點加和；ehomo是分子的最高占據分子軌道能量；gats3m是質量加權的lag3geary自相關系數，表征原子特性的相關性；saaac是aaac電狀態(tài)的數量；narno2是分子含有的連接在芳香集團上的硝基數量；mpc05是5階分子路徑計數；ic1是信息含量指數(1階鄰域對稱性)；spmax_ea(dm)是偶極矩加權的邊界連接指數標準化后的主特征值。

所得模型的r²adj＝0.819，表明模型具有良好的擬合能力，q²loo＝0.805，q²boot＝0.787，表明模型具有很強的穩(wěn)健性，r²ext＝0.786，q²ext＝0.784，表明模型具有良好的預測能力，能夠成功應用到訓練集以外的化合物中，預測值與實驗值擬合關系見附圖1。從圖中可見預測值和實驗值擬合較好，說明模型預測能力良好，可以應用在訓練集以外的化合物。

用標準殘差s對杠桿值h的williams圖對模型的應用域進行表征(見圖2)，一般認為，s的絕對值大于3.0時，該化合物是離群點。當杠桿值h大于警戒值h^*時，表明該化合物結構與其他化合物結構有顯著差異。

圖2中h^*為警戒杠桿值h^*＝3(m+1)/n＝3(8+1)/158＝0.171,其中m為描述符個數，n為訓練集化合物個數。由圖可見，共七種化合物的杠桿值h高于警戒杠桿值h^*，表明這七種化合物的結構與數據集中其他化合物結構存在較大差異，但標準殘差的絕對值均小于3，表明本模型可以推廣到這幾種化合物中。全部化合物的標準殘差均在(-3，+3)范圍內，沒有離群點的出現，這表明模型具有較理想的外推能力。七種杠桿值超過警戒杠桿值的化合物分別為乙酸(cas號：64-19-7)、萘普生(cas號：22204-53-1)、色氨酸(cas號：73-22-3)、2-甲硫基苯并噻唑(cas號：615-22-5)、硫唑嘌呤(cas號：446-86-6)、苯并噻唑(cas號：95-16-9)以及1-異丙基氨基-3-(萘-1-氧基)丙-2-醇(cas號：525-66-6)。其中除乙酸外，其他化合物均含有兩個以上的環(huán)結構，這可能是它們與其它化合物結構存在較大差異的原因。

本發(fā)明的有益效果：采用本發(fā)明方法可以有效地預測有機化合物的數值。該方法具有預測準確、成本較低、操作簡便、易于程序化、人力物力消耗較小等優(yōu)勢，本模型為現有的涵蓋化合物種類和數量最多的模型，擬合優(yōu)度、穩(wěn)健性與預測能力良好，并進行了應用域表征，建立和驗證的各個過程嚴格遵守了oecd模型構建及使用導則。因此，使用通過本模型得到的有機化合物的數值，可以為化學品風險評價與管理提供數據支持，對污染物的生態(tài)風險評價具有重要意義。

本發(fā)明的有益效果：

1.模型可預測的化合物范圍較廣，主要包括含有碳碳雙鍵、碳碳三鍵、醇羥基、酚羥基、羰基、醛基、羧基、酯基、酰胺基、硝基、氨基、氰基、醚鍵、二硫鍵、氟原子、氯原子、溴原子、碘原子、砷原子等基團的化合物，相比前人工作增加了含有酰胺基、硝基、氰基、氟原子、碘原子、砷原子等基團的化合物，針對這些化合物數值的預測，將有助于污水處理中判斷目標物是否適合利用硫酸自由基氧化降解，選擇合適的高級氧化工藝、估算藥物的投加量與投加速率，也有助于評價該類污染物在環(huán)境中的持久性與生態(tài)風險；

2.本發(fā)明采用多元線性回歸方法，共采用了8個dragon描述符建立模型，所得模型形式簡潔，透明性好，易于程序化推廣應用；

3.本發(fā)明所得模型建立和驗證的各個過程嚴格遵守了oecd模型構建及使用導則，所得模型擬合優(yōu)度高，穩(wěn)健性好，預測能力強，有明確的機理解釋和應用域。

附圖說明

圖1為訓練集和驗證集的實測值與預測值的擬合圖。

圖2為模型的williams圖。

具體實施方式

以下結合附圖和技術方案，進一步說明本發(fā)明的具體實施方式。

實施例1

隨機給定羧酸類化合物丁二酸(cas號為110-15-6)，預測其數值。首先優(yōu)化丁二酸的分子結構，然后基于優(yōu)化的分子結構，計算出8種分子描述符avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)的數值，分別為2.79,-11.506,1.293,0,0,1.609,2.522,0.334。由matlab計算的h＝0.036<0.171＝h^*，所以該化合物在應用域內。將描述符值代入所建模型，計算結果為6.72，實驗值為6.85，預測結果良好。

實施例2

隨機給定含碳碳雙鍵、氰基的化合物丙烯腈(cas號為107-13-1)，預測其數值。

首先優(yōu)化丙烯腈的分子結構，然后基于優(yōu)化的分子結構，計算出8種分子描述符avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)的數值，分別為3.081,-10.982,0.903,0,0,0,2.128,0.678。由matlab計算的h＝0.131<0.171＝h^*，所以該化合物在模型的應用域內。將描述符值代入所建模型，計算結果為7.90，實驗值為7.91，預測結果良好。

實施例3

隨機給定含氨基、硫原子的抗生素類化合物磺胺噻唑(cas號為72-14-0)，預測其數值。

首先優(yōu)化磺胺噻唑的分子結構，然后基于優(yōu)化的分子結構，計算出8種分子描述符avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)的數值，分別為3.53,-8.886,1.084,0,0,3.555,3.433,0。由matlab計算的h＝0.032<0.171＝h^*，所以該化合物在模型的應用域內。將描述符值代入所建模型，計算結果為10.69，實驗值為10.44，預測結果良好。

實施例4

隨機給定含醇羥基、酚羥基的荷爾蒙類化合物雌二醇(cas號為50-28-2)，預測其數值。

首先優(yōu)化雌二醇的分子結構，然后基于優(yōu)化的分子結構，計算出8種分子描述符avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)的數值，分別為3.427,-9.235,1.065,0,0,4.585,3.009,0。由matlab計算的h＝0.025<0.171＝h^*，所以該化合物在模型的應用域內。將描述符值代入所建模型，計算結果為9.72，實驗值為9.08，預測結果良好。

實施例5

隨機給定含氨基的含砷有機物對氨基苯胂酸(cas號為98-50-0)，預測其數值。

首先優(yōu)化對氨基苯胂酸的分子結構，然后基于優(yōu)化的分子結構，計算出8種分子描述符avs_b(p)，ehomo，gats3m，saaac，narno2，mpc05，ic1和spmax_ea(dm)的數值，分別為3.35,-8.583,1.07,0,0,2.944,3.09,0。由matlab計算的h＝0.028<0.171＝h^*,所以該化合物在模型的應用域內。將描述符值代入所建模型，計算結果為10.48，實驗值為10.29，預測結果良好。

表1建立模型所用化合物信息