通过分子对接筛选14-3-3σ天然产物抑制剂

Screening of 14-3-3σ Natural Product Inhibitors by Molecular Docking

摘要: 14-3-3s通过蛋白–蛋白相互作用来调控多种生物功能，如信号传导、新陈代谢、细胞生长及细胞增殖等。本研究选取14-3-3 σ 作为靶标来筛选蛋白–蛋白相互作用的天然产物抑制剂。通过分子对接方法从台湾天然产物数据库中筛选出可能的抑制剂，并通过分子动力学方法模拟抑制剂与14-3-3 σ 的复合物的结合特性。由MMGBSA方法计算结合自由能定量探讨抑制剂与4-3-3 σ 的结合强度。研究表明存在两种潜在的14-3-3 σ 抑制剂。

Abstract: 14-3-3s regulates a variety of biological functions through protein-protein interactions, such as signaling, metabolism, cell growth, and cell proliferation. In this study, 14-3-3 σ was selected as a target to screen for natural product inhibitors of protein-protein interactions. The possible inhib-itors were screened from the Taiwan natural product database by molecular docking and the binding models between the inhibitor and 14-3-3σ were understood by molecular dynamics sim-ulations. The binding free energy of inhibitor/4-3-3σ was quantitatively investigated by the MM-GBSA method. This study has shown that there are two potential 14-3-3 σ inhibitors.

1. 引言

生物体通过蛋白–蛋白相互作用来调控生理过程中的多种生物功能，如信号传导、新陈代谢、细胞生长及细胞增殖等 [1]。通过调节蛋白–蛋白相互作用可以实现对相关生理活动的人为调控，所以蛋白–蛋白相互作用作为新的药物靶标获得医药研究者的广泛关注 [2] [3] [4]。14-3-3s是研究蛋白–蛋白相互作用的模型 [5] [6]。

14-3-3s是广泛存在于真核生物中的信号蛋白，人类主要有七种亚型(a/β, ε, η, γ, τ/θ, σ和ζ) [7] ，这些不同的亚型在人体分布不均 [8]。14-3-3s各亚型之间除能形成同源或异源二聚体之外，还能与上千种的配体蛋白形成蛋白–蛋白相互作用 [9]。研究表明14-3-3s主要与磷酸丝氨酸或苏氨酸蛋白相互作用，主要结合模式有RSXpSXP (I型结合模式)和RXXXpSXP (II型结合模式) [10]。还有一种motif-pS/pT(X_1–2)-COOH结合模式(III结合模式) [11] [12]。此外，14-3-3s也与非磷酸化的蛋白作用，称为IV型结合模式，如14-3-3s与ExoS蛋白相互作用 [13] [14]。这些结合模式保证了14-3-3s蛋白在人体中的广泛作用。

14-3-3s涉及多种疾病，如肺癌 [15] ，CNS疾病 [16] 以及神经系统疾病 [17]。14-3-3s作为治疗这些疾病的新靶标引起了药物研究者的广泛兴趣 [1] [16] [18]。14-3-3s单体由九个a螺旋构成一个杯状两亲性的高度保守的活性空腔 [19] [20]。针对14-3-3s靶标的研究主要集中在两方面。一方面设计小分子来稳定蛋白–蛋白相互作用，如fusicoccin稳定14-3-3s与PMA2之间的相互作用 [21] [22] [23]。pyrrolidone 1和epibestatin稳定14-3-3s与PMA2之间的相互作用 [24]。Cotylenin A能有效稳定14-3-3s与H⁺-ATPase复合物 [23] [25]。另一方面设计抑制剂分子破坏蛋白–蛋白相互作用 [5] [6] [26] [27]。如Thiel等人就根据14-3-3s与蛋白结合模式在ZINC数据库中筛选了含有磷酸根的分子作为抑制剂破坏蛋白–蛋白相互作用 [28]。

天然产物作为药物分子有着其独特的生理相容性，已经广泛应用于药物开发。对于14-3-3s而言，天然产物的稳定剂在实验阶段获得了成功 [21] [22] [23] [29] [30] ，而抑制剂的研究有待进一步研究。

