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

1. 前言

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

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

敲黑板处，请注意！！！在模拟开始前，须事先准备好两者的复合物的初始结构，复合物结构大多来自分子对接的结果。分子对接流程请查看相关教程。 （含共晶配体蛋白与小分子的分子对接（DOCK 6.9） | 殷赋的博客)

2. 流程图

3. 结构信息

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

配体分子：共晶配体6DC

4. 平台操作步骤

4.1 准备受体三维结构

1. 打开殷赋云平台（殷赋云 - 科研计算快易准)【处理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号残基的位置约束，生产阶段放开约束。