分子动力学模拟研究活化剂对葡萄糖激酶活性的影响及活化机理

姚雪霞

（南京农业大学 工学院，江苏 南京210031）

分子动力学模拟研究活化剂对葡萄糖激酶活性的影响及活化机理

姚雪霞

（南京农业大学 工学院，江苏 南京210031）

葡萄糖激酶（GK）催化葡萄糖转变为6－磷酸葡萄糖，是糖代谢的第一步；GK的活性改变对葡萄糖的代谢影响很大.通过分子动力学模拟研究了葡萄糖激酶活化剂（GKA）对GK活性的影响，探讨了GKA的活化机理.包结自由能分析表明，范德华作用是GK和GKA体系的主要驱动力；与此同时，葡萄糖使得GK和GKA的结合增强.构象分析结果表明，GKA能够在一段时间内限制GK的构象，但不能长久维持GK的活化构象不变，这和以往的实验结果相一致.相关理论分析结果不仅有利于从原子水平解释GKA的活化机理，也可为设计新的用于糖尿病治疗的GKA提供一定的理论参考.

葡萄糖激酶；活化剂；活性；活化机理；分子动力学模拟

葡萄糖激酶（GK）催化葡萄糖转变为6－磷酸葡萄糖，是糖代谢极为关键的一步，所以GK活性异常在糖代谢紊乱的发生发展中起着重要作用［1－2］.文献［3］对GK的晶体结构和催化过程进行了研究，认为GK存在三种构象，即：关闭型、开启型和超开启型.其中，前两种为活性构象，第3种为非活性构象.其催化过程如下：当血糖浓度较低时，GK主要以稳定的超开启构象存在，当葡萄糖分子与之结合后，其构象缓慢转变为开启型构象，然后转变为关闭型构象，在ATP存在下同时发生酶催化反应——葡萄糖磷酸化.反应结束后GK回复到开启型构象，同时释放出6－磷酸葡萄糖（G6P）和ADP.若有新的葡萄糖分子结合则继续发生磷酸化，进入快循环；若没有新的葡萄糖分子结合，GK则转变为超开启构象，进入慢循环.若GK的氨基酸残基发生失活突变，则影响相应的催化过程，导致青年型早发糖尿病（MODY 2）［4］发生.为了治疗MODY 2型糖尿病，以GK为靶点的小分子活化剂（GKA）的研究非常活跃，并取得了很大的进展［5－9］.

观察GK的晶体结构［8］（如图1所示），发现GK由448个氨基酸残基构成，折叠成大小不等的两个结构域，中间凹槽为连接域，葡萄糖分子磷酸化的活性位点位于该连接域内，而GKA则位于GK的异构点，距离葡萄糖底物大约2nm左右.虽然实验已经证实GKA可以提高GK的活性，但是还是有很多基本问题没有解决.如，异构点距离葡萄糖底物的距离为2nm左右，GKA如何提高GK的活性？文献［10］虽然讨论了其中一个GKA对GK的影响，但是还有很多细节问题没有解决.如GKA是否一直对GK的构象起着约束作用？另外葡萄糖的存在反过来对GKA有什么影响？以及其他GKA对GK的影响是否相同？为了解决这些疑问，很有必要对其他GKA对GK的活化机理做进一步的理论研究.

以往我们曾分析GK的活化突变对GK活性的影响［11－12］，发现它们主要通过约束GK的构象来改变GK的活性.活化剂是否采用同样的方式来改变GK的活性呢？为了研究GKA的活化机理，作者通过分子动力学（molecular dynamical，MD）模拟和 MM－PBSA（molecular mechanics／Poisson Boltzmann surface area）方法［13－14］，对GKA、GK以及葡萄糖的相互作用进行了理论研究.

1 理论与方法

1.1 蛋白质构建

由于第一个包含GK、GKA和葡萄糖复合物的晶体结构［3］缺少实验数据EC50值，因此本文作者选用了一个新研制出来的晶体结构［8］作为研究对象.该体系和第一个体系相比，只有GKA不同，GKA的分子结构如图2所示.