在本项工作中，我们使用了分子对接从台湾中医天然产物数据库TCM [31] 筛选14-3-3s抑制剂。人体中七种14-3-3s亚型在一级结构和三级结构都是高度保守的，因此我们选择14-3-3s作为14-3-3s的模型蛋白。通过分子动力学定性研究了抑制剂与14-3-3s的结合模式，结合自由能定量研究结合强度。通过计算研究发现了两种潜在的14-3-3s抑制剂。

2. 计算方法

2.1. 分子对接

台湾中医化合物数据库TCM作为全新的先导化合物数据库已经得到了广泛的使用 [32] [33]。该数据库主要包含了从453种中药材中提取的超过20,000种的小分子化合物。

本研究选取14-3-3s晶体结构中的单体作为初始结构和受体蛋白(PDB编码为1YZ5) [34]。14-3-3s单体中缺失的氨基酸残基通过手动的方式补全，通过AutoDockTools软件 [35] 生成受体分子的pdbqt文件。在对接中设置14-3-3s的K122的C_原子作为盒子中心，盒子大小设置为30 Å ´ 30 Å ´ 30 Å。每个分子都对接十个合理的结合构象，并选取能量最低的构象作为排序构象。所有对接都采用AutoDock Vina软件 [36] 实现。

从对接结果中选出能量排前15个的化合物，根据结构的相似性对化合物归类。归类规则如下：结构相似手性不同的配体分子归为一类；基本结构单元一致不同取代位置和不同取代基的化合物归为一类。根据分类规则将前15个天然产物分子归为五类，在每一类中选能量最低的化合物代表该类分子。

2.2. 分子动力学

通过分子对接构建五种抑制剂与14-3-3s的复合物结构。将复合物体系加TIP3P [37] 水盒子，水盒子边缘距离复合物最近的距离为12 Å。通过加抗衡离子将系统的整体电荷保持为电中性。溶剂化后的体系包含了232个氨基酸和一个配体小分子，并有约15,350个水分子，构成79 Å ´ 85 Å ´ 91 Å的系统。14-3-3s使用ff14SB力场 [38] ，配体小分子没有标准的力场参数，使用GAFF力场 [39] 来生成相关的力场。使用Gaussian 09软件 [40] 在B3LYP/6-311G水平 [41] [42] 下优化配体分子。随后使用RESP法计算配体分子的原子电荷，所有配体分子的力场参数均通过Antechamber程序 [43] 实现。

溶剂化后的体系可能局部不稳定，因此优先对体系进行优化。首先使用10,000步的最陡梯度法优化水盒子，接着在固定蛋白主链的情况下使用9000步的最陡下降法和1000步的共轭梯度法来优化非蛋白主链部分。最后使用10,000步的共轭梯度法优化整个体系。前两步优化过程中添加的限制条件为200 kcal/mol∙Ų，第三步优化不加任何限制条件。优化过的体系通过200 ps将温度从0 K升到300 K。在300 K使用200 ps改变体系的体积维持体系的压强为1 atm。在300 K和1 atm下200 ps平衡让体系整体达到稳定。最后运行300 ns动力学用于数据分析。在动力学过程中使用PME方法 [44] 描述长程相互作用。为了有效减少计算量使用SHAKE方法 [45] 处理共价连接的氢原子。在计算过程中步长设置为0.02 fs，并以1 ps间隔保存动力学轨迹。所有的计算都基于AMBER14 Tools程序 [46] 实现。

2.3. 结合自由能计算

通过结合自由能的计算评价抑制剂与14-3-3s结合强度。这项工作中选用了隐性水溶剂模型的MMGBSA方法来计算。关于MMGBSA方法在文献中 [47] - [52] 已经进行了详细的介绍，这里简单的描述如下：

$\Delta G_{binding}=\Delta G_{14-3-3\sigma /inhibitor}-\Delta G_{14-3-3\sigma }-\Delta G_{inhibitor}$ (1)

$G=E_{gas}+G_{sol}-TS$ (2)

$E_{gas}=E_{\mathrm{int}}+E_{vdW}+E_{ele}$ (3)

$G_{sol}=G_{GB}+G_{np}$ (4)

$G_{np}=\gamma \cdot SA+b$ (5)

