基于分子對接和QSAR方法預測B-Raf II型抑制劑活性

劉海春 盧 帥 冉 挺 張艷敏 徐金星 熊 瀟徐安陽 陸 濤,2,* 陳亞東,*

(1中國藥科大學理學院, 南京 211198; 2中國藥科大學, 天然藥物活性組分與藥效國家重點實驗室, 南京 210009)

基于分子對接和QSAR方法預測B-Raf II型抑制劑活性

劉海春1盧 帥1冉 挺1張艷敏1徐金星1熊 瀟1徐安陽1陸 濤1,2,*陳亞東1,*

(1中國藥科大學理學院, 南京 211198;2中國藥科大學, 天然藥物活性組分與藥效國家重點實驗室, 南京 210009)

B-Raf激酶在促分裂素原活化蛋白激酶(MAPK)信號轉(zhuǎn)導通路中起著重要作用, 已被確定為癌癥治療非常有吸引力的靶標. 新型高效B-Raf抑制劑的開發(fā)成為癌癥治療的一個熱門研究領(lǐng)域. 本文以結(jié)構(gòu)多樣的B-Raf II型抑制劑為研究對象, 聯(lián)合應用分子對接和定量構(gòu)效關(guān)系(QSAR)模型研究其定量構(gòu)效關(guān)系去探討抑制活性的起源. 兩個主題作為研究重點: 生物活性構(gòu)象和描述符. 首先對分子對接方法(Glide、Gold、LigandFit、Cdocker和Libdock)進行準確性評價, 后將研究的對象分子對接到B-Raf活性位點并獲得生物活性構(gòu)象. 基于準確的對接結(jié)果, 計算得到16個打分評價函數(shù)和21個能量描述符, 以此構(gòu)建定量構(gòu)效關(guān)系模型. QSAR結(jié)果表明模型具有高度精確的擬合和強的預測能力(模型同時探討了對抑制活性有重要影響的描述符, 結(jié)果表明打分評價函數(shù)(G_Score, -ECD, Dock_Score, PMF)與能量描述符(S(hb_ext), DE(int), Emodel)對抑制活性影響非常大. 通過虛擬篩選和QSAR模型理論預測, 一些新的具有潛在抑制活性的化合物作為B-Raf II型抑制劑被獲得. 上述信息對于進一步設(shè)計新穎高效的B-Raf II型抑制劑提供了有用的指導.

B-Raf II型抑制劑; 分子對接; 打分評價函數(shù); 對接能量描述符; 定量構(gòu)效關(guān)系模型

1 引 言

促分裂素原活化蛋白激酶(MAPK)信號轉(zhuǎn)導通路在細胞的生長、分化和增殖中起著重要作用.1Raf是MAPK信號轉(zhuǎn)導通路中一個重要的絲氨酸-蘇氨酸蛋白激酶家族, 其激活突變能夠促進細胞的轉(zhuǎn)化與增殖, 因此與腫瘤發(fā)生密切相關(guān). Raf激酶家族有3 個亞型: A-Raf、B-Raf和C-Raf (Raf-1). 與ARaf和C-Raf相比, B-Raf有更高的激酶活性. 目前已確定B-Raf被激活后的體細胞突變率在黑色素瘤中高達50%–70%, 卵巢癌中達35%, 甲狀腺癌中達30%, 結(jié)腸癌中達10%.2最常見的B-Raf突變發(fā)生在臨近保守的DFG (Asp-Phe-Gly) 基序的600位殘基,即谷氨酸取代了纈氨酸(V600E). B-RafV600E具有癌基因的特點, V600E突變的B-Raf的催化活性比野生型B-Raf (B-RafWT)提高了500倍.3–5

通過B-Raf激酶抑制劑與其蛋白結(jié)合的晶體復合物分析, 根據(jù)DFG序列及α-C螺旋位置的不同, BRaf激酶抑制劑可分為I型、II型和I1/2型.6目前已有4個B-Raf激酶抑制劑上市, 包括II型抑制劑Sorafenib和Regorafenib,7–10以及I1/2型抑制劑Vemurafenib和Dabrafenib11–16(圖1). 據(jù)報道I型和I1/2型抑制劑可以異常激活B-RafWT或者MAPK通路中Raf-1, 引起腫瘤的二次增長. 而II型抑制劑因能夠抑制多種激酶的活性, 包括Raf激酶、VEGFR-2、VEGFR-3、PDGF-β、C-KIT和Ret等, 可以逆轉(zhuǎn)這種耐藥性.17–19

圖1 已上市的B-Raf激酶抑制劑的結(jié)構(gòu)Fig.1 Structures of the approved B-Raf kinase inhibitors

近年計算機輔助藥物分子設(shè)計(CADD)在新藥開發(fā)中得到了越來越多的應用, 定量構(gòu)效關(guān)系(QSAR)和對接(DOCK)是CADD的兩個主要工具.分子對接方法是將小分子化合物放入大分子蛋白靶標活性口袋中, 預測小分子與靶標的結(jié)合模式,同時應用打分函數(shù)評估其與大分子的結(jié)合自由能.分子對接目前已相對成熟, 但絕非完美, 仍有大量的問題有待進一步研究. 目前打分函數(shù)只能粗略評估結(jié)合自由能及預測生物活性.20–23幸運的是, QSAR方法具有準確可靠評估結(jié)合自由能, 預測新的先導物生物活性的能力. QSAR通過理論計算已知結(jié)構(gòu)的先導物的物理化學參數(shù)進行統(tǒng)計分析, 找出結(jié)構(gòu)與其生物效應(如生物利用度、活性或毒性等)的定量關(guān)系, 建立結(jié)構(gòu)活性模型并以此指導藥物分子設(shè)計. QSAR方法應用日益廣泛, 具有較高的預測準確度, 但也遭遇到了不少問題. 目前選擇和設(shè)計合適的描述符參數(shù)是QSAR研究的一大熱點. 隨著計算能力的提高, 計算化學也將提供更為準確的參數(shù), 為進行準確的QSAR研究奠定基礎(chǔ).24,25

為準確預測分子結(jié)合自由能提高虛擬篩選打分排序, 本研究結(jié)合分子對接和QSAR方法的技術(shù)力量建立以對接產(chǎn)生的打分函數(shù)和能量描述子為自變量的QSAR模型. 先是分子對接準確性評價, 對Glide、Gold、Cdocker、LigandFit和Libdock26–28對接策略進行可靠性驗證. 在準確的對接基礎(chǔ)上, 對結(jié)構(gòu)多樣的B-Raf II型抑制劑進行對接建立正確的結(jié)合模式, 產(chǎn)生分子對接評價參數(shù), 包括16個打分評價函數(shù)和21個對接產(chǎn)生的能量描述符. 利用一元線性及多重線性(MLR)回歸方法對打分評價函數(shù)及能量描述符和實驗活性之間的相關(guān)性進行統(tǒng)計分析, 建立QSAR模型.29,30MLR-QSAR模型獲得了預測活性和實驗pIC50值之間的高相關(guān)性系數(shù)并具有良好的預測能力. 應用具有良好預測能力的MLRQSAR模型得到多個具有較好的預測活性且骨架新穎的化合物可作為潛在的B-Raf先導化合物. 分子對接聯(lián)合QSAR模擬方法為合理設(shè)計和優(yōu)化新的具有潛在活性的B-Raf II型抑制劑提供了新的策略.

2 實驗材料和方法

2.1 B-Raf晶體復合物結(jié)構(gòu)準備

從RCSB蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)庫(PDB)選擇最具代表性也是最早發(fā)表的B-Raf晶體復合物(PDB ID: 1UWH, 分辨率為0.295 nm, 配體為Sorafenib)作為對接模板.3復合物結(jié)構(gòu)用Discovery Studio (DS)28進行預處理, A鏈及其配體被選定保留為對接目的. 將Sorafenib從復合物中提取出來, 除去所有水分子, 使用Sybyl31將配體和受體分別加氫加電荷進行結(jié)構(gòu)能量最小化.

