基于迭代自組織數據分析算法與蟻群算法建立有機物黏度的QSPR模型

時靜潔 陳利平 陳網樺

(南京理工大學化工學院安全工程系,南京210094)

1 引言

結構決定性能是化學中的一條基本規律.定量構效關系(QSPR)方法可以將化合物的結構信息、物理化學參數與化合物的性質進行分析計算,得到表征分子結構的描述符,建立合理的數學模型,發現和確定對有機物的物理化學性質起決定作用的結構因素.1-3目前,已有研究者4-8運用QSPR方法對預測黏度進行了一定的研究.如張鶯4采用分子電性距離矢量(MEDV)表征酯分子的結構,分別應用多元線性回歸(MLR)和偏最小二乘(PLS)回歸進行了線性建模,結果較好,但只適用于酯類液體,涉及到的化合物較少,使用受到限制;而且并未對樣本集進行訓練集和測試集劃分,所得結果有待驗證.Gharagheizi5對聚合物溶液黏度進行了QSPR研究,分別建立了遺傳算法與多元線性回歸(GA-MLR)模型和徑向基函數神經網絡(RBFNN)模型,其結果令人滿意,但該研究并未對所建立的模型進行全面的評價與驗證,且模型僅適用于聚合物溶液,使用范圍較窄.

針對以上一些局限性,本文進行了一些改進.在對樣本集進行初步分類的基礎上,再對訓練集和測試集進行劃分,這也是本文的創新點之一.由于樣本中一些化合物結構比較復雜,且有些化合物同時屬于幾個類別,很難將其分類,因此采用了具有試探性和人機交互功能的動態聚類算法——迭代自組織的數據分析法(ISODATA)對其進行分類.本文采用蟻群優化算法(ACO)對分子描述符進行篩選,且分別與多元線性回歸法和支持向量機法組合建立預測黏度的模型,并對所建立的模型進行了機理解釋.此外,本文還對所建立的模型進行了比較全面的評價與驗證,對模型的應用域做了一定的研究.

2 樣本與方法

2.1 分子描述符的計算

分子描述符是建立QSPR模型的基礎.分子描述符是用邏輯和數學程序將分子結構中的化學信息轉化成數值的符號表征.本文首先借助化學軟件Hyperchem7.5進行分子結構的輸入,隨后采用分子力學方法MM+進行初步優化,用量子化學半經驗方法PM3進行進一步幾何優化,所有計算都限制在Hartree-Fock能級,采用Polar-Ribiere方法直至RMS梯度達到0.41868 kJ·mol-1,最后獲得能量最低的穩定構型.將優化好的分子結構導入DRAGON2.1軟件計算分子描述符,其中包括組成描述符、拓撲描述符、幾何描述符、電性描述符等18類共1481種分子描述符.

2.2 ISODATA算法

ISODATA算法是一種常用的動態聚類劃分算法.目前,ISODATA算法在油水層判別問題,9原油開發中判斷漏失層的存在以及漏失層的位置,10復雜體制雷達信號分選,11網絡路由選擇算法12等方面得到了很好的應用.鑒于ISODATA算法的上述優點,本文提出基于ISODATA算法對樣本集進行初步分類.

ISODATA主要算法思想是首先根據最小距離準則獲得初始聚類,然后判斷初始聚類結果是否符合要求.若不符,則將聚類集進行分裂和合并處理,得到新的聚類中心,再判斷聚類結果是否符合要求.如此反復迭代直到完成聚類劃分操作.ISODATA算法步驟13如下:(1)設置聚類分析控制參數是ISODATA的關鍵步驟,主要包括:期望得到的聚類數K;一個聚類中的最少樣本數θN,如小于此數就不作為一個獨立的聚類;一個聚類域中樣本距離分布的標準差θS;θC是兩聚類中心之間的最小距離,若小于此數,合并兩個聚類;一次迭代運算中可以合并的聚類中心的最多對數L;允許迭代的代數為I.(2)初始分類,讀入準備分類的N個模式樣本(Xi,i=1,2,…,N),預選Nc個樣品作為初始聚類中心,根據最小距離準則將各樣本分類.(3-5)按控制參數給定的要求,將前一次獲得的聚類集進行分裂和合并處理,以獲得新的聚類中心和分類集(其中(4)為分裂處理,(5)為合并處理).(6)再次迭代運算,重新計算各項指標,判別聚類結果是否符合要求,以此反復經過多次迭代運算,直至得到理想的聚類結果.