其中 $\Delta G_{14-3-3\sigma /inhibitor}$ ， $\Delta G_{14-3-3\sigma }$ 和 $\Delta G_{inhibitor}$ 分别表示抑制剂与14-3-3s复合物，14-3-3s和抑制剂的能量变化。其中式(3)需要的各项能量( $E_{\mathrm{int}},E_{vdW}$ 和 $E_{ele}$ )可以直接通过提取动力学轨迹中结构并分析各项能量获取。而极性溶剂化自由能贡献部分( $G_{GB}$ )可以通过Generalized Born (GB)方程 [49] 来计算。非极性溶剂化能( $G_{np}$ )可以通过式(5)来计算，其中γ = 0.0072 kcal/Ų，b = 0.0 kcal/mol。SA表示可接触表面积，可以通过程序MSMS [53] 来计算。对于抑制剂与14-3-3s形成的复合物，从动力学轨迹中提取出1000帧用来计算焓贡献部分。

平动，转动和振动相关的熵贡献( $\Delta S$ )是质量和分子内部的偶极矩的函数，可以通过统计平均来计算。在该工作中，通过normal mode分析法 [54] 计算。由于系统大小合适，计算采用100帧结构来计算熵变化。结合自由能中焓和熵计算使用的构象都是从动力学最后200 ns轨迹中平均采样获得。最后获取的能量值是将每一帧的各项能量统计平均得到，从而计算得到了结合自由能。

3. 结果与讨论

3.1. 筛选结果分析

使用AutoDock Vina程序 [36] 可以快速实现小分子与靶标蛋白的筛选。该项工作中筛选出对接结合自由能排名为前15的化合物，这些化合物的对接结合自由能都低于−10.7 kcal/mol，表明这些化合物都能与14-3-3s稳定结合。当这些抑制剂占据14-3-3s的活性空腔之后就可以竞争性的破坏14-3-3s与配体蛋白的相互作用，从而影响相关生理过程。

为了有效降低后续的计算时间，将前15个化合物手动归为五类，并从这五类中分别选取对接结合自由能最低的化合物继续研究。这五个化合物的结构如图1所示。化合物85531185和95911592的结构中都存在一个大环。而化合物85593830中包含有很多含氮基团，化合物85530919和85532194则包含有较多的含氧基团，这些基团可能有助于与氨基酸残基之间形成氢键从而稳定彼此间的结合。筛选出的五种化合物都包含有大量的苯环，有助于与14-3-3s两亲性空腔之间形成强的疏水性相互作用。

Figure 1. The structures for the selected five natural compounds

图1. 虚拟筛选选择的天然产物结构

化合物85531185含有一个十六元的大环，而95911592则是一个二十一元的大环，另外三种化合物虽然没有这样的大环，不过有着直接相连的环状单元，这些环状单元导致五种化合物的构象能保持相对刚性。化合物85593830，85530919和85532194的可旋转键主要存在于环状单元之间的连接，这样就导致每个环状单元都能相对的适应14-3-3s的活性空腔。化合物中可旋转键有利于化合物与受体蛋白更好的结合。

3.2. 动力学稳定性

通过分子对接得到了五种天然产物与14-3-3s的复合物结构。为了考察配体小分子在14-3-3s空腔中的结合稳定性，分析了整个动力学过程中14-3-3s蛋白的主链重原子的均方根偏差(RMSD)，图2所示。由图2可知除了85530919之外的四个体系RMSD都基本保持稳定。85532194，85593830和95911592三个体系的稳定性非常好，在整个动力学过程中偏差不超过0.6 Å，RMSD值分别为2.87 ± 0.40 Å，1.98 ± 0.30 Å和4.28 ± 0.55 Å。其中95911592体系的RMSD值在初始阶段有小幅度的变化，在100 ns之后基本维持稳定。体系85532194和85593830中14-3-3s结构相对于初始结构变化较小，特别是85593830体系不到2 Å的偏差，这表明化合物85532194和85593830两个化合物对14-3-3s结构的影响小。而95911592体系稳定存在，却对14-3-3s的结构带来了较大的变化。85532194，85593830和95911592三种配体分子导致14-3-3s结构发生不同变化的原因可能是这三种配体分子体积上的差异。95911592可能是由于该受

Figure 2. The RMSD value although the MD times

图2. 动力学过程中RMSD值

