基于中草藥活性成分的GSK-3β抑制劑的分子模擬

唐文強，張靜曉，張麗雷

1陜西國際商貿(mào)學院醫(yī)藥學院，西安 712046；2湖北民族學院化學與環(huán)境工程學院，恩施 471934； 3陜西省中藥綠色制造技術協(xié)同創(chuàng)新中心，西安 712046

糖原合成酶激酶-3(glycogen synthase kinase-3，GSK-3)是一種組成型活性蛋白激酶，通過對100多種底物的磷酸化參與眾多細胞功能的調(diào)節(jié)[1]。GSK-3通過β-連環(huán)蛋白[2]和c-Myc[3]控制細胞的增殖，并且通過TNFα，IL-1β和IL-6等細胞因子表達炎癥[4]。因此，GSK-3的過量表達會導致癌癥和炎癥等疾病[5,6]。GSK-3抑制劑的研究最初集中于神經(jīng)性疾病，如阿爾茨海默病的治療，如今GSK-3抑制劑作為抗癌藥物也已經(jīng)成為研究的熱點，并開發(fā)了許多GSK-3抑制劑[1,7-9]。

GSK-3具有GSK-3α和GSK-3β兩種亞型，GSK3α主要參與糖原代謝過程，GSK3β主要調(diào)控癌細胞基因轉錄、加速細胞周期、參與腫瘤細胞侵襲與轉移及凋亡過程，是腫瘤治療的重要靶點。GSK-3β具有折疊的N端和α螺旋的C端兩個結構域，ATP結合口袋位于兩個結構域之間，主要的氨基酸殘基包括Lys85、Aspl33和Vall35等，另外，其活性區(qū)域包括T-loop區(qū)和正調(diào)節(jié)因子Tyr2l6，主要的氨基酸包括Arg96、Argl80和Lys205等。目前研制的GSK-3β抑制劑主要包括蒎酮類、靛玉紅類和馬來酰胺類等多種類型[10]，然而，多數(shù)抑制劑的作用機理是結合于GSK-3β的ATP結合口袋部位，從而不僅阻斷目標底物蛋白的磷酸化，同時阻止了所有底物的磷酸化，導致Wnt信號傳導通路中的β連環(huán)蛋白濃度增加，長期使用具有引發(fā)癌癥的危險[11,12]。如果設計作用靶標是T-loop區(qū)的抑制劑，可以阻礙底物與T-loop區(qū)的結合，從而阻斷底物的磷酸化，但不會影響axin,β-連環(huán)蛋白，因這些蛋白不需要預磷酸化，可以直接與活性部位結合[13]。因此，此類藥物不會擊中Wnt信號路徑，與ATP的競爭抑制劑相比，減少癌癥發(fā)生的可能性。故識別與T-loop區(qū)的氨基酸的相互作用是這類抑制劑開發(fā)的可行方法。

課題組前期利用系統(tǒng)藥理學方法建立了一個中藥有效成分篩選的模擬體系，通過ADME篩選，靶點預測和驗證等方法，從多味中草藥中篩選出了與GSK-3β靶點有相互作用的22個有效成分[14]。本文在前期研究基礎上，采用分子對接和分子動力學方法研究了這22個活性成分分別與GSK-3β靶點的ATP結合口袋和T-loop區(qū)域的結合方式，揭示其與GSK-3β靶點的相互作用機理，為中草藥的藥用機理研究和GSK-3β抑制劑的開發(fā)提供依據(jù)。

1 材料與方法

1.1 材料的準備

