TNNI3K抑制剂3D-QSAR的研究及虚拟筛选

徐闻曦，蔡晓然，郑小娇，刘根炎，巨修练

武汉工程大学化工与制药学院，湖北 武汉 430205

肌钙蛋白I相关激酶（TNNI3K）是一种心肌特异性激酶［1］。它包含三个可识别的结构域：N-末端的十个锚蛋白重复序列，一个中心蛋白激酶结构域和C末端富丝氨酸结构域［2］。TNNI3K是促进分裂原活化蛋白激酶（MAP3K）的超家族成员［2］，然而目前仅有一个TNNI3K上游靶点被阐明，其余靶点尚不明确［3］。对TNNI3K在心脏生物学中的作用知之甚少。据Vagnozzi等［4］研究发现TNNI3K会加重缺血/再灌注损伤，氧化应激和心肌细胞死亡。TNNI3K介导的损伤通过产生更多线粒体超氧化物及破坏线粒体功能发生，并且在很大程度上依赖于p38丝裂原活化蛋白激酶（MAPK）激活。TNNI3K被抑制能够减少线粒体超氧化物产生，保留心脏功能并限制慢性不良重塑［4］。这项研究结果表明TNNI3K调节缺血性心脏再灌注损伤，是急性冠状动脉综合症（ACS）的治疗新靶点。并且TNNI3K仅在心脏中表达，所以靶向TNNI3K的抑制剂可能相对比较安全。

本文通过Brain等［5］研究得到的TNNI3K抑制剂构建CoMFA及TopomerCoMFA模型对TNNI3K抑制剂进行三维定量构效关系研究。并通过已构建的 TopomerCoMFA 模型结合 Topomer search［6］技术对ZINC数据库进行基于R基团的虚拟筛选，最终整合得到25个分子，这些分子预测活性均高于建模分子中活性最高者。将这些分子对接到TN⁃NI3K蛋白上，筛选得到11个分子可以作为潜在的TNNI3K抑制剂，为TNNI3K的优化设计提供了理论依据和进一步研究的基础。

1 实验部分

1.1 数据来源及其预处理

本研究所涉及的所有计算工作都在WIN7系统SYBYL-X 2.1软件上完成。除特别指出外，计算参数均为缺省值。本研究所用53个TNNI3K抑制剂均来自Brain等［5］文献，其中化合物的生物活性需先将文献中所给IC50（nmol/L）值转化为pIC50（mol/L）值。化合物结构与生物活性值如表1所示。

表1 化合物结构和生物活性数据Tab.1 Structure and biological activity data of compounds

续表1

续表1

续表1

按随机化原则，从53个TNNI3K抑制剂中选出43个分子作为训练集用于构建CoMFA及TopomerCoMFA模型，其余10个分子作为测试集用于检测模型的预测能力。通过Sketch Molecule模块绘制出全部53个分子结构，并通过Minimize模块利用Tripos力场进行能量最小化计算实现结构优化。其中指定电荷类型为Gasteiger-Huckel，将能量计算终止标准值缩小为0.005，最大重复次数增大到1 000。保存优化后分子用于后续计算。

1.2 CoMFA模型的构建

传统的CoMFA是3D-QSAR中应用最广泛的方法之一。本研究通过Align database模块，选用生物活性最高的第52号化合物为模板，通过43个分子的部分公共结构进行叠合（图1）。叠合完毕后程序自动计算所有分子的静电场及立体场参数，并将其作为自变量，以对应的pIC50值作为因变量，通过偏最小二乘法分析（Partial least squares）拟合分子场与生物活性间的关系生成模型［7］。进而通过抽一法交叉验证［8］检测模型的内部预测能力，并得到最佳主成分数N以及交叉验证相关系数q2，通过最佳主成分数运用非交叉验证法计算得到非交叉验证相关系数r2，显著性检验值F，以及标准误差SEE。并通过计算测试集外部预测相关系数（r2pred）预测模型预测能力［9-10］。

图1 分子叠合图Fig.1 Superposition of molecules

1.3 TopomerCoMFA模型构建

传统的CoMFA存在一些缺陷，其中最大的是其对输入的要求，一方面需要提供每个配体的3D结构，另一方面需要通过考虑所有其他配体结构来选取配体分子的某一构象进行适当的“叠合”［11］。而第二代CoMFA技术，即TopomerCoMFA，可自动创建用于预测化合物生物活性或性质的模型［12］。TopomerCoMFA模型可以在几分钟内创建，使用起来更方便，并且结果通常与传统的CoMFA结果相当且重复性高。本研究以生物活性最高的52号分子为模板。运用TopomerCoMFA模块将训练集中分子切割两次得到R1和R2两个R基团，并自动得到其三维构象，切割方式如图2所示。与CoMFA相同，切割完毕后通过偏最小二乘法分析得到TopomerCoMFA模型。