2.2 訓練集和測試集的準備

B-Raf抑制劑數(shù)據(jù)來源于Binding Database (http://www.bindingdb.org/bind/index.jsp), 收集有BRaf抑制活性的化合物, 刪除重復的I型和I1/2型結(jié)構(gòu)及IC50≥ 10 μmolL–1的化合物, 最終選擇70個具有明顯活性變化梯度(IC50最小2.5 nmolL–1, 最大8800nmolL–1)和結(jié)構(gòu)多樣性的II型化合物用于模型構(gòu)建和評價. 通常數(shù)據(jù)集隨機被分為訓練集和測試集.本文為提高QSAR模型可靠性和預測能力, 采用化學指紋圖譜相似度聚類分析方法合理劃分數(shù)據(jù)集.32,33化學指紋圖譜聚類分析是一種將數(shù)據(jù)集劃分為小區(qū)域的集群定義劃分方法. 通過化學指紋圖譜最大不同, 一些有代表性的分子被選為簇的中心,剩余分子被分配到最近的聚類中心找到相應的集群. 通過DS計算化學指紋圖譜, Tanimoto系數(shù)為描述符描述分子相似度, 同時考慮活性梯度選擇, 數(shù)據(jù)集被分為訓練集38分子和測試集32分子(表中帶星號化合物). 圖2顯示劃分的訓練集分子和測試集分子(藍色小球)在化學空間的投影, 進一步說明數(shù)據(jù)集劃分的合理性. 表1顯示所有化合物結(jié)構(gòu)和IC50值.

2.3 分子對接、打分函數(shù)和描述符

分子對接程序Glide、Gold、LigandFit、Cdocker和Libdock被應用, 16個打分函數(shù)和21個對接描述符被選擇進行QSAR研究. 16個打分函數(shù)包括LigandFit的十一個打分函數(shù)(Dock_Score、PMF、PLP1、PLP2、Jain、Ligscore-1、Ligscore-2、PMF04、Ludi1、Ludi2、Ludi3), Gold的二個打分函數(shù)(Glod_Score和Chem_Score), 以及Glide的G_Score、Cdcoker的-Cdocker Interaction Energy和Libdock的Libdockscore; 對接描述符分別來源于不同的對接打分功能, 十一個來源于G_Score, 四個來源于Gold_Score, 六個來源于Chem_Score (表2).

2.3.1 Glide對接

Glide是近似為完全系統(tǒng)搜索配體的構(gòu)象、取向和空間位置的對接軟件. Glide對接過程中, 蛋白質(zhì)的結(jié)構(gòu)是剛性的, 而小分子的構(gòu)象是柔性的. 蛋白準備模塊Protein Preparation Wizard在OPLS-2005分子力場下準備受體; 以晶體復合物中的配體位置為對接中心, Receptor Grid Generation生成格點文件; Ligpre準備配體及其抑制劑. Ligand Docking進行分子對接, 選擇標準精度(SP)柔性對接, 每個配體保留10個構(gòu)象, 其它參數(shù)保持默認值. 選取打分最優(yōu)的構(gòu)象用于進一步研究, G_Score和十一個對接描述符被獲得.

表1 B-Raf抑制劑的訓練集和測試集分子結(jié)構(gòu), 實驗活性值與QSAR模型(模型M1及M2)預測活性值(Pre.)及殘差(Res.)Table1 Molecular structures, actual activities vs predicted activities (Pre.) and residuals (Res.) of QSAR models (M1 and M2) for the training set and the test set of B-Raf inhibitors

*compounds in the test set

2.3.2 Gold對接

Gold (Genetic Optimisat ion for Ligand Docking)程序是應用遺傳算法在局部柔性蛋白中尋找柔性配體構(gòu)象的自動柔性對接程序. 可以充分搜索構(gòu)象空間, 找到相互作用的最優(yōu)模式. 本文采用Gold_ Score和Chem_Score打分函數(shù)得到分子相互作用的擬合度得分(Fitness score, 函數(shù)值大小用來判斷分子對接情況的好壞, 函數(shù)值越大則分子間相互作用越強, 形成的包合物能量越低, 越穩(wěn)定)數(shù)值, 同時計算得到多種對接能量描述符. 結(jié)合位點用共結(jié)晶配體坐標位置來定義. 各參數(shù)保持缺省值, 每個配體保留10個構(gòu)象.

表2 分子對接方法, 打分評分函數(shù)及生成的各種描述符Table2 Molecular docking methods, scoring functions, and the various descriptors used for generating models

2.3.3 Cdocker對接

Cdocker是基于CHARMm力場的柔性對接程序, 設(shè)定蛋白質(zhì)構(gòu)象保持不變, 配體構(gòu)象發(fā)生改變,計算配體分子和受體之間的非鍵相互作用, 以-Cdocker Interaction Energy(-ECD)表示受體和配體間結(jié)合自由能. 分子對接選取PDB晶體結(jié)構(gòu)1UWH為受體模型, 通過加氫及CHARMm力場對其進行優(yōu)化. 以配體原子來定義結(jié)合位點, 其范圍半徑都采用自動生成值, 將配體分子置入結(jié)合位點球中, 變更配體構(gòu)象, 進行多構(gòu)象分子對接. 將這些不同構(gòu)象的配體分子與蛋白質(zhì)視為一個整體復合物,采用CHARMm力場和退火模擬算法分別優(yōu)化這些復合物的能量. 每個化合物對接后都得到10個構(gòu)象,綜合考慮其能量打分值(-ECD)及結(jié)合狀態(tài)選取最優(yōu)構(gòu)象進行分析.

2.3.4 LigandFit 對接

LigandFit對接具備了受體分子活性位點的自動尋找和確認、構(gòu)象柔性的多配體對接以及基于力場的相互作用打分評價的功能. LigandFit設(shè)定活性位點內(nèi)關(guān)鍵的氨基酸上的基團相互作用作為篩選標準, 要求被對接的小分子發(fā)生相應的基團匹配,形成氫鍵, 疏水相互作用等等. 研究基于B-Raf晶體復合物中的配體Sorafenib定義對接結(jié)合位點, 應用蒙特卡洛算法對配體的構(gòu)象空間進行取樣, 以基于力場參數(shù)的方法進行受體和配體分子相互作用的評價. 使用Ligscore對虛擬篩選結(jié)果進一步評價, 包括LigScore-1、LigScore-2、PLP1、PLP2、Jain、PMF等在內(nèi)的多種不同算法的打分函數(shù).

2.3.5 LibDock 對接

LibDock是基于特征的對接(Site Feature Docking), 在將小分子放入到結(jié)合腔內(nèi)之前需對結(jié)合腔內(nèi)部的特征(熱點區(qū)域“hot spot”, 即小分子配體要滿足的化學特征)進行計算. LibDock預先生成化合物的低能量構(gòu)象, 對接過程中不再生成新的分子構(gòu)象. LibDock有自己獨特的評分系統(tǒng), 也可以使用PLP進行評分. 本研究以B-Raf晶體蛋白1UWH為受體, ATP結(jié)合位點為對接中心, 半徑為1.5 nm, 每次保留30個構(gòu)象, 配體構(gòu)象產(chǎn)生方法為best, Hotspots為100, 其他參數(shù)默認, 對接結(jié)果以Libdockscore排序.

2.4 線性回歸模型

本研究以不同的對接打分函數(shù)及對接產(chǎn)生的能量描述子作為獨立變量, 利用一元線性和多元線性回歸(MLR)算法構(gòu)建強大的QSAR模型. 綜合考慮對接打分函數(shù)和化合物與激酶的結(jié)合狀態(tài)選取其最優(yōu)構(gòu)象. 對化合物生物活性值與其最優(yōu)構(gòu)象對接打分函數(shù)及能量描述子先進行一元線性回歸, 分別得到一元線性回歸方程. 同時運用TSAR軟件34對打分函數(shù)及描述符進行MLR分析, 選出影響化合物抗癌活性的重要獨立變量, 建立MLR模型. 一些統(tǒng)計指標, 如復相關(guān)性系數(shù)(r2)、交互驗證相關(guān)系數(shù)檢驗值(F)、均方根誤差(S)、測試集復相關(guān)系數(shù)評估模型的質(zhì)量, 模型的擬合能力檢驗包括復相關(guān)系數(shù)和均方根誤差, 模型的內(nèi)部預測能力驗證采用留一交叉驗證法的CV), 模型的外部預測能力驗證通過系數(shù)來評估.35,36

