楊 蕾，王 鵬，王 虹，3，蔣益林

(1.哈爾濱工業(yè)大學 基礎交叉與科學技術研究院，150001哈爾濱，leiyang84@vip.sina.com;2.哈爾濱工業(yè)大學市政環(huán)境工程學院，150090哈爾濱;3.黑龍江大學化學化工與材料學院，150080哈爾濱)

在環(huán)境化學領域，QSAR是進行有毒化學品生態(tài)風險評價的重要手段之一.目前，QSAR研究中常用的方法有多元線性回歸分析(MLR)、人工神經網絡(ANN)等，后者在處理非線性關系方面有著非常強大的功能，而化合物的結構和毒性之間大多存在非線性的關系，使ANN成為QSAR研究的熱點[1－3].

由輸入變量ANN便能預測輸出變量，但網絡內部的作用機理往往被忽略，因而被認為是個黑箱模型.近年來，研究者提出了許多方法來描述變量在神經網絡QSAR模型中所起的作用，然而大多數方法被用來消除不相關變量，因而被稱為修剪方法[4－6].簡單地說，修剪算法就是由具有高度聯(lián)結的網絡(i.e.神經元之間有許多連接)開始，逐步移除弱連接，或當連接移除時，網絡誤差無明顯變化的連接.事實上，在QSAR建模中，不僅需要好的預測能力，還要了解每個輸入變量對輸出的相對貢獻大小.本文回顧和比較了用ANN建模并解釋QSAR模型的6種方法，這些方法被用來確定每個輸入變量對輸出的貢獻，因而不是修剪算法.

以35種硝基化合物對黑呆頭魚96 h的生物毒性為研究對象，探討了ANN中引入變量選擇方法后，QSAR模型的解釋能力.結果表明，偏微分方法分析所建模型能得出最為全面準確的結果，模型具有良好的預測和解釋能力;其次為分布圖方法.擾動法和權重法對輸入參數能實現(xiàn)較好的分類，但過于簡化且方法不穩(wěn)定;傳統(tǒng)的逐步回歸法所得結果最差.

1 建模數據庫及方法

以文獻[7]報道的35種硝基芳烴化合物對黑呆頭魚的96 h半致死濃度cL50(mmol/L)數據作為研究對象，來討論和比較用于分析和解釋人工神經網絡QSAR模型的不同方法.該硝基芳烴主要由具有不同硝基取代基的甲苯、苯胺和苯酚、鹵代苯組成，具體結構和活性見表1.數據lg1/cL50見文獻[7].

表1 硝基芳烴化合物結構及其毒性

在QSAR研究中，用于描述化合物的結構參數有很多，包括拓撲的、量子的、實驗值等[7－9].本文在 Hall[7]和黃慶國等[10]研究硝基芳烴類化合物的基礎上，利用HyperChem 6.03軟件和自編的C++軟件分別計算了7種量子化學參數和5種自相關拓撲指數來表征化合物的結構，具體見表2(為表述方便，以下選擇變量方法中都以表中的編號來代替變量).

表2 硝基芳烴化合物的量子化學參數和拓撲指數

2 方法

2.1 人工神經網絡結構

目前，在QSAR建模中，多層前饋性人工神經網絡結構得到了最為廣泛的應用，其依據誤差反向傳播算法訓練而得.本文采用較為普遍的3層網絡結構(其中輸入層12個神經元，隱含層5個神經元和輸出層1個神經元)，具體結構見圖1.

圖1 人工神經網絡結構

建模過程主要分兩步:

1)隨機選擇75%的化合物作為訓練集，25%的化合物作為測試集，利用訓練數據集來訓練模型，測試集來驗證模型，反復多次來確定最佳的網絡結構[12];

2)在整個數據集上，利用第一步所獲得的網絡最佳結構進行QSAR建模，采用不同方法研究輸入變量對網絡輸出，即化合物的生物毒性的相對貢獻大小，并分析解釋QSAR模型.

2.2 偏微分法

由該法可以獲得兩種結果:一是每個輸入變量的微小變化導致網絡輸出變化的偏微分圖;二是每個輸入變量相對輸出的貢獻大小排序.

為了獲得輸入變量的微小變化導致輸出變化的偏微分圖，計算輸出對輸入變量的偏微分.以具有ni個輸入節(jié)點、nh個隱含節(jié)點和1個輸出節(jié)點，第j個樣本的輸出yj關于輸入xj(j=1…N，N為樣本總數)的偏微分為

其中:Sj為輸出神經元對其輸入的偏微分;Ihj為第h個神經元的響應;who和wih分別為輸出與第h個隱含層神經元、第i個輸入神經元與第h個隱含神經元之間的連接權重.

