燈盞花揮發性組分結構與保留時間關系研究

李建鳳，廖立敏*

1內江師范學院化學化工學院；2“果類廢棄物資源化”四川省高等學校重點實驗室，內江 641100

燈盞花是菊科植物短葶飛蓬(Erigeronbreviscapus)的全草，又名燈盞細辛、東菊[1]。燈盞花主要生長在我國云南、廣西、四川、貴州、西藏等高海拔地區。燈盞花性寒、微苦、甘溫辛，具有散寒解表、祛風除濕、活血化瘀、通經活絡、消炎止痛等功效[2]。目前，燈盞花在臨床上主要用于心腦血管疾病、糖尿病、腎病和眩暈等疾病的治療。Xu等[3]利用固相微萃取、氣相色譜-質譜聯用技術對燈盞花中的揮發性成分進行了定性分析。化合物結構與色譜保留時間關系研究，對于分析化合物的色譜保留行為、輔助鑒定化合物等具有重要的意義[4]。化合物結構參數化表征是建立化合物結構與性質關系的關鍵步驟之一，目前研究者在這方面做過許多工作，如He等[5,6]利用分子連接性指數和分子拓撲指數對果酒、香水百合頭中部分揮發性成分進行了結構-保留關系研究，Zhang等[7]利用分子連接性指數對燃料油中有機硫化物進行了結構-保留關系研究，均取得較好的結果。分子連接性指數和分子拓撲指數都屬于分子二維結構描述符，它的優點就是計算快速、簡便、易懂，缺點是不能區分順反異構、光學異構等現象；又如Jing等[8]利用三維全息原子場作用矢量對部分脂肪酸進行了結構-保留關系研究，Liu等[9]利用化合物三維結構計算結構參數值對飽和醇類化合物進行了結構-保留關系研究，均取得較好的結果。分子三維結構描述符，它的優點就是基于化合物分子三維立體結構計算，更加接近化合物分子實際存在的狀態，可以區分順反異構、光學異構等現象，缺點是需要進行分子結構優化、計算復雜、難懂、工作量大，有時還可能涉及到探針選取、網格劃分、分子重疊等不確定因素[10]。本文在前人的基礎上構建新的二維結構描述符，對燈盞花中出現的多類有機化合物進行結構表征，進而運用多元線性回歸(MLR)和偏最小二乘回歸(PLS)兩種方法建立化合物結構與保留時間關系模型，分析影響化合物色譜保留時間的結構因素，為揮發性化合物結構-性質關系研究提供參考。

1 材料和方法

1.1 實驗材料

以燈盞花中的64種揮發性有機化合物為研究樣本，化合物保留時間是以HP-5 ms(30 m×250 μm×0.25 μm)為色譜柱分離得到，柱箱升溫程序為：初始溫度50 ℃，保持1 min，以5 ℃/min升至220 ℃保持6 min，再以10 ℃/min升至260 ℃保持15 min。進樣口溫度280 ℃。采用脈沖不分流進樣模式。載氣使用高純氦氣，流速為1.0 mL/min。按照色譜保留時間(tR)[3]大小順序列于表1，取序號個位數為“0”和“5”的化合物(用“*”標注，共12個)為測試樣本集，測試樣本集不參與建模，僅用于評價模型的外部預測能力，其余52個化合物為訓練樣本集用于建立模型。

表1 64種揮發性有機化合物及其色譜保留時間Table 1 64 volatile organic compounds and their chromatographic retention times

續表1(Continued Tab.1)

續表1(Continued Tab.1)

多元線性回歸使用的是SPSS13.0軟件、偏最小二乘回歸使用的是SIMCA-P 11.5軟件。

1.2 實驗方法

1.2.1 化合物分子結構表征

構建化合物結構與性質之間的關系模型，化合物結構參數化表征是關鍵步驟之一。認為化合物中的非氫原子以及它們之間的關系對化合物色譜保留時間產生影響，而氫原子影響通常被忽略。不同類型的非氫原子以及不同類型非氫原子之間的關系對化合物色譜保留時間影響不同，參照文獻[11-13]方法將化合物中的非氫原子按照其連接的其它非氫原子數分為4類，與k個其它非氫原子相連的非氫原子屬于第k類非氫原子，例如與2個其它非氫原子相連的仲碳原子屬于第2類非氫原子。在參考文獻[14-16]的基礎上，對化合物中的非氫原子進行參數化染色，見式(1)。

(1)

式中i為非氫原子在分子中的編碼；ni為非氫原子i的原子核外電子層數；mi為價電子數，di為非氫原子i的成鍵電子數，也就是氧化數；Li為非氫原子i與相鄰非氫原子直接連接的化學鍵數(δ鍵取值為1，π鍵取值為0.5)。非氫原子i的ni越大，其半徑就越大，相應的原子體積就越大，相應的Zi值也就越大。