表3 Top3對接構(gòu)象均方根偏差(RMSD)值驗證對接可靠性Table3 Root mean square deviation (RMSD) values of Top3 configurations obtained through docking for docking reliability validated

3 結(jié)果與討論

3.1 分子對接研究

3.1.1 對接可靠性驗證

本工作的主要目的是利用不同對接程序產(chǎn)生的評分函數(shù)及描述符產(chǎn)生基于結(jié)構(gòu)的QSAR模型,以對確定新的有效的B-Raf抑制劑起指導作用. 生成模型之前, 驗證分子對接方法的可靠性是非常重要和必要的. 采用B-Raf及其配體Sorafenib復合物的X射線結(jié)構(gòu)1UWH執(zhí)行此驗證. 通過Glide、Gold、LigandFit、Cdocker和Libdock將配體Sorafenib還原對接到激酶結(jié)合口袋, 驗證對接方法能否再現(xiàn)結(jié)合底物在晶體結(jié)構(gòu)中的構(gòu)象. 對接后的分子構(gòu)象與晶體結(jié)構(gòu)中配體分子構(gòu)象的均方根偏差(RMSD)作為評價標準, 排名前Top3最佳對接(相應最高打分或最低能量)構(gòu)象進行驗證考慮, 一般RMSD不超過0.20nm即認為對接還原成功. 可接受的對接結(jié)果(表3顯示RMSD在0.0292–0.2154 nm)顯示對接構(gòu)象在結(jié)合位點內(nèi)具有相似的結(jié)合位置和方向, 與晶體構(gòu)象類似(圖3顯示Top1對接結(jié)果), 表明上述對接方法能成功在B-Raf晶體結(jié)構(gòu)中還原再現(xiàn)結(jié)合底物的構(gòu)象.