從蛋白質(zhì)晶體數(shù)據(jù)庫(RCSB PDB)中獲取GSK-3β的晶體結構模型(PDB ID:5K5N)，除去結晶水和配體PF-04802367，并補全氫原子，記為GSK-3β，作為分子模擬的GSK-3β蛋白模型。采用該晶體結構的原因是它是迄今為止分辨率最高的GSK-3β蛋白晶體模型。小分子化合物的分子結構從化學結構數(shù)據(jù)庫(http://www.chemspider.com/)中獲取，作為配體的初始構象。

1.2 分子對接

分子對接采用AutoDock Vina程序[15]進行，使用半柔性對接方法將配體分子與GSK-3β進行對接，其中，GSK-3β被視為一個剛體，與配體分子所有的可旋轉鍵進行對接。結合位點分別設定為GSK-3β的ATP結合口袋和T-loop區(qū)域，口袋大小為20×20×20 ?3的盒子，網(wǎng)格間距為1.0 ?，中心分別為ATP結合口袋和T-loop區(qū)域的中心。根據(jù)在AutoDock Vina的評分函數(shù)計算出的結合親和能篩選出10個最有可能的模型。

1.3 分子動力學模擬

分子動力學(MD)模擬采用GROMACS(版本2016.4)程序[16]，將分子對接得到的復合物結構作為MD模擬的初始構象，GSK-3β部分采用amber99sb-ildn力場，小分子配體部分采用GAFF力場和AM1-BCC電荷。在復合物結構的周圍建立可將其完全包圍且延伸8 ?的立方體水模型的周期性結構，并加入5個Cl-離子，使體系保持電中性，體系原子的總數(shù)目大約為66 000個。在進行MD模擬之前，首先進行500步的能量極小化動力學，再進行100 ps的限制性動力學，令溶劑弛豫。優(yōu)化結束后，采取緩慢升溫的方法，將體系溫度在200 ps內(nèi)由0 K緩慢加熱至310 K，之后進行20 ns的NPT模擬，時間步長2 fs。在動力學模擬過程中，運用Velocity-rescale方法控制體系溫度，將蛋白和配體作為一組，其他原子作為另外一組，分別進行控溫。設置所有和氫原子相連的鍵被認為是不振動的，使用PME(particle mesh ewald)方法計算長程靜電相互作用，范德華相互作用的截斷值為12?。

1.4 結合自由能的計算

采用連續(xù)介質(zhì)模型MM-PBSA(molecular mechanics Poisson-Boltzmann surface area)方法計算配體與蛋白復合物的結合自由能，此方法采用分子力學和連續(xù)介質(zhì)模型估算復合物的結合自由能。

MM-PBSA方法的計算方程式：

G=Emm-TSmm+Gsolvatc

Gsolvate=Gpolar+Gnonpolar

式中EMM為氣相內(nèi)能項，TSMM為氣相熵項，Gsolvate為溶劑能量項，Gpolar和Gnonpolar分別為極性項和非極性項。計算采用g_MMPBSA程序[17]進行，非極性項Gnonpolar采用溶劑可及表面(SASA)方法進行計算:

Gnonpolar=γSSASA+β

其中 γ=2.2 kJ·mol-1·nm-2，β=3.84 kJ·mol-1。每個復合物采用51個結合構象，在全部20 ns模擬中，選取最后1 ns為平衡狀態(tài)，即從19 ns開始(包含第19 ns)取樣，間隔為20 ps，選取一個結構，至20 ns結束，共51個構象，采用這51個構象的極性和非極性項的平均值作為計算值。

2 結果與討論

2.1 分子對接

首先考察了對接方法的合理性，對已知晶體結構(PDBID:5K5N)[18]中的配體PF-04802367與重新對接后的PF-04802367的結構進行了疊合，結果如圖1所示。圖中黑色部分是原始蛋白晶體中配體PF-04802367的結構，白色部分是重新對接后配體PF-04802367的結構，從圖中可見，對接后配體的結構與原晶體中的結構基本重合，能夠很好地重現(xiàn)晶體構象，表明對接方案是合理的。

圖1 配體PF-04802367在原晶體中的結構與重新對接后的結構疊合圖Fig.1 The overlapped structures of ligand PF-04802367 in the crystal and re-docking systems

表1 分子對接結果Table 1 The results of molecular docking

將篩選得到的22個中草藥有效成分作為配體，與GSK-3β的ATP結合口袋以及T-loop區(qū)域進行對接，對接結合能量的結果如表1所示。從結果可見，所有的結合能量均為負值，表明所篩選得到的化合物均可以與GSK-3β的ATP結合口袋以及T-loop區(qū)域結合。另外，絕大部分配體結合于ATP結合口袋的能量高于結合于T-loop區(qū)域的能量，表明這些配體傾向于結合在GSK-3β的ATP結合口袋區(qū)域，從而與ATP競爭結合達到抑制的目的。

其中，蘆丁(M01)、楊酶酮(M16)和二氫丹參酮I (M17)等成分的結合能量均高于配體PF-04802367的結合能量，表明這些化合物有可能獲得更好的抑制活性，并且已經(jīng)有研究表明，蘆丁(M01)[19]、楊酶酮(M16)[20]和二氫丹參酮I (M17)[10]具有良好的抗腫瘤活性，并且對GSK-3β具有一定的抑制作用，與分子對接結果一致。另外，人參皂苷Rb1(M20)在GSK-3β的T-loop區(qū)的結合能高于在ATP結合口袋的結合能，表明人參皂苷Rb1(M20)可能更易結合于GSK-3β的T-loop區(qū)域，從而達到抑制作用。

2.2 分子動力學模擬

圖2 GSK-3β復合物模型分子動力學模擬期間RMSD值隨時間的變化Fig.2 The values of RMSD during molecular dynamics simulation of GSK-3β complexes

為了進一步研究所篩選得到的配體與GSK-3β的相互作用機理，選取蘆丁(M01)、楊酶酮(M16)、二氫丹參酮I(M17)和人參皂苷Rb1(M20)分別結合于GSK-3β的ATP結合口袋區(qū)域和T-loop區(qū)域的構象為初始結構，進行了20 ns的分子動力學模擬。選取這些成分的原因是蘆丁(M01)、楊酶酮(M16)、二氫丹參酮I(M17)結合于GSK-3β的ATP結合口袋區(qū)域的結合能均大于原晶體結構中配體PF-04802367的結合能，并且是所篩選配體中結合能最高的化合物，可能具有較優(yōu)的抑制活性。另外，人參皂苷Rb1(M20)結合于GSK-3β的T-loop區(qū)域的結合能高于在ATP結合口袋區(qū)域的結合能，可能是可結合于T-loop區(qū)域的GSK-3β有效抑制劑。

通過分子模擬過程中構象的均方根偏差(RMSD)驗證結構是否達到平衡，結果如圖2所示。從圖中可見，經(jīng)過約2 ns的模擬后，構象中蛋白和配體的RMSD值上下波動范圍均穩(wěn)定在1 ?的范圍之內(nèi)，表明所研究構象經(jīng)過平衡后均達到穩(wěn)定結構。

從最后1ns的運動軌跡中，取出最低能量結構對配體與蛋白之間的相互作用進行了分析。使用軟件Ligplot +計算得到蘆丁(M01)、楊酶酮(M16)、二氫丹參酮I (M17)和人參皂苷Rb1(M20)與GSK-3β的ATP結合口袋區(qū)域和T-loop區(qū)域的相互作用網(wǎng)絡。

配體和蛋白之間的相互作用主要包括疏水作用和氫鍵作用。蘆丁(M01)、楊酶酮(M16)、二氫丹參酮I(M17)和人參皂苷Rb1(M20)分別與GSK-3β的ATP結合口袋區(qū)域的關鍵氨基酸如Val135、Phe67和Asp133等，以及T-loop區(qū)域的關鍵氨基酸Leu88、Asp181和Cys218等形成了若干強度不同的疏水作用和氫鍵作用。

表2 小分子配體與GSK-3β殘基間的氫鍵相互作用Table 2 The hydrogen bond interactions between small molecule ligands and GSK-3β residues

氫鍵是維系蛋白質(zhì)與配體分子結合穩(wěn)定性的重要作用力。分別對蘆丁(M01)、楊酶酮(M16)、二氫丹參酮I(M17)和人參皂苷Rb1(M20)與GSK-3β形成的氫鍵進行了分析，并對模擬過程中10 ns至20 ns的運動軌跡中的氫鍵統(tǒng)計了存活率。判斷氫鍵采用的幾何依據(jù)為鍵長不大于3.5 ?，鍵角不大于30o。結果如表2所示。

從結果可見，形成的氫鍵均為配體中的含氧和含氮基團與蛋白的殘基之間的氫鍵，表明配體的含氧和含氮基團為形成氫鍵的關鍵部分。蘆丁(M01)與GSK-3β的ATP結合口袋區(qū)域的殘基Asp133和Gln185 形成了3個主要氫鍵，其中兩個氫鍵的存活率達到95%以上，是其能夠與GSK-3β的ATP結合區(qū)域形成穩(wěn)定相互作用的主要原因，另外，蘆丁(M01)與GSK-3β的T-loop區(qū)域的殘基Leu88、Gly202和Glu125形成了4個主要氫鍵，然而其存活率均僅在40%左右，氫鍵的穩(wěn)定性小于在ATP結合口袋區(qū)域的氫鍵，這也是蘆丁(M01)更傾向于與GSK-3β的ATP結合區(qū)域結合的原因。楊酶酮(M16)與GSK-3β的ATP結合口袋區(qū)域的殘基Val135形成了一個存活率為65%的氫鍵，相較于與T-loop區(qū)域的殘基Asp181和Asp200形成的存活率分別為20%和5%的兩個氫鍵更穩(wěn)定，是該化合物更傾向于與GSK-3β的ATP結合區(qū)域結合的原因。二氫丹參酮I(M17)與GSK-3β的ATP結合區(qū)域的殘基Lys85形成存活率為35%的氫鍵，相較于與T-loop區(qū)域的殘基Arg96形成的存活率4%左右的三個氫鍵更加穩(wěn)定，因此，二氫丹參酮I(M17)更傾向結合于GSK-3β的ATP結合區(qū)域。

另外，人參皂苷Rb1(M20)與GSK-3β的ATP結合口袋區(qū)域的殘基Asn186和Lys183形成兩個氫鍵，存活率分別為8%和5%，氫鍵并不穩(wěn)定，與T-loop區(qū)域的殘基Lys183、Asp181、Tyr222、Ser219和Ile217形成了6個氫鍵，存活率分別為17%、5%、8%、11%、16%和28%，盡管存活率不高，但是由于形成氫鍵數(shù)目較多，并且存活率高于與ATP結合口袋區(qū)域形成的氫鍵，故人參皂苷Rb1(M20)更易結合于T-loop區(qū)域。

綜上可見，所研究的4個配體與GSK-3β的殘基形成氫鍵的數(shù)目和穩(wěn)定性，與其和ATP結合區(qū)域和T-loop區(qū)域的結合能力具有相同的趨勢，是影響化合物與GSK-3β結合能力的重要影響因素。

2.3 結合自由能分析

由于對接得到的結合能量只考慮了分子內(nèi)能和分子間的能量，是較為粗糙的，并且也沒有考慮溶劑對配體結合的影響，為了更準確反應配體和蛋白在水溶劑作用下的結合能量，采用MM-PBSA(Molecular Mechanics-Poisson Bolzmann Surface Area,分子力學泊松玻爾茲曼表面積)方法計算了配體與GSK-3β的結合自由能，結果如表3所示。

將結合自由能分為范德華力、靜電力、極性溶劑化能和非極性溶劑化能四個部分。從表3中可見，蘆丁(M01)、楊酶酮(M16)和二氫丹參酮I(M17)在ATP結合口袋內(nèi)的結合能量均大于在T-loop區(qū)域的能量，這與前述分子對接能量分析和分子動力學氫鍵分析的結果一致，這主要是由于小分子配體與GSK-3β的ATP結合口袋區(qū)域形成了更為穩(wěn)定的氫鍵的原因，并且從結合自由能的分解可見，范德華力，靜電力和非極性溶劑化能有利于配體的結合，而極性溶劑化能對配體的結合是不利的。另外，人參皂苷Rb1(M20)與GSK-3β的T-loop區(qū)域的結合自由能大于ATP結合口袋區(qū)域的結合自由能。這表明人參皂苷Rb1(M20)更易結合于T-loop區(qū)域從而起到抑制作用，與前述結果一致，研究類似與人參皂苷Rb1(M20)的相關結構化合物，特別是其與GSK-3β的T-loop區(qū)域形成氫鍵部分的結構，有可能獲得與GSK-3β的T-loop區(qū)域具有相互作用的抑制劑。

表3 配體與GSK-3β相互作用的結合自由能分析結果Table 3 Binding free-energy analysis of the interaction between ligands and GSK-3β

3 結論

以一類中草藥的活性成分為研究對象，通過分子對接發(fā)現(xiàn)，蘆丁、楊酶酮、二氫丹參酮I和人參皂苷Rb1能夠與GSK-3β良好地結合，其中，蘆丁、楊酶酮、二氫丹參酮I主要結合于GSK-3β的ATP結合口袋，人參皂苷Rb1主要結合于GSK-3β的T-loop區(qū)域。通過分子動力學模擬獲取了配體與蛋白復合物的穩(wěn)定構象，通過氫鍵分析發(fā)現(xiàn)，配體和蛋白形成的氫鍵的數(shù)目和存活率是影響抑制能力的主要因素，氫鍵的形成主要取決于配體的含氧和含氮基團與蛋白殘基的相互作用。基于人參皂苷Rb1的結構，特別是其與GSK-3β的T-loop區(qū)域殘基形成氫鍵部分的結構，通過分子設計有望獲得結合于GSK-3β的T-loop區(qū)域的有效抑制劑。