离子型GABA受体新型拮抗剂的虚拟筛选

赵 锦，王雯洁，王晨晨，巨修练

(1.湖北轻工职业技术学院，湖北 武汉 430070；2.辉瑞(武汉)研究开发有限责任公司，湖北 武汉 430000；3.武汉工程大学化工与制药学院，湖北 武汉 430205)

γ-氨基丁酸(γ-aminobutyric acid，GABA)是中枢神经系统中重要的抑制性递质，可作用于细胞膜上的GABA门控氯离子通道(GABA-gated chloride channel，GABACl)，使氯离子内流导致细胞膜超极化，从而抑制生物体的神经活动[1-3]。GABACl是一种五聚体蛋白，又被称为离子型GABA受体(以下简称GABA受体)，广泛存在于昆虫的神经和肌肉细胞中。GABA受体的拮抗剂已有很多被开发成农用杀虫剂，如氟虫腈等，但随着该类杀虫剂的大量使用，出现了害虫抗药性和环境不友好等问题[1-3]。近年来报道的两类新型GABA受体拮抗剂：异噁唑啉类化合物(如Fluralaner)[3-5]和间二酰胺类化合物(如BPB-1)[6-7]，它们可以作用于GABA受体的一个新位点，同时选择性高，对哺乳动物低毒甚至无毒，较传统杀虫剂更安全[8]。因此，基于GABA受体的这一新作用位点研究Fluralaner等异噁唑啉类化合物的作用机制，可以为合理设计新型安全高效杀虫剂提供重要的理论依据。

探寻与Fluralaner具有相似作用机制和活性结构的新化学实体具有重要意义。作者采用分子模拟和计算机辅助药物设计方法，以序列一致性较高的人GABA受体为模板，同源模建黑腹果蝇(Drosophilamelanogaster)GABA受体的三维结构，探索Fluralaner的作用模式，并通过药效团模型和分子对接的虚拟筛选探寻GABA受体拮抗剂的新型苗头化合物。

1 实验

1.1 药效团模型的构建和验证

选择10个结构多样且活性高的Fluralaner衍生物(表1)用于药效团模型的构建[9]，运用SYBYL-X 2.1软件Sketch模块构建这10个化合物，并通过Compute命令对化合物的分子结构进行能量优化，获取合理的分子构象。具体方式如下：设定Tripos立场，添加Gasteiger-Huckel电荷，终止时的梯度差设为0.005kcal·mol-1·Å-1，最大重复次数设为1 000，其余参数保持不变。药效团模型的构建通过SYBYL-X 2.1软件Galahad模块完成，参数设置为系统默认。

表1 10个Fluralaner衍生物的化学结构及其对草地贪夜蛾的杀虫活性(LC50)

为了验证构建的药效团模型是否具有良好的识别活性化合物的能力，采用诱饵集方法计算阳性分子富集因子(enrichment factor，EF)来评价模型的可靠性，计算公式如下：

式中：TP为命中阳性分子的数量；n为所有命中分子的数量；A为所有阳性分子的数量；N为数据库的分子数量。

诱饵集方法的实验步骤如下：从ZINC数据库中选择338个无活性的化合物分子，同时选出12个已知活性但未用于构建药效团模型的分子共同作为阳性参照分子，然后将这350个分子进行能量优化后构建一个小型数据库用于虚拟筛选。通过SYBYL-X 2.1软件Unity模块将构建的药效团模型生成提问式，对该数据库进行搜索，计算富集因子。

1.2 同源模建

1.2.1 模型构建

从UniProt数据库中下载黑腹果蝇GABA受体β亚基序列(编号：P25123)，通过序列比对发现人GABA受体(PDB ID：4COF)β亚基与其序列一致性最高，为46.28%，故选用人GABA受体蛋白4COF作为模板。在同源模建前删去模板蛋白小分子配体，通过Swiss-model网站自动执行黑腹果蝇GABA受体的同源模建。

1.2.2 蛋白优化

采用GROMACS 2016.5软件对模建的黑腹果蝇GABA受体蛋白进行3 ns的自由模拟以消除不合理的键长、键角和扭转角。利用pdb2gmx命令在AMBER99SB立场下生成蛋白的拓扑文件[10]，然后将模型置于填充有TIP3P水模型和氯离子的立方体盒子中。为消除不合理的范德华力，通过多次迭代计算使模型能量收敛到10 kJ·mol-1，接着使整个体系在300 K恒温条件下达到NVT和NPT平衡，最后进行3 ns时长的自由模拟[3，11]。