非氫原子自身對化合物色譜保留時間的影響，按照式(2)進行分類累加。

(2)

式中，k表示非氫原子i的原子類型，Zi按式(1)計算。根據非氫原子的分類，對于一個有機化合物分子中最多含有4類非氫原子，因此最終可得到4個非氫原子自身對化合物保留時間的影響項，用x1、x2、x3和x4表示。

4類非氫原子之間可以有10種不同的組合，即m11、m12、…、m44，簡寫為x5、x6、…、x14，m13表示第1類非氫原子與第3類非氫原子之間的關系，以此類推。非氫原子之間的關系并非某種具體的相互作用，而是要反映出相關程度隨著非氫原子間的距離的增大而減小，隨著非氫原子自身某種性質的增大而增大，式(3)可以滿足此要求。

(n=1，2，3，4；n≤l≤4)

(3)

Z按式(1)計算；dij為非氫原子i、j之間的相對距離(即鍵長之和與碳碳單鍵鍵長的比值，如果i、j之間有多條路徑，則以最短的為準)；n和l為原子所屬類型，α=0.5。這樣一來，對于一個化合物經參數化表達后最多可得14個變量(結構描述符)。

1.2.2 建模與評價

(4)

(5)

2 結果與討論

2.1 分子結構表征

化合物經分子結構表征后得到14個變量，由于變量數較多，在此沒有全部列出。部分用到的變量列于表2。

表2 化合物的結構參數化表征結果Table 2 Structural parameterized characterization results of compounds

續表2(Continued Tab.2)

2.2 化合物結構與色譜保留關系

2.2.1 多元線性回歸模型

本研究的訓練集樣本數為52個，而變量數多達14個，不符合N/n≥5的經驗規則，因而在確定最優模型之前應該對變量進行篩選，將與化合物色譜保留時間(tR)相關性不大的變量進行剔除。逐步回歸(SMR)是篩選變量的常用方法，建模前運用逐步回歸依據變量顯著性大小順序逐步引入候選變量，根據每一步所得模型的相關系數(R2)及標準偏差(SD)進行綜合考量后選出最佳變量組合進行回歸建模。對選定的模型變量間的共線性進行評價，計算方差膨脹因子(VIF)，如果某變量的方差膨脹因子(VIF)大于或等于20，則應減少變量建模。在逐步回歸中，相關系數(R2)和標準偏差(SD)隨變量的引入而發生變化，為便于直觀分析將其變化情況繪圖于圖1及圖2中。

圖1 R2在逐步回歸中的變化情況Fig.1 Changes of R2 in stepwise regression

圖2 SD在逐步回歸中的變化情況Fig.2 Changes of SD in stepwise regression

從圖1中可以發現，相關系數(R2)隨著變量的引入而增大，起初增大的趨勢明顯，當逐步回歸進行到第7步時相關系數(R2)接近最大值，之后相關系數(R2)增大趨勢放緩。同樣，圖2中可以發現，標準偏差(SD)隨著變量的引入而逐漸減小，起初減小的趨勢明顯，當逐步回歸進行到第7步時標準偏差(SD)達到最小值，之后標準偏差(SD)略有增大。對7變量模型進行回歸診斷，發現變量x8的方差膨脹因子(VIF)最大，且值僅為13.045 2，變量之間沒有明顯的共線性。當變量繼續引入，逐步回歸進行到第8步，雖然相關系數(R2)略有增大，然而此時最大變量方差膨脹因子(VIF)遠大于20，說明8變量模型不可靠。綜合各方面的考量后，認為應該選擇逐步回歸第7步所得變量組合(x1、x2、x3、x4、x5、x7、x8，列入表2)進行建模(M1)，如式(6)。

tR= -18.629 3+0.479 4×x1+2.434 4×x2+1.640 6×x3

-1.601 5×x4-8.006 6×x5+1.755 8×x7+4.322 7×x8(6)

2.2.2 偏最小二乘回歸模型

偏最小二乘回歸(PLS)也是化合物結構和性質關系研究中常用的方法之一，特別適合于樣本數較少而變量數較多的情況下建模。變量x1、x2、x3、...x14作為自變量X，化合物色譜保留時間(tR)作為因變量Y，建立偏最小二乘回歸模型(M2)。相關系數(R2)及交互檢驗的相關系數(RCV2)隨主成分數(A)的變化情況見圖3。圖3中可以看出當主成分數(A)達到6個時，相關系數(R2)接近最大值，交互檢驗的相關系數(RCV2)達到最大值，應該選擇6個主成分進行建模(M2)。52個訓練集樣本在PLS前2個主成分得分空間散點分布繪于圖4，圖4中可以發現全部樣本得分點都落在95%置信度的橢圓置信圈內，沒有出現一個異常點，反映出構建的結構描述符能較好地表現揮發性有機化合物的分子結構特征，并在統計模型中得到正確的表現。

