抑制剂RN-486与BTK靶标作用机理的分子动力学研究

Insight into Interaction Mechanism of Inhibitor RN-486 to BTK Target by Using Molecular Dynamics Simulation

摘要: 布鲁顿式酪氨酸激酶(Bruton’s tyrosine kinase)是近来治疗B细胞恶性肿瘤的明星靶点，在很多疾病和信号通路中起着关键作用，现已成为药物的研发热点。搞清楚BTK-抑制剂之前的相互作用及关键作用残基对药物的研发是十分必要的，本文选择罗氏药业公司开发的BTK抑制剂(研发代号为RN-486)，使用MM/GBSA方法以及正则模式熵来计算靶标与抑制剂的结合自由，同时采用基于残基的自由能分解方法评估主要残基对抑制剂结合的贡献。结果发现，RN-486与靶标结合的能量主要以疏水相互作用和静电相互作用为主，其中，LEU408和VAL416为两个重要的疏水残基，与RN-486发生范德华相互作用，而LYS430、TYR476和ASP539三个带电残基与RN-486发生静电相互作用。结合计算数据，我们找到了抑制剂与靶标的主要结合位点与作用形式，这为后期继续设计和开发BTK靶标的抑制剂提供根据指向性的理论指导。

Abstract: Bruton’s tyrosine kinase (BTK) is becoming an increasingly attractive target for drug discovery due to its critical role in multiple pathways and a variety of diseases. Understanding the binding mech-anism and key residues in BTK-inhibitor interaction is crucial for effective drug discovery. In this study, a BTK inhibitor developed by Roche (RN-486) was selected. The MM/GBSA method and the normal mode entropy were used to calculate the binding free energy between the target and the inhibitor. Meanwhile, Residue-based free energy decomposition method was used to reveal the contribution of the major residues. The results showed that the binding mainly from hydrophobic interaction and electrostatic interaction. Leu408 and Val416 were two important hydrophobic res-idues, which had van der Waals interaction with RN-486. The charged residues Lys430, Tyr476 and Asp539 interact with RN-486 electrically. The results indicate that the main binding sites and in-teraction mechanism between inhibitor and target, which are expected to provide theoretical guidance for the subsequent design and development of inhibitors targeting BTK.

1. 引言

布鲁顿氏酪氨酸激酶(Bruton’s tyrosine kinase)是近期药物研发的热点，在治疗B细胞恶性肿瘤上非常关键，是非受体酪氨酸激酶TEC激酶家族的成员之一，几乎在所有的血液细胞中都会发挥其调控作用 [1]。

BTK是B细胞抗原受体(BCR)信号通路的重要组成部分，也是与B细胞正常生长、成熟紧密关联的重要分子。B细胞抗原受体(BCR)及其前体(Pre-BCR)可调节与控制B细胞的稳态、分化及其功能 [2]。BTK结构包括5个区域，从N端到C端依次为：PH、TH、SH2、SH3、Kinase [3]。2007年，从Celera基因组学的科学家报告的中可知，一些小分子可以通过不可逆地结合到其ATP结合域中的Cys481残基而使BTK失活 [4]。首个被批准上市的BTK抑制剂Ibrutinib就是通过选择性地与Cys481位点发生有效且不可逆的结合，阻止BTK的磷酸化与活化，从而抑制BTK的活性。Ibrutinib对治疗与B细胞恶性增殖相关的恶性肿瘤如慢性淋巴细胞白血病、套细胞淋巴瘤、弥漫性大B细胞淋巴瘤等也有较好的疗效 [5]。罗氏药业公司开发的小分子抑制剂RN-486与Ibrutinib不同的是，Ibrutinib为共价抑制剂，而RN-486以非共价形式与BTK靶标结合 [6]。临床发现Cys481的突变后，抑制剂与靶标蛋白的共价键结合变弱，这样会导致药物产生耐药性 [7] [8] [9]，RN-486这类的非共价结合的抑制剂可以抵抗这种耐药性，那么它与BTK靶标的相互结合作用机制就十分值得探讨了，这也为后续开发和设计BTK靶向抑制剂提供更好的理论支持。

