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

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

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

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

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

1 建模數(shù)據(jù)庫及方法

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

表1 硝基芳烴化合物結(jié)構(gòu)及其毒性

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

表2 硝基芳烴化合物的量子化學(xué)參數(shù)和拓撲指數(shù)

2 方法

2.1 人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

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

圖1 人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

建模過程主要分兩步:

1)隨機選擇75%的化合物作為訓(xùn)練集，25%的化合物作為測試集，利用訓(xùn)練數(shù)據(jù)集來訓(xùn)練模型，測試集來驗證模型，反復(fù)多次來確定最佳的網(wǎng)絡(luò)結(jié)構(gòu)[12];

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

2.2 偏微分法

由該法可以獲得兩種結(jié)果:一是每個輸入變量的微小變化導(dǎo)致網(wǎng)絡(luò)輸出變化的偏微分圖;二是每個輸入變量相對輸出的貢獻大小排序.

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

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

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

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

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

2.3 擾動法

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

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

2.4 權(quán)重方法

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

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

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

2.5 輪廓圖法

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

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

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

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

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

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

3 結(jié)果與討論

3.1 ANN－QSAR模型的預(yù)測能力

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

3.2 不同方法選擇變量結(jié)果及其比較

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網(wǎng)絡(luò)輸出lg1/cL50對變量QNO2的偏微分圖

3.2.2 輪廓圖法

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

3.2.3 權(quán)重法和擾動法

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

圖3 12個參數(shù)變量對網(wǎng)絡(luò)輸出lg1/cL50的輪廓圖

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

3.2.4 逐步回歸法

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

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

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

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

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

圖4 不同方法獲得的12個結(jié)構(gòu)參數(shù)對ANN輸出的相對貢獻分布圖

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

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

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

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

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

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

4 結(jié)論

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

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]陳國華，陸瑤，陳虹.基于逐步回歸所得變量集的遺傳反向傳播神經(jīng)網(wǎng)絡(luò)的QSAR研究[J].計算機與應(yīng)用化學(xué)，2010，27(9):1257 －1262.

[6]杜雨靜，范英芳.人工神經(jīng)網(wǎng)絡(luò)用于三苯基丙烯腈衍生物的定量結(jié)構(gòu) －活性關(guān)系模型[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]李鳴建，馮長君.取代苯甲酸對植物生長調(diào)節(jié)活性的拓撲QSAR研究[J].哈爾濱工業(yè)大學(xué)學(xué)報，2009，41(5):195－197.

[9]陳炫，聶長明，蔣司同，等.量子拓撲方法對硫醇的定量構(gòu)效關(guān)系研究[J].南華大學(xué)學(xué)報，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]于秀娟，王鵬，龍明策，等.有機化學(xué)品點價自相關(guān)拓撲指數(shù)與生物降解性的定量關(guān)系[J].環(huán)境科學(xué)學(xué)報，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)境科學(xué)與技術(shù)，2011，34(2):17 －21.