分子表面隨機采樣分析用于苯乙基噻唑硫脲衍生物抗艾滋病藥物的3D－QSAR研究

仝建波，趙 翔，鐘 黎，常 佳，李夢龍

（1．陜西科技大學化學與化工學院，西安710021；2．教育部輕化工助劑化學與技術重點實驗室，西安710021；3．四川大學化學學院，成都610064）

1 引 言

由人類免疫缺陷病毒（HIV）引起的艾滋病（即獲得性免疫缺陷綜合癥，AIDS）自1981年被發現以來，在全球迅速蔓延，其中速度最快的地區是非洲和亞洲，中國也屬于艾滋病高發地區．艾滋病的迅速蔓延引起了各國政府和科研機構的高度重視，紛紛投入大量的人力物力進行艾滋病的預防和治療研究．自1987年第一個抗艾滋病藥物齊多夫定（AZT）問世以來，已經有30余種不同類型的藥物問世．

非核苷類逆轉錄酶抑制劑（NNRTIs）是一類與核苷結構、作用機理迥異的特異性抑制HIV－1 RT的化合物，該類抑制劑具有結構多樣、高效低毒以及可與其他藥物協同作用等特性，而且與核苷類逆轉錄酶相比，還具有許多獨特的優點，比如它們的活性形式不需要細胞內的磷酸化，因而在靜止和活化的細胞或不同的細胞系中都具有廣泛的活性；它們亦可與細胞外的逆轉錄酶結合，減少游離病毒顆粒的逆轉錄，從而降低病毒的傳染性[1－2]．因而對NNRTIs的研究一直是發現新型抗艾滋病藥物的有效途徑和重要方向之一．

噻唑硫脲化合物N－（2－苯乙基）－N｀－（2－噻唑基）－硫脲（PETT）及其衍生物[2]作為一種典型的HIV NNRTIs，如果能通過進一步結構優化提高活性，有望成為一類新的高效抗艾滋病藥物．而定量結構活性關系（QSAR）研究是藥物設計的一種重要方法，它對于設計和篩選生物活性顯著的藥物以及闡述藥物的作用機理等具有指導作用[3－5]．因而構建這些化合物分子結構與生物活性之間定量相關模型，對于研究、設計和開發出高效抗艾滋病藥物具有重要意義．

QSAR為藥物設計和篩選提供了相對簡便的途徑，并受到了廣泛的重視．其中分子結構表征（MSC）是QSAR的一個關鍵環節，當前流行的MSC技術大致可以分為基于二維結構的2D結構描述子和基于三維空間構型的3D結構描述子兩大類．由于2D描述子很難再現配基與受體的空間契合方式，因此3D結構描述子正成為QSAR研究的主流．采用本實驗室新提出的一種3D分子結構描述方法—分子表面隨機采樣分析（RASMS），得到用于表征藥物分子結構特征的三維矢量描述子，將其用于PETT類抗艾滋病藥物的QSAR研究，取得了較好的結果．

2 原理和方法

2．1 RASMS中的幾個概念

2．1．1 蛋白質受體原子探針

由于大多數藥物直接作用靶標是具有一定生物活性功能的肽和蛋白質分子，因此RASMS法使用20種標準天然氨基酸中的各類原子作為探針．考慮到具有不同雜化狀態的原子往往因所處基團和區域不同而使其對活性貢獻表現出一定的差異，進而以此得到8個探針原子．為了反映這些探針的性質，本文分別給其賦予了平均電荷指數、范德瓦耳斯指數和平均疏水指數的概念．

平均電荷指數（Mean charge index，MCI）：探針原子電性特征取其出現于氨基酸中的平均電荷數．具體計算如下：首先使用Chemoffice 8．0自帶數據庫生成20個天然氨基酸的初始分子立體結構；并利用分子模擬軟件HyperChem 7．5進行分子力學構象優化；其結果進一步應用Gaussian 98W量子化學計算軟件在密度泛函水平基于廣義梯度近似法最終優化得到分子三維結構，并采用Mulliken布居分析法以單點形式計算出原子的凈電荷數量[6－8]；自編程序計算每種探針原子在氨基酸中出現的平均電荷數作為其MCI值．