3.1.2 B-Raf II型抑制劑對接結(jié)果

70個具有生物活性B-Raf II型抑制劑分子(表1)被選擇進行分子對接和特殊的QSAR模型研究.對所有化合物進行幾何優(yōu)化, 用上述對接方法對接到B-Raf (PDB 1UWH)晶體結(jié)構(gòu)模板的結(jié)合位點. 基于對接打分評價函數(shù)及同B-Raf激酶之間的相互作用確定化合物的最佳對接構(gòu)象. II型B-Raf激酶抑制劑作用模式表現(xiàn)為: (1)與DFG基序天冬氨酸(Asp)和α-C螺旋谷氨酸(Glu)形成一對關(guān)鍵性的氫鍵; (2)與鉸鏈區(qū)半胱氨酸(Cys)上的NH和羰基形成1–2個氫鍵; (3)與變構(gòu)位點附近的疏水口袋形成疏水作用; (4)鉸鏈區(qū)附近有一個溶劑可及區(qū)(見圖4).根據(jù)對接結(jié)果研究, 所有的化合物都能夠充分準確對接, 基本滿足上述作用模式, 整體取向一致, 展示對接模型良好的質(zhì)量.

圖3 Sorafenib共晶體構(gòu)象(綠色)與Top1最佳對接構(gòu)象(其他顏色)在B-Raf活性部位疊加圖Fig.3 Superimposition of co-crystallized Sorafenib (green) and redocked Top1 conformation by multiple docking (the other colors) in the active site of B-Raf

3.2 相關(guān)模型的計算和驗證

以Glide、Gold、LigandFit、Cdocker和Libdock對接方法計算得到的打分評價函數(shù)和能量描述子(表2)作為QSAR建模的自變量函數(shù), 抑制劑分子的–lgIC50(pIC50)作為因變量函數(shù)用來構(gòu)建對接打分函數(shù)/能量描述符和B-Raf II型抑制劑pIC50的相關(guān)性模型.

3.2.1 一元線性回歸模型

首先研究實驗活性與單一對接打分評價函數(shù)及其能量描述子的定量構(gòu)效關(guān)系, 建立一元線性回歸模型, 計算其相關(guān)性系數(shù)r2. 結(jié)果列于表4, r2在0.002到0.622之間. 對于單一打分評價函數(shù), G_Score給出最好的相關(guān)系數(shù)(r2= 0.622), Libdockscore給最差的相關(guān)系數(shù)(r2= 0.027). -ECD、Dock_Score、LigScore-1和LigScore-2等也具有相對較好的相關(guān)性系數(shù), r2在0.50左右. 這一結(jié)果表明, QSAR模型基于柔性對接的打分函數(shù)比剛性對接打分函數(shù)更能獲得高的相關(guān)系數(shù)r2. 因為通過柔性對接方法可以得到更準確的結(jié)合模式, 而基于正確的結(jié)合模式基礎(chǔ)上的活性預測被認為是合理的. 進一步說明對接方法和蛋白質(zhì)模板的驗證是必要的. 對于單一能量描述符, Emodel給出最好的相關(guān)系數(shù)(r2= 0.471), H(rot)給最差的相關(guān)系數(shù)(r2= 0.002). Rewards和ENERGY的r2分別為0.261和0.272, Lipo, Hbond, LE, S(vdw_int)的r2也高于0.20. 說明對接產(chǎn)生的能量描述符也可能對估計準確的活性值是有用的. 一元線性回歸模型顯示每一個獨立變量同實驗活性之間具有一定的相關(guān)性. 然而其相關(guān)性較低,大部分相關(guān)系數(shù)r2小于0.50 (表4). 因此, 單一使用對接打分函數(shù)/能量描述符對于確定有效的抑制劑和預測準確的抑制活性具有明顯的局限性.