GK－GKA－glucose复合体系的初始结构取自于文献（PDB code：3FR0）［8］.目前单独的GK－glucose和单独的GK－GKA晶体结构还没有文献报道，因此单独的GK－glucose和GK－GKA只能从复合体系中去掉另一个小分子作为初始结构.另外本文中所指的GK－free是指直接从GK－GKA－glucose复合体系中去掉GKA和葡萄糖两个小分子得到的GK结构，并不是真正的晶体结构GK－free.

图1 GK－活化剂－葡萄糖复合物的晶体结构Fig.1 The crystal structure of GK in complex with GKA and glucose

图2 GKA的分子结构式Fig.2 Structural formulas of the GKA

1.2 MD模拟

MD模拟的初始结构如上所述.所有的MD模拟均采用AMBER 10程序包［15］.AMBER03力场用于蛋白质分子，GAFF力场用于葡萄糖和GKA小分子.葡萄糖和GKA小分子都采用Gaussian 03程序包［16］的HF／6－31G＊基组优化，采用RESP程序计算其电荷.程序LeaP用于中和体系，使其总电荷为零.

对于所有体系，均在溶质外围加上1nm的水分子层.采用周期性边界条件和particle mesh ewald（PME）方法用于处理长程静电作用［17］.非键截距为1.2nm，使用SHAKE选项用于约束和氢原子连接的键.

在运行MD之前，先进行两步能量优化：第一步约束溶质，用最速下降法优化5 000步，再用共轭梯度法优化5 000步；第二步去约束后再优化整个体系，用最速下降法优化2 500步，再用共轭梯度法优化2 500步，收敛条件为能量梯度小于4.18×10－4kJ·mol－1·nm－1.

MD模拟分为两步：首先进行200ps的约束溶质分子MD模拟，即通过200ps预热，将体系从0K加热到300K；然后进行10ns无约束恒温MD模拟，温度为300K，压力为1×105Pa（即NPT系综）.系统的温度和压力则通过Langevin piston和Hoover’s thermostat实现控制.模拟中采用的积分步长为2fs，非键连接表每50步更新一次.每4ps保存一次轨迹结构用于后期数据处理.

2 结果与讨论

2.1 MD轨迹稳定性的分析

体系在模拟过程中产生的结构相对于晶体结构随时间变化的均方根偏差（RMSD）是衡量体系是否稳定和达到收敛的重要依据.为了检测MD轨迹的稳定性，分别对GK－free，GK－GKA，GK－glucose和 GK－GKA－glucose四个轨迹中GK的Cα原子进行RMSD分析，结果如图3所示.分析图3，发现相比于其他轨迹，GK－free的轨迹RMSD值最大，最不稳定.而GK－glucose和GK－GKA－glucose的RMSD值比较小，轨迹相对稳定一些.但是仔细分析，发现GK－GKA－glucose的轨迹最稳定.这表明GKA的存在，可以使得GK－glucose体系更加稳定.

2.2 GKA和GK之间静电作用的分析

分析其晶体结构（如图1所示），可以看到GKA处于GK背部的空腔内，不仅和GK有非键的范德华作用，而且还有氢键作用.晶体结构中出现的氢键，在动态的MD过程中会有什么变化？下面分别讨论两种MD情况下GK和GKA之间的氢键随时间的变化：（1）GK和GKA单配体作用.从图4（A）中可以看到其中的两个氢键非常稳定，只有第三个氢键不太稳定，在MD运行到3ns附近发生断裂，在6ns左右又重新形成；（2）GK和GKA、葡萄糖双配体作用.分析三个氢键距离随模拟时间的变化，可以发现在整个MD过程中，三个氢键都非常稳定.由于以上两种情况下MD模拟参数都完全相同，唯一的区别在于第二种情况下含有葡萄糖小分子.因此我们可以得出这样的结论，即葡萄糖的存在使GK和GKA的结合更加稳定.这和文献报道［3］一致，即GKA可以提高GK对葡萄糖的活性，同时葡萄糖的存在，反过来又可以促进GK和GKA的结合.迄今为止，还没有发现在没有葡萄糖存在的条件下，GK和GKA单配体结合的晶体结构，这也说明GK单独和GKA在一起是不能长期稳定存在的.