范德瓦耳斯指數（Van der Waals index，VWI）：通常文獻給出的是孤立原子Van der Waals半徑，但實際原子由于其所處分子的化學微環境和自身雜化狀態不同該半徑有所變化，因此本文對其進行了調整，即使用經校正后的原子Van der Waals半徑作為探針原子的范德瓦耳斯指數[9－10]，VWI＝Ch×RVDW＊．

平均疏水指數（Mean hydrophobic index，MHI）：類似于MCI，MHI取每類探針的疏水性在天然氨基酸中出現的平均值．這里以文獻[11]所定義的原子溶解參數（atomic solvation parameter，ASP）作為疏水性量度．

2．1．2 虛擬受體可及表面

本文提出了虛擬受體可及表面（pseudo－receptor accessible surface，PRAS）的概念．當作為藥物作用靶點的生物分子（蛋白質、核酸、糖類等）中所包含的原子可以抵達的該藥物分子表面稱為分子虛擬受體可及表面（pseudo－receptor accessible surface of molecule，PRASM），如果以上述劃分蛋白質8類探針原子中的氫作為受體探針，其定義如下：利用單個氫原子球體在藥物分子Van der Waals表面滾動其球心所經歷的曲面則稱為虛擬受體氫原子可及表面（H－PRASM）．同理可以計算其余7種虛擬受體探針原子可及表面（圖1）．依照上述PRASM的計算方法可以定義孤立原子的原子虛擬受體可及表面（pseudo－receptor accessible surface of atom，PRASA），顯然該表面是一個球面（圖2），其半徑為該原子與探針半徑之和．可以看到藥物分子中每個原子的PRASA可能有一部分參與形成該分子的PRASM．

表1 8個探針原子在20種天然氨基酸中出現的頻數及其MCI，VWI和MHI取值Table 1 The appearing frequencies of eight types of atomic probes in twenty types of natural amino acids as well as the values of its MCI，VWI and MHI

圖1 分子虛擬受體可及表面Fig．1 Pseudo-receptor accessible surface of molecule

圖2 原子虛擬受體可及表面Fig．2 Pseudo-receptor accessible surface of atom

2．2 藥物分子中常見原子分類

有機化合物中常見的原子包括氫、碳、氮、磷、氧、硫、氟、氯、溴、碘，它們分別屬于元素周期表的IA、IVA、VA、VIA、VIIA共計5個主族．基于“具有相同化學性質的原子應屬于同一類”的觀點，本文很自然的按照這些原子所處元素周期表的族將其劃分為5大類．為了更好的表現分子細微結構特征，并考慮到同一族原子處于不同雜化狀態時化學性質也有較大差異，繼而在上述分類的基礎上再進一步把不同族中的原子按其雜化狀態細分為10類（表2）．

表2 RASMS中按周期表的族和原子雜化類型劃分的10類原子Table 1 Ten types of atoms according to element periodic tables

2．3 探針原子與藥物配體的作用模式

經典藥學理論認為藥物在抵達受體并與之發生作用絕大多數都是暫時的、可逆的非鍵效應，其表現為靜電、立體、疏水、氫鍵、電荷轉移等多種因素．本文考慮到靜電、立體、疏水效果幾乎包含了大部分這類信息，故在RASMS中主要計算這3種作用類型．對于氫鍵、電荷轉移等可以看成是靜電和立體效應的特殊表現形式．

靜電作用（electrostatic interaction）是一類重要的非鍵效應，經典的點電荷作用方式服從Coulomb定理（式1）．其中rij是探針到配體原子間的Euclid距離；e為單位電荷電量1．6021892×10－19C；ε0為真空中的介電常數8．85418782×10－12C2／J·m；Z為配體原子凈電荷數；p和l分別為探針和配體原子所屬類型．

立體作用（steric interaction）是空間原子間存在的非偶極－偶極或偶極誘導作用，這里采用Lennard－Jones方程來描述這種作用方式（式2）．該式中εij＝（εii·εjj）1／2為探針－受體原子勢能阱深，取文獻值[12]；D為經驗推導的原子間作用能修正常數，取0．01；Rij＊＝（VWIi＋Ch·Rj＊）／2，為經過校正后的探針－受體原子間van der Waals碰撞半徑，其校正因子Ch，當sp3雜化時取1．00，sp2雜化取0．95，sp雜化取0．90[13]．

