蛋白与小分子的分子动力学模拟（隐式溶剂）

分子动力学
分子动力学方案
发布日期:   2023-06-25
更新日期:   2024-03-07
文章字数:   3.2k
阅读时长:   11 分
阅读次数:  

1. 前言

分子动力学（Molecular Dynamics, MD）模拟技术的一个广泛应用是研究蛋白与小分子之间的相互作用。本教程讲述如何使用殷赋云平台的大方案和小工具，实现蛋白与小分子复合物的隐式溶剂分子动力学模拟。

**在模拟开始前，须事先准备好两者的复合物的初始结构。在实际研究中，复合物结构大多来自分子对接的结果。**为方便读者学习，本文采用晶体结构作为示例，分子对接流程请查看相关教程。同时，由于MD分析内容繁多驳杂，因项目而异，本文只介绍简单常规的分析，更多分析内容详见相关教程。

敲黑板处，请注意！！！在模拟开始前，须事先准备好两者的复合物的初始结构，复合物结构大多来自分子对接的结果。分子对接流程请查看相关教程。（http://blog.yinfotek.com/posts/10101/)

2. 流程图

graph LR
A(准备受体结构) --> B(准备Amber文件)
B(准备Amber文件) --> C(开展动力学模拟)
C(开展动力学模拟) --> D(分析结果)
E(准备配体结构) --> F(添加原子电荷)
F(添加原子电荷) --> B

3. 结构信息

受体蛋白：人Serine/threonine-protein kinase B-raf蛋白，uniprot AC号：P15056，PDB编号：5ITA

配体分子：共晶配体6DC

4. 平台操作步骤

4.1 准备受体三维结构

  1. 打开殷赋云平台（https://cloud.yinfotek.com/)【处理PDB结构（进阶版）】小工具；

  2. 输入PDB编号5ITA，点击【确定】；

  3. 删除多余结构（水分子），保留蛋白和共晶配体，点击【下一步】；

    该晶体结构为同源二聚体，含有A、B两条链，各包含一个共晶配体。我们保留A链蛋白和配体，删除B链结构和所有水分子，如下图所示。

  1. 进入修复PDB页面，查看折叠卡片中给出的蛋白问题；

    • 序列缺口

    当序列出现缺口时，通常使用NME和ACE封端。但当缺口两端残基编号仅差1-2号，要特别留意“非标准残基”卡片中的残基了。

    例如，缺口两段残基为Arg10和Ser12，而在非标准残基中发现磷酸化的Phe11，那么，这里其实并无缺口，只是因为非标准残基而不被识别罢了。因此，不应封端，以免NME和ACE与Phe11重叠。

    • 非标准残基

      除标准氨基酸、标准核苷酸、水分子、金属离子以外，其他残基分子都被认为是非标准残基，通常是修饰后的残基或有机小分子。若残基修饰对计算模拟不重要，则应替换为标准残基。

  2. 点击【生成文件】，下载文件。

    动力学模拟会给残基编号重新编号，为了知道新旧编号的对应关系，可勾选对残基重新编号，从renumber.csv文件查看。

4.2 准备配体结构

  1. 打开殷赋云平台【准备化合物结构】小工具；

  2. 上传分子文件，不要勾选【在MMFF94力场中能量最小化】；

    配体分子的构象是适应蛋白口袋的，若进行能量最小化，将改变构象，会与蛋白原子过近甚至重叠，这是应当避免的。

  3. 勾选【质子化/去质子化】，点击【准备】，下载文件。

    一般来说，蛋白-小分子所处环境是水溶液。如果结构中带有氨基、羧基等可被质子化或去质子化的基团，小分子会呈离子状态，因此勾选质子化/去质子化。若明确知道其分子状态，应按实际情况处理。

    建议：下载compounds.mol2文件，改为意义更明确的名称，本例改为ligand.mol2。

4.3 添加原子电荷

经过上述步骤，分子结构已带有MMFF94电荷，但在MD模拟中，推荐使用精度更高的电荷：AM1-BCC或RESP电荷。前者操作简单，精度接近后者，但不适合于较大的分子（分子量 > 1000）；后者需要先进行量化计算，相对繁琐，但是公认最好的电荷之一，强烈推荐使用。

