一种靶向作用于蛋白激酶非活性构象的化合物筛选方法

摘要

本发明属于蛋白质结构预测和药物分子虚拟筛选技术领域，具体为一种靶向作用于蛋白激酶非活性构象的化合物筛选方法。本发明包括蛋白激酶活性链段构象预测方法，从激酶的DFG-in活性构象来产生相应的DFG-out非活性构象；还包含II型抑制剂对接后结合构象的挑选方法，用于构象预测和虚拟筛选中小分子的挑选。本发明已经在7种已知非活性构象的蛋白激酶上进行了计算验证，成功率接近96%。本发明方法已用于预测结核杆菌的PknB蛋白激酶的非活性构象，并虚拟筛选了PknB可能的II型抑制剂，经抑菌实验验证，已发现2种小分子已证明具有抑菌作用。

说明书

技术领域

本发明属于蛋白质结构预测和药物分子虚拟筛选技术领域，具体涉及蛋白激酶的非活性构象的预测与筛选方法，蛋白激酶II型抑制剂对接构象的预测与虚拟筛选方法。

背景技术

随着人类基因组测序的完成，人们发现基因组中约有518种不同的基因编码了各类激酶，占了总基因数的2%左右。这个庞大的蛋白激酶家族调节着细胞内几乎方方面面的活动，包括信号的转导、DNA的转录以及细胞代谢、细胞周期的运行，进而控制着细胞的分化和凋亡。所以激酶的活化与失活关系重大，一旦产生异常便会导致各种疾病的发生，包括癌症、中枢神经系统紊乱、心血管疾病以及自体免疫失调等。正因为如此，长期以来蛋白激酶都被认为是最重要的药物靶点之一。

尽管激酶家族的成员众多，但是X射线衍射得到的晶体结构表明所有蛋白激酶的催化结构域都十分的相似——一般由一个较小的N-lobe和一个较大的C-lobe组成，在2个lobe之间存在一道裂缝，ATP结合位点以及对于磷酸基团转移重要的活性链段（activation loop）都位于那里。几乎所有蛋白激酶在ATP结合位点附近都有一段3个氨基酸残基的序列――Asp-Phe-Gly（DFG），而DFG又同时处于活性链段的N端。这段基序的构象状态已经被证明对于激酶的活性具有决定性的作用。在活性状态下，Phe位于ATP结合位点附近的疏水腔中，而Asp位于腔外活性链段的另一侧（DFG-in构象）；然而在一个被称为DFG-flip的变化发生后，激酶就会处于非活性状态。在这个状态下，这2个残基的位置恰恰相反：Asp位于疏水腔侧，而Phe被旋转到外侧（DFG-out构象）。

由于处于DFG-out构象的蛋白激酶是非活性的，这启发了人们去寻找结合这种结构的药物小分子来抑制激酶。这种小分子被称为二型（II型）抑制剂，它们能够诱导蛋白激酶形成非活性的DFG-out构象并占据由此产生的疏水腔。研究证明它们较之靶向激酶活性结构的type-I抑制剂具有更好的特异性和有效性。抗癌药物伊马替尼（imatinib）就是其中的一个例子。它能够特异性地结合BCR-ABL、c-Abl、 c-Kit和PDGFR等激酶，从而治疗慢性骨髓白血病和其他癌症。

基于受体结构的药物设计是现今新药开发中的一个重要手段。但是至今只有小部分的DFG-out晶体结构被人们得到，所以在结构数据上存在的空缺大大延缓了发现新II型抑制剂的进度。在2008年Protein Data Bank的哺乳动物激酶数据中，70%以上的结构都是DFG-in构象的，22%处于过渡构象，还有约3%的结构是II型抑制剂不相容的DFG-out结构。因此，目前已知的激酶结构数据绝大部分都无法用于基于受体结构的II型抑制剂的开发研究。为了解决这个问题，发展出一种能从已有的大量激酶DFG-in构象得到相应的DFG-out构象的计算模拟方法至关重要。

目前提出一个可能的DFG-flip的机制分为3个步骤：活性构象：DFG-in/αC-in ?过渡构象：DFG-in/αC-out? 非活性构象：DFG-out/αC-in。