圖4 B-Raf抑制劑作用模式示意圖Fig.4 Schematic diagram of interaction patterns of B-Raf kinase inhibitor

3.2.2 多重線性回歸模型

一元線性回歸模型的結(jié)果促使我們開發(fā)新的預測模型, 因此采用TSAR軟件, 通過MLR方法對訓練集分子對接得到的16個打分函數(shù)及21個對接描述符分別進行特定的QSAR研究, 得到相關(guān)的回歸方程M1和M2. 對回歸方程進行了測試和驗證, 以評估模型的質(zhì)量. 最終得到的多元線性回歸模型與一元線性回歸模型相比, 活性相關(guān)性及預測能力大大提高. M1和M2的相應統(tǒng)計參數(shù)如表5所示.

QSAR模型M1使用所有16個打分評價函數(shù)生成, 回歸方程:

在QSAR建模中, 訓練集樣本直接參與模型的建立, 模型的擬合能力由相關(guān)系數(shù)(r2)、標準偏差(SD或S)來評價. r2主要體現(xiàn)模型對訓練集樣本數(shù)據(jù)的擬合能力, 通常認為r2越大越接近1說明模型的相關(guān)性越好, 擬合效果越好, QSAR模型質(zhì)量越高. S表示數(shù)據(jù)與平均數(shù)據(jù)之間差值的平均數(shù). S的大小表明數(shù)據(jù)的分散性, 同樣表示模型的擬合能力, 越接近于零說明模型擬合效果越好. 一個模型對內(nèi)部訓練集樣本有較好的擬合能力并不代表模型具有較好的預測能力. QSAR模型最主要目的是對未知化合物進行活性預測, 基于這個意義層面考慮, 模型的預測能力顯得比擬合能力要重要. QSAR模型的穩(wěn)健性也是必要的, 否則模型受到異常值的影響會出現(xiàn)抖動或者模型自身會存在機會相關(guān)問題. 為檢驗所建立QSAR模型的穩(wěn)健性和內(nèi)部預測能力, 采用“留一法”(LOO) 計算交叉驗證(CV)系數(shù)評價模型穩(wěn)健性和預測能力. 當大于0.5時, 通常認為所建立的模型具有較好的預測能力, 越接近于1說明穩(wěn)定性越好, 預測能力越強.

表4 實驗pIC50值與評分函數(shù)及描述符一元線性回歸模型相關(guān)性系數(shù)Table4 Correlations of individual scoring function and descriptors with the experimental pIC50values

表5 QSAR模型M1和M2的相應統(tǒng)計參數(shù)Table5 QSAR modes M1, M2 and the corresponding statistics

綜合考慮模型的擬合能力、穩(wěn)健性及預測能力三個要素, Alexander Tropsha曾提出, 一個較好的QSAR模型需要滿足下面兩個條件: r2> 0.600,所得打分函數(shù)QSAR模型M1復相關(guān)系數(shù)r2、交互驗證相關(guān)系數(shù)分別為0.852和 0.790, 說明模型具有良好的擬合、穩(wěn)健性及內(nèi)部預測能力, 此B-Raf II型抑制劑QSAR模型M1是可靠的.

21個對接產(chǎn)生的表示配體和受體之間的靜電力、范德華力、氫鍵以內(nèi)部能量相互作用的描述符也被視為獨立的變量生成QSAR模型M2, 回歸方程:

模型對訓練集和測試集中化合物的活性分別進行預測, 活性實驗值和預測值以及殘差見表1, 線性關(guān)系如圖5和圖6所示. 表1顯示模型M1和M2訓練集化合物的預測值與實驗值絕大部分相差都很小,說明該模型有較強的擬合能力. 圖5(a)和圖6(a)所示為訓練集實驗值和預測值線性回歸關(guān)系圖, 由圖可見實驗值與預測值有較好的相關(guān)性, 進一步顯示模型具有較好的擬合及內(nèi)部預測能力. 同時表1顯示模型M1和M2測試集化合物預測值與實驗值相差也都很小, 說明該模型有較強的外部預測能力. 圖5(b)和圖6(b)所示為測試集的實驗值和預測值線性回歸關(guān)系圖, 得到同訓練集一樣的結(jié)果, 說明模型具有好的外部預測能力可用來準確預測新抑制劑的生物活性.

圖5 基于打分函數(shù)的QSAR模型M1實驗活性值與預測活性值關(guān)系圖Fig.5 Actual activity vs predicted activity using QSAR model M1 based on scoring functions

圖6 基于打分函數(shù)的QSAR模型M2實驗活性值與預測活性值關(guān)系圖Fig.6 Actual activity vs predicted activity using QSAR model M2 based on docking descriptors

3.3 模型參數(shù)解釋

