基于Gromacs的蛋白分子识别研究

研究内容：本文采用分子对接和分子动力学模拟方法研究了APOE3、APOE4与VAMP2蛋白的分子识别过程。

1. 研究背景

APOE蛋白由317个氨基酸残基组成，其中1-18位的氨基酸为信号肽，经过翻译修饰后被切掉。具有生物活性的APOE蛋白由299个氨基酸残基组成。有三个主要的APOE等位基因变体，分别为APOE2,APOE3和APOE4。三个相应的基因等位基因编码的APOE同工型仅在位置112和158位彼此不同（APOE2中的Cys112和Cys158; APOE3中的Cys112和Arg158; APOE4中的Arg112和Arg158）。

VAMP2蛋白的主要功能是将突触小泡靶向和/或融合至靶膜。实验研究发现，VAMP2可以与三种APOE酶进行结合，但是结合的亲和力大小不同。因此为了深入研究VAMP2与APOE3、APOE4蛋白亲和力差异的主要原因，本文采用分子动力学模拟（Molecular dynamics simulation）方法对VAMP2与APOE蛋白的分子识别展开了研究。

1. 模拟方法

• APOE与Vamp2晶体结构的获取

从蛋白质数据库(Protein data bank, PDB)获取完整的APOE3蛋白晶体结构，PDB编号为2L7B[1]。APOE4与APOE3仅112位的氨基酸存在差异，APOE2与APOE3158位的氨基酸残基存在差异。因此先用PyMOL 2.1[2]将APOE3蛋白112位的Cys残基突变为Arg作为APOE4蛋白的初始结构，APOE3蛋白的158位的Arg残基突变为Cys作为APOE2蛋白的初始结构，随后用Gromacs软件分别对APOE4和APOE2蛋白进行20 ns的分子动力学模拟，以最终平衡时候的结构用于与Vamp2蛋白对接。

由于Vamp2蛋白没有完整的晶体结构，因此需要同源建模来构建其完整结构。用软件Modeller 9.24[3]进行单模板建模，Swiss-model[4]进行模板搜索，选择PDB编号为3HD7[5]的蛋白作为建模的模板，同源性为100%，分辨率为3.4 Å。并采用Rosetta软件中的FastRelax模块进行结构优化[6]，构象生成数为1000。并用线上模型质量评价服务器SAVES v5.0[7]对建模后的蛋白进行质量评价。

• 蛋白-蛋白对接

线上服务器Swarmdock Server[8]分别对接APOE2、APOE3、APOE4与Vamp2。对接类型为盲对接，产生模型数量设置为20。最后根据对接打分，选择最优复合物构象作为分子动力学模拟的初始结构。

• 分子动力学模拟

MD模拟采用Gromacs 2018.4程序[9]，在恒温恒压以及周期性边界条件下进行。应用Amber99sb-ildn全原子力场，TIP3P水模型[10]。在MD模拟过程中，所有涉及氢键采用LINCS算法[11]进行约束，积分步长为2 fs。静电相互作用采用（Particle-mesh Ewald）PME方法[12]计算，截断值设为1.2 nm。非键相互作用截断值设为10 Å，每10步更新一次。采用V-rescale [13]温度耦合方法控制模拟温度为300 K，采用Berendsen方法控制压力为1 bar。首先，采用最陡下降法对两个体系进行能量最小化，以消除原子间过近的接触；然后，在300 K进行100 ps的NVT平衡模拟；最后，对APOE3-VAMP2和APOE4-VAMP2两个复合物体系不同体系分别进行100 ns的MD模拟，每10 ps保存一次构象，模拟结果可视化采用Gromacs内嵌程序和VMD完成。并用g_mmpbsa程序分别计算APOE3和APOE4与Vamp2蛋白的结合自由能。

1. 结果与分析

• APOE4蛋白结构收敛性分析

均方根偏差(Root mean square deviation, RMSD)表示某一时刻的构象与目标构象所有原子偏差的加和，是衡量体系是否稳定的重要依据。从图1可以看出，APOE4经过20 ns的动力学模拟，从10 ns开始APOE4蛋白的RMSD值趋于稳定，平均值为0.5238±0.0212 nm。APOE2经过20 ns的动力学模拟，从5 ns开始APOE2蛋白的RMSD值趋于稳定，平均值为0.3919±0.0129 nm因此两者可以用于后续的蛋白相互作用研究。