2.3 樣本集的確定

由于黏度受溫度等因素影響較大,因此本文選用了310種有機物在25°C的液體黏度作為研究樣本,并將所有的黏度值取自然對數(用lnη表示)作為模型的因變量.可靠的QSPR預測模型必須以可靠的實驗樣本為基礎,因此為了最大程度地消除實驗數據之間的差異可能給預測模型所帶來的影響,本研究選用權威數據庫《有機化合物實驗物性數據手冊》14作為實驗樣本來源,以確保實驗數據的準確性和權威性.

隨后,本文并未對樣本集進行直接劃分.應用ISODATA算法對計算出的所有化合物描述符進行聚類分析,再在每個類別中隨機選擇樣本作為訓練集和測試集,訓練集和測試集中均涵蓋了每個類別的化合物,避免訓練集中缺失測試集中的部分信息,使得模型預測能力受到影響.本文選擇253個物質作為訓練集,用于變量選擇和建立模型,57個物質作為測試集,用于外部驗證.

2.4 分子描述符的確定

大量的分子描述符中存在冗余信息,故首先利用DRAGON2.1軟件的篩選功能進行初步篩選,篩除取值為常數或者近似常數的描述符以及相關系數達到0.95以上的描述符.經過預先篩選,描述符的數目得到了一定程度的降低,減少至628個分子描述符(見表S1,Supporting Information),但仍無法滿足QSPR建模的需要.蟻群優化算法(ACO)作為一種新興的種群智能優化算法,已經在離散的組合優化等問題中被證明是一種高效的變量選擇方法.而參數選擇問題正屬于離散的組合優化問題.因此,我們將蟻群優化算法用于分子描述符參數的變量篩選,以期取得良好的效果.作為一種全局搜索算法,蟻群算法能夠有效地避免局部極優.15因此采用蟻群算法進一步篩選,最終確定5個描述符作為模型的輸入.

3 結果與討論

3.1 樣本集的劃分

運用MATLAB語言來實現ISODATA算法,其主要思路是,把類的分裂合并操作看成是一種二維數組中行矢量位置移動的過程,每一個樣本作為數組中的一個行矢量,而每一行的每一列都是樣本的屬性值,使用MATLAB的矩陣運算可以完成對樣本位置的調整,從而模擬對類的調整,最終達到聚類分析的結果.16本文向量的行數310為化合物數,列數628為屬性值,即預篩選后的分子描述符數.ISODATA算法部分參數設置為:K=8,θN=3,θS=1,θC=2,L=10,I=20.基于ISODATA算法聚類分析對化合物的分類如表1所示,類中心間的距離如表2所示.

表1和表2中的距離均為歐氏距離,歐氏距離是最易于理解的一種距離計算方法,源自歐氏空間中兩點間距離公式.兩個n維向量a(x11,x12,…,x1n)與b(x21,x22,…,x2n)間的歐氏距離:由表1可以看出,每一類化合物與相應的類中心距離比較小,說明每一類化合物的相似度比較高.由表2可以看出,類中心間的距離均比較大,說明類與類之間差異度比較大.完成對化合物的初步分類之后再從每個類中隨機選擇樣本劃分訓練集和測試集,其結果見表S2(Supporting Information).

3.2 描述符的篩選結果