图3 GK Cα原子在10ns MD模拟过程中的RMSD轨迹Fig.3 The RMSD trajectories of Cαatoms of GK during 10ns MD simulations

图4 GK和GKA之间的氢键键长随时间的变化Fig.4 Time evolution of three hydrogen bond lengths between GK and GKA

2.3 GKA和GK之间包结自由能的分析

通过MM－PBSA方法［11］计算了两个模拟体系 （GK－GKA和GK－GKA－glucose）中GKA和GK之间的相互作用能，结果如表1所示.分析计算结果，发现不论是哪一种模拟体系，GK和GKA之间的主要作用力都是范德华力.尽管GK和GKA之间真空下的静电作用能为负值，但是极性的溶剂化能却是一个更大的正值，导致它们之间总的静电作用能为正值，不利于GK和GKA的稳定结合.通过能量计算表明，GKA主要是通过范德华力和GK结合在一起.

表1 GK和配体的包结自由能 （kJ·mol－1）aTable 1 Binding free energy of the ligands and GK

另外分析GK－GKA和GK－GKA－glucose两个体系的能量差别，发现葡萄糖的存在，使得GK和GKA之间的能量变得更低.具体分析各个能量项，发现主要的差别不是范德华力，而是静电作用力.葡萄糖的存在使得GK－GKA之间的静电作用能变低，最终导致总的包结能降低.以上的计算结果和文献报道［3］一致，即葡萄糖可以促进GK和GKA的稳定结合.

2.4 GKA对GK构象影响的分析

前面我们分析了GK和GKA之间的主要作用力.下面将分析GKA对GK构象的影响.首先分析GKA对GK夹角的影响（如图5所示）.从图5中可以看到对于GK－free，大小域的夹角在开始的1ns变化很大，从15°迅速增加到30°左右.随后，GK－free夹角的值趋于稳定，说明GK达到一种相对较稳定的构象.GK－free夹角的这种变化趋势和其RMSD的变化趋势相一致，从图3中可以看到，其RMSD值在开始的1 ns上升非常快，迅速达到3nm后趋于稳定.和GK－free相比，加了GKA以后，GK大小域的夹角和初始结构夹角相比，变化较小.从图5中可以看到，在整个5ns时间内，夹角的值比较稳定，在18°左右波动.因此可以得出GKA对GK的构象变化确实起到了一种门闩的作用，约束着GK大小域的开合.但是5ns以后，我们发现GKGKA体系的夹角从18°增加到22°.以上结果表明GKA可以在一段时间内约束GK的构象，使之处于活化状态，但不能一直保持GK的构象不变，这和以往的实验结果是一致的［3］.

为了进一步研究GKA对GK构象的影响，分别用动态相关性矩阵（Dynamic cross－correlation matrices，DCCM）［18－19］分析了四种MD情况下GK 448个氨基酸残基内部的相关性，如图6所示.DCCM最吸引人的地方在于预测一些非相邻的氨基酸之间的动态相关性，也就是我们通常所说的蝴蝶效应，这是静态的晶体结构无法提供的.因此DCCM已被广泛用于蛋白质动态相关性分析［18－19］，并取得了很好的预测效果.

图5 GK－free和GK－GKA的夹角随时间的进展图Fig.5 Cleft－angle profiles along the GK－free and GK－GKA trajectories

图6 GK的DCCM动态相关性分析Fig.6 DCCM illustrating the correlation of motion between residues in GK

图6 （A）为GK－free的DCCM 图，从图6（A）可以看到，大部分氨基酸残基的相关性很小，颜色以浅蓝色和绿色为主，很小的黄色区域表示该区域的氨基酸残基之间存在正的相关性.另外还存在一些很小的蓝色区域，表示该部分残基之间的负相关性较强.