体分子中存在一个较大的环，而环状一般都能保持相对的刚性，14-3-3s为了更好的包合该配体分子就需要更大的活性空腔。而85532194和85593830两者虽然也包含了很多的五元环，六元环和螺环，但是并没有导致体积急剧增加的大环，同时这两种化合物含有一定的可旋转键，这样就可以与受体活性空腔更好的诱导契合达到局部构型调整以保持较好的包合。85531185体系的RMSD值保持起伏不定，可以看到RMSD值为3.48 ± 0.76 Å。说明该配体与受体的结合不太稳定，配体在受体的结合位点有较大的移动导致受体分子的结构也有适当的调整，但是整体上来说还是维持在一定的区间。85530919体系随着动力学的进行，在150 ns之后RMSD快速升高，在200 ns左右达最大值，可能该配体分子与受体的结合不是太稳定。

为了获取动力学过程中天然产物与14-3-3s包合物的结构信息，我们从分子动力学轨迹最后的200 ns中间隔50 ns均匀提取五个结构。由图3可知85530919体系随着动力学的进行，天然产物分子从最初的包埋于受体活性空腔到最后逃离了该空腔。这个与RMSD反应的情况一致，也表明该受体分子不能与14-3-3s形成稳定的包合物。85531185小分子并没有处于受体分子活性空腔的中心，而是处于活性空腔中靠近a9，这种有限的结合模式可能导致配体分子在空腔中适当的移动，这与RMSD的起伏相吻合。85531185体系中抑制剂处于14-3-3s活性空腔中，这就表明该天然产物能作为14-3-3s的抑制剂。与该体系类似的还有95911592体系，由提取出来的结构可以发现抑制剂处于受体分子活性空腔中，并与a7和a9形成较好的包合作用，配体分子在空腔的位置保持稳定，这也解释了为何该复合物中14-3-3s蛋白的RMSD值基本维持稳定。85593830化合物在14-3-3s空腔中的包合不完全，可以明显观察到配体分子有部分处于活性空腔的外部，这样的包合模式不能有效参与竞争性抑制。而85532194化合物在动力学过程中与14-3-3s形成了稳定的包合物，而且该配体分子处于受体分子的活性空腔的底部，与形成空腔底部的a5和a7形成稳定的相互作用。

Figure 3. The frame structures for the five selected natural products binding with the 14-3-3 protein with the MD times

图3. 动力学过程中提取的复合物构象

通过14-3-3s蛋白主链重原子的RMSD值分析和动力学过程复合物结构的分析可以看到85531185和95911592两种天然产物中存在大环能被14-3-3s的空腔所包合。85593830配体分子的长度较长，很难被整个的活性空腔所包合。而85530919与85532194在结构上有相似性，但是与受体的包合作用差别较大，其中85530919化合物脱离了受体的活性空腔，而85532194化合物能与受体形成很好的复合物。导致这两者之间差别的原因可能是85532194分子中的氢键受体和供体比较分散。

通过分子对接和分子动力学可以发现五个分子的结合位点在14-3-3s活性空腔的靠近小口端这边，与之前实验发现含有磷酸基团的化合物都处于大口端不同 [28]。这表明可能磷酸类化合物的结合位点与其它化合物在于14-3-3s形成包合物过程有不同的驱动力。

3.3. 结合自由能

通过对选出的五个天然产物与14-3-3s形成的复合物体系的动力学研究发现，85593830和85530919两个体系的结合都不稳定，所以在接下来的计算过程中排除了这两个体系。通过分子动力学定性的研究外，我们也需要对结合强度进行定量的评估，以找到更加合理的抑制剂。这里对85532194，85531185和95911592三个复合物体系通过MMGBSA方法计算了配体分子与受体之间的结合自由能，表1~3。

Table 1. Binding free energy for 14-3-3s/85532194 complex and decomposition to electrostatic interaction, van der Walls interaction, solvation free energies and entropy. Energies are in kcal/mol

表1. 14-3-3s/85532194复合物的结合自由能和能量分解为静电、范德华、溶剂化和熵效应。能量单位为kcal/mol

Table 2. Binding free energy for 14-3-3s/85531185 complex and decomposition to electrostatic interaction, van der Walls interaction, solvation free energies and entropy. Energies are in kcal/mol