由式(1)可獲得一系列輸出相對輸入變量的偏微分圖，能直接評價每個輸入對輸出的影響.例如，偏微分為負，對于研究變量，輸出隨輸入的增大而減小.

另外，對于整個數據集，由偏微分方法可得到ANN輸出對每個輸入變量的相對貢獻大小，具體計算如下:

其中:SSDi為第i個變量對所有化合物毒性網絡輸出的偏微分平方之和，SSD值越大，其對網絡輸出，即對化合物毒性的影響最大.

2.3 擾動法

該法旨在評價每個輸入的微小變化對ANN輸出的影響.算法首先調整一個變量的值，而保持其他變量不變，同時記下每個輸入對輸出的響應.輸入變量變化對輸出影響最大的變量被視為對網絡輸出影響最大[13]，為最重要的變量.

基本思想如下:假定xi為選定變量，δ為變化量，則xi的變化可表示為xi=xi+δ.一般δ可選定變量值的10%～50%不等，這樣便可獲得按重要性排序的輸入變量分類.

2.4 權重方法

該法通過分割網絡連接權重來確定輸入變量的相對重要性，是由 Garson[14]首先提出的.方法主要涉及兩部分:一是按隱含層節(jié)點分割隱含－輸出層間連接權重;二是按輸入層節(jié)點劃分輸入－隱含層間連接權重.本文對此方法進行了簡化，而所得結果一致，具體如下:

1)對于隱含神經元h，用其輸入－隱含層連接權重的絕對值除以所有輸入－隱含層間的連接權重絕對值之和，即

2)對于輸入神經元i，用每個隱含神經元與其連接的輸入所獲得的Qih之和，除以所有隱含神經元與其連接的輸入所獲得的Qih之和，再乘以100便可獲得每個輸入變量對所有樣本輸出，權重分布貢獻的相對重要性(Relative Importance)，即

2.5 輪廓圖法

Lek等[15]首先提出了該方法，其主要思想是構建隸屬于所有輸入變量范圍的假定矩陣，同一時刻固定其他變量的值，在假定矩陣范圍內連續(xù)變化某個輸入變量來觀察網絡輸出的變化.詳細地說，就是每個輸入變量在區(qū)間范圍內被分成等間距的一系列值，該間距被稱為標度.其他變量被依次固定在不同倍數的標度上，一般取5個點，分別是最小值、1/4區(qū)間、1/2區(qū)間、3/4區(qū)間和最大值上.對于每個輸入變量，根據不同的取值，便可得到輸出變量的分布圖.由分布曲線圖(見圖3)可以直觀地看到隨著輸入變量的遞增，網絡輸出變量的變化趨勢和垂直波動范圍，波動范圍越大，表明該變量越重要.

本文在利用輪廓圖法研究輸入變量對輸出貢獻的過程中，分別將輸入變量的最大值和最小值區(qū)間范圍分成12、24、48、96、144 和 192 標度.圖 3代表了24標度的輪廓圖.事實上，不管采用什么標度，該方法相當穩(wěn)定，不同標度下變量的輪廓圖具有相似的形狀，唯一不同的是標度越大，變量的輪廓圖越精細.

2.6 傳統(tǒng)逐步回歸法

該法主要包括逐步地增加或刪除一個輸入變量來考察對輸出結果的影響，根據網絡輸出均方差(MSE)的變化，輸入變量便能按照重要性進行排序.例如在逐步減少輸入變量個數的過程中，引起均方差最大程度增大的變量，便是問題空間最重要的變量[16];反之，在逐步增加輸入變量的過程中，引起均方差最大程度減小的變量，便是問題空間最重要的變量.本文利用這兩種逐步回歸建模方法來評價12個輸入變量的影響，分別可以獲得變量之間的相對重要性排序.

1)前進法:首先產生12個模型，每個模型僅包含一個輸入變量，產生最小誤差的變量x最為重要，并參與下一步建模;接著產生11個模型，每個模型由x和剩余變量中的任意一個組成，這個過程反復進行，直到所有的變量都進入模型.網絡模型中輸入變量的出現(xiàn)排序即為它們對網絡輸出的相對重要性關系;

2)后退法:首先產生12個模型，每個模型由11個變量組成，如果模型中不包含變量x引起網絡輸出的最大誤差，則該變量x最為重要;接著產生11個模型，每個ANN模型由10個輸入變量組成.這個過程反復進行，直到模型中剩下一個變量為止.網絡中輸入變量的消除順序即為它們對網絡輸出的重要性排序.

3 結果與討論

3.1 ANN－QSAR模型的預測能力