下面分别介绍两种电荷的添加方法：

4.3.1 AM1-BCC电荷

  1. 打开殷赋云平台【添加原子电荷】小工具；

  2. 上传分子文件，填写净电荷，选择AM1-BCC电荷方法，点击【生成】；

    净电荷是指分子结构中不能抵消的正电荷或负电荷（整数）。中性分子净电荷为0，带电分子的净电荷即是其所带电荷的代数和，例如：乙酸净电荷为0，去质子化（脱去羧基氢）后为-1，十八烷基胺的净电荷为0，质子化后为1。

    不填写时，平台会自动计算。但对于特殊结构，可能会计算出错。必要时，请用户根据自己的化学知识判断填写。

  3. 检查电荷分布，下载文件。

    注意：因显示问题，下图的苯环只显示单键，没有双键，实为共轭大π键。

4.3.2 RESP电荷

  1. 打开殷赋云平台【量化计算（Gaussian 09）】大方案，创建任务，进入提交任务页面；

  2. 上传分子结构，设置净电荷，其他保持默认；

  3. 按下图设置计算步骤，点击【提交】任务；

  4. 待任务完成，点击【查看】，进入分析页面；

  5. 选择log格式，点击selected下载文件或从文件列表中找到*.log.zip文件。

  6. 打开殷赋云平台【添加原子电荷】小工具；

  7. 上传分子文件，填写净电荷，选择RESP电荷方法，上传量化计算文件，选择程序，点击【生成】，下载ligand.mol2文件；

  8. 检查电荷分布，下载文件。

    可见，RESP电荷与AM1-BCC电荷是有区别的。

4.4 准备Amber文件

  1. 打开殷赋云平台【准备Amber文件】小工具；

  2. 上传蛋白和小分子的结构文件；

    • 蛋白、核酸、多肽、水、金属离子等结构采用pdb格式文件，有机小分子采用mol2格式文件；
    • 建议修改有机小分子的残基名为有意义的名称（3个字符，仅英文大写字母或阿拉伯数字），不同分子应具有唯一名称，且避免采用标准氨基酸、核苷酸的缩写名称。

  3. 设置【环境】，不勾选【添加溶剂】，则不添加水分子，点击【生成】；

    尽管不添加溶剂，但若体系净电荷不为0，本工具仍会添加抗衡离子以中和体系。

  4. 查看结构，下载amber_files.zip文件；

  5. 查看residue_numbers.csv文件，了解残基编号（后面会用到）。

    该体系由蛋白、小分子和氯离子构成，残基编号分别为：1-250、251和252-257。

4.5 开展动力学模拟

  1. 打开殷赋云平台《分子动力学（Amber20）》大方案，创建任务，进入提交任务页面；

  2. 上传拓扑文件和坐标文件；

  3. 设置模拟环境；

    非键截断值（cutoff）不能大于生成Amber文件设置的水盒子边距。

  4. 设置计算步骤；

    简单解释上图各计算步骤的目的：

    1. Minimization：优化整个体系，释放（晶体结构中的）张力或过近的原子接触（close contacts）或碰撞（clashes）；

    2. Heat：将体系（从0）升温至300 K；b

    3. Equilibration：让体系达到平衡；

      • 需要预先估计该体系总费用的用户，请参考以下指导完成分步操作。

        • 完成第1到第3步操作后，先不设置第4步，点击【提交】，等计算完毕，到 用户中心-费用账单-账单明细查看费用，估算出跑1ns的费用。

        • 如需完成后续模拟，可以进入我的项目，点击查看，从文件列表中下载prmtop（拓扑）和rst（坐标）文件（下拉文件列表至末尾，选取计算步骤3生成的rst文件）作为后续模拟的起点。

        • 新建分子动力学任务，上传 prmtop 和 rst 这两个文件，参考步骤4Production的时长分段建议设置production时长。（步骤1-3的信息已在prmtop和rst文件中，无需设置重跑）

      • 不需要估计该体系总费用的用户，可直接进行第4步骤。

    4. Production：开始长时间采样。

      • 通常所说跑多少ns的动力学，是指生产（Production）阶段的总时长；这里设置一个50 ns，当然可以继续增加Production阶段以便达到预期时长；

      • 建议将长时间模拟分成若干小段，例如，100 ns的动力学可设置为5个连续的20 ns Production，200 ns的动力学可以设置成4个50 ns。这样好处很多：可避免意外事故发生时毁掉整段轨迹（只要从最后一段完整的轨迹接着算），可方便后期分析时选择理想时段而不用进行轨迹分割处理，可了解进度并在合适时候提前终止任务……

  5. 点击【提交】。

