取代苯酚類化合物生物毒性的模擬預測

廖立敏

內江師范學院化學化工學院，內江 641100

苯酚類化合物是重要的精細化學品或中間體，在農藥、醫藥、炸藥和染料等方面有廣泛應用。苯酚類化合物大多具有毒性，個別還具有較強的致癌作用。隨著現代工農業的發展，越來越多的苯酚類化合物隨著廢水排放到環境中，給動植物造成嚴重的危害。因此，全面評估苯酚類化合物的毒性，對于規范其生產、管理、使用和排放具有重要的意義。定量構效關系(QSAR)利用計算機模擬預測有機污染物性質參數，簡便、省時、節約實驗資源且降低費用，在全面獲取有機污染物各種參數方面顯示出了明顯的優勢[1-4]。QSAR研究中，化合物結構參數化是關鍵步驟之一，目前應用廣泛的主要有二維結構表征法[5-6]和三維結構表征法[7-8]。二維結構表征法的優勢是計算簡單、快速，但無法分辨順反異構、光學異構等現象；三維結構表征法的優勢是基于化合物三維立體結構計算，對各種異構體能很好地區分，但計算量大、復雜難懂。本文基于化合物分子二維結構構建一種簡易的結構描述符對化合物結構參數化表達，通過多元線性回歸(MLR)和偏最小二乘回歸(PLS)進行建模，分析影響化合物毒性的結構因素，為苯酚類化合物結構-性質研究提供參考。

1 材料與方法(Materials and methods)

1.1 實驗材料

以25個取代苯酚類化合物為研究樣本，化合物毒性以其誘發浮萍萎黃的活性值的負對數(pC)表示，其實驗值取自文獻[9]，結果如表1所示。pC值越大表示化合物對浮萍的毒性越強，反之越弱。

表1 化合物及其毒性活性值的負對數(pC)Table 1 Compounds and their negative logarithm of toxic activity (pC)

1.2 實驗方法

化合物結構參數化是建立化合物結構與毒性關系模型的關鍵一環，認為化合物中的非氫原子及非氫原子之間的關系對化合物的生物毒性有影響，而氫原子僅僅影響與其直接相連的非氫原子的取值。非氫原子的取值除與其直接相連的氫原子有關外，還與非氫原子自身的電子結構相關，首先在參考文獻[10-11]的基礎上，利用式(1)對化合物中的非氫原子進行賦值。

Zi=[(ni-hi)/(mi-ni-1)]1/2

(1)

式中：i為原子在分子中的編碼，ni為非氫原子i的價電子數，mi為原子核外電子總數，hi為與其直接連接的氫原子數。

不同類型的非氫原子對化合物的毒性值影響不盡相同，而相同類型的非氫原子對化合物的毒性值影響具有加和性，因而需要將化合物中的非氫原子進行分類。參閱文獻[12-14]方法將化合物中的非氫原子分為4類，與k(≤4)個非氫原子直接相連的非氫原子為第k類非氫原子，如與2個非氫原子直接相連的仲碳原子為第2類非氫原子，以此類推。不同類型的非氫原子自身對化合物生物毒性的影響按式(2)進行累加。

(2)

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

不同類型的非氫原子的關系對化合物的生物毒性的影響可能不同，相同類型的非氫原子之間的關系對化合物毒性的影響具有加和性。非氫原子之間的關系并非某種具體的相互作用，這種關系隨著非氫原子之間的距離的增加而減弱，隨著非氫原子的某種性質的增加而加強，式(3)可以滿足這一要求。

(3)

式中：Z按式(1)計算；dij為非氫原子i、j之間的相對距離(即鍵長之和與碳碳單鍵鍵長的比值，如果i、j之間有多條路徑，則以最短的為準)；n和l為原子所屬類型，α=0.5。4類非氫原子之間可以有10種不同的組合，即m11、m12、…、m44，簡寫為x5、x6、…、x14，m13表示第1類非氫原子與第3類非氫原子之間的關系，以此類推。這樣一來，對于一個化合物經參數化表達后最多可得14個變量(結構描述符)，暫將其稱為化合物非氫原子及其關系(non-hydrogen atoms and their relationship, NATR)。

運用多元線性回歸(MLR)及偏最小二乘回歸(PLS)進行建模，“留一法”對模型進行交互檢驗。

2 結果與討論(Results and discussion)

將化合物進行結構參數化表征后得到14個變量，但由于研究樣本中不具有第4類非氫原子，因而與第4類非氫原子相關的x4、x8、x11、x13、x14這5個變量全為“0”，其余9個非全“0”變量用于建模分析。

變量數達到9個，而樣本數僅有25個，在建模之前有必要對變量進行篩選。剔除與化合物生物毒性相關性不大的變量，尋找最佳變量組合進行建模。運用逐步回歸對變量進行篩選，依據變量的顯著性大小依次將變量引入模型。逐步回歸中建模相關系數(r)及標準偏差(SD)的變化情況如圖1和圖2所示。

由圖1可知，建模的r隨著變量的引入逐漸變大，起初增大的趨勢非常明顯，當變量數達到3個時，建模的r接近最大值，此后r略有增大，但增大趨勢不明顯。同樣，由圖2可知，SD隨著變量的引入逐漸變小，起初變小的趨勢非常明顯，當變量數達到3個時，SD接近最小值，此后SD略減小，但減小趨勢不明顯。綜上，應該選擇3個變量進行建模，此時符合樣本數(N)/變量數(n)≥5的經驗規則，三變量模型如式(4)計算。