疏水作用（hydrophobic interaction）是影響藥物分子與生物體結合的重要因素，由于其往往表現于體系熵的改變，因此很難用一個統一的公式來描述．對于有關疏水作用的研究已有眾多文獻報道，考慮到RASMS要求深入到配基分子內部原子與其表面受體探針相互作用，我們使用Kellogg等人提出的 HINT方法[14－18]來表達該類勢能形式．在HINT中定義了一個簡單的計算兩個原子之間的疏水相互作用表達式（式3），該式中S為原子溶劑可及面積（solvent accessible surface area of atom，SASA），是水分子在原子表面滾動其球心形成的表面面積[19]；a為原子疏水性常數，這里同樣使用上文所提到的ASP作為其表達值；T是作用形式的二值判別函數，以表明不同類型原子疏水作用的熵效應變化方向．

2．4 RASMS實現過程

使用Chemoffice 8．0分子圖形構建軟件包自動生成初始藥物分子立體結構；同時應用MOPAC 6．0[20－21]半經驗量子化學計算軟件在 AM1（Austin model 1）[22]水平上采用共扼梯度法完成幾何優化，并在最終分子結構上進行單點計算以求得每個原子所帶電荷的Mulliken布居數；繼而利用自編程序Sampling－tool進行RASMS分析：當以氫原子為探針時，首先在藥物分子虛擬受體氫原子可及表面（H－PRASM）上進行隨機取點作為該探針的探點，每一次采樣都計算該點探針與藥物分子中10類原子的相互作用情況，而每種作用又分為靜電、立體和疏水3種效應，這樣可以得到30個作用項，它包含了分子表面該點的勢場分布狀況．可以看到，經過大規模重復上述隨機采樣，探點將幾乎完全均勻的覆蓋分子周圍，從而得到整個表面的場能分布信息．完成所有采樣之后將多次探測得到的30項作用每項對應加和，并取其平均值作為該分子表面勢場分布密度．其中第一項表示HPRASM上藥物分子中第1類原子（氫原子）與探針平均作用情況，第二項表示第2類原子（sp3雜化的碳原子）與探針平均作用情況，以此類推．這樣通過使用8個探針對一個藥物分子表面進行采樣計算將得到8×30＝240個作用分量，即為RASMS法對該分子的最終分析結果（圖3）．

圖3 RASMS法對于每一個藥物分子得到的240個作用分量示意圖Fig．3 240types of interaction with RASMS method of every drug molecule

3 結果與討論

3．1 數據集的選取

61個PETT類衍生物的結構及活性數據來自文獻[23]（表3）．活性標度為化合物對HIV－1感染細胞的半數有效抑制濃度IC50的負對數pIC50：

3．2 模型驗證

為了深入研究RASMS方法與此類抗艾滋病藥物活性的內在聯系，采用多元線性回歸（MLR）這種典型的線性方法進行建模．另外對所建模型的外部預測能力和真實有效性進行驗證是定量構效關系中非常重要的一個部分，用來表征模型穩定性與預測能力的參數有rm、r0和k[24－27]等，檢測方法在文獻[24]中提到．為此，本文將數據集分成了51個內部樣本和10個外部樣本，10個外部樣本將用于模型檢驗．

3．3 結構表征

使用Chemoffice 8．0自動搭建61個PETT類抗艾滋病藥物分子的立體結構，用Chem3D自帶的MOPAC半經驗量子化學軟件在AM1水平上最終優化得到分子結構（截斷值0．0001kJ／mol），并采用Mulliken布居分析法以單點形式計算出原子的凈電荷數量，將分子中每個原子的空間位置及電荷分別以笛卡爾坐標和凈電子數目的形式輸入實驗室自編的應用程序Sampling－tool加以處理，得到分子的RASMS描述子．

表3 苯乙基噻唑硫脲衍生物母體結構及其pIC50Table 3 Structure and pIC50activity of phenylethylthiazolylthiourea derivatives

（continued）

3．4 多元線性回歸建模