图2 52号分子的切割方式Fig.2 Cutting ways of molecule 52

1.4 运用Topomer Search进行虚拟筛选

Topomer search只是简单的从一个“标准”的R基团拓扑结构相似性搜索来定义命中结构，并运用TopomerCoMFA技术来计算每一个结构足够相似的R基团的PIC50值（当然也需要计算R基团的原子电荷和静电场）［13］。本研究通过TopomerCoM⁃FA构建的模型在ZINC（2012）中的Leads Now数据库对R1及R2基团进行虚拟筛选，该数据库中共包括3 687 621个分子，Topomer距离设置为150。

1.5 分子对接

通过Docking suite模块对新设计化合物与TNNI3K蛋白进行对接研究。受体晶体结构来自PDB数据库（PDB：4YHT），将其导入SYBYL中并通过BIOPOLYMER模块对其进行修补侧链、主链末端处理、加氢、设定残基质子化状态、给配体指定AMBER7-FF99的原子类型等蛋白准备工作。准备筛选得到的新化合物。设置对接模式为Surflex-Dock GeomX（SFXC），输出构象个数为 20，并选择进行CScore计算。一般认为输出构象的总打分函数Total Score大于6时构象被判定为较好的输出构象。当总打分函数和其他四种经验打分函数均满足要求时，CScore打分为5（总分为5）。一般认为CScore为5时构象具有良好的选择性。

2 结果与讨论

2.1 CoMFA

如表2，训练集分子所构建CoMFA模型的交叉验证相关系数为0.622，最佳组成分数n为6，非交叉验证相关系数r2为0.952，标准偏差SEE为0.211，F值为62.931。这些数据说明模型有较好的稳定性和内部预测能力。外部预测相关系数为0.823，说明模型同样具有良好的外部预测能力。如图3，对实验活性和预测活性做线性相关分析发现，测试集与训练集分子实测活性值与预测活性值具有较好的相关性（见图3），这进一步说明模型具有较好的预测能力。另外，值得关注的是，立体场贡献值为93%，而静场场贡献值仅为7%，说明立体场在该模型中发挥主导作用。

表2 CoMFA模型的统计参数Tab.2 Statistical parameters of CoMFA model

图3 采用（a）CoMFA和（b）Topomer CoMFA模型得到训练集和测试集分子实验值与预测值的线性回归图Fig.3 Linear regression between actual and predicted values of training and test sets by（a）CoMFA model，（b）Topomer CoMFA model

PLS分析的结果可以转换为X个变量的回归系数，它可以用于生物活性值的计算和预测。由于有大量的回归系数，所以直接解释相应的方程是不可能实现的。所以将结果以色块图的方式显示［8］。分子周围静电场显示为红蓝色块，红色代表增加此区域电负性可以使生物活性增加，而蓝色代表减少该区域电负性可以使生物活性增加。立体场显示为黄绿色块，绿色代表增加该区域基团体积可以增加活性，而黄色代表减小该区域基团体积可以增加活性。

以生物活性最高的52号化合物为模板分析CoMFA模型的三维等势图（见图4）。图4（a）为静电场等势图，如图所示苯环6号位有一蓝色区域，此处连有电负性弱的基团时分子活性更高，如氟、氯两个原子的电负性大小为F＞O＞Cl，所以苯环6号位连有氟原子的9号分子生物活性（pIC50=6.6 mol/L）小于连有羟基的20号分子生物活性（pIC50=7.4 mol/L），小于连有氯原子的10号分子生物活性（pIC50=7.7 mol/L）。蓝色色块附近有一小的红色色块，说明此处如连有多个基团时，为了避免空间位阻的影响，两个基团不应同时过大。另外，20号分子生物活性较好可能是由于与相邻苯胺上的N形成分子内氢键使得与TNNI3K结合所需的共面构象更稳定造成的。

图4（b）为立体场等势图，从图中可以看到在苯环6号位附近有一大的绿色色块，说明此处取代基体积越大活性越强，如34、35、36三个分子，随着苯环6位四氢吡咯上连有甲基数目依次增大，化合物体积依次增大，生物活性依次增大。苯环五位氨基左侧有多个黄绿色块，说明此处立体场对活性的影响较为复杂。1，3，5-哒嗪环左侧有多个黄色色块，说明此处连有体积较小基团时活性较好，如喹唑啉上连有氨基的42号分子（pIC50=7.5 mol/L）比连有氨甲基的41号分子（pIC50=7.3 mol/L）和连有氧甲基的43号分子（pIC50=7.1 mol/L）生物活性都要高。另外友谊绿色色块横跨52号分子左侧末端苯环，这也解释了53号分子（pIC50=7.3 mol/L）活性高于51号分子（pIC50=7.1 mol/L）的原因。