由2.1給出的建模過程，最終確定最佳網絡結構為12－5－1(見圖1).對于化合物學習樣本集，步驟 1(見 2.1)的結果為 R2=0.923(P＜0.01);對于測試樣本集，其結果為 R2=0.930(P＜0.01).這表明該網絡結構可以應用于步驟2(見2.1)，即分析結構參數對所有化合物毒性的相對重要性.所有樣本參與建模，結果為R2=0.938(P＜0.01)，驗證了該網絡模型的預測能力.

3.2 不同方法選擇變量結果及其比較

3.2.1 偏微分方法

由偏微分方法可以獲得一系列輸入變量對輸出的偏微分圖.圖2給出了QNO2對硝基芳烴化合物毒性lg1/cL50的偏微分圖.可以看出，其偏微分值都為正，且隨著QNO2的增大，偏微分接近于0，表明隨著QNO2的增大，lg1/cL50也隨之增大并最終達到一個穩(wěn)定值，類似情況的還有變量QC、QN、FH、－ ELUMO、A2，其中 QN無明顯的規(guī)律性;此外變量 μ、－ EHOMO、A1、B1、B2、C1對硝基芳烴化合物毒性lg1/CL50的偏微分值大多為負，其中μ、B1無明顯的規(guī)律性.

圖2 ANN網絡輸出lg1/cL50對變量QNO2的偏微分圖

3.2.2 輪廓圖法

圖3代表了24標度的輪廓圖，可以看出，網絡輸出lg1/cL50隨QNO2、B2和 － EHOMO的增大有明顯變化，其中QNO2在整個取值范圍內對網絡輸出影響最大，是最重要的變量.另外，lg1/cL50隨QNO2的增大而增大，且當QNO2增大到一定程度時，lg1/cL50保持恒定，而B2和－EHOMO的增大將導致lg1/cL50減小，這與偏微分方法所得結果一致.由輪廓圖法獲得的變量間的相對重要性關系見表3.

3.2.3 權重法和擾動法

圖4(a)給出了由偏微分圖得到的輸入變量對輸出的相對貢獻圖，可以看出，該方法非常的穩(wěn)定且有較小的置信度區(qū)間，QNO2是化合物結構變量中對生物毒性貢獻最大的變量(45.2%)，其次是－ELUMO(16.1%)和B2(11.3%).由權重方法獲得輸入變量對輸出的相對貢獻率見圖4(b).與偏微分方法比較，其置信度區(qū)間更大，因而穩(wěn)定性較差.由圖可見，QNO2變量對網絡輸出的貢獻率最大，其次為 － ELUMO、QC、FH和 B2，而其他變量貢獻率差異不大.

圖3 12個參數變量對網絡輸出lg1/cL50的輪廓圖

圖4(c)給出了利用擾動法(δ=50%)獲得的輸入對輸出的相對貢獻分布圖.由圖可見，QNO2變量對網絡輸出的貢獻率最大，其次為 QC、－ELUMO和B2.該方法同樣也不夠穩(wěn)定，因為有些變量如QN和A2、－EHOMO和B1等之間對輸出的貢獻差異并不明顯.

3.2.4 逐步回歸法

逐步回歸方法分為前進法和后退法，獲得的變量之間的相對重要性排序結果見表3，可以看出，除了最重要的變量QNO2，兩種方法對變量重要性分類結果不盡相同.根據前進法，QNO2之后依次為FH、QC、－ EHOMO，而后退法依次為QC、－ ELUMO、QN.

3.3 ANN－QSAR模型的解釋性能力

本文采用12－5－1的3層神經網絡結構，對硝基芳烴對黑呆頭魚生物毒性進行QSAR建模，并將各種選擇變量方法作用于模型，來研究不同變量對網絡輸出，即生物毒性的相對重要性，進而來闡釋硝基芳烴化合物的作用機理，提高網絡模型的解釋能力.

早期的研究表明[17]，硝基芳烴化合物是一類重要的遺傳毒性化合物，其致毒機理為:苯環(huán)上硝基N原子的親電中心與生物組織中作為親核中心的DNA分子相互反應引起的.

結構參數QNO2是用來表征苯環(huán)上所有硝基中N原子的最大凈正電荷數，由表3可以看出，所有變量選擇方法都得出QNO2是影響化合物毒性的最重要參數，這正好驗證了文獻[17]所報道的該類化合物的致毒機理.

圖4 不同方法獲得的12個結構參數對ANN輸出的相對貢獻分布圖

表3 采用不同方法對輸入變量相對重要性分類結果

比較表3和圖4可以看出，－ELUMO、QC、B2是另外3個影響硝基芳烴毒性的重要結構參數.

其中－ELUMO表示分子最低未占用軌道能量，其值愈大，接受電子的能力越強，化合物對黑呆頭魚毒性也越大.這表明－ELUMO與生物毒性之間正相關，這與偏微分法和輪廓圖法所得的結論一致.可以認為，當化合物的親電中心與DNA分子的親核中心發(fā)生反應時，－ELUMO越大則化合物越容易接受電子發(fā)生反應，因而化合物的生物毒性越強.