分子动力学(molecular dynamics, MD)模拟和结合自由能计算是现阶段研究抑制剂与靶标蛋白结合作用机制的主要工具，本文采用基于经验方程的MM/PB (GB) SA [10] [11] 方法以及基于残基的自由能分解方法来计算RN-486与BTK靶标的结合自由能，图1展示了小分子抑制剂的平面分子结构以及与BTK靶标的复合体结构。同时，希望此项工作能够为接下来设计非共价抑制剂提供有价值的结合位点信息和理论指导。

Figure 1. RN-486 with BTK complex structure chart and plane structure of RN-486

图1. RN-486与BTK的复合物结构图以及RN-486平面分子结构图

2. 理论和方法

2.1. 分子动力学模拟

在蛋白质库(PDB)中获取复合物晶体结构，PDBID为5P9G [12] 作为做分子动力学模拟的初始构象，本文使用的是AMBER力场进行分子动力学模拟，蛋白质和水分子力场参数由ff14SB产生，小分子抑制剂的力场参数由gaff力场生成，小分子的电荷处理以及参数生成都使用Amber Tools中Antechamber程序，并利用Tleap模块对复合物模型进行处理。复合体被溶解在显性的TIP3P水盒子中，水盒子边缘与复合体最近原子的距离设定为12 Å。在进行分子动力学模拟之前先约束溶质重原子，优化氢原子和水，在约束溶质主链原子，优化侧链和水，最后对体系进行整体优化、升温和一段时间的模拟平衡。体系在NVT系综中从0 K缓慢升温到300 K，用时1 ns；随后在NPT系综中进行平衡模拟，用时2 ns，模拟温度为300 K，压强为101 kPa；最后再在NPT系综中进行10 ns的分子动力学模拟，时间步长为2.0 fs，每100 fs保存一帧。在分子动力学模拟过程中采用SHAKE [13] 算法控制含氢化学键的伸缩，PME [14] 方法计算长程静电相互作用，非键相互作用的截断值为10 Å。

2.2. 结合自由能的计算

MM/GBSA (分子力学广义玻尔兹曼比表面积)方法，是端点自由能方法的一种，以分子动力学模拟采样和连续介质模型方法为基础去计算两个分子间的结合自由能 [15]。结合自由能焓和熵贡献的计算分别由AMBER16中的perl脚本mm_pbsa.pl和nmode模块完成，GB模型的“igb”设置为2，从最后2 ns中的轨迹中均匀的取100帧用于MM/GBSA的计算。计算过程中去掉构象中的水分子和离子，采用下面的公式计算小分子抑制剂和和靶标蛋白的结合自由能：

$\Delta G=\Delta E_{ele}+\Delta E_{vdw}+\Delta G_{pol}+\Delta G_{nonpol}-T\Delta S$ (1)

式中， $\Delta E_{ele}$ 和 $\Delta E_{vdw}$ 分别表示气相中小分子抑制剂与靶标蛋白的静电相互作用和范德华作用，第三项 $\Delta G_{pol}$ 表示极性溶解自由能对结合自由能的贡献，第四项 $\Delta G_{nonpol}$ 表示非极性溶剂对结合自由能的贡献，最后一项 $-T\Delta S$ 表示熵变对结合自由能的贡献，由传统热力学计算中的正则模式算出。第四项 $\Delta G_{nonpol}$ 非极性项由溶剂可及表面积(SASA)得到，有公式：

$\text{Δ}{G}_{nonpolar}=\gamma \Delta \text{SASA}+\beta$ (2)

式中，γ和β是常数，γ = 0.00542 kcal/(mol∙Å)，β = 0.92 kcal/mol。

3. 结果和讨论

3.1. 体系的稳定性与自由能计算