在最近以c-Abl为对象的研究中发现，DFG中Asp的质子化对产生DFG-out构象起着关键的作用，Shan et al. 用分子动力学（MD）的方法成功模拟出了从DFG-in到DFG-out构象的变化 （Shan et al., Proc Natl. Acad. Sci. U.S.A, 2007,  Vol106: 139-144），发现对于Asp和Phe来说最重要的是它们所处的环境的改变。在DFG-in构象中，Asp处于一个极性、充满电荷的水性坏境中，而Phe处于一个疏水的环境中。在发生DFG-flip之后，情况正好相反，Phe处于极性的环境而Asp处在一个几乎疏水的环境。电荷在疏水的环境中要以极大的自由能为代价，所以Asp的质子化是十分重要的。同时，在DFG-in的构象中Asp骨架的扭转角是一个高能状态，而Phe此时就处在疏水环境来使能量达到平衡。Shan et al.和过去的研究都暗示着DFG-in可能并不是一个能量最低的构象而是一个在功能上最优化的结构。DFG-flip可能是用来在磷酸化后促进ADP的释放，而这在激酶催化反应中是一个重要的限速步骤。因此蛋白激酶的DFG-in和DFG-out构象的能量很可能是相近的，它们都是最低的能量构象并在实际溶液中处于动态平衡。

另一方面，Kufareva et al.尝试使用所谓的DOLPHIN（deletion-of-loop Asp-Phe-Gly-in）模型来筛选药物小分子，即直接将DFG-in构象的激酶 中包含DFG的6个残基删除后的模型来预测II型抑制剂的结合构象以及虚拟筛选，取得了意想不到的效果 （Kufareva et al.，Journal of Medicinal Chemistry， 2009， Vol51: 7921-7932）。这说明除了DFG以外，激酶的其他部分也存在与II型抑制剂结合有关的决定性结构信息。

目前，虽然对药物敏感的肺结核是可以治愈的，但是需要很长的疗程并使用多种药物。不仅如此，对于结核分枝杆菌（Mycobacterium tuberculosis）来说，出现了越来越多的抗药性菌种，而且如何有效地控制肺结核的传染也面临严峻的挑战。另外，大量出现的多药抗性的菌种等（至少抗利福平和异烟肼）都表明发现新的治疗肺结核的药物已经刻不容缓。

另一方面，在结核分枝杆菌中，Ser/Thr蛋白激酶已发现可磷酸化几百种蛋白质底物。在11种结核杆菌的Ser/Thr蛋白激酶中有2种是可溶的；还有9种被推测为跨膜受体——N端是类似真核激酶的催化结构域，通过单一跨膜螺旋与胞外感受器结构域连接。以结核杆菌的跨膜受体激酶PknB为例，它对结核分枝杆菌的生长至关重要。其磷酸化底物涉及肽聚糖的合成、细胞分裂、应激反应、转录、代谢调控以及其他Ser/Thr蛋白激酶。这使得结核杆菌中的蛋白激酶成为一个潜在的药物靶点。

PknB的催化结构域是细菌激酶中第一个得到结构的，其活性构象显示出与原核生物的Ser/Thr激酶相似的特征，这一结构在其他细菌的激酶家族成员中也相当的保守。

发明内容

本发明的第一个目的是为了能够既快速又准确地从大量已有的DFG-in构象来获得全原子的DFG-out构象模型，供II型抑制剂的筛选和分子对接使用。

本发明的第二个目的是为了虚拟筛选出能抑制结核杆菌PknB激酶活性的小分子，从而找到新的抗结核病药物。

本发明利用活性链段重塑的构象预测技术，对于每个蛋白激酶大约可产生31-55%的DFG-out预测构象模型。这些DFG-out模型与其相应的X射线衍射晶体结构有很高的相似性。

本发明还提供了蛋白激酶活性状态的分类方法，按照DFG中Asp和Phe残基相对于活性链段所处的位置可以从预测之后的结果中挑选出DFG-out构象的模型。

本发明还提供了分子对接后小分子构象聚类和筛选的方法。本发明预测得到的蛋白激酶DFG-out非活性构象能够很好的预测已知II型抑制剂的结合构象并从分子库中准确地识别出已知的II型抑制剂。

具体而言，本发明提供一种预测蛋白激酶非活性构象的计算模拟方法，其步骤为：