模型M1包括4個獨立的自變量, X1是G_Score, X2是-ECD, X4是Dock_Score, X8是PMF. G_Score基于經(jīng)驗的打分函數(shù), 結(jié)合配體和受體之間的庫侖和范德華相互作用能, 重點考慮氫鍵相互作用, 同時引入溶劑化作用. G_Score負的回歸系數(shù)表明對親和力有不利的貢獻(G_Score打分為負值, 越小越有利). -ECD描述配體與受體的相互作用能, 包括范德華能與靜電能量. Dock_Score是基于力場的打分函數(shù), 包括配體的內(nèi)能以及配體與受體的相互作用能,主要考慮靜電和立體互補性對于結(jié)合能的貢獻.PMF是基于知識的打分函數(shù), 總結(jié)成對相互作用包括所有在受體-配體原子間相互作用連同金屬離子相互作用和鹵素的電極電位, 將其轉(zhuǎn)換成距離依賴性Helmholtz自由能. -ECD、Dock_Score和PMF正的回歸系數(shù)表明其與親和力正相關(guān), 對親和力有積極有利的貢獻. 模型M2中, 三個對接能量描述符S(hb_ext)、DE(int)、Emodel對活性預測有重要意義. S(hb_ext)為蛋白質(zhì)-配體的氫鍵相互作用能, 有助于提高回歸系數(shù)與其正相關(guān), 對親和力有有利的影響. DE(int)為內(nèi)部扭轉(zhuǎn)應變張力. Emodel是Glide對接得到的能量描述符,是對G_Score的拓展. Emodel結(jié)合了G_Score, 非鍵相互作用能和內(nèi)部應變能. DE(int)和Emodel描述符與回歸系數(shù)負相關(guān),對親和力有負貢獻, 數(shù)值增大導致親和力降低. QSAR模型M1和M2表明在QSAR研究中選擇正確有效的描述符對建立良好模型的重要性. 同時結(jié)果表明通過改變某些重要描述符的大小可增強抑制活性.

表6 苗頭化合物化學結(jié)構(gòu)及活性預測值Table6 Chemical structures and predicted activity values of the hit compounds

3.4 新化合物的篩選與活性預測

基于MLR模型, 采用分子對接對ChemDiv數(shù)據(jù)庫進行虛擬篩選, 用QSAR模型進行活性預測. 應用Lipinski五規(guī)則進行初始過濾以減少數(shù)據(jù)庫的樣本數(shù)量. 首先使用一元線性回歸得到最高相關(guān)系數(shù)的Glide進行對接虛擬篩選, 打分評價函數(shù)G_Score大于10.0且符合II抑制劑結(jié)合模式的化合物被選擇進行深入研究. 對這些化合物進一步進行Gold, Cdocker、LigandFit、LibDock對接, 以得到其他打分評價函數(shù)和能量描述子進行QSAR模型活性預測.應用MLR-QSAR模型預測進行活性預測, 最終8個具有較好的預測pIC50值(pIC50> 6.5)且骨架新穎的化合物被選擇作為潛在的先導化合物. 8個命中的苗頭化合物化學結(jié)構(gòu)和預測pIC50值列于表6, 仍有待進一步的實驗確認和評價.

4 結(jié) 論

研究以定量構(gòu)效關(guān)系和分子對接方法來預測B-Raf II型抑制劑的活性. 通過對接還原配體結(jié)晶構(gòu)象驗證分子對接方法(Glide、Gold、LigandFit、Cdocker和Libdock)的準確性. 結(jié)果表明, 對接方法可以預測正確的結(jié)合模式, 對接研究對象分子到BRaf的活性位點同樣可以獲得生物活性構(gòu)象. 基于正確的結(jié)合模式, 計算得到對接打分函數(shù)及其能量相關(guān)的描述符. 利用一元線性及其多重線性回歸方法對打分函數(shù)及能量描述符和實驗活性之間的相關(guān)性進行統(tǒng)計分析, 找出結(jié)構(gòu)與活性的定量關(guān)系, 建立QSAR模型. 一元線性模型顯示描述符與實驗活性的相關(guān)性較低. MLR-QSAR模型獲得了高相關(guān)性系數(shù)及良好的內(nèi)外部預測能力, 同時探討了對抑制活性有重要影響的描述符. 應用MLR-QSAR模型的8個具有較好的預測活性且骨架新穎的化合物被選擇作為潛在的先導化合物, 對于進一步設(shè)計新穎高效的B-Raf II型抑制劑提供了有用的指導. 分子對接聯(lián)合QSAR模擬方法為合理設(shè)計和優(yōu)化新的具有潛在活性的B-Raf II型抑制劑提供了新的策略. 該策略是一種很有前途的方法, 可用于新靶標的分子設(shè)計.

(1)Peyssonnaux, C.; Eychene, A. Biol. Cell 2001, 93, 53. doi: 10.1016/S0248-4900(01)01125-X

(2)Davies, H.; Bignell, G. R.; Cox, C.; Stephens, P.; Edkins, S.; Clegg, S.; Teague, J.; Woffendin, H.; Garnett, M. J.; Bottomley, W.; Davis, N.; Dicks, E.; Ewing, R.; Floyd, Y.; Gray, K.; Hall, S.; Hawes, R.; Hughes, J.; Kosmidou, V.; Menzies, A.; Mould, C.; Parker, A.; Stevens, C.; Watt, S.; Hooper, S.; Wilson, R.; Jayatilake, H.; Gusterson, B. A.; Cooper, C.; Shipley, J.; Hargrave, D.; Pritchard-Jones, K.; Maitland, N.; Chenevix-Trench, G.; Riggins, G. J.; Bigner, D. D.; Palmieri, G.; Cossu, A.; Flanagan, A.; Nicholson, A.; Ho, J. W. C.; Leung, S. Y.; Yuen, S. T.; Weber, B. L.; Seigler, H. F.; Darrow, T. L.; Paterson, H.; Marais, R.; Marshall, C. J.; Wooster, R.; Stratton, M. R.; Futreal, P. A. Nature 2002, 417, 949. doi: 10.1038/nature00766