为了评估复合物体系动力学过程的稳定性以及结合自由能取样是否合理，我们以初始构象作为参考对象，计算了BTK蛋白主链的均方根偏差(root-mean-square displacement, RMSD)，并绘制了RMSD随时间变化的图像(图2)。如图2所示，经过10 ns的MD模拟，体系的RMSD最终稳定在RMSD的值波动范围在1.2 Å左右，方差和标准差分别为0.0138和0.1175，说明所有体系都处于相对平衡的状态，因此我们认为对这段轨迹采样计算蛋白配体结合自由能是合理的。

Figure 2. Molecular dynamics simulation of the root mean square deviation of 5P9G main chain atoms

图2. 分子动力学模拟5P9G主链原子的均方根偏差

3.2. 计算结果与分析

我们从最后2 ns中的轨迹中均匀的取100帧用于MM/GBSA的结合自由能计算，如表1所示，结合自由能包含5个部分，静电相互作用( $\Delta E_{ele}$ )、范德华相互作用( $\Delta E_{vdw}$ )、极性溶剂自由能( $\Delta G_{pol}$ )、非极性溶剂自由能( $\Delta G_{nonpol}$ )和熵变对结合自由能的贡献( $-T\Delta S$ )。其中，范德华相互作用、静电相互作用和非极性溶剂自由能都为负值，有利于小分子抑制剂和BTK靶标的结合，而极性溶剂自由能和熵变为正值，即不利于小分子抑制剂与靶标的结合。而气相的静电作用和液相的极性溶剂作用总体看来不利于结合自由能的增加，范德华相互作用比起非极性溶剂自由能又大得多，因此，范德华相互作用在小分子抑制剂与靶标的结合中起着更重要的作用。

Table 1. Binding Free Energy of RN-486 and BTK Target Calculated by MM/GBSA Method (in Kcal/mol)

表1. MM/GBSA方法计算的RN-486与BTK靶标的结合自由能(单位kcal/mol)

如表2所示，通过基于残基的自由能分解方法，将结合自由能大于1 kcal/mol的残基列出，共有11个热点残基对结合自由能贡献较大，将它们从大到小进行排列。LYS430在与rn486结合时自由能最大，可以看出发挥主要作用的是范德瓦尔斯作用。如图3所示，赖氨酸LYS430提供了含有氢原子的基团，与rn486结构中的电负性原子F产生了较强的氢键相互作用。除此之外，LYS430与rn486也产生了π-烷基相互作用，在作用部位，二者的距离分别为3.3 Å和5.2 Å。极性的LYS430、TYR476、ASP539与rn486相互作用所得溶剂能中的极化作用能对结合自由能有更多贡献，非极性的MET471、LEU408、VAL416与rn486相互作用所得的极化作用能与非极性能的大小则没有那么明显的差别。以典型的非极性氨基酸亮氨酸LEU408为例，如图4所示，它与rn486形成了两个较强的 -烷基相互作用，在作用部位，二者之间距离分别为3.9 Å和4.8 Å。

Table 2. Interactions of RN-486 with key residuals in BTK target (in kcal/mol)

表2. RN-486与BTK靶标中主要残基的相互作用(单位kcal/mol)

Figure 3. The hydrogen bond and π-alkyl interaction of Lys430 with RN-486

图3. LYS430与RN-486的氢键与π-烷基相互作用

Figure 4. The hydrophobic interaction of Lys430 with RN-486

图4. LEU408与RN-486的疏水相互作用

4. 结论

本文通过分子动力学模拟及结合自由能计算的方法研究了RN-486与BTK靶标的相互作用机制，计算得到结合自由能为36.85 kcal/mol，其中，范德华相互作用能和静电相互作用的数值最大，由此可知BTK与rn486相互作用的过程中，范德华相互作用和静电力相互作用占有主导地位。通过残基自由能分解的方法得到LYS430、MET471、TYR476和ASP539更易与rn486发生静电相互作用，LEU408和VAL416则更易与rn486发生范德华相互作用。明确rn486与热点残基占主导地位的相互作用后，可以此为依据修饰rn486结构，使其更易于与对应的热点残基相互作用，最大程度地释放结合自由能。

基金项目

内蒙古医科大学青年创新基金(YKD2018QNCX002)。