1.2.3 蛋白评价

通过SAVEs在线蛋白检测服务器(http://services.mbi.ucla.edu/saves/)和PROSA程序(https://prosa.services.came.sbg.ac.at/prosa.php)来评价蛋白模型[3]。利用PROCHECK方法进行拉氏图评价，同时结合Z-score分布图来验证蛋白模型的质量。

1.3 分子对接

运用SYBYL-X 2.1软件Surflex-Dock Geom模块进行分子对接实验。在对接之前需要进行蛋白质准备，在对动力学优化后的蛋白质进行结构分析之后，对其进行侧链修复、末端处理、加氢、加电荷等优化。同源模建的蛋白质三维结构中不包含共晶配体，因此使用残基模式来设定结合口袋。文献报道，Fluralaner等异噁唑啉类化合物与黑腹果蝇GABA受体相互作用的关键氨基酸可能是Ile218、Leu222、Gly277、Phe280、Val281[12-13]，因此选择黑腹果蝇GABA受体模型A链上的Ile218、Leu222和B链上的Gly277、Phe280、Val281附近3 Å以内的氨基酸来定义结合口袋[3]。

1.4 虚拟筛选

通过基于配体的药效团筛选和基于受体的分子对接对ZINC数据库中可购买的22 724 825个化合物进行虚拟筛选。首先将构建的药效团模型生成提问式，采用Flex Search对ZINC数据库进行初步筛选，然后将初筛所得的分子对接到模建蛋白的结合口袋中，结合Total_Score和Cscore值进行进一步筛选。观察经过两次筛选后的小分子与GABA受体的结合模式，优先考虑能与结合口袋中关键氨基酸形成相互作用的化合物。

2 结果与讨论

2.1 药效团模型

最优的药效团模型分子叠合图如图1所示。

：疏水中心 ：氢键受体 ：正氮原子

该模型具有3个氢键受体(AA_1、AA_6、AA_7)、5个疏水中心(HY_3、HY_9、HY_4、HY_2、HY_8)和1个正氮原子(NP_5)，代表了异噁唑啉类化合物的关键药效结构特征。计算表明，该模型具有相对高的Specificity值(5.394)和相对低的Energy值(-0.51 kcal·mol-1)，富集因子为15.4，远大于1，证明该药效团模型具有较好的活性分子识别能力，故利用药效团模型进行后续的虚拟筛选具有合理性。

2.2 同源模建

2.2.1 序列比对

同源模建是预测目标蛋白质合理三维结构的有效方法，模板的选择和序列比对的精确度是决定模建蛋白质质量的重要因素。黑腹果蝇GABA受体β亚基和人GABA受体β亚基的序列一致性为46.28%，符合建模要求。此外，在解析人GABA受体4COF的三维结构时，在TM3和TM4的细胞内环区域氨基酸Arg357～Ser564被序列SQPARAA替换[14]，故黑腹果蝇的GABA受体β亚基序列的相同位置也用SQPARAA替换。黑腹果蝇GABA受体β亚基和人GABA受体β亚基的序列比对如图2所示，整个跨膜区域(TM1～TM4)氨基酸残基的一致性和相似性较高，因此选择人GABA受体作为模板具有合理性。

图2 黑腹果蝇与人GABA受体β亚基的序列比对Fig.2 Sequence alignment of GABA receptor β subunits of Drosophila melanogaster and human

2.2.2 蛋白优化

同源模建后的蛋白质通过3 ns分子动力学模拟进行优化，在模拟过程中，黑腹果蝇GABA受体结构的势能变化和骨架碳原子(Cα)的均方根偏差(root mean square deviation，RMSD)如图3所示。

由图3a可以看出，当能量收敛到设置范围10 kJ·mol-1时，优化步数达到2 718步，此时整个蛋白质体系达到能量最小化。由图3b可以看出，在1 600 ps后，黑腹果蝇GABA受体RMSD值达到收敛，稳定在0.16 nm左右，这表明所得模型已基本达到稳定状态，可用于分子对接。

图3 黑腹果蝇GABA受体模型的3 ns分子动力学模拟过程中的势能(a)和骨架碳原子RMSD(b)变化Fig.3 Potential energy(a) and carbon skeleton RMSD(b) of Drosophila melanogaster GABA receptor model during 3 ns molecular dynamics simulation

2.2.3 蛋白评价

模建得到黑腹果蝇GABA受体的PROCHECK结果(图4左)表明，模型中落在最佳和允许区域的氨基酸残基总和为99.6%。因此，经过动力学优化后的黑腹果蝇GABA受体的二面角是合理的。描述模型质量的Z-score分布图(图4右)显示，黑腹果蝇GABA受体模型的Z-score值为-3.16，落在了已知结构蛋白的Z-score值形成的分布区内，这表明模型能量也是合理的。上述评价结果表明，构建的黑腹果蝇GABA受体模型结构和能量均较为合理，可运用于分子对接研究。

图4 黑腹果蝇GABA受体模建蛋白的拉氏图(左)和Z-score图(右)Fig.4 Ramachandran plot(left) and Z-score plot(right) of optimized GABA receptor model of Drosophila melanogaster

2.3 分子对接

Fluralaner与黑腹果蝇GABA受体的对接模式如图5所示。对接打分Total_Score值为5.47，Cscore值为5，表明其与黑腹果蝇GABA受体有较强的结合能力[15-16]。

从图5可以看出，Fluralaner垂直地插入到黑腹果蝇GABA受体模型TM1和TM3之间靠近细胞外膜的间隙区域(图5a)，其含有三氟甲基的脂肪尾链嵌入2个亚基之间，并通过与Ile218、Leu222、Leu250、Met256形成疏水作用来维持自身与蛋白质结合的稳定性[3]。而Fluralaner空间位阻较大的异噁唑啉环(B)和苯环(C)未能完全嵌入该口袋中(图5b)。此外，Fluralaner异噁唑啉环上的O原子和N原子均能与残基Ala285形成氢键作用(表2)，同时苯环(B)可与残基Phe280形成π-π相互作用(图5b)，这些作用可能是Fluralaner与GABA受体稳定结合的原因。

图5 Fluralaner与黑腹果蝇GABA受体的对接模式Fig.5 Docking mode of Fluralaner with Drosophila melanogaster GABA receptor

表2 Fluralaner与黑腹果蝇GABA受体结合形成的氢键

2.4 虚拟筛选

根据所建立的药效团模型和分子对接方法，对ZINC可购买小分子数据库进行3轮虚拟筛选。在第1轮筛选中，利用Unity模块将最优的药效团模型转化成提问式，对ZINC可购买化合物数据库进行搜索，得到18 886个化合物，然后以QFIT>50为标准筛选出115个化合物；在第2轮筛选中，通过SYBYL-X 2.1软件Surflex-Dock Geom模块将这115个化合物对接到黑腹果蝇GABA受体的结合口袋中，综合考虑候选化合物的Total_Score值和Cscore值，选择排名前20且打分高于Fluralaner的化合物作为备选；在第3轮筛选中，通过观察候选化合物与氨基酸残基Leu222、Met256、Phe280、Ala285的相互作用来判断对接模式是否具有合理性，从而剔除14个化合物，最终得到6个苗头化合物，其对接打分见表3,化学结构见图6。

图6 筛选得到的6个苗头化合物的化学结构Fig.6 Chemical structures of screened six hit compounds

表3 筛选得到的6个苗头化合物的对接打分

化合物N1与黑腹果蝇GABA受体的对接模式

如图7所示。

图7 化合物N1与黑腹果蝇GABA受体的对接模式Fig.7 Docking mode of compound N1 with Drosophila melanogaster GABA receptor

由图7可知，与Fluralaner类似，化合物N1垂直地插入黑腹果蝇GABA受体模型TM1和TM3之间(图7a)，但深度不及Fluralaner，这可能是因为N1苯环支链的体积比Fluralaner脂肪长链更大，容易与氨基酸发生碰撞，不利于深入结合口袋[3]。此外，化合物N1可以与Phe280形成π-π堆积作用，与氨基酸残基Ile218、Pro219、Met256、Asp273、Leu276、Val281形成疏水作用。化合物N1与Ala285并未形成氢键作用，但却可以与Gln215、Met256、Asn260等氨基酸形成氢键作用。

3 结论

针对黑腹果蝇离子型GABA受体靶标及其新型拮抗剂进行了系统分子模拟，构建了异噁唑啉类化合物Fluralaner的药效团模型，通过同源模建和分子对接探索了Fluralaner与黑腹果蝇GABA受体的可能作用模式。随后基于药效团以及分子对接的虚拟筛选，发现了6个潜在的黑腹果蝇GABA受体新型拮抗剂，为靶向GABA受体新型杀虫剂的设计和研究提供了重要理论依据。