基于中草药活性成分的GSK-3β抑制剂的分子模拟

唐文强，张静晓，张丽雷

1陕西国际商贸学院医药学院，西安 712046；2湖北民族学院化学与环境工程学院，恩施 471934； 3陕西省中药绿色制造技术协同创新中心，西安 712046

糖原合成酶激酶-3(glycogen synthase kinase-3，GSK-3)是一种组成型活性蛋白激酶，通过对100多种底物的磷酸化参与众多细胞功能的调节[1]。GSK-3通过β-连环蛋白[2]和c-Myc[3]控制细胞的增殖，并且通过TNFα，IL-1β和IL-6等细胞因子表达炎症[4]。因此，GSK-3的过量表达会导致癌症和炎症等疾病[5,6]。GSK-3抑制剂的研究最初集中于神经性疾病，如阿尔茨海默病的治疗，如今GSK-3抑制剂作为抗癌药物也已经成为研究的热点，并开发了许多GSK-3抑制剂[1,7-9]。

GSK-3具有GSK-3α和GSK-3β两种亚型，GSK3α主要参与糖原代谢过程，GSK3β主要调控癌细胞基因转录、加速细胞周期、参与肿瘤细胞侵袭与转移及凋亡过程，是肿瘤治疗的重要靶点。GSK-3β具有折叠的N端和α螺旋的C端两个结构域，ATP结合口袋位于两个结构域之间，主要的氨基酸残基包括Lys85、Aspl33和Vall35等，另外，其活性区域包括T-loop区和正调节因子Tyr2l6，主要的氨基酸包括Arg96、Argl80和Lys205等。目前研制的GSK-3β抑制剂主要包括蒎酮类、靛玉红类和马来酰胺类等多种类型[10]，然而，多数抑制剂的作用机理是结合于GSK-3β的ATP结合口袋部位，从而不仅阻断目标底物蛋白的磷酸化，同时阻止了所有底物的磷酸化，导致Wnt信号传导通路中的β连环蛋白浓度增加，长期使用具有引发癌症的危险[11,12]。如果设计作用靶标是T-loop区的抑制剂，可以阻碍底物与T-loop区的结合，从而阻断底物的磷酸化，但不会影响axin,β-连环蛋白，因这些蛋白不需要预磷酸化，可以直接与活性部位结合[13]。因此，此类药物不会击中Wnt信号路径，与ATP的竞争抑制剂相比，减少癌症发生的可能性。故识别与T-loop区的氨基酸的相互作用是这类抑制剂开发的可行方法。

课题组前期利用系统药理学方法建立了一个中药有效成分筛选的模拟体系，通过ADME筛选，靶点预测和验证等方法，从多味中草药中筛选出了与GSK-3β靶点有相互作用的22个有效成分[14]。本文在前期研究基础上，采用分子对接和分子动力学方法研究了这22个活性成分分别与GSK-3β靶点的ATP结合口袋和T-loop区域的结合方式，揭示其与GSK-3β靶点的相互作用机理，为中草药的药用机理研究和GSK-3β抑制剂的开发提供依据。

1 材料与方法

1.1 材料的准备