圖1 相關系數(r)在逐步回歸中的變化情況Fig. 1 Changes of correlation coefficient (r) in stepwise regression

圖2 標準偏差(SD)在逐步回歸中的變化情況Fig. 2 Changes of standard deviation (SD) in stepwise regression

pC=6.385-3.853×x1+2.401×x3-0.293×x7

(4)

建模：N=25，r=0.984，SD=0.216，F=219.822；交互檢驗：N=25，rCV=0.979，SDCV=0.248，FCV=165.810。

建模的r高達0.984，說明模型具有良好的擬合能力；“留一法”交互檢驗相關系數(rCV)高達0.979，并且比建模的r略小，說明模型具有良好的預測能力與穩定性。判斷模型質量好壞的另一個指標是SD，優良的模型通常要求SD與數值范圍之比處于10%范圍內。本文所建模型的SD為0.216，而化合物毒性值最大值-最小值=6.20-1.80=4.40，0.216/4.40=4.91%，遠<10%，說明模型對化合物的生物毒性值預測誤差在可接受的范圍內。

為進一步分析影響化合物生物毒性的結構因素，采用偏最小二乘回歸(PLS)建模和分析。將化合物結構描述符作為自變量X，化合物誘發浮萍萎黃的活性參數pC作為因變量Y建立PLS模型(M2)。PLS模型(M2)的主成分數(A)為6，r為0.988，SD為0.174，rCV為0.980。從r、rCV及SD值來看，模型(M2)也具有良好的擬合效果、預測能力和穩定性，并且建模效果還優于模型(M1)，可以用于預測取代苯酚類化合物的毒性值和分析影響化合物毒性的結構因素。

25個樣本在PLS前2個主成分得分空間散點分布如圖3所示，由圖3可知，所有的樣本得分點都落在95%置信度的橢圓置信圈內，說明本文構建的化合物結構描述符能恰當表達該類化合物的結構特征，并且在PLS模型中得到反映。

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

變量重要性投影(VIP)可以反映各個變量與化合物的生物毒性相關性大小，通常認為VIP值>1的變量與所研究的性質相關性大，而其他變量與所研究的性質相關性可能較小。VIP如圖4所示，由圖4可知，x5、x3和x2的VIP值>1，說明這3個變量與化合物的毒性值相關性較大。進一步分析發現，VIP值最大的x5為第1類非氫原子之間的關系對化合物毒性值的影響，第1類非氫原子對應的是化合物取代基末端原子；x3為第3類非氫原子自身對化合物毒性值的影響，第3類非氫原子對應的大部分是苯環上被取代后的碳原子，x5和x3這2個變量都與化合物的取代基數量有關。綜上，化合物取代基越多，該化合物的生物毒性值可能就越大，這與表1中的數據特征基本相符。化合物產生毒性通常需要透過生物脂質膜，而化合物的體積與疏水性與其透過脂質膜的能力密切相關，當化合物分子都比較小時，化合物的疏水性起主導作用。本文研究樣本中的取代基都為疏水性基團，取代基越多，該化合物的疏水性就越強，越易透過脂質膜而產生毒性。

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

2個模型對化合物的毒性值進行了預測，預測值分別為表1的Cal.1和Cal.2，Err.1和Err.2分別為其誤差。模型預測值與實驗值相關圖如圖5所示。

由圖5可知，所有的樣本點都落在正方形的45°對角線上及附近，說明模型預測值與實驗值高度相關，兩者之間相差性小。2個模型對化合物毒性值預測誤差的絕對值分布如圖6所示。由圖6可知，絕大部分樣本的預測誤差的絕對值都小于模型的2倍標準偏差(2SD)。模型(M1)中，有2個樣本的預測誤差的絕對值超出了2SD的范圍，而模型(M2)中所有樣本的預測誤差的絕對值均處于2SD范圍內，說明模型(M2)對化合物毒性值的預測比模型(M1)更為準確，這與上述分析結果一致。

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

圖6 模型對樣本毒性的預測誤差注：SD表示標準偏差。Fig. 6 Model prediction error of samples’ toxicityNote: SD stands for standard deviation.

文獻[15]采用DFT-B3LYP計算描述符，運用多元線性回歸(MLR)進行了建模研究，將文獻結果與本文結果進行了對比(表2)。

由表2可知，本文所建模型的r及rCV明顯大于文獻結果，而SD明顯小于文獻結果，說明本文所建模型優于文獻結果。另外，文獻采用的分子結構表征方法DFT-B3LYP需經過分子結構優化，結算量大、復雜，而本文構建的結構描述符是基于分子二維結構計算得到，無需結構優化，計算量小，簡單易懂。

表2 模型比較Table 2 Model comparison

以化合物中的非氫原子及非氫原子之間的關系構建結構描述符，對部分取代苯酚類化合物結構進行了參數化表征。通過MLR和PLS方法建立了化合物結構與毒性關系模型，模型具有良好的擬合能力與預測能力。分析了影響化合物毒性的結構參數，并將結果與文獻進行了對比，化合物苯環上的取代基越多，該化合物可能就有較大的毒性值，本文取得的結果明顯優于文獻結果。本文所構建的化合物結構描述符簡單、易懂、計算方便，對于環境中的有機污染物的結構與生物毒性關系研究具有一定的參考價值。但結構描述符還存在不能區分順反異構、光學異構等現象，這些在后續研究中都應予以考慮。

◆