图4 CoMFA等势图：（a）静电场，（b）立体场Fig.4 Contour map of CoMFA：（a）electrostatic，（b）steric fields

2.2 Topomer CoMFA

TopomerCoMFA模型的最佳组成分数为13，交叉验证相关系数q2为0.768，q2的标准偏差为0.38，非交叉验证相关系数r2为0.981，r2的标准偏差为0.11，截距为7.8，计算TopomerCoMFA模型外部预测相关系数R2pred为0.754。作为以上参数的补充，图3反映了训练集及测试集分子生物活性的预测值与理论值间的差异及散点离散程度，图中各点基本分布在直线两侧。模型参数及线性回归图综合分析表明模型具有良好的稳定性及内部和外部预测能力。

2.3 新化合物筛选

基于得到的TopomerCoMFA模型，对ZINC数据库中Leads Now数据库（共包含36 876 219个分子）进行基于R基的虚拟筛选，筛选分为两个阶段进行。第一阶段，数据库中分子被切割成为R基团并通过Topomer相似性与TopomerCoMFA模型所基于的训练集R基团进行比较。第二阶段，对满足Topmer相似性要求的R基团通过它们对预测值的贡献来打分。再将筛选得到的R1及R2基团通过生物活性最高的52及33号分子进行过滤，得到贡献值最高的5个R1及5个R2。根据排列组合原理，依次替换模板分子52号分子中的R1及R2基团，最终得到25个新分子。用与训练集及测试集中分子相同的方式对这些分子进行结构优化构建新表单，放入TopomerCoMA模型中预测活性，发现设计得到的新分子预测活性均大于活性最高的模板分子（PIC50值均大于9.5）。

将这25个新设计的分子对接到TNNI3K蛋白上，由于分子对接和TopomerCoMFA两种方法评价原理差异，与TopomerCoMFA不同，分子对接会考虑疏水、极性、排斥、熵、溶剂化的综合作用，所以对接打分与TopomerCoMFA模型预测活性之间并不存在线性关系［14］。为提高筛选的成功率，本文选择CScore打分为5的分子作为对接评价优秀的分子，并选择Total Score最高的构象为每个分子的最佳构象。25个新化合物中共有11个化合物有CScore打分为5的构象（见表3），从这11个化合物中选择对接打分最高的2号化合物进行具体的对接研究，2号化合物与TNNI3K蛋白口袋对接如图5所示。从图中可以看到，Lys482、Val470、Leu513、Ala480、Phe582、Trp530形成疏水口袋与配体形成疏水作用，Trp53与配体形成π-π堆积作用。另外，值得注意的是对接打分最高的化合物2与4YHT晶体中原配体与TNNI3K形成的氢键相同，与Cys531、Asp593、Thr528形成氢键作用。研究这些化合物与TNNI3K蛋白的作用发现，11个化合物中有9个化合物与Asp593形成分子间氢键，并且两个不能和此蛋白形成此氢键的化合物对接打分较低，小于7（11个分子中有8个分子对接打分大于8），说明这3个氨基酸特别是Asp593在蛋白与配体结合过程中作用显著。

表3 新设计分子的结构，预测活性值以及对接打分Tab.3 Structure predicted pIC50and surflex-dock scores of new designed molecules

图5 2号分子与TNNI3K活性位点对接图Fig.5 Docking results of compound 2 with TNNI3K

3 结 语

本研究通过53个TNNI3K抑制剂构建了其CoMFA及TopomerCoMFA模型，模型参数及线性回归综合分析可知模型具有良好的稳定性和预测能力。将TopomerCoMFA模型作为Topomer Search提问式搜索zinc数据库，并将搜索得到的R1、R2基团组合共得到25个新化合物，结合分子对接最终筛选得到11个潜在的TNNI3K抑制剂，并且观察蛋白受体与配体相互作用发现，Cys531、Asp593、Thr528可能是配体与TNNI3K蛋白结合过程中产生关键氢键的氨基酸。对3D-QSAR模型的研究可以方便理解肌钙蛋白I相关激酶抑制剂结构与活性的关系，为设计新的具有更高活性的抑制剂分子提供信息。而TopomerCoMFA与分子对接技术的结合，可以使我们从大的小分子数据库中筛选得到潜在的抑制剂分子，提高抑制剂分子设计的效率。