運用蟻群算法對DRAGON2.1軟件初步篩選后的628個分子描述符進行進一步篩選,其過程是在VC++6.0中采用C語言編程實現.將螞蟻數量設為500,揮發率設為0.9.進行不斷迭代直至收斂,最后篩選出5個分子描述符,分別為 G3u、Hy、H1p、R5m、X3sol.G3u屬于WHIM(權重整體不變分子指數)描述符,G3u定義為第三成分對稱定向WHIM指數,主要表征分子的結構對稱性.Hy是一個與化合物親水性相關的簡單經驗指數,它是通過有關親水性官能團的函數獲得,主要反映該分子的親水性特征.H1p屬于GETAWAY(Geometry,Topology and Atom-Weights Assembly,幾何、拓撲和原子權重組合)描述符,以原子極化率作為加權性質計算得到的參數,用來反映相鄰原子之間原子極化率的相互關系.R5m也是GETAWAY描述符,它根據杠桿矩陣單元計算而得,該杠桿矩陣通過面心原子坐標獲得.該描述符主要表征原子量的分布.X3sol為拓撲描述符,由分子圖論獲得,是一種溶解連接性指數.175個分子描述符相關性分析如表3所示.

相關系數大于0稱為正相關,否則為負相關.相關系數絕對值越大,則相關性越強.從表3中可看出其相關系數均小于0.8,這說明參數之間不存在共線性問題,可作為建模的輸入參數.18

3.3 ACO-MLR模型

針對訓練集樣本,以ACO篩選的5個描述符作為輸入參數,應用SPSS軟件中的多元線性回歸(MLR)模塊,得到了MLR預測模型,結果如下:

式中lnη為黏度自然對數值.MLR模型的性能參數見表4,表中R2為復相關系數,為“留一法”交互內部驗證的復相關系數,為交互外部驗證的復相關系數,為模型外部預測能力的驗證系數,一致性相關系數(CCC)為模型調和系數,RMSE為均方根誤差,n為樣本數.方程中,回歸系數為正,表明此描述符與黏度值正相關,反之,回歸系數為負則為負相關,且系數越大,相關程度越高.由此得出,在此模型中,各描述符與黏度的相關性大小順序依次為:Hy、R5m、H1p、X3sol、G3u.根據各描述符回歸系數的正負號和大小可知,Hy越大,即分子中所含親水基團越多,親水性越強,則黏度越大;R5m越大,即分子量越大,則黏度也越大;H1p越大,即相鄰原子之間原子極化率越強,則黏度也越大;X3sol越大,即分子中的溶解連接性指數越大,則黏度也越大;而G3u指的是分子的結構對稱性,其前面的系數最小,相關性最低,說明它對黏度的影響最小.ACO-MLR模型對訓練集和測試集的實驗值和預測值的結果見表S2.

表1 基于ISODATA算法聚類分析對化合物的分類Table 1 Classification of compounds based on ISODATAclustering analysis

表2 類中心間的距離Table 2 Distances among cluster centroids

表3 分子描述符相關系數Table3 Correlation coefficients of molecular descriptors

表4 ACO-MLR和ACO-SVM模型的主要性能參數對比Table4 Performance parameter comparison between ACO-MLR andACO-SVM models

隨后,對測試集樣本的黏度進行預測,以驗證模型的外部預測能力.ACO-MLR模型所得的黏度預測值與實驗值的相關系數為0.934,其預測值與實驗值的比較見圖1.

3.4 ACO-SVM模型

隨后,為了進一步對黏度與分子結構間可能存在的非線性關系進行研究,并獲得性能更優的預測模型,本文以MLR模型所選擇的分子描述符作為輸入變量,以黏度自然對數值作為輸出變量,應用非線性的支持向量機(SVM)方法建立新的QSPR模型.用SVM做預測時,相關參數(主要是懲罰參數C和核函數參數γ)的選擇是個難點,參數選擇不好,將會嚴重影響預測的精度和準確率.

圖1 ACO-MLR模型對測試集所得預測值與實驗值的比較Fig.1 Comparison between the predicted and experimental values byACO-MLR for test set

