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

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

1. 前言

分子动力学（Molecular Dynamics, MD）模拟技术的一个广泛应用是研究蛋白与小分子之间的相互作用。本教程讲述如何使用殷赋云平台的大方案和小工具，实现蛋白与小分子复合物的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链结构和所有水分子，如下图所示。

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

    • 序列缺口

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

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

    • 非标准残基

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

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

    动力学模拟会给残基编号重新编号，为了知道新旧编号的对应关系，可勾选对残基重新编号，从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. 设置盒子形状和盒子边距，点击【生成】；

    由于该体系整体是接近球形的，因此使用Truncated Octahedron盒子比较合适；若呈长条形，则应选择Cuboid盒子。这样能最大程度匹配体系形状，减少水分子个数，从而提高MD运算速度，节省费用。

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

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

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

4.5 开展动力学模拟

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

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

  3. 设置模拟环境；

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

  4. 设置计算步骤，添加约束条件；（继续下拉可以看见计算步骤的分步解释）

    ​

    ​

    约束条件中的残基编号可从residue_numbers.csv文件中查看。如下图，残基编号总是从1开始，第1步约束蛋白和小分子需要选取1-251号进行约束，第2步约束蛋白骨架需要选取1-250号进行约束。

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

    1. Minimization：约束蛋白和小分子的原子，优化水分子、盐离子；

    2. Minimization：约束蛋白骨架原子（也可加上小分子的重原子），优化其他原子；

    3. Minimization：不加约束，优化整个体系；

    4. Heat：将体系（从0）升温至300 K，为避免剧烈运动破坏体系，与第1步Minimization一样设置约束；

    5. Equilibration：与第2步Minimization一样设置约束，使用NPT系综来控制压强，让体系调整自身密度；

    6. Equilibration：体系密度已达平衡，不再有大变动，可使用NVT系综控制体积，并放开所有约束，让体系继续平衡；

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

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

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

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

      • 不需要估计该体系总费用的用户，可继续第7步骤。

    7. Production：与上一步Equilibration的条件相同，从此开始长时间采样。

      • 通常所说跑多少ns的动力学，是指生产（Production）阶段的总时长；这里设置一个20 ns；

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

  5. 点击【提交】；

4.6 分析结果

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

  2. 查看热力学性质；

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

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

  3. 查看RMSD；

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

    下图所示，在20 ns的动力学模拟中，蛋白和全部溶质重原子的变化趋势一致，逐步上升并在18 ns附近达到顶峰，然后有回落的趋势；小分子重原子则除了2-7 ns间出现了两个短暂跃升的平台外，其余时间均保持平稳。但如果条件允许，建议延长模拟时间。

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

  4. 分析氢键；

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

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

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

  5. 下载各阶段的结构文件，使用自己熟悉的软件（比如 PyMOL）查看结构。

    观察发现，上述分析揭示的氢键与晶体结构观察到的存在差异：在晶体结构中，配体磺胺N原子（NAS）与ASP147（原始编号549）骨架酰胺N原子形成氢键作用，但在模拟中却是它与ASP147羧基O原子（OD1）形成氢键。经过分析，我们认为晶体结构有误。原因是，配体N原子上有H原子，是氢键供体，ASP147酰胺N原子也是，两者是相互排斥的，不应形成氢键，MD模拟充分显示了这一点。

    在晶体结构中，ASN134（原始编号581）的侧链朝向配体，并与之形成水桥，但在MD模拟中该水桥消失。通过观察各个阶段的*-dry.pdb文件，发现是从加热后的第一个平衡阶段开始，ASN134的侧链逐渐远离配体（想想为什么加热阶段没有出现这种现象？请把答案写在评论区）。我们认为该水桥是可能存在的，为确认这一点，考虑重新模拟：在平衡阶段将相关残基进行位置约束，生产阶段再放开约束，来观察水桥能否保持。

    • 文件名带有-dry的是删除水和离子后的结构；
    • 将来会上线专门的结构可视化功能。

  6. 下载rst文件作为后续模拟的起点。

    根据上述第5点分析，考虑从加热后开始重新模拟。从文件列表里下载complex.prmtop（拓扑）和4.heat.rst（坐标）文件，新建分子动力学任务，上传这两个文件，设置计算步骤；其中，两段平衡阶段均增加对134号残基的位置约束，生产阶段放开约束。

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