第一步，构建初始结构模型；

第二步，通过蛋白激酶的活性链段重塑，得到其各种可能的预测构象；

第三步，将得到预测构象按几何分类方法，进行DFG-in和DFG-out分类；

第四步，按活性位点的空间体积大小挑选预测出的DFG-out构象，构成可用于分子对接的非活性构象系综。

进一步，本发明提供一种挑选II型抑制剂对接构象的方法，其步骤为：

第一步，按照前面所述的方法来预测蛋白激酶的非活性构象，与小分子进行分子对接；

第二步，通过分子对接构象的分析，排除一部分没有正确对接到活性位点的预测结果；

第三步，通过构象聚类和对接能量排序来实现II型抑制剂对接构象的筛选。

本发明还获得2个从Calbiochem筛选出来的对结核杆菌生长起抑制作用的小分子。抑菌实验证明它们具有明显的抑制结核杆菌生长的作用。

下面对本发明内容进行具体描述。

一、预测蛋白激酶非活性构象的计算模拟方法

从PDB（Protein Data Bank）数据库（www.rcsb.org）中下载附表1中这些既有DFG-in又有DFG-out构象的蛋白激酶PDB文件。对所得的PDB文件进行预处理，去除其中的水分子或其他分子；如果是多聚体则取其中的一个单体。如果晶体结构中含有磷酸基团的话，用MODELLER同源模建程序将其突变为原始残基。对于缺少残基的蛋白激酶也同样采用MODELLER程序以自身为模板进行补齐。然后将所得的PDB文件用Rosetta软件的relax功能进行优化得到初步的起始模型。

将优化后的模型与DFG-out晶体结构叠加后发现，两者之间在N-lobe之间往往存在结构上的不同。如附图 1A所示，蛋白激酶LCK的DFG-out晶体结构较之优化模型有明显的向外旋转的变化，使得它在II型抑制剂的活性位点上有更大的空间。对于蛋白激酶来说，DFG-flip的发生以及II型抑制剂的结合都需要N-lobe的外旋；而在DFG-in的晶体结构中，N-lobe，特别是其中的αC螺旋会更多的内旋，形成更紧密的结构，使得活性位点变得狭小。

这个构象变化可以被看作是N-lobe绕一个轴旋转所致。如附图1B所示，这个轴被定义为垂直于一个包含催化域坐标中心和αC螺旋的首尾残基在内的平面，并且穿过连接N-lobe和C-lobe的铰链区中点残基的一条直线。然后DFG-out构象的N-lobe相对于相应的DFG-in构象向外旋转了大约5-15度。为了定量地描述由于N-lobe旋转引起的激酶活性位点的大小，需要使用4个保守残基两两之间的4对距离之和。这4个保守残基是形成盐桥的Lys273和Glu288以及DFG中的Asp、Phe（残基编号取自PDB：3LCK）；这4对距离是Lys273到Asp和Phe的距离以及Glu288到Asp和Phe的距离。如果这4个距离之和小于30 ?，则认为活性位点过小；如果和大于32 ?，则认为活性位点足够大了；介于30-32 ?的情况则认为需要略微的放大。基于上述的观察，将N-lobe需要旋转的情况分为以下3类：上述的距离之和大于32 ?的N-lobe不需要旋转；介于30-32 ?的将其N-lobe外旋5度；小于30 ?的需要外旋15度 （附表 2）。在旋转之后，初始结构模型的N-lobe与晶体DFG-out结构的就非常相似了（附图 1C）。

为了预测激酶的DFG-out构象，我们提出的活性链段重塑方法是通过Rosetta 3软件的的“loop_relax”功能来进行具体实施的。Rosetta 3是一款获得公认的蛋白质结构预测软件。它的loopmodeling功能可以设定目标蛋白质的若干片断为对象，对其进行重新的结构预测，运用模拟退火等算法寻找目标片断的最低能量构象。

为了进行活性链段重塑，我们首先需要从蛋白激酶的FASTA序列得到其fragments文件，其中包含了蛋白质片段的结构信息；然后设定需要重塑的链段在蛋白质中的位置，其中设定的链段对象具体是以DFG基序之前第二个残基为起始，蛋白激酶活性链段的C末端为结尾的；以这段序列中位于中点附近最亲水的残基为剪切断点；最后设定文件所在的路径以及默认的参数之后递交计算机进行计算。在活性链段重塑之后，对结果的侧链进行全原子优化来获得最终的全原子最低能量构象。流程如附图2。