為了得到適合的相關參數,選擇了應用最為廣泛的徑向基函數(RBF)作為核函數,采用格點搜索(GS)的方法來選擇最佳的參數組合.懲罰系數C和RBF核函數的寬度γ的搜索范圍為[2-8,28],步長為1,并根據對訓練集進行“留一法”交互驗證所得的來確定最佳的模型參數.最終選擇的最優參數為:懲罰系數C=64,核函數的寬度γ=1.4142,不敏感損失函數ε=0.1.ACO-MLR模型對訓練集和測試集的實驗值和預測值的結果見附表S2.所得的SVM模型的主要性能參數見表4,ACO-SVM模型所得的黏度預測值與實驗值的相關系數為0.950,其預測值與實驗值的比較見圖2.

3.5 結果比較與分析

從圖1和圖2中可以看出,ACO-MLR模型和ACO-SVM模型對測試集中57個樣本的預測值均與實驗值有較好的一致性,預測精度令人滿意.比較模型中訓練集和測試集的預測結果發現,各子集的復相關系數均比較高,預測誤差較低,而且比較接近.一般認為,若均大于0.6,CCC大于0.85,19,20則說明所建立的模型不但比較穩定,而且具備較強的預測能力和泛化推廣性能.

充分評價一個模型的預測能力,針對測試集必須要滿足以下三個條件:21(1)復相關系數R2大于0.6;(2)為實驗值與預測值的相關系數的平方;(3)通過原點的回歸線斜率k和k′在0.85和1.15之間.由表4可知,測試集的復相關系數R2均大于0.6.ACO-MLR模型中,測試集的0.8724,模型中,測試集的兩模型的值均遠遠小于0.1.通過擬合可得ACOMLR模型的k=0.9313,k′=0.9445,ACO-SVM模型的k=1.0646,k′=0.8513.計算結果表明,兩模型通過原點的回歸線斜率k和k'均在0.85和1.15之間.

圖2 ACO-SVM模型對測試集所得預測值與實驗值的比較Fig.2 Comparison between the predicted and experimental value byACO-SVM for test set

式中,r2是實驗值與預測值的回歸線斜率,為交換橫縱坐標后實驗值與預測值的回歸線斜率,其回歸線均需通過原點.是實驗值為縱坐標時所計算的值,而是實驗值為橫坐標時所計算的值.計算式如式(5)所示：

隨后,本文對ACO-MLR模型和ACO-SVM模型的樣本集進行了殘差分析,討論模型在建立過程中是否有系統誤差產生,兩個模型的殘差圖分別見圖3和圖4.由圖3和圖4可以看出,各模型的計算殘差均隨機分布于基準線的兩側,不存在明顯的規律性.由此可以推斷,兩個預測模型在建立過程中未產生系統誤差.

為檢驗所建模型是否存在機會相關,本文將“Y-隨機性檢驗”方法針對ACO-MLR和ACO-SVM模型分別運行50次.對于ACO-MLR模型所得最大R2為0.3360,ACO-SVM模型為0.7196,其結果均不如原始模型.由此可見,本文所建立的預測模型不存在“偶然相關”現象,具備較強的穩定性.

用Williams圖表征模型的應用域(AD).當化合物的標準殘差落在(-3,+3)以外時,認為其實驗值為離群點;當化合物的臂比值(Hat(hi))大于警戒值(h*)時,認為化合物顯著影響模型的回歸效果.27

圖3 ACO-MLR模型預測值殘差與黏度實驗值lnη關系圖Fig.3 Plot of the residuals of prediction values versus the experimental lnη values forACO-MLR model

圖4 ACO-SVM模型預測值殘差與黏度實驗值lnη關系圖Fig.4 Plot of the residuals of prediction values versus the experimental lnη values forACO-SVM model