4.6 分析结果

  1. 待任务完成，点击【查看】，进入分析页面；

  2. 查看热力学性质；

    在平衡（Equilibration）或生产（Production）阶段，能量、温度、压强、体积等曲线应当平稳。如果发现异常，应检查前后阶段的结构（见下面第5点）。

    详情见《分子动力学模拟（Amber 20）》文章。

  3. 查看RMSD；

    通过RMSD分析，了解体系的构象变化，体系是否达到平衡；详见《分析Amber轨迹（一）》。

    下图所示，在50 ns的动力学模拟中，蛋白骨架、溶质重原子和小分子重原子的变化比较趋势一致。在经历短暂的跃升后，体系迅速达到平稳，并保持至采样最后。尽管某些时刻存在一些相对剧烈的波动，但整体振幅保持在2 Å以下。这表明，该体系已经达到平稳，可做进一步分析。

    如果出现明显的上升趋势，或者波动十分剧烈，应当查看该阶段的始终状态的结构（见下面第5点），分析原因，必要时延迟模拟时间，或者调整条件重新模拟。

  4. 分析氢键（可选）；

    下载氢键数据文件hbond.xlsx，查看统计数据，确认关键的氢键；详见《分子动力学模拟（Amber 20）》。

    下图整理了配体分子（6DC）相关的氢键作用（删除Fraction < 1%的氢键），可见所有氢键的比例均不到50%，表明这些氢键并不稳定（起码达到70%才算稳定）。

    • 注意：由于使用隐式溶剂模型，没有加入水分子，因此不存在水桥作用的统计数据。

    • 影响Fraction的因素有很多，可能是采样时间太长，把前期波动较厉害的阶段纳入统计范围，自然拉低了平均值，也可能是小分子的结合模式有问题，应查看结构分析原因（详见第5点分析）。

  5. 下载各阶段的结构文件，使用熟悉的软件（比如 PyMOL或平台【处理PDB结构】小工具）查看结构。

    观察发现，上述分析揭示的氢键与晶体结构观察到的存在一定差异：

    • 在晶体结构中，配体磺胺N原子（NAS）与ASP147（原始编号549）骨架酰胺N原子形成氢键作用，但在模拟中却是它与ASP147羧基O原子（OD1）形成氢键。经过分析，我们认为晶体结构有误。原因是，配体N原子上有H原子，是氢键供体，ASP147酰胺N原子也是，两者是相互排斥的，不应形成氢键，MD模拟充分显示了这一点。
    • 在晶体结构中，磺胺O原子（OAF）与Phe148（原始编号550）形成氢键，但经过动力学模拟，它却与Asp147的骨架N原子形成氢键，同时，磺胺另一O原子（OAE）也与之形成氢键。这样导致磺胺围绕N-S键发生了旋转，末端丙基的结合位置也有了很大变化，因而使得Gly149（原始编号551）失去了与磺胺的氢键作用，并远离配体。
    • 在晶体结构中，ASN134（原始编号581）的侧链朝向配体，并与之形成水桥，但在模拟中，由于使用了隐式溶剂模型，没有水分子，观察不到水桥，Asn134侧链的构象也因此发生了很大改变。
    • 文件名带有-dry的是删除离子后的结构；
    • 将来会上线专门的结构可视化功能。

文章作者: 殷赋量子氢
文章链接: http://blog.yinfotek.com/posts/10406/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 殷赋量子氢 !
分子动力学
  目录