圖3 相關系數隨主成分數的變化情況Fig.3 Correlation coefficient change with the number of principal components

圖4 樣本在前2個主成分得分分布Fig.4 Distribution of the top 2 principal component scores of the sample

為驗證模型結果是否為偶然所得，對模型進行20次Y向量隨機排序驗證。以Y原始向量和經過隨機排序的Y向量相關系數對模型的R2和RCV2作圖于圖5。根據Andersson等[20]提出的判斷標準，R2和RCV2在縱軸上的截距分別不得超過0.300和0.050。從圖5中可以發現所建PLS模型的R2和RCV2的截距分別為0.105和-0.882，因而認為模型的優良結果并非偶然，偏最小二乘回歸模型(M2)可以用于分析化合物結構對色譜保留時間的影響。

圖5 Y向量隨機排序驗證結果Fig.5 Y vector random sorting verification results

變量重要性可以反映出變量與Y之間的相關程度，通常認為變量重要性投影(VIP)值大于1的變量與化合物色譜保留時間(tR)相關性大。變量重要性投影見圖6，圖6可以看出變量x2、x1、x9這3個變量的VIP值大于1，說明這3個變量與化合物色譜保留時間(tR)相關性大。x1、x2分別為第1、2類非氫原子自身對化合物色譜保留時間的影響，x9為第2類非氫原子之間的關系對化合物色譜保留時間的影響。第1類非氫原子對應取代基的末端原子，而第2類非氫原子的多少決定了鏈的長度，說明化合物支鏈越多、鏈越長，即伯、仲碳原子越多，化合物可能具有較大的色譜保留時間(tR)值，這與表1中的數據特征基本吻合。

圖6 變量重要性投影圖Fig.6 Projection of variable importance

2.2.3 模型結果與比較

多元線性回歸(MLR)模型(M1)和偏最小二乘回歸(PLS)模型(M2)對訓練集化合物的色譜保留時間進行了計算，對不參與建模的預測集化合物色譜保留時間進行了預測，M1和M2對化合物的計算值及預測值分別列于表1的Cal.1和Cal.2列，Err.1和Err.2分別為誤差，Err.1%、Err.2%分別為百分誤差。兩模型對化合物保留時間的計算值與實驗值相關性，見圖7。圖7中容易看出所有樣本點都落在正方形45°對角線附近，說明兩模型對化合物保留時間計算值與實驗值相關性好，兩者間誤差不大。另外，圖上可以看出Cal.1的樣本點與Cal.2的樣本點分布大體相似，說明模型(M1)與模型(M2)的質量大體相當。

圖7 模型預測值與實驗值的相關圖Fig.7 Correlation diagram between model predicted values and experimental values

模型對化合物保留時間計算值的誤差可以反映模型預測的準確性，兩模型計算誤差分布見圖8。圖8中可以發現大部分樣本的預測誤差值都處于2倍標準偏差(±2SD)范圍以內，對于模型(M1)、模型(M2)均僅有2個化合物(3.13%)超出±2SD1、±2SD2，說明兩模型對化合物保留時間預測較為準確，誤差都處于可以接受的范圍內。另外，絕大多數樣本的Err.1小于Err.2，±2SD1線處于±2SD2線內側并更加靠近中的“0”線，說明模型(M1)的質量略優于模型(M2)。

圖8 模型對樣本保留時間預測誤差Fig.8 Model prediction error of samples′ tR

本文所研究的化合物結構跨度大，包括了33種烴類、9種酮類、8種芳香族類、4種醇類、3種醛類、3種醚類、2種酚類和2種酸類等類型化合物，含有氧、硫、氮等雜原子，含有單鍵、雙鍵、三元環、四元環、五元環、六元環、七元環、八元環、十一元環等結構，說明構建的結構描述符普適性強。對于復雜樣本體系的定量結構與保留時間關系研究已有一些報道，為便于對比，現將本文結果與部分文獻結果列于表3，可以發現本文所取得的結果明顯優于文獻。

表3 模型比較Table 3 Comparisons among different QSRR models

3 結論

通過構建非氫原子之間的關系得到新的結構描述符，對燈盞花中的64種揮發性有機化合物結構進行了參數化表征，運用多元線性回歸(MLR)和偏最小二乘回歸(PLS)方法構建了化合物結構與色譜保留時間的關系模型。經檢驗模型具有通良好的擬合能力、穩定性和外部預測能力。分析了影響化合物色譜保留時間的結構因素，結果表明化合物的伯、仲碳原子越多，色譜保留時間(tR)值就越大。構建的化合物結構描述符為二維結構描述符，計算簡便、快速、易懂、通用性強，對于化合物結構與性質關系研究具有一定的參考價值。