表2. 14-3-3s/85531185复合物的结合自由能和能量分解为静电、范德华、溶剂化和熵效应。能量单位为kcal/mol

Table 3. Binding free energy for 14-3-3s/95911592 complex and decomposition to electrostatic interaction, van der Walls interaction, solvation free energies and entropy. Energies are in kcal/mol

表3. 14-3-3s/95911592复合物的结合自由能和能量分解为静电、范德华、溶剂化和熵效应。能量单位为kcal/mol

通过计算可以看到85532194体系的结合自由能为−2.26 kcal/mol，而85531185体系的结合自由能为−10.71 kcal/mol，是前者的4.7倍，这充分显示了后者与14-3-3s结合非常强。结合最强的为95911592体系，能达到−15.10 kcal/mol。85532194化合物结合自由能显示其不能很好的作为抑制剂被使用。相反另外两个大环的天然产物能与受体很好的结合，这与之前的研究相一致，即破坏蛋白–蛋白相互作用需要使用大大分子量的化合物。实验和理论研究表明植物14-3-3与模拟H+-ATPases的磷酸化五肽之间的结合强度为−8.3 kcal/mol [22] ，由此可以看到天然产物85531185和95911592可能会对该结合起到抑制作用。通过结构可以看到这里的两种可能的抑制剂分子也是存在于磷酸化五肽的结合位点，而不是稳定剂fusicoccin的结合位点 [25]。

结合自由能分解表明溶剂化作用对这些配体分子结合到14-3-3s两亲性空腔影响较大，这表明受体的两亲性活性空腔中的水分子被抑制剂分子竞争性排斥到空腔外相对比较难。将结合能量分解为极性部分和非极性部分可以看到这两体系中极性相互作用起着负作用，如在85531185和85532194两个体系做分别为12.92和22.51 kcal/mol，而非极性相互作用为−47.27和−51.68 kcal/mol，表明在这些配体分子结合过程中非极性相互作用是主要的贡献项。在后续的优化过程中需要考虑在一定程度上增加非极性的基团来增强疏水相互作用，而降低极性相互作用。

对比85532194和85531185体系可以发现两者最终的结合自由能差别非常大，可是对于熵效应贡献部分差别只有3.28 kcal/mol，而95911592体系的熵变化和85531185体系相差约为0.2 kcal/mol，这表明在这三个体系之间的不同主要是由于焓的不同引起。熵之间小的差异可以通过配体结构的刚性来解释，如85531185和95911592两个天然产物都是大环化合物，结构保持相对的刚性，而85532194中并联的四个环和支链的其它环也导致其结构刚性。85532194，85531185和95911592三个体系的焓贡献部分分别为−29.18l，−34.35和−38.82 kcal/mol，焓的绝对值大于熵的绝对值表明在结合过程中焓起到了主导作用，这些结合为焓驱动的结合。

通过结合自由能的计算表明85531185和95911592两个大环的天然产物可以作为潜在的14-3-3s抑制剂。而能量的分解也表明在结合中范德华相互作用起主要的作用，而抑制剂与受体的结合是焓驱动的反应。由于这些抑制剂都是刚性分子，在进行优化的时候可以将两个大环化合物适当的破坏，组成新的化合物保证更好的结合。

4. 结论

14-3-3s广泛存在于人体中，并且与癌症、神经系统紊乱等疾病相关。为了有效的治疗这些疾病，14-3-3s作为一种全新的靶标得到了广泛的研究，不过目前还没有相关的药物进入临床应用，所以急需相关的研究获取针对该靶标的抑制剂。本文通过分子对接从TCM数据库中来寻找天然产物抑制剂。通过分子对接筛选得到五个具有代表性的天然产物。为了研究这些天然产物化合物与14-3-3s的结合特性，我们进行了动力学研究。根据结合特征发现其中有两种分子无法稳定的结合在14-3-3s的两亲性空腔中。使用MMGBSA方法定量的研究三个配体与受体的结合强度，根据计算可以看到其中95911592的结合最强为−15.10 kcal/mol，表明该化合物可以作为潜在的14-3-3s的竞争性抑制剂分子先导化合物使用。通过我们的计算发现了85531185和95911592两种可能的14-3-3s抑制剂，不过这里没有进行相关的生物活性的测试，为了获得更加可靠的数据还需要更多的实验研究。