对于每个激酶来说，将独立地产生200个全原子最低能量模型。在一台普通Intel Pentium IV电脑上产生一个这样的模型需要3-4个CPU时。对于附表1中的每个蛋白激酶来说，大约有31-55%的预测模型为DFG-out构象。这些DFG-out模型与其相应的X射线衍射晶体结构用TM-align程序进行了对比，对比产生的Cα-RMSD和TM-score值就列于附表3。附表 3的结果显示预测模型相对于晶体结构的RMSD平均值都低于3 ?而TM-score的平均值都大于0.85，表明预测出的DFG-out模型与晶体结构具有很高的相似性。

和DFG-out分类

对于获得的模型，首先要将它们按DFG-in、DFG-out和过渡构象进行分类。

根据DFG-flip的定义，Asp和Phe的位置几乎翻转了180°，也就是说在DFG-flip的前后Asp和Phe分别位于主链的两侧，如附图 3所示。这样的话，通过将得到的模型与初始的DFG-in结构进行叠加，如果模型的Asp和Phe残基和相应的DFG-in结构中的这2个残基都位于同侧的话，那么就认为它也是DFG-in构象的；如果情况正好相反的话，那么就把它归为DFG-out那类；其它所有的模型都认为是过渡构象。而残基的位置就通过构建向量的方法来进行判别。

最终如附图 3B所示取出模型和初始DFG-in结构中各4个点的坐标后，即Asp的C_γ原子、Asp的C_α原子、Phe382的C_α原子和Phe382的C_γ原子，分别标为R1、R2、R3、R4和R1'、R2'、R3'、R4'。然后我们就能得到8个向量，即

r₂₁=R1-R2；r₂₃=R3-R2；r₃₂=R2-R3；r₃₄=R4-R3；

r₂₁'=R1'-R2'；r₂₃'=R3'-R2'；r₃₂'=R2'-R3'；r₃₄'=R4'-R3'。

以r₂₁、r₂₃、r₂₁'、r₂₃'为例，我们注意到r₂₁×r₂₃后所得向量的方向是向下的，而r₂₁'×r₂₃'后所得向量的方向是向上的，所以叉乘得到的2个向量夹角必然为一钝角，即这两个向量的点乘小于0。因此我们还需要将上述8个向量叉乘后得到4个向量：

p₁=r₂₁×r₂₃    p₂=r₃₄×r₃₂    p₁'=r₂₁'×r₂₃'     p₂'=r₃₄'×r₃₂'

最终，将所有模型和已有的DFG-out构象叠加，并取出上述8个点，计算各个向量并规定：若p₁·p₁'< 0且p₂·p₂' < 0，那么就认为模型是一个DFG-out构象；若p₁·p₁' > 0且p₂·p₂' > 0，那么就认为模型是一个DFG-in构象；其它的就认为是过渡构象。

模型的选择

为了进一步分析已知DFG-out结构中活性位点的疏水腔，使用PASS（Putative Active Sites with Spheres)程序来探测活性位点附近可能的结合位点。PASS的结果显示在DFG-out结构的活性位点存在有3个结合口袋。例如，如附图 1D所示，ABL1的DFG-out结构的三个结合口袋中心被分别显示为B₁、B₂和B₃。这三个口袋也同时被Seeliger et al. 分别定义为疏水区I、腺嘌呤口袋以及异构位点。为了定量的描述这3个口袋，使用LIGSITE程序来计算口袋中相距1 ?的格点数。LIGSITE所得格点数的多少直接与腔的大小有关。为了有个直观的认识，简单地将格点数转化为1 g?ml^-1密度的水分子个数，即将活性口袋中的格点放入一个足够大的水盒子中，并且认为任何离格点最近距离大于1.6 ?的水分子是未被占据的。一些已知的DFG-out结构中占据上述活性口袋的水分子个数列于附表 4（一些在附表 1中列出的DFG-out晶体结构由于其活性链段上有缺失的残基，所以它们活性口袋中占据的水分子数未在附表 4中列出）。水分子数据同样说明一个足够大的活性位点疏水腔对于II型抑制剂的结合来说是至关重要的。经过几次的测试，确定使用活性腔中能包含大于20个水分子为标准来选择DFG-out模型，构成可用于分子对接的非活性构象系综。并用于后续的抑制剂结合构象和虚拟筛选的试验。