图1 APOE4(a), APOE2(b)蛋白的RMSD值随MD模拟时间的变化

• 同源建模模型质量评价

本文分别用PROCHECK程序对优化后的蛋白模型进行了评价。Ramachandran plot用于阐述蛋白质或肽立体结构中肽键内α碳原子和羰基碳原子间的键的旋转度对α碳原子和氮原子间的键的旋转度，主要用来指明蛋白质或肽类中氨基酸残基的允许和不允许的构象。Ramachandran plot主要分为三个区域，允许区（红色区域），最大允许区（黄色区域），不允许区（空白区域）。图2a给出了APE蛋白模型的Ramachandran plot结果，位于允许区的氨基酸占比98%，没有氨基酸残基位于扭转角禁止区域，因此蛋白的所有氨基酸残基的二面角均在合理范围内，符合立体化学能量规则。

图2 VAMP2蛋白模型的Ramachandran plot

• MD模拟收敛参数分析

均方根偏差(Root mean square deviation, RMSD)表示某一时刻的构象与目标构象所有原子偏差的加和，是衡量体系是否稳定的重要依据，可用于量化蛋白质三级结构的稳定性[21]。图3a给出了APOE3-VAMP2，APOE4-VAMP2和APOE2-VAMP2复合物体系骨架原子的RMSD值随时间的变化情况。从图3a中可以看出两个体系在模拟初期(0-20 ns)的RMSD变化趋势及变化值相近，呈现逐步上升的趋势，这主要是由于模拟刚开始蛋白与周围水溶剂之间的相互作用造成蛋白结构的较大波动。整体来讲，三个体系的RMSD值从50 ns后趋于稳定，APOE4-Vamp2体系的RMSD值为2.0738±0.0241 nm，波动幅度为1.16%，APOE3-Vamp2体系的RMSD值为2.4434±0.0993 nm，波动幅度为4.06%。APOE2-Vamp2体系的RMSD值为2.0574±0.0453 nm，波动幅度为2.21%。因此，APOE2-Vamp2，POE4-VAMP2体系在模拟过程中的RMSD值相对更小，表明AOPE2,APOE4与VAMP2形成的复合物体系比更加稳定。

回转半径（Radius of Gyration, Rg）可以用来描述整体结构的变化情况，可用于表征蛋白质结构的紧密程度，Rg变化越大表明体系越膨胀。从图3b可以看出三个体系的Rg值在模拟的前20 ns呈逐渐降低的趋势，主要是由于VAMP2较长的α螺旋在模拟过程中逐渐靠近APOE蛋白，二者的结构变得更加紧密。从60 ns开始AOPE3-Vamp2，AOPE4-Vamp2复合体系的Rg趋于稳定，AOPE3-Vamp2，AOPE4-Vamp2的Rg值分别为2.500±0.013nm，2.382±0.030 nm， APOE4-VAMP2体系的Rg值相对更小，表明VAMP2与APOE4形成 额 复合物结构比APOE3更加紧密，因此推测二者的相互作用可能更强。AOPE2-Vamp2在60ns处有一个波动，从80 ns开始处于一个稳定的状态，Rg值为2.465±0.031 nm，表明VAMP2与APOE2形成复合物结构比APOE3更加紧密，但弱于APOE4与Vamp2的复合物。

图3 APOE3-VAMP2，APOE4-VAMP2和APOE2-VAMP2体系的RMSD（a）和Rg（b）随MD模拟时间的变化

• APOE蛋白柔性分析

均方根涨落（Root mean square fluctuation, RMSF）可以表示蛋白质中氨基酸残基的柔性大小。图4分别给出了APOE2，APOE3和APOE4蛋白的氨基酸残基RMSF分布。从图中可以看出，APOE3和APOE4二者差别较大的区域主要有Ala160~Val190和Arg191~Glu205。从图4b中可以看出，APOE4蛋白的Ala160~Val190区域更靠近VAMP2，受到VAMP2中氨基酸残基的相互作用影响，因此其柔性相比APOE3更大；而APOE3结构中的Arg191~Glu205区域与VAMP2结构相互作用比APOE4中对应区域更强，因此其柔性相对更大。整体来讲，APOE3和APOE4在模拟过程中蛋白的整体柔性分布基本一致，变化较大的区域主要原因是受到VAMP2结合的影响。APOE2与APOE4二者差别较大的区域主要有Ala237~Lys261和Arg191~Gly211。蛋白之间的相互作用可能导致这部分氨基酸残基构象发生变化较其他区域剧烈，因此柔性更大。