和GK－free相比，GK－GKA的氨基酸残基负相关性非常明显，部分区域的颜色显示为深蓝色（如图6（B）所示）.仔细分析，我们发现GK 300～400的氨基酸残基和100～150的氨基酸残基之间有很明显的负相关性.进一步分析晶体结构（如图7所示），发现300～400的氨基酸残基位于GK的大域，而100～150的氨基酸残基位于GK的小域，它们之间相距甚远.在没有GKA存在的情况下，它们之间并没有负相关性运动，但是在GKA存在的情况下，这种负相关运动则非常明显.表明GKA可以影响GK大小域的开合运动.这和我们前面分析GKA对GK夹角的影响是一致的.

GK－glucose的DCCM如图6（C）所示.分析图6（C）发现，在葡萄糖存在的情况下，GK的正负相关性运动都非常不明显，可忽略不计.我们推测是因为GK和葡萄糖的强烈静电作用，导致GK的构象非常稳定，所以观察不到GK的相关性运动.

最后一个DCCM图用来说明GK－GKA－glucose的氨基酸残基的相关性运动（见图6（D））.分析图6（D）发现，GK－GKA－glucose的负相关性介于GK－GKA和GK－glucose之间.我们推测是因为葡萄糖和GK的相互作用，导致GKA对GK的负相关性影响减小.

以上DCCM分析和作者以前所研究的Y214C［11］和M197V［12］活性突变的DCCM结果基本相同.表明该活化剂的作用类似于GK本身的一些活性突变，都是通过约束GK的构象来提高GK的活性.

图7 GK（绿色），残基ID：100～150（红色）和残基ID：300～400（紫色）Fig.7 GK （green），RES ID：100～150（red）and RES ID：300～400（purple）

3 结论

通过分子动力学模拟分析了GKA和GK的相互作用.包结自由能分析表明，范德华作用是GK和GKA体系的主要驱动力.同时葡萄糖的存在使得GK和GKA结合得更加稳定.构象分析表明，GKA可以在一段时间约束GK的构象，但是不能一直维持GK的活化构象不变，这和以往的实验结果相一致.本研究从原子水平解释了GKA的活化机制，为今后设计新的GKA小分子提供了一定的理论参考.

致谢：南京大学理论与计算化学研究所提供计算机服务器用于计算.

［1］FERRÉT，RIU E，BOSCH F，et al.Evidence from transgenic mice that glucokinase is rate limiting for glucose utilization in the liver［J］.FASEB J，1996，10（10）：1213－1218.

［2］NOGUEIRA F N，SANTOS M F，NICOLAU J.Influence of streptozotocin－induced diabetes on hexokinase activity of rat salivary glands［J］.J Physiol Biochem，2005，61（3）：421－428.

［3］KAMATA K，MITSUYA M，NISHIMURA T，et al.Structural basis for allosteric regulation of themonomeric allosteric enzyme human glucokinase［J］.Structure，2004，12（3）：429－438.

［4］VIONNET N，STOFFEL M，TAKEDA J，et al.Nonsense mutation in the glucokinase gene causes early－onset non－insulin－dependent diabetes mellitus［J］.Nature，1992，356（6371）：721－722.

［5］GRIMSBY J，SARABU R，CORBETT W L，et al.Allosteric activators of glucokinase：potential role in diabetes therapy［J］.Science，2003，301（5631）：370－373.

［6］FUTAMURA M，HOSAKA H，KADOTANI A，et al.An allosteric activator of glucokinase impairs the interaction of glucokinase and glucokinase regulatory protein and regulates glucose metabolism［J］.J Biol Chem，2006，281（49）：37668－37674.

［7］IINO T，TSUKAHARA D，KAMATA K，et al.Discovery of potent and orally active 3－alkoxy－5－phenoxy－N－thiazolyl benzamides as novel allosteric glucokinase activators［J］.Bioorg Med Chem，2009，17（7）：2733－2743.

［8］NISHIMURA T，IINO T，MITSUYA M，et al.Identification of novel and potent 2－amino benzamide derivatives as allosteric glucokinase activators［J］.J Bioorg Med Chem Lett，2009，19（5）：1357－1360.

［9］MATSCHINSKY F M.Assessing the potential of glucokinase activators in diabetes therapy：nature reviews drug discovery［J］.Nat Rev Drug Discov，2009，8（5）：399－416.

［10］ZHANG Jian，LI Chen Jing，CHEN Kai Xian，et al.Conformational transition pathway in the allosteric process of human glucokinase［J］.Proc Natl Acad Sci USA，2006，103（36）：13368－13373.

［11］姚雪霞.葡萄糖激酶Y214C活性突变的理论研究 ［J］.化学研究，2010，21（6）：71－76.

［12］姚雪霞，徐素明.葡萄糖激酶 M197V突变的活化机理研究 ［J］.化学学报，2011，69（4）：405－410.

［13］WANG Jun Mei，MORIN P，WANG Wei，et al.Use of MM－PBSA in reproducing the binding free Energies to HIV－1 RT of TIBO derivatives and predicting the binding mode to HIV－1RT of efavirenz by docking and MM－PBSA ［J］.J Am Chem Soc，2001，123（22）：5221－5230.

［14］GOUDA H，KUNTZ I D，CASE D A，et al.Free energy calculations for theophylline binding to an RNA aptamer：MMPBSA and comparison of thcrmodynamic integration methods［J］.Biopolymers，2003，68（1）：16－34.

［15］CASE D A，DARDEN T A，CHEATHAM T E，et al.AMBER 10［CP］.University of California：San Francisco，2008.

［16］FRISCH M J，TRUCKS G W，SCHLEGEL H B，et al.Gaussian 03，Revision C.02［CP］.Gaussian，Inc.：Wallingford CT，2004.

［17］DARDEN T，YORK D，PEDERSEN L.Particle mesh Ewald：An N.log（N）method for Ewald sums in large systems［J］.J Chem Phy，1993，98（12）：10089－10092.

［18］LI Li Wei，UVERSKY V N，DUNKER A K，et al.A computational investigation of allostery in the catabolite activator protein［J］.J Am Chem Soc，2007，129（50）：15668－15676.

［19］GHOSH A，VISHVESHWARA S.A study of communication pathways in methionylt RNA synthetase by molecular dynamics simulations and structure network analysis［J］.Proc Natl Acad Sci USA，2007，104（40）：15711－15716.

Effect of activator on activity of glucokinase and activation mechanism of the activator—study by molecular dynamics simulation

YAO Xuexia
（College of Engineering，Nanjing Agricultural University，Nanjing210031，Jiangsu，China）

Glucokinase（GK）is a glycolic enzyme that catalyzes the phosphorylation of glucose to glucose－6－phosphate in the first step of glycolysis，and the change in GK activity plays an important role in abnormal glycemia.Thus molecular dynamics simulation was performed to study the effect of glucokinase activator（GKA）on the activity of GK，and the activation mechanism of GKA was examined.The binding free－energy analysis suggests that van der Walls forces are the main driving forces for GK－GKA complex，and glucose can enhance the binding between GK and GKA.Moreover，conformation analysis shows that GKA can constrain the conformation of GK in a period of time，but it cannot keep the active state of GK in an extended period，which is in agreement with previously reported experiment results.These observations are beneficial to understanding the activation mechanism of GKA at an atomic level and to guiding the design of novel GKA for the treatment of diabetes.

glucokinase；activator；activity；activation mechanism；molecular dynamic simulation

O 629.8

A

1008－1011（2012）04－0091－06

2011－12－20.

江苏省农机局科研启动基金资助项目（gxz10008）.

姚雪霞（1975－），女，讲师，研究方向为物理化学.E－mail：yaoxx0623＠sina.com.