从蛋白质晶体数据库(RCSB PDB)中获取GSK-3β的晶体结构模型(PDB ID:5K5N)，除去结晶水和配体PF-04802367，并补全氢原子，记为GSK-3β，作为分子模拟的GSK-3β蛋白模型。采用该晶体结构的原因是它是迄今为止分辨率最高的GSK-3β蛋白晶体模型。小分子化合物的分子结构从化学结构数据库(http://www.chemspider.com/)中获取，作为配体的初始构象。

1.2 分子对接

分子对接采用AutoDock Vina程序[15]进行，使用半柔性对接方法将配体分子与GSK-3β进行对接，其中，GSK-3β被视为一个刚体，与配体分子所有的可旋转键进行对接。结合位点分别设定为GSK-3β的ATP结合口袋和T-loop区域，口袋大小为20×20×20 Å3的盒子，网格间距为1.0 Å，中心分别为ATP结合口袋和T-loop区域的中心。根据在AutoDock Vina的评分函数计算出的结合亲和能筛选出10个最有可能的模型。

1.3 分子动力学模拟

分子动力学(MD)模拟采用GROMACS(版本2016.4)程序[16]，将分子对接得到的复合物结构作为MD模拟的初始构象，GSK-3β部分采用amber99sb-ildn力场，小分子配体部分采用GAFF力场和AM1-BCC电荷。在复合物结构的周围建立可将其完全包围且延伸8 Å的立方体水模型的周期性结构，并加入5个Cl-离子，使体系保持电中性，体系原子的总数目大约为66 000个。在进行MD模拟之前，首先进行500步的能量极小化动力学，再进行100 ps的限制性动力学，令溶剂弛豫。优化结束后，采取缓慢升温的方法，将体系温度在200 ps内由0 K缓慢加热至310 K，之后进行20 ns的NPT模拟，时间步长2 fs。在动力学模拟过程中，运用Velocity-rescale方法控制体系温度，将蛋白和配体作为一组，其他原子作为另外一组，分别进行控温。设置所有和氢原子相连的键被认为是不振动的，使用PME(particle mesh ewald)方法计算长程静电相互作用，范德华相互作用的截断值为12Å。

1.4 结合自由能的计算

采用连续介质模型MM-PBSA(molecular mechanics Poisson-Boltzmann surface area)方法计算配体与蛋白复合物的结合自由能，此方法采用分子力学和连续介质模型估算复合物的结合自由能。

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

G=Emm-TSmm+Gsolvatc

Gsolvate=Gpolar+Gnonpolar

式中EMM为气相内能项，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为平衡状态，即从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的结构，从图中可见，对接后配体的结构与原晶体中的结构基本重合，能够很好地重现晶体构象，表明对接方案是合理的。

图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区域进行对接，对接结合能量的结果如表1所示。从结果可见，所有的结合能量均为负值，表明所筛选得到的化合物均可以与GSK-3β的ATP结合口袋以及T-loop区域结合。另外，绝大部分配体结合于ATP结合口袋的能量高于结合于T-loop区域的能量，表明这些配体倾向于结合在GSK-3β的ATP结合口袋区域，从而与ATP竞争结合达到抑制的目的。

其中，芦丁(M01)、杨酶酮(M16)和二氢丹参酮I (M17)等成分的结合能量均高于配体PF-04802367的结合能量，表明这些化合物有可能获得更好的抑制活性，并且已经有研究表明，芦丁(M01)[19]、杨酶酮(M16)[20]和二氢丹参酮I (M17)[10]具有良好的抗肿瘤活性，并且对GSK-3β具有一定的抑制作用，与分子对接结果一致。另外，人参皂苷Rb1(M20)在GSK-3β的T-loop区的结合能高于在ATP结合口袋的结合能，表明人参皂苷Rb1(M20)可能更易结合于GSK-3β的T-loop区域，从而达到抑制作用。

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结合口袋区域和T-loop区域的构象为初始结构，进行了20 ns的分子动力学模拟。选取这些成分的原因是芦丁(M01)、杨酶酮(M16)、二氢丹参酮I(M17)结合于GSK-3β的ATP结合口袋区域的结合能均大于原晶体结构中配体PF-04802367的结合能，并且是所筛选配体中结合能最高的化合物，可能具有较优的抑制活性。另外，人参皂苷Rb1(M20)结合于GSK-3β的T-loop区域的结合能高于在ATP结合口袋区域的结合能，可能是可结合于T-loop区域的GSK-3β有效抑制剂。

通过分子模拟过程中构象的均方根偏差(RMSD)验证结构是否达到平衡，结果如图2所示。从图中可见，经过约2 ns的模拟后，构象中蛋白和配体的RMSD值上下波动范围均稳定在1 Å的范围之内，表明所研究构象经过平衡后均达到稳定结构。

从最后1ns的运动轨迹中，取出最低能量结构对配体与蛋白之间的相互作用进行了分析。使用软件Ligplot +计算得到芦丁(M01)、杨酶酮(M16)、二氢丹参酮I (M17)和人参皂苷Rb1(M20)与GSK-3β的ATP结合口袋区域和T-loop区域的相互作用网络。

配体和蛋白之间的相互作用主要包括疏水作用和氢键作用。芦丁(M01)、杨酶酮(M16)、二氢丹参酮I(M17)和人参皂苷Rb1(M20)分别与GSK-3β的ATP结合口袋区域的关键氨基酸如Val135、Phe67和Asp133等，以及T-loop区域的关键氨基酸Leu88、Asp181和Cys218等形成了若干强度不同的疏水作用和氢键作用。

表2 小分子配体与GSK-3β残基间的氢键相互作用Table 2 The hydrogen bond interactions between small molecule ligands and GSK-3β residues

氢键是维系蛋白质与配体分子结合稳定性的重要作用力。分别对芦丁(M01)、杨酶酮(M16)、二氢丹参酮I(M17)和人参皂苷Rb1(M20)与GSK-3β形成的氢键进行了分析，并对模拟过程中10 ns至20 ns的运动轨迹中的氢键统计了存活率。判断氢键采用的几何依据为键长不大于3.5 Å，键角不大于30o。结果如表2所示。

从结果可见，形成的氢键均为配体中的含氧和含氮基团与蛋白的残基之间的氢键，表明配体的含氧和含氮基团为形成氢键的关键部分。芦丁(M01)与GSK-3β的ATP结合口袋区域的残基Asp133和Gln185 形成了3个主要氢键，其中两个氢键的存活率达到95%以上，是其能够与GSK-3β的ATP结合区域形成稳定相互作用的主要原因，另外，芦丁(M01)与GSK-3β的T-loop区域的残基Leu88、Gly202和Glu125形成了4个主要氢键，然而其存活率均仅在40%左右，氢键的稳定性小于在ATP结合口袋区域的氢键，这也是芦丁(M01)更倾向于与GSK-3β的ATP结合区域结合的原因。杨酶酮(M16)与GSK-3β的ATP结合口袋区域的残基Val135形成了一个存活率为65%的氢键，相较于与T-loop区域的残基Asp181和Asp200形成的存活率分别为20%和5%的两个氢键更稳定，是该化合物更倾向于与GSK-3β的ATP结合区域结合的原因。二氢丹参酮I(M17)与GSK-3β的ATP结合区域的残基Lys85形成存活率为35%的氢键，相较于与T-loop区域的残基Arg96形成的存活率4%左右的三个氢键更加稳定，因此，二氢丹参酮I(M17)更倾向结合于GSK-3β的ATP结合区域。

另外，人参皂苷Rb1(M20)与GSK-3β的ATP结合口袋区域的残基Asn186和Lys183形成两个氢键，存活率分别为8%和5%，氢键并不稳定，与T-loop区域的残基Lys183、Asp181、Tyr222、Ser219和Ile217形成了6个氢键，存活率分别为17%、5%、8%、11%、16%和28%，尽管存活率不高，但是由于形成氢键数目较多，并且存活率高于与ATP结合口袋区域形成的氢键，故人参皂苷Rb1(M20)更易结合于T-loop区域。

综上可见，所研究的4个配体与GSK-3β的残基形成氢键的数目和稳定性，与其和ATP结合区域和T-loop区域的结合能力具有相同的趋势，是影响化合物与GSK-3β结合能力的重要影响因素。

2.3 结合自由能分析

由于对接得到的结合能量只考虑了分子内能和分子间的能量，是较为粗糙的，并且也没有考虑溶剂对配体结合的影响，为了更准确反应配体和蛋白在水溶剂作用下的结合能量，采用MM-PBSA(Molecular Mechanics-Poisson Bolzmann Surface Area,分子力学泊松玻尔兹曼表面积)方法计算了配体与GSK-3β的结合自由能，结果如表3所示。

将结合自由能分为范德华力、静电力、极性溶剂化能和非极性溶剂化能四个部分。从表3中可见，芦丁(M01)、杨酶酮(M16)和二氢丹参酮I(M17)在ATP结合口袋内的结合能量均大于在T-loop区域的能量，这与前述分子对接能量分析和分子动力学氢键分析的结果一致，这主要是由于小分子配体与GSK-3β的ATP结合口袋区域形成了更为稳定的氢键的原因，并且从结合自由能的分解可见，范德华力，静电力和非极性溶剂化能有利于配体的结合，而极性溶剂化能对配体的结合是不利的。另外，人参皂苷Rb1(M20)与GSK-3β的T-loop区域的结合自由能大于ATP结合口袋区域的结合自由能。这表明人参皂苷Rb1(M20)更易结合于T-loop区域从而起到抑制作用，与前述结果一致，研究类似与人参皂苷Rb1(M20)的相关结构化合物，特别是其与GSK-3β的T-loop区域形成氢键部分的结构，有可能获得与GSK-3β的T-loop区域具有相互作用的抑制剂。

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

3 结论

以一类中草药的活性成分为研究对象，通过分子对接发现，芦丁、杨酶酮、二氢丹参酮I和人参皂苷Rb1能够与GSK-3β良好地结合，其中，芦丁、杨酶酮、二氢丹参酮I主要结合于GSK-3β的ATP结合口袋，人参皂苷Rb1主要结合于GSK-3β的T-loop区域。通过分子动力学模拟获取了配体与蛋白复合物的稳定构象，通过氢键分析发现，配体和蛋白形成的氢键的数目和存活率是影响抑制能力的主要因素，氢键的形成主要取决于配体的含氧和含氮基团与蛋白残基的相互作用。基于人参皂苷Rb1的结构，特别是其与GSK-3β的T-loop区域残基形成氢键部分的结构，通过分子设计有望获得结合于GSK-3β的T-loop区域的有效抑制剂。