图4 APOE3和APOE4蛋白氨基酸残基的柔性分布（a）以及各个区域在复合物结构中的位置（b）

• VAMP2与APOE蛋白相互作用分析

氢键是生物大分子维持结构稳定性的重要作用力之一，也是蛋白-蛋白相互作用的重要形式。因此，为了进一步分析VAMP2与APOE2，APOE3和APOE4蛋白的相互作用，本文对MD模拟过程中VAMP2与APOE蛋白之间的氢键数量进行了统计分析，如图5a所示。从图中可以看出，模拟过程中APOE4与VAMP2之间的氢键数量一直比APOE3，APOE2更多，平均值分别为16.74和12.20，15.42.

氢键作用是可以增强VAMP2与APOE蛋白之间的亲和力，为了进一步评价蛋白之间亲和力的大小差异，图5b给出了APOE2，APOE3和APOE4蛋白与VAMP2在模拟过程中的结合自由能随模拟时间的变化。从图中可以看出，APOE蛋白与VAMP2的结合能在60 ns之后基本稳定，APOE2，APOE3和APOE4与VAMP2的结合能平均值分别为-548.851±116.578 kJ/mol，481.21±110.52 kJ/mol和-568.57±140.47 kJ/mol。

通过分析APOE蛋白与VAMP2之间的氢键作用和结合能，可以发现VAMP2与APOE4之间的氢键作用更强，进而使得二者的亲和力也更强。

图5 APOE蛋白与VAMP2之间氢键数量（a）和结合能（b）随MD模拟时间的变化

• APOE3和APOE4与VAMP2结合的差异

为了进一步研究VAMP2与APOE蛋白结合的差异，图6给出了VAMP2蛋白在APOE表面的结合模式。从图中可以看出，VAMP2蛋白主要依靠结构中间的α螺旋与APOE蛋白进行结合，主要结合在APOE蛋白的亲水性凹槽中，因此可以推测亲水作用是二者的分子识别过程中的关键作用力之一。另外VAMP2蛋白与APOE蛋白结合表面附近也存在一些疏水区域，疏水作用也可以进一步增强蛋白的结合。图6中还给出了VAMP2结构中可以与Snap25、Syntaxin蛋白形成SNARE复合物的区域（31-91），从图中可以看出，SNARE区域可以与APOE蛋白完全结合，进而影响VAMP2与其他蛋白质的作用。

图6 VAMP2在APOE3和APOE4蛋白表面的结合模式（cartoon模型表示VAMP2，红色区域表示与Snap25、Syntaxin蛋白形成SNARE复合物的区域，APOE3和APOE4用surface模型表示，蛋白表面蓝色和橙色区域分别表示亲水和疏水区域）

前面通过分析模拟过程中APOE蛋白与VAMP2分子识别各种参数的差异，发现APOE4与VAMP2形成的复合物结构更加稳定，亲和力更强。APOE3和APOE4的主要结构差异在112位氨基酸的变化，因此为了深入研究112位氨基酸突变所造成的与VAMP2结合的差异，本文对MD模拟之后APOE蛋白112位附近的氨基酸构象进行了分析，如图7所示。从图7a中可以看出，Cys112侧链在APOE3和VAMP2的结合中没有起到明显的作用，APOE3和VAMP2的分子识别主要依靠Asp110和Lys91，Val185和Asn92之间的氢键作用，以及Val116、Val185（APOE3）与Leu99、Ile98、Met96、Met95、Leu93（VAMP2）之间的疏水作用。而在APOE4与VAMP2的结合模式中（图7b），Arg112的侧链翻转朝向VAMP2蛋白，与Glu109形成了氢键作用，并且与VAMP2结构中的芳香性氨基酸Tyr88形成了Cation-π作用，Arg112与芳香性苯环质心距离为4.22 Å。另外APOE4中的Arg189也与VAMP2的Trp89形成了Cation-π作用，Arg189的主链和Asn92、Glu238和Lys91之间也存在氢键作用，进一步增强了APOE4与VAMP2的亲和力。