由圖5和圖6可以看出,在ACO-MLR和ACOSVM模型中,測試集中沒有出現“異常值”物質.在訓練集中,ACO-MLR模型和ACO-SVM模型分別有四個物質和七個物質的標準殘差落在了(-3,+3)以外.ACO-MLR模型中的四個物質分別是叔丁醇、二丙酮醇、甲醇、1,4-丁二醇.ACO-SVM模型中的七個物質分別是1,1,1-三氟乙烷、1-癸醇、1-壬醇、3-甲酚、1,4-丁二醇、1-辛醇、三溴甲烷.這些物質的實驗值可能出現了比較大的偏差,被稱為“異常值”.除此以外,兩模型中十氟丁烷、甘油、鄰苯二甲酸雙(2-乙基己基)酯)、鄰苯二甲酸二辛酯,這四個物質的臂比值hi均超過了警戒值h*,可能有兩種原因造成.第一是這些分子的某些結構特征并未被所篩選出的分子描述符很好表征;第二是這些分子的某些結構對于整個樣本集來說比較特別.

十氟丁烷、甘油、鄰苯二甲酸二辛酯是訓練集樣本中的,這些物質影響變量篩選,進而影響所建立模型的質量.測試集中只有鄰苯二甲酸雙(2-乙基己基)酯的hi值超過了警戒值h*,說明模型對此化合物的預測結果可能是不可靠的,它的預測結果可以被認為是由模型外推出來的.

圖5 ACO-MLR模型訓練集和測試集的Williams圖Fig.5 Williams plot of theACO-MLR model for the train and test sets

圖6 ACO-SVM模型訓練集和測試集的Williams圖Fig.6 Williams plot of theACO-SVM model for the train and test sets

4 結論

運用迭代自組織數據分析技術(ISODATA)對樣本集進行初步分類,以避免模型中部分類別的缺失.隨后本文以蟻群算法(ACO)作為分子描述符篩選方法,分別與多元線性回歸(MLR)和支持向量機(SVM)進行組合建立了蟻群-多元線性回歸模型(ACO-MLR)與蟻群-支持向量機模型(ACO-SVM),對310個有機化合物的黏度進行了QSPR研究.其中,以ACO-SVM模型為最佳,揭示了有機化合物黏度與其分子結構間可能有較強的非線性關系,也說明了將SVM用于研究QSPR的優越性,它能有效地解決小樣本、非線性、過擬合、維數災難和局部極小等問題,泛化推廣的性能也較強.

本文所建立的2個黏度QSPR模型的預測值與實驗值非常接近,其結果令人滿意,模型性能參數均在可接受范圍之內,所建模型穩定且可用于預測黏度.本文的研究為預測有機化合物的黏度提供了一種新的有效的方法,對于化工安全設計和風險評價研究具有重要的意義.

Supporting Information:Molecular descriptors after a preliminary screening have been included in Table S1 and experimental and predicted values of the compound viscosity by ACO-MLR and ACO-SVM have been included in Table S2.This information is available free of chargeviathe internet at http://www.whxb.pku.edu.cn.

(1) Shi,J.J.;Chen,L.P.;Chen,W.H.;Shi,N.;Yang,H.;Xu,W.Acta Phys.-Chim.Sin.2012,28,2790.[時靜潔,陳利平,陳網樺,石 寧,楊 惠,徐 偉.物理化學學報,2012,28,2790.]doi:10.3866/PKU.WHXB201209273

(2) Shi,J.J.;Chen,L.P.;Chen,W.H.J.Chemom.2013,27,251.

(3)Han,C.;Yu,G.R.;Wen,L.;Zhao,D.C.;Asumana,C.;Chen,X.C.Fluid Phase Equilib.2011,300,95.doi:10.1016/j.fluid.2010.10.021

(4)Zhang,Y.Comput.Appl.Chem.2013,30,195.[張 鶯.計算機與應用化學,2013,30,195.]

(5) Gharagheizi,F.Comput.Mater.Sci.2007,40,159.doi:10.1016/j.commatsci.2006.11.010

(6)Chen,B.K.;Liang,M.J.;Wu,T.Y.;Wang,H.P.Fluid Phase Equilib.2013,350,37.