(3)Wan, P. T.; Garnett, M. J.; Roe, S. M.; Lee, S.; Niculescu-Duvaz, D.; Good, V. M.; Jones, C. M.; Marshall, C. J.; Springer, C. J.; Barford, D.; Marais, R. Cell 2004, 116, 855. doi: 10.1016/S0092-8674(04)00215-6

(4)Samowitz, W. S.; Sweeney, C.; Herrick, J.; Albertsen, H.; Levin, T. R.; Murtaugh, M. A.; Wolff, R. K.; Slattery, M. L. Cancer Res. 2005, 65, 6063. doi: 10.1158/0008-5472.CAN-05-0404

(5)Riesco-Eizaguirre, G.; Gutiérrez-Martínez, P.; García-Cabezas, M. A.; Nistal, M.; Santisteban, P. Endocr-Relat. Cancer 2006, 13, 257. doi: 10.1677/erc.1.01119

(6)Wang, X.; Kim, J. J. Med. Chem. 2012, 55 (17), 7332. doi: 10.1021/jm300613w

(7)Wilhelm, S. M.; Carter, C.; Tang, L.; Wilkie, D.; McNabola, A.; Rong, H.; Chen, C.; Zhang, X.; Vincent, P.; McHugh, M.; Cao, Y.; Shujath, J.; Gawlak, S.; Eveleigh, D.; Rowley, B.; Liu, L.; Adnane, L.; Lynch, M.; Auclair, D.; Taylor, I.; Gedrich, R.; Voznesensky, A.; Riedl, B.; Post, L. E.; Bollag, G.; Trail, P. A. Cancer Res. 2004, 64, 7099. doi: 10.1158/0008-5472.CAN-04-1443

(8)Wilhelm, S.; Carter, C.; Lynch, M.; Lowinger, T.; Dumas, J.; Smith, R. A.; Schwartz, B.; Simantov, R.; Kelley, S. Nat. Rev. Drug Discov. 2006, 5 (10), 835. doi: 10.1038/nrd2130

(9)Wilhelm, S. M.; Dumas, J.; Adnane, L.; Lynch, M.; Carter, C. A.; Schütz, G.; Thierauch, K. H.; Zopf, D. Int. J. Cancer 2011, 129 (1), 245. doi: 10.1002/ijc.v129.1

(10)Strumberg, D.; Schultheis, B. Expert Opin. Investig. Drugs 2012, 21 (6), 879. doi: 10.1517/13543784.2012.684752

(11)Chapman, P. B.; Hauschild, A.; Robert, C.; Haanen, J. B.; Ascierto, P.; Larkin, J.; Dummer, R.; Garbe, C.; Testori, A.; Maio, M.; Hogg, D.; Lorigan, P.; Lebbe, C.; Jouary, T.; Schadendorf, D.; Ribas, A.; O'Day, S. J.; Sosman, J. A.; Kirkwood, J. M.; Eggermont, A. M. M.; Dreno, B.; Nolop, K.; Li, J.; Nelson, B.; Hou, J.; Lee, R. J.; Flaherty, K. T.; McArthur, G. A. N. Engl. J. Med. 2011, 364, 2507. doi: 10.1056/NEJMoa1103782

(12)Bollag, G.; Tsai, J.; Zhang, J.; Zhang, C.; Ibrahim, P.; Nolop, K.; Hirth, P. Nature. Rev. Drug. Discov. 2012, 11, 873. doi: 10.1038/nrd3847

(13)Rheault, T. R.; Stellwagen, J. C.; Adjabeng, G. M.; Hornberger, K. R.; Petro, K. G.; Waterson, A. G.; Dickerson, S. H.; Mook, R. A., Jr.; Laquerre, S. G.; King, A. J.; Rossanese, O. W.; Arnone, M. R.; Smitheman, K. N.; Kane-Carson, L. S.; Han, C.; Moorthy, G. S.; Moss, K. G.; Uehling, D. E. ACS Med. Chem. Lett. 2013, 4, 358. doi: 10.1021/ml4000063

(14)Kim, D. H.; Sim, T. Arch. Pharm. Res. 2012, 35, 605. doi: 10.1007/s12272-012-0403-5

(15)Holderfield, M.; Deuker, M. M.; McCormick, F.; McMahon, M. Nat. Rev. Cancer 2014, 14, 455. doi: 10.1038/nrc3760

(16)Ackerman, A.; Klein, O.; McDermott, D. F.; Wang, W.; Ibrahim,N.; Lawrence, D. P.; Gunturi, A.; Flaherty, K. T.; Hodi, F. S.; Kefford, R.; Menzies, A. M.; Atkins, M. B.; Long, G. V.; Sullivan, R. J. Cancer 2014, 120, 1695. doi: 10.1002/cncr.28620

(17)Hartsough, E.; Shao, Y.; Aplin, A. E. J. Invest. Dermatol. 2014, 134, 319. doi: 10.1038/jid.2013.358

(18)Bucheit A. D.; Davies M. A. Biochem. Pharmacol. 2014, 87 (3), 381. doi: 10.1016/j.bcp.2013.11.013

(19)Zhang, F.; Lu, T.; Tang, W. F. Prog. Pharm. Sci. 2014, 38, 31. [張 帆, 陸 濤, 唐偉方. 藥學進展, 2014, 38, 31.]

(20)Enyedy, I. J.; Egan, W. J. J. Comput. Aided Mol. Des. 2008, 22, 161. doi: 10.1007/s10822-007-9165-4

(21)Kroemer, R. T. Curr. Protein Pept. Sci. 2007, 8, 312. doi: 10.2174/138920307781369382

(22)Lin, J.; Li, Z. G.; Zou, J. W.; Lu, S. Y. Acta Chim. Sin. 2012, 11, 1309. [林 軍, 李祖光, 鄒建衛(wèi), 陸紹永. 化學學報, 2012, 11, 1309.]

(23)Li, X. D.; Hou, T. J.; Xu, X. J. Acta Phys. -Chim. Sin. 2005, 21, 504. [李旭東, 侯廷軍, 徐筱杰. 物理化學學報, 2005, 21, 504.] doi: 10.3866/PKU.WHXB20050509

(24)Li, M. Y.; Xia, L. Journal of China Pharmaceutical University 2003, 34 (6), 586. [李敏勇, 夏 霖. 中國藥科大學學報, 2003, 34 (6), 586.]

(25)Patel, H. M.; Noolvi, M. N.; Sharma, P.; Jaiswal, V.; Bansal, S.; Lohan, S.; Kumar, S. S.; Abbot, V.; Dhiman, S.; Bhardwaj, V. Med. Chem. Res. 2014, 23, 4991. doi: 10.1007/s00044-014-1072-3

(26)Maestro, Version 9.0; Schr?dinger, L. L. C.: New York, 2011

(27)GOLD, Version 4.0; Astex Technology: Cambridge, 2001

(28)Discovery Studio Client, Version 2.5; Accelrys Inc.: San Diego, 2008

(29)Chen, L.; Chen, X. J. Mol. Graph. Model. 2012, 33, 35. doi: 10.1016/j.jmgm.2011.11.003

(30)Fr?czek, T.; Siwek, A.; Paneth, P. J. Chem. Inf. Model. 2013, 53, 3326. doi: 10.1021/ci400427a

(31)SYBYL, Version 7.1; Tripos Inc.: St. Louis, Missouri, 63144, U. S. A.

(32)Hassan, M.; Bielawski, J.; Hempel, J.; Waldman, M. Mol. Divers. 1996, 2 (1–2), 64. doi: 10.1007/BF01718702

(33)Matter, H.; Potter, T. J. Chem. Inf. Comput. Sci. 1999, 39 (6), 1211. doi: 10.1021/ci980185h

(34)TSAR Software, Version 3. 3; Accelrys Inc.: Oxford, England

(35)Roy, K. Expert Opin. Drug Discov. 2007, 2, 1567. doi: 10.1517/edc.2007.2.issue-12

(36)Qin, L. T.; Liu, S. S.; Xiao, Q. F.; Wu, Q. S. Environ. Chem. 2013, 32, 1205. [覃禮堂, 劉樹深, 肖乾芬, 吳慶生. 環(huán)境化學, 2013, 32, 1205.]

(37)Tropsha, A.; Gramatica, P.; Gombar, V. K. QSAR Comb. Sci. 2003, 22, 69.

Accurate Activity Predictions of B-Raf Type II Inhibitors via Molecular Docking and QSAR Methods

LIU Hai-Chun1LU Shuai1RAN Ting1ZHANG Yan-Min1XU Jin-Xing1XIONG Xiao1XU An-Yang1LU Tao1,2,*CHEN Ya-Dong1,*
(1School of Basic Science, China Pharmaceutical University, Nanjing 211198, P. R. China;2State Key Laboratory of Natural Medcines, China Pharmaceutical University, Nanjing 210009, P. R. China)

B-Raf kinase plays an important role in the mitogen-activated protein kinase (MAPK) signaling transmission pathway and has been identified as an attractive target for cancer therapy. The exploitation of novel and efficient B-Raf inhibitors has become a hot research topic. In this study, we investigated quantitative structure–activity relationship (QSAR) to probe the origins of the inhibitory activities of B-Raf Type II inhibitors. We used structurally diverse B-Raf Type II inhibitors and an integrated docking and QSAR extended method. We focused mainly on two themes: bioactive conformations and descriptors. First, various molecular docking methods (Glide, Gold, LigandFit, Cdocker, and Libdock) were evaluated, and then all molecules were docked into the B-Raf active site to obtain the bioactive conformations. Secondly, based on the docking results, 16 scoring functions and 21 docking-generated energy-based descriptors were calculated to construct regressionmodels. The results gave highly accurate fitting and had strong predictive abilitiesThe important descriptors were also explored to elucidate the main factors influencing the inhibition activities. The models suggested that the scoring functions (G_Score, -ECD, Dock_Score, and PMF) and docking-generated energy-based descriptors (S(hb_ext), DE(int), and Emodel) were significant. Some new compounds that are potential B-Raf inhibitors were obtained through virtual screening and theoretical predictions using the established models. Such information is useful in guiding the design of novel and robust B-Raf Type II inhibitors.

B-Raf Type II inhibitor; Molecular docking; Scoring function; Docking-generated energybased descriptor; Quantitative structure–activity relationship model

O641

10.3866/PKU.WHXB201510134

Received: July 9, 2015; Revised: October 12, 2015; Published on Web: October 13, 2015.

*Corresponding authors. LU Tao, Email: lutao@cpu.edu.cn; Tel: +86-25-86185179. CHEN Ya-Dong, Email: ydchen@cpu.edu.cn; Tel: +86-25-86185163.

The project was supported by the National Natural Science Foundation of China (21102181, 81302634).

國家自然科學基金(21102181, 81302634)資助項目

?Editorial office of Acta Physico-Chimica Sinica