图7 APOE3（a）和APOE4（b）APOE3（c与VAMP2的结合模式（绿色和红色虚线分别表示氢键作用和cation-π作用）

1. 结论

本文通过分子动力学模拟方法分别研究了Vamp2与APOE2，APOE3和APOE4的分子识别过程。通过分析模拟过程中的收敛参数变化、氢键以及结合能，发现APOE4-Vamp2复合物体系在模拟过程中比APOE3-Vamp2体系更加稳定，Vamp2与APOE4之间的亲和力相对更强。最后通过分析112位突变点附近的氨基酸构象差异，发现Cys112突变为Arg112之后，不仅可以影响周围氨基酸残基的构象发生变化，形成更强的氢键作用，并且Arg112侧链还可以与Vamp2上的Tyr88形成了Cation-π作用，Tyr88位于Vamp2的SNARE区域上，在进一步增强APOE4与VAMP2的相互作用的同时，还可以影响Vamp2与Snap25、Syntaxin蛋白的结合。Arg158突变为Cys后，APOE2蛋白与Vamp2在Cys附近的结合也会有所差别,APOE2蛋白的Glu109氨基酸残基与Vamp2的Leu93氨基酸残基之间形成氢键，距离为3.3A。另外APOE2蛋白的Gly120氨基酸残基与Vamp2的Gly100氨基酸残基之间形成氢键, 距离为4.2A。

最后，有相关需求欢迎通过微信公众号联系我们。

微信公众号：320科技工作室。

参考文献

1. Chen J, Li Q, Wang J. Topology of human apolipoprotein E3 uniquely regulates its diverse biological functions. Proc Natl Acad Sci U S A. 2011;108(36):14813-14818.

2. Delano W.L. The PyMOL Molecular Graphics System[J]. My Publications, 2014, 30:442-454.

3. Benjamin, Webb, Andrej, et al. Comparative protein structure modeling using modeller.[J]. Current Protocols in Protein Science, 2016.

4. Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., Heer, F.T., de Beer, T.A.P., Rempfer, C., Bordoli, L., Lepore, R., Schwede, T. SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 46(W1), W296-W303 (2018).

5. Stein A, Weber G, Wahl MC, Jahn R. Helical extension of the neuronal SNARE complex into the membrane. Nature. 2009;460(7254):525-528.

6. Conway P , Tyka M D , Dimaio F , et al. Relaxation of backbone bond geometry improves protein energy landscape modeling.[J]. Protn ence, 2013, 23(1):47-55.

7. Pontius J, Richelle J, Wodak SJ. Deviations from standard atomic volumes as a quality measure for protein crystal structures. J Mol Biol. 1996 Nov 22;264(1):121-36.

8. M. Torchala, I. H. Moal, R. A. G. Chaleil, J. Fernandez-Recio, P. A. Bates, 'SwarmDock: a server for flexible protein-protein docking', Bioinformatics 29, 807-809 (2013).

9. Van d S D , Lindahl E , Hess B , et al. GROMACS: Fast, flexible, and free. Journal of Computational Chemistry, 2005, 26(16):1701-1718.

10. Jorgensen W L , Chandrasekhar J , Madura J D , et al. Comparison of simple potential functions for simulating liquid water. Journal of Chemical Physics, 1983, 79(2):926-935.

11. Hess B , Bekker H , Berendsen H J C , et al. LINCS: A linear constraint solver for molecular simulations. Journal of Chemical Theory & Computation, 1997, 4(1):1463–1472.

12. Darden T A , York D M , Pedersen L G . Particle mesh Ewald-an N.log(N) method for Ewald sums in large systems. The Journal of Chemical Physics, 1992, 98:10089-10092.

13. Berendsen H J C , Postma J P M , Van Gunsteren W F , et al. Molecular dynamics with coupling to an external bath. Journal of Chemical Physics, 1984, 81(8):3684-3690.