QC代表苯環(huán)上與硝基相連的C原子的凈正電荷數，其值越大，則與之相鄰的硝基N原子親電中心越強，越容易與DNA分子反應，因而QC與化合物毒性之間正相關，這與偏微分法的結論一致.

自相關拓撲指數B2代表分子中間位原子電子信息總和，可以認為取代基電子相互作用同樣影響了DNA分子的反應活性.供電基團，如NH2、CH3、OH(見表1)可能離域了硝基上N原子的正電荷，提高了反應的活化能，因而化合物的毒性與B2負相關，這與偏微分和輪廓圖方法所得結論一致.

利用不同的變量選擇方法，剩余變量的相對貢獻大小排序有較大的出入，這主要是由方法的局限性引起的.不管怎樣，ANN中引入選擇方法有助于識別影響問題空間的主因子，提高模型的解釋能力.

4 結論

1)在ANN中引入不同的變量選擇方法，可大大增強QSAR模型的解釋能力，其中偏微分方法能得出最為全面準確的結果，其次為輪廓圖方法.擾動法和權重法對輸入參數能實現(xiàn)較好的分類，但過于簡化且方法不穩(wěn)定;而傳統(tǒng)的逐步回歸法結果最差.

2)硝基芳烴對黑呆頭魚毒性的QSAR模型中，選擇方法識別的重要變量包括QNO2、－ELUMO、QC和B2，它們能準確揭示化合物的致毒機理，從而證明了變量選擇方法的有效性.

[1]WU J H，MEI J，WEN S X，et al.A self－adaptive genetic algorithm－artificial neural network algorithm with leave－one－out cross validation for descriptor selection in QSAR study[J].Journal of Computational Chemistry，2010，31(10):1956 －1968.

[2]JAGDISH C P，ONKAR S.Artificial neural networksbased approach to design ARIs using QSAR for diabetes mellitus[J].Journal of Computational Chemistry，2009，30(15):2494－2508.

[3]JAGDISH C P，BOON H C.Artificial neural networkbased drug design for diabetes mellitus using flavonoids[J].Journal of Computational Chemistry，2011，32(4):555－567.

[4]APILAK W，CHANIN N，THANAKORN N，et al.Modeling the activity of furin inhibitors using artificial neural network[J].European Journal of Medicinal Chemistry，2008，44:1664 －1673.

[5]陳國華，陸瑤，陳虹.基于逐步回歸所得變量集的遺傳反向傳播神經網絡的QSAR研究[J].計算機與應用化學，2010，27(9):1257 －1262.

[6]杜雨靜，范英芳.人工神經網絡用于三苯基丙烯腈衍生物的定量結構 －活性關系模型[J].化工進展，2010，29(1):25 －28.

[7]KIER L B，HALL L H.Molecular connectivity in structure － activity analysis[M].[S.l.]:Research Studies Press，1987:232 －256.

[8]李鳴建，馮長君.取代苯甲酸對植物生長調節(jié)活性的拓撲QSAR研究[J].哈爾濱工業(yè)大學學報，2009，41(5):195－197.

[9]陳炫，聶長明，蔣司同，等.量子拓撲方法對硫醇的定量構效關系研究[J].南華大學學報，2009，23(4):84－87.

[10]HUANG Qingguo，LIU Yongbin.Genotoxicity of substituted nitro benzenes and the quantitative structureactivity relationship[J].Journal of Environmental Sciences，1996，8:103 －109.

[11]于秀娟，王鵬，龍明策，等.有機化學品點價自相關拓撲指數與生物降解性的定量關系[J].環(huán)境科學學報，2000，20(增刊):93 －96.

[12]GEMAN S，BIENENSTOCK E，DOURSAT R.Neural networks and the bias/valance dilemma[J].Neural Computation，1992，4(3):51 －58.

[13]SCARDI M，HARDING L W.Developing an empirical model of phytoplankton primary production:a neural network case study[J].Ecological Modelling，1999，120(2/3):213－223.

[14]GARSON G D.Interpreting neural－network connection weight[J].Artificial Intelligence，2001，6(8):47 －51.

[15]LEK S，DELACOSTE M，BARAN P，et al.Application of neural networks to modelling nonlinear relationships in ecology[J].Ecological Modelling，1996，90(32):39－52.

[16]SUNG A H.Ranking importance of input parameters of neural networks[J].Expert Systems with Applications，1998，15(12):405 －411.

[17]沈洪艷，張國霞，劉寶友，等.地表水體中常見硝基芳烴對鯉魚的聯(lián)合毒性[J].環(huán)境科學與技術，2011，34(2):17 －21.