上述步骤1--4的流程图见附图 4。

型抑制剂对接构象的方法

第一步，DFG-out模型系综分子对接

本发明使用AutoDock 4.2来进行抑制剂结合构象预测和虚拟筛选中的分子对接。将之前的B₂位点定义为分子对接的中心，对接的格点盒子大小为60×60×60 ?。对于每个蛋白激酶的每个DFG-out模型与其对应的II型抑制剂小分子计算20次对接结果，其他关于分子对接算法的参数都按照默认的进行设置。最终，每个配体-DFG-out构象对有20个对接构象产生。它们都要被集中起来进行下一步的构象选择分析。

由于分子对接软件将蛋白受体作为刚性分子来处理，且对于每个蛋白激酶来说，预测出的抑制剂结合构象可能是非常多的，所以单纯通过AutoDock的打分来判断难以得到正确的结果。首先要排除一部分没有正确对接到活性位点的预测结果。对于一个抑制剂的预测构象来说，如果其中没有一个原子与之前定义 的B₂或B₃距离在3 ?以内，则它就被认为是位于活性位点之外而被排除。剩下的构象用聚类的方法，即使用一个较大的RMSD（均方根偏差，是常用来表示两个结构间差异的参数）值阈将其进行大致的分类。每个类群中能量最低的构象就会被认为是这个类群的代表构象，而与它结合的DFG-out模型就被认为是相应的DFG-out代表模型。而这些代表构象最后再按照对接能量打分重新排序——拥有最低能量代表构象的类群就被指定为第一类群，以此类推。具体流程见附图 5。

对于附表1中的每个激酶-抑制剂对，分子对接的结果经过筛选和聚类会得到了1个、2个或3个类群及其代表构象。其中类群的个数取决于抑制剂在活性位点对接的取向：如果所有的构象都只有一个相近的取向，那么它们通常被归类在一起；如果有2个或3个类群存在的话，那么就表示存在2个或3个不同的对接取向。特别在2个类群的情况中，它们的取向往往是正好相反的。然后，将所有的代表构象与晶体复合物中抑制剂的结构进行了对比，其中重原子RMSD的结果列于附表 5中。由于对于大部分的配体来说，它们RMSD最小的结果就是第一类群的代表构象，那么其它类群代表构象的RMSD值就不再列出了。在附表 5所列的结果中，大部分代表构象相对于晶体结构的RMSD值都低于2 ?或在2 ?附近。对于一些抑制剂，存在2个或3个代表构象，其中一个与X射线衍射晶体结构相似，而有一个较小的RMSD值；其它代表构象与晶体结构的取向不同或相反而有一个较大的RMSD值。然而，在所有的结果中，只有一个激酶-抑制剂对（ABL1-PRC）的最低RMSD构象是第二类群中的代表构象。从这个角度来说，这些抑制剂的预测构象其取向都与晶体结构能很好的符合，准确率接近96%。这点也可以从附图 6抑制剂预测构象与晶体结构叠加图中看出。

第三步，药物的虚拟筛选

将小分子库与蛋白激酶预测出来的非活性构象系综进行分子对接，将结果经上述步骤挑选和聚类，并将第一类群的代表构象作为该DFG-out预测构象-小分子对的对接结果，按代表构象的对接能量排序来确定可能的II型抑制剂。

方法验证时所用的小分子库含有750种抑制剂的Calbiohem分子库和已知的II型抑制剂，根据它们在最后排序表中的位置做相应的ROC（receiver operating characteristic curve，受试者工作特征曲线）图，ROC曲线下的面积即为AUC（Area Under the ROC Curve，即ROC图曲线下面积）值。最后所得的AUC值全部都大于0.9，说明预测得到的DFG-out模型对它们的II型抑制剂有很好的选择性，见附图 7。

附表说明

表 1  本发明中用于验证方法可靠性的7种蛋白激酶，其活性与非活性构象的PDB ID以及相应的II型抑制剂。

表 2  DFG-in初始模型所需作的3类外旋，按照4个保守残基的4对距离之和来确定其外旋的程度。

表 3  本发明方法预测所得7中蛋白激酶的非活性模型的统计及其与已知结构的叠加对比所得C_α-RMSD和TM-score值。

表 4  具有type-II抑制剂的DFG-out晶体复合物其活性位点所能容纳的水分子数目。

表 5  II型抑制剂的预测最低能量的代表构象与晶体结构对比的RMSD值。

表1

。

表2

。

表3

。

表 4

。

表 5

。

附图说明

图 1中， A：将蓝色的LCK DFG-out晶体结构与其DFG-in初始模型叠加后发现N-lobe部分存在明显外旋；标出的保守残基之间的距离可以用于描述II型抑制剂结合的疏水腔的大小；B：N-lobe绕一垂直于既定平面的轴旋转的示意图；C：旋转15°之后的初始DFG-in模型的N-lobe与DFG-out晶体结构叠加图，可以看出它们现在重合的很好；D：用PASS软件分析ABL1的一个DFG-out晶体结构后得到在活性位点附近的3个可结合位点，用圆球表示其中心。

图 2  活性链段重塑中Rosetta 3 Loopmodelling设置流程。

图 3中，  A：蛋白激酶DFG-in和DFG-out构象叠加之后的俯视图。深色的是DFG-in构象，浅色的是DFG-out构象。从这张图中可清晰地看出在2种构象中Asp和Phe分别位于主链的两侧。B，C：向量法图示。其中2种构象各取4个原子的坐标，即R1：Asp的C_γ原子；R2：Asp的C_α原子；R3：Phe的C_α原子；R4：Phe的C_γ原子。点1'、2'、3'、4'为另一构象中相应的原子。

图 4  蛋白激酶DFG-out模型预测流程图。

图 5  分子对接之后小分子构象处理流程图。

图 6  II型抑制剂预测得到的结合构象与晶体结构的叠加比较图。

图 7  本发明所提出方法得到的非活性构象模型用于虚拟筛选已知II型抑制剂所作的ROC图及AUC值。

图 8  预测得到的9种潜在的结核杆菌PknB抑制剂对细菌生长抑制作用的实验结果。其中第1、10列为对照组；A、B、C…为小分子的浓度梯度，按50, 50/2, 50/4, …, 50/128 μg/ml。

具体实施方式

本发明的方法应用于结核杆菌蛋白激酶PknB的计算流程：

第一步，激酶非活性构象系综的构建：

a)  下载PknB蛋白激酶的PDB文件1O6Y，去除其中的杂原子，补齐后进行初步优化；

b)  按前所述，提取4个原子的坐标计算所需距离之和，得到其值为29.64 ?，说明活性链段之前1O6Y的N-lobe需要向外旋转15度；

c)  旋转之后得到所需初始构象，设定各参数文件后进行活性链段重塑的计算；

d) 在得到的200个预测模型中有163个经过向量法判断为DFG-out构象；

e)  将这些构象再做活性位点空间大小的检测，其中有8个构象符合要求，可以用于PknB可能II型抑制剂的虚拟筛选；

第二步，潜在PknB II型抑制剂的虚拟筛选：

a)  将这8个预测模型与Calbiochem抑制剂数据库中的小分子一一做分子对接，每个构象-小分子对产生20个结果；

b) 排除所有这些结果中没有对接到活性位点的构象，然后作聚类与排序，得到所有结果代表构象的对接能量排序表；

c) 拿出最低能量的20个小分子，仔细查看后确定其中9个为实验对象，做随后的抑菌效能实验；

第三步，检验结核杆菌PknB蛋白激酶潜在抑制剂的抑菌作用

a) 培养结核杆菌至对数生长期；

b) 将小分子化合物溶于DMSO，浓度为10mg/ml；

c) 取结核杆菌接种至细菌培养板，其中每孔接种10⁵个结核杆菌；

d) 将化合物按浓度梯度稀释，形成终浓度50, 50/2, 50/4, …, 50/128 μg/ml；

e) 将培养板放至37℃培养箱，生长3-4周，在此期间，每周观察结果，记录生长情况；

f)  根据生长情况，确定最小抑菌浓度。

最后我们发现其中有2种存在明显的抑菌作用，见附图 8。