(7) Gharagheizi,F.;Mirkhani,S.A.;Keshavarz,M.H.;Farahani,N.;Tumba,K.J.Taiwan Inst.Chem.E2013,44,359.doi:10.1016/j.jtice.2012.12.015

(8)Tochigi,K.;Yamamoto,H.J.Phys.Chem.C2007,111,15989.doi:10.1021/jp073839a

(9)Feng,Y.L.;Yi,S.Q.;Feng,Z.L.;Yu,Z.G.;Xu,S.H.Pet.Geol.Oilfield Dev.in Daqing2005,24,87.[馮亞麗,伊三泉,馮卓利,于志剛,許少華.大慶石油地質與開發,2005,24,87.]

(10) Wang,S.L.;Jiang,H.Q.Transp.Porous.Med.2011,86,483.doi:10.1007/s11242-010-9634-4

(11) Zhang,H.L.;Yang,C.Z.Electro.Warfare Technol.2010,25,34.[張洪亮,楊承志.電子信息對抗技術,2010,25,34.]

(12) Ma,Y.;Tan,Z.H.;Chang,G.R.;Wang,X.Y.Procedia Eng.2011,15,2966.doi:10.1016/j.proeng.2011.08.558

(13)Yang,X.M.;Luo,Y.Mining Technol.2006,6,67.[楊小明,羅 云.采礦技術,2006,6,67.]

(14) Ma,P.S.Data Manual for the Physical Property of Organic Compounds;Chemical Industry Press:Beijing,2006.[馬沛生.有機化合物實驗物性數據手冊.北京:化學工業出版社,2006.]

(15) Li,S.Y.;Chen,Y.Q.;Li,Y.Ant Colony Algorithms with Application;Institute of Technology Press:Harbin,2004.[李士勇,陳永強,李 妍.蟻群算法及其應用.哈爾濱:哈爾濱工業大學出版社,2004.]

(16) Li,Y.P.;Li,Y.L.Mining Res.Dev.2005,25,79.[李元萍,李元良.礦業研究與開發,2005,25,79.]

(17) Todeschini,R.;Consonni,V.Handbook of Molecular Descriptors;Wileyt-VCH:Weinheim,2000.

(18)Zhang,Y.M.;Yang,X.S.;Sun,C.;Wang,L.S.Sci.China Chem.2011,41,893.[張一鳴,楊旭曙,孫 成,王連生.中國科學:化學,2011,41,893.]

(19) Chirico,N.;Gramatica,P.J.Chem.Inf.Model.2011,51,2320.doi:10.1021/ci200211n

(20) Chirico,N.;Gramatica P.J.Chem.Inf.Model.2012,52,2044.doi:10.1021/ci300084j

(21)Tropsha,A.;Gramatica,P.;Gombar,V.K.QSAR Comb.Sci.2003,22,69.

(22) Roy,K.;Matra,I.;Kar,S.;Ojha,P.K.;Das,R.N.;Kabir,H.J.Chem.Inf.Model.2011,52,396.

(23)Ojha,O.K.;Mitra,I.;Das,R.N.;Roy,K.Chemom.Intell.Lab.Syst.2011,107,194.doi:10.1016/j.chemolab.2011.03.011

(24) Roy,K.;Chakraborty,P.;Mitra,I.;Ojha,P.K.;Kar,S.;Das,R.N.J.Comput.Chem.2013,34,1071.doi:10.1002/jcc.23231

(25) Mitra,I.;Saha,A.;Roy,K.Mol.Simul.2010,36,1067.doi:10.1080/08927022.2010.503326

(26) Roy,P.P.;Paul,S.;Mitra,I.;Roy,K.Molecules2009,14,1660.doi:10.3390/molecules14051660

(27) Gu,Y.L.;Fei,Z.H.;Zhang,Y.Y.J.Anal.Sci.2012,28,333.[顧云蘭,費正皓,張玉瑩.分析科學學報,2012,28,333.]doi:10.2116/analsci.28.333