RASMS含有240個描述子，所研究51個PETT內部樣本通過計算最終得到192個非零矢量，使用SPSS 18．0對這192個非零矢量進行逐步線性回歸（stepwise multiple regression，SMR），并采取有進有出的原則，SMR按Fisher顯著性檢驗依次引入變量，選出15個顯著性矢量；同時將SMR每一步得到的原始變量多元線性回歸（Multiple Linear Regression，MLR）建模，并使用留一法交互校驗法計算模型的QCV大小，以其達到最大值時變量數目來確定最終模型．表征模型穩定性與預測能力的參數有rm、r0和k，定義如下：

結合交互校驗的復相關系數（QCV）變化情況，最終確定選用11變量回歸模型，按照文獻[22]的方法分別輸入內部樣本和外部樣本前11個描述子和實驗值即可得到：Rcum、QCV、r0（test）、rm（test）、rm（all）和k分別為0．919、0．856、0．936、0．848、0．837和0．967（0．85＜k＜1．15），模型的標準化形式如下：

其中，E表示靜電作用，S表示立體作用，H表示疏水作用，如，E1－3表示第1類探針與第3類原子之間的靜電作用項，S1－10表示第1類探針與第10類原子之間的立體作用項，H8－9表示第8類探針與第9類原子之間的疏水作用項，以此類推．樣本活性實驗觀測值與模型計算相關情況如圖4，幾乎所有樣本都均勻分布于過原點45°直線周圍，無明顯異常點，結果表明RASMS能夠較好地表征香豆素類化合物的活性與結構之間的關系．

圖4 苯乙基噻唑硫脲衍生物活性實驗值與計算值關系圖Fig．4 Plot of calculated versus observed pIC50for PETT derivatives with MLR

4 結 論

本文所建模型包含了靜電、立體、疏水3種經典作用項，說明RASMS能夠較好地表征有機化合物活性與環境間復雜相互作用，所建模型以化合物本身的立體參數與電性參數為基礎．與當今大多3D藥物設計技術相比，RASMS在應用過程中所建的模型物化意義明確，可解釋性強，參數簡單易得．本文使用RASMS對61個PETT類化合物進行了建模和定量構效關系研究，結果表明RASMS描述子與分子生物活性呈線性相關，具有一定的指導意義，可應用于潛在和新型抗艾滋病藥物的設計和開發，值得進一步深入研究．

[1]Zhan P．Design，synthesis and anti－HIV evaluation of novel arylazolyl（azinyl）thioacetanilides as potent NNRTIs[D]．China Shandong University，2010（in Chinese）[展鵬．新型芳唑（嗪）巰乙酰胺類HIV－1 NNRTI的設計、合成與活性研究 [D]．中國：山東大學，2010]

[2]Tong J B，Wang P，Che T，etal．Studies of the quantitative structure activity relationship for phenylethylthiazolythiourea of anti-HIV drug using new three-dimensional structure descriptors[J]．J．At．Mol．Phys．，2012，29（6）：960（in Chinese）[仝建波，王平，車挺，等．苯乙基噻唑硫脲衍生物抗艾滋病藥物3D－QSAR研究[J]．原子與分子物理學報，2012，29（6）：960]

[3]Zhang X，Yang L Y，Zheng Y T．Advances of screening methods in vitro for HIV－1intergrase in-hibitors[J]．Chin．Pharmacol．Bull．，2013，29（1）：14（in Chinese）[張璇，楊柳萌，鄭永唐．HIV－1整合酶抑制劑體外篩選方法研究進展[J]．中國藥理學通報，2013，29（1）：14]

[4]Xu Y S，Zhong R G，Zeng C C，etal．The Research progress of diketo Acid HIV Integrase Inhibitors[J]．J．BeijingUniv．Technol．，2005，31（6）：635（in Chinese）[徐義生，鐘儒剛，曾程初，等．二酮酸類HIV整合酶抑制劑研究進展[J]．北京工業大學學報，2005，31（6）：635]

[5]Zhang N，Zhang H B，Huang S．The research progress of anti-HIV integrase inhibitors[J]．Pharm．Biotechnol．，2009，16（3）：269（in Chinese）[張楠，張惠斌，黃山．抗HIV整合酶抑制劑研究進展[J]．藥物生物技術，2009，16（3）：269]

[6]Hohenberg P，Kohn W．Inhomogeneous electron gas[J]．Phys．Rev．，1964，136：864．

[7]Kohn W，Sham L J．Self-consistent equations including exchange and correlation effects [J]．Phys．Rev．，1965，140：1133．

[8]Mulliken R S．Electronic population analysis on LCAO-MO molecular wave functions．I [J]．J．Chem．Phys．，1955，23：1833．

[9]Hahn M．Receptor surface models．1．definition and construction [J]．J．Med．Chem．，1995，38：2080．

[10]Bondi A．Van der Waals volumes and radii[J]．J．Chem．Phys．，1964，68：441．

[11]Pei J，Wang Q，Zhou J，etal．Estimating proteinligand binding free energy：Atomic solvation parameters for partition coefficient and solvation free energy calculation[J]．Proteins：Structure，Functionand Genetics，2004，57：651．

[12]Hasel W，Hendrikson T F，Still W C．A rapid approximation to the solvent accessible surface areas of atoms[J]．TetrahedronComp．Method，1988，1：103．

[13]Zhou P，Tong J B，Tian F F，etal．A novel comparative molecule／pseudo receptor interaction analysis[J]．Chin．Sci．Bull．，2006，51（15）：1824．

[14]Kellogg G E，Semus S F，Abraham D J．HINT-a new method of empirical hydrophobic field calculation for CoMFA [J]．J．Comput．AidedMol．Des．，1991，5：545．

[15]Wireko F C，Kellogg G E，Abraham D J．Allosteric modifiers of hemoglobin．2．crystallographically determined binding sites and hydrophobic binding／interaction analysis of novel hemoglobin oxygen effectors[J]．J．Med．Chem．，1991，34：758．

[16]Kellogg G E，Joshi G S，Abraham D J．New tools for modeling and understanding hydrophobicity and hydrophobic interactions[J]．Med．Chem．Res．，1992，1：444．

[17]Kellogg G E，Abraham D J．KEY，LOCK，and LOCKSMITH．Complementary hydrophobicity map predictions of drug structure from a known receptor／receptor structure from known drugs[J]．J．Mol．Graphics，1992，10：212．

[18]Nayak V R，Kellogg G E．Cyclodextrin-barbiturate inclusion complexes：a CoMFA／HINT 3－D QSAR study[J]．Med．Chem．Res．，1994，3：491．

[19]Hasel W，Hendrikson T F，Still W C．A rapid approximation to the solvent accessible surface areas of atoms[J]．TetrahedronComp．Method，1988，1：103．

[20]Stewart J J P．MOPAC：A semiempirical molecular orbital program [J]．J．Comput．AidedMol．Des．，1990，4（1）：1．

[21]Tong J B，Wang P，Che T，etal．Quantitative structure activity relationship studies of benzoxazinone of anti－HIV drug using new three－dimensional structure descriptors．[J]．J．At．Mol．Phys．，2012，29（5）：777（in Chinese）[仝建波，王平，車挺，等．苯并噁嗪酮衍生物類抗艾滋病藥物3D－QSAR研究[J]．原子與分子物理學報，2012，29（5）：777]

[22]Dewar M J S，Zoebisch E G，Healy E F，etal．AM1：A new general purpose quantum mechanical molecular model[J]．J．Am．Chem．Soc．，1985，107：3902．

[23]Razieh S，Afshin F，Behzad M．QSAR study of PETT derivatives as potent HIV－1reverse transcriptase inhibitors[J]．J．Mol．GraphicsModell．，2009，28：146．

[24]Roy K，Chakraborty P，Mitra I，etal．Some case studies on application of“r2m♂”metrics for quality of quantitative structure-activity relationship predictions：emphasis on scaling of response data[J]．J．Comput．Chem．，2013，34：1071．

[25]Mitra I，Saha A，Roy K．Exploring quantitative structure-activity relationship studies of antioxidant phenolic compounds obtained from traditional Chinese medicinal plants[J]．Mol．Simul．，2010，13：1067．

[26]Roy K，Mitra I，Kar S，etal．Comparative studies on some metrics for external validation of QSPR models[J]．J．Chem．Inf．Model．，2012，52：396．

[27]Roy K，Mitra I．On various metrics used for validation of predictive QSAR models with applications in virtual screening and focused library design [J]．Comb．Chem．HighThroughputScreening，2011，14：450．