文章/答案/技术大牛

发布

社区首页 >专栏 >分子动力学--蛋白配体模拟

分子动力学--蛋白配体模拟

原创

追风少年i

发布于 2025-09-18 12:49:43

9900

代码可运行

运行总次数：0

代码可运行

作者，Evil Genius

大家的下班空余时间，或者科研的空闲时间多多扩展自己的知识内容，视野越广，看问题越全面。

很多人说自己失业了，很焦虑，没关系，多学点东西，我也失业过，失业的时间坚持把客户的课题都做完了，现在行情不好，失业是很正常的事情，正好失业有大把的时间，不如多学习学习。

今天我们继续分子动力学，其实在单细胞空间基因组培训上说过突变对蛋白的影响，当时想达到的目的是通过计算的方法计算出某个点突变对蛋白结构、酶活的影响，当时由于知识面比较少，以为没有，后来有个做分子动力学模拟的同学告诉我是可以的，所以我们要扩展这部分内容。

大家有空也学一学，万一我学成了准备开个培训班，大家那就不可避免要付费了啊。

今天我们来分享蛋白配体结合的分子动力学模拟（教程）。

当然了，我们先学会用，深层次的内容需要一步一步扩展了。

准备的PDB文件：T4溶菌酶L99A/M102Q

下载结构后，可通过VMD、Chimera、PyMOL等可视化程序进行观察，我一般通过PyMOL查看。

去除晶体水、PO4和BME分子。注意这种操作并非普遍适用（例如结合活性位点水分子的情况）。

分两步进行

1、使用pdb2gmx制备蛋白质拓扑

2、使用外部工具制备配体拓扑

鉴于需要分别制备这两种拓扑，必须将蛋白质和JZ4配体分别保存至独立的坐标文件中。

首先保存保存JZ4坐标

grep JZ4 3HTB_clean.pdb > jz4.pdb

从3HTB_clean.pdb文件中删除JZ4相关行。

教程将使用CHARMM36力场，使用pdb2gmx为T4溶菌酶生成拓扑文件

gmx pdb2gmx -f 3HTB_clean.pdb -o 3HTB_processed.gro -ter

三项选择

力场选择

水模型选择

末端类型选择

选择CHARMM36力场（选项1）

选择默认水模型（CHARMM改良版TIP3P），随后为末端基团选择"NH3+"和"COO-"。

为JZ4添加氢原子

使用jz4.mol2文件前需进行若干修正。用文本编辑器打开jz4.mol2文件

@<TRIPOS>MOLECULE
*****
 22 22 0 0 0
SMALL
GASTEIGER

@<TRIPOS>ATOM
      1 C4         24.2940  -24.1240   -0.0710 C.3   167  JZ4167     -0.0650
      2 C7         21.5530  -27.2140   -4.1120 C.ar  167  JZ4167     -0.0613
      3 C8         22.0680  -26.7470   -5.3310 C.ar  167  JZ4167     -0.0583
      4 C9         22.6710  -25.5120   -5.4480 C.ar  167  JZ4167     -0.0199
      5 C10        22.7690  -24.7300   -4.2950 C.ar  167  JZ4167      0.1200
      6 C11        21.6930  -26.4590   -2.9540 C.ar  167  JZ4167     -0.0551
      7 C12        22.2940  -25.1870   -3.0750 C.ar  167  JZ4167     -0.0060
      8 C13        22.4630  -24.4140   -1.8080 C.3   167  JZ4167     -0.0245
      9 C14        23.9250  -24.7040   -1.3940 C.3   167  JZ4167     -0.0518
     10 OAB        23.4120  -23.5360   -4.3420 O.3   167  JZ4167     -0.5065
     11 H          25.3133  -24.3619    0.1509 H       1  UNL1        0.0230
     12 H          23.6591  -24.5327    0.6872 H       1  UNL1        0.0230
     13 H          24.1744  -23.0611   -0.1016 H       1  UNL1        0.0230
     14 H          21.0673  -28.1238   -4.0754 H       1  UNL1        0.0618
     15 H          21.9931  -27.3472   -6.1672 H       1  UNL1        0.0619
     16 H          23.0361  -25.1783   -6.3537 H       1  UNL1        0.0654
     17 H          21.3701  -26.8143   -2.0405 H       1  UNL1        0.0621
     18 H          21.7794  -24.7551   -1.0588 H       1  UNL1        0.0314
     19 H          22.2659  -23.3694   -1.9301 H       1  UNL1        0.0314
     20 H          24.5755  -24.2929   -2.1375 H       1  UNL1        0.0266
     21 H          24.0241  -25.7662   -1.3110 H       1  UNL1        0.0266
     22 H          23.7394  -23.2120   -5.1580 H       1  UNL1        0.2921
@<TRIPOS>BOND
     1     4     3   ar
     2     4     5   ar
     3     3     2   ar
     4    10     5    1
     5     5     7   ar
     6     2     6   ar
     7     7     6   ar
     8     7     8    1
     9     8     9    1
    10     9     1    1
    11     1    11    1
    12     1    12    1
    13     1    13    1
    14     2    14    1
    15     3    15    1
    16     4    16    1
    17     6    17    1
    18     8    18    1
    19     8    19    1
    20     9    20    1
    21     9    21    1
    22    10    22    1

需要修改的第一个地方是MOLECULE标题行。将"*****"替换为"JZ4"

@<TRIPOS>MOLECULE
JZ4

接下来需要统一修正残基名称和编号。由于该.mol2文件仅包含一个分子，因此应只指定一个残基名称和编号。编辑完成后，jz4.mol2文件中修正后的ATOM部分应显示为：

@<TRIPOS>ATOM
      1 C4         24.2940  -24.1240   -0.0710 C.3     1  JZ4        -0.0650
      2 C7         21.5530  -27.2140   -4.1120 C.ar    1  JZ4        -0.0613
      3 C8         22.0680  -26.7470   -5.3310 C.ar    1  JZ4        -0.0583
      4 C9         22.6710  -25.5120   -5.4480 C.ar    1  JZ4        -0.0199
      5 C10        22.7690  -24.7300   -4.2950 C.ar    1  JZ4         0.1200
      6 C11        21.6930  -26.4590   -2.9540 C.ar    1  JZ4        -0.0551
      7 C12        22.2940  -25.1870   -3.0750 C.ar    1  JZ4        -0.0060
      8 C13        22.4630  -24.4140   -1.8080 C.3     1  JZ4        -0.0245
      9 C14        23.9250  -24.7040   -1.3940 C.3     1  JZ4        -0.0518
     10 OAB        23.4120  -23.5360   -4.3420 O.3     1  JZ4        -0.5065
     11 H          25.3133  -24.3619    0.1509 H       1  JZ4         0.0230
     12 H          23.6591  -24.5327    0.6872 H       1  JZ4         0.0230
     13 H          24.1744  -23.0611   -0.1016 H       1  JZ4         0.0230
     14 H          21.0673  -28.1238   -4.0754 H       1  JZ4         0.0618
     15 H          21.9931  -27.3472   -6.1672 H       1  JZ4         0.0619
     16 H          23.0361  -25.1783   -6.3537 H       1  JZ4         0.0654
     17 H          21.3701  -26.8143   -2.0405 H       1  JZ4         0.0621
     18 H          21.7794  -24.7551   -1.0588 H       1  JZ4         0.0314
     19 H          22.2659  -23.3694   -1.9301 H       1  JZ4         0.0314
     20 H          24.5755  -24.2929   -2.1375 H       1  JZ4         0.0266
     21 H          24.0241  -25.7662   -1.3110 H       1  JZ4         0.0266
     22 H          23.7394  -23.2120   -5.1580 H       1  JZ4         0.2921

最后请注意@<TRIPOS>BOND部分中的异常键级顺序。所有程序似乎都有自己生成该列表的方法，但效果参差不齐。如果化学键未按升序列出，在构建具有匹配坐标的正确拓扑时会出现问题。要解决此问题，请下载sort_mol2_bonds.pl脚本并执行：

perl sort_mol2_bonds.pl jz4.mol2 jz4_fix.mol2

使用CGenFF生成JZ4拓扑参数

现在jz4_fix.mol2文件已可用于生成拓扑。访问CGenFF服务器，登录账户，点击页面顶部的"Upload molecule"上传jz4_fix.mol2文件。CGenFF服务器将快速返回CHARMM"流文件"格式（扩展名为.str）的拓扑参数。将网页内容保存为"jz4.str"文件

CHARMM流文件包含所有拓扑信息——原子类型、电荷和键连接关系。其中还包含通过类比生成的新键合参数章节，用于补充力场未覆盖的内相互作用参数。CGenFF同时提供每个参数的惩罚值评分，用于评估所分配参数的可靠性。低于10分的参数可直接使用；10-50分之间的参数需要一定验证；超过50分的参数通常需要手动重新参数化。这种惩罚值评分是CGenFF服务器最重要的特性之一。其他多数服务器生成拓扑时如同"黑盒"，用户只能盲目信任。将整个研究项目寄托于未经验证的自动化程序非常危险——劣质配体拓扑会导致大量时间浪费和不可靠结果。请务必验证新参数化分子的拓扑！至少需要对照力场中现有基团检查电荷量和分配的原子类型。

虽然CHARMM格式的JZ4拓扑很好，但要在GROMACS中运行模拟还需转换格式。从GitHub网站下载相应版本的cgenff_charmm2gmx.py脚本。执行转换命令：

python cgenff_charmm2gmx.py JZ4 jz4_fix.mol2 jz4.str charmm36-jul2022.ff

成功转换后请注意屏幕输出

============ DONE ============
Conversion complete.
The molecule topology has been written to jz4.itp
Additional parameters needed by the molecule are written to jz4.prm, which needs to be included in the system .top
============ DONE ============

构建复合物

通过pdb2gmx我们获得了包含处理后蛋白质结构的"3HTB_processed.gro"文件。同时从cgenff_charmm2gmx.py获得的"jz4_ini.pdb"包含所有必需氢原子且与JZ4拓扑中的原子名称匹配。使用editconf将此.pdb文件转换为.gro格式：

gmx editconf -f jz4_ini.pdb -o jz4.gro

将3HTB_processed.gro复制为新文件（如命名为"complex.gro"），通过添加JZ4坐标构建蛋白质-配体复合物。只需将jz4.gro的坐标部分复制到complex.gro中蛋白质原子末行与盒子矢量之间：

  163ASN      C 1691   0.621  -0.740  -0.126
  163ASN     O1 1692   0.624  -0.616  -0.140
  163ASN     O2 1693   0.683  -0.703  -0.011

    1JZ4     C4    1   2.429  -2.412  -0.007
    1JZ4     C7    2   2.155  -2.721  -0.411
    1JZ4     C8    3   2.207  -2.675  -0.533
    1JZ4     C9    4   2.267  -2.551  -0.545
    1JZ4    C10    5   2.277  -2.473  -0.430
    1JZ4    C11    6   2.169  -2.646  -0.295
    1JZ4    C12    7   2.229  -2.519  -0.308
    1JZ4    C13    8   2.246  -2.441  -0.181
    1JZ4    C14    9   2.392  -2.470  -0.139
    1JZ4    OAB   10   2.341  -2.354  -0.434
    1JZ4     H1   11   2.531  -2.436   0.015
    1JZ4     H2   12   2.366  -2.453   0.069
    1JZ4     H3   13   2.417  -2.306  -0.010
    1JZ4     H4   14   2.107  -2.812  -0.407
    1JZ4     H5   15   2.199  -2.735  -0.617
    1JZ4     H6   16   2.304  -2.518  -0.635
    1JZ4     H7   17   2.137  -2.681  -0.204
    1JZ4     H8   18   2.178  -2.476  -0.106
    1JZ4     H9   19   2.227  -2.337  -0.193
    1JZ4    H10   20   2.458  -2.429  -0.214
    1JZ4    H11   21   2.402  -2.577  -0.131
    1JZ4    H12   22   2.374  -2.321  -0.516

   5.99500   5.19182   9.66100   0.00000   0.00000  -2.99750   0.00000   0.00000   0.00000

向.gro文件添加了22个原子，需要相应增加complex.gro第二行的原子总数。当前坐标文件应包含2636个原子。

构建拓扑体系

在系统拓扑中包含JZ4配体参数非常简单：在topol.top文件包含位置限制文件后插入#include "jz4.itp"语句。位置限制文件的包含标志着"Protein"分子类型部门的结束。

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "./charmm36-jul2022.ff/tip3p.itp"

变成

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Include water topology
#include "./charmm36-jul2022.ff/tip3p.itp"

配体引入的新二面角参数已被cgenff_charmm2gmx.py脚本写入"jz4.prm"。在topol.top顶部添加包含语句：

; Include forcefield parameters
#include "./charmm36-jul2022.ff/forcefield.itp"

; Include ligand parameters
#include "jz4.prm"

[ moleculetype ]
; Name            nrexcl
Protein_chain_A     3

此#include语句的位置至关重要——必须出现在所有[ moleculetype ]条目之前，因为所有参数都必须在构建分子前定义；同时必须出现在父力场#include语句之后，因为所有键合参数都需要在原子类型已知的前提下引入。

最后需要在[ molecules ]指令中添加新分子：

[ molecules ]
; Compound        #mols
Protein_chain_A     1
JZ4                 1

至此，拓扑文件与坐标文件在系统组成方面已达到一致。

第三步：定义晶胞与添加溶剂

define the unit cell and fill it with water.

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

GROMACS程序始终采用数值效率最高的坐标表示方式，即将所有分子重包装到三斜晶胞中。mdrun执行的物理计算在不同坐标包装方式下具有等效性，因此采用最高效的方式更为可取。在生成.tpr文件后，可以随时恢复所需的晶胞形状。

第四步：添加离子

已经获得了包含带电蛋白质的溶剂化体系。pdb2gmx的输出信息显示该蛋白质净电荷为+6e（根据其氨基酸组成计算）。由于生物体系不能以净电荷形式存在，我们必须向系统中添加离子。

使用grompp程序组装.tpr文件时，可采用任意.mdp文件。建议使用能量最小化的.mdp文件，因为这类文件所需参数最少，最容易维护。

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

We now pass our .tpr file to genion:

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

[ molecules ]指令现在应显示为

[ molecules ]
; Compound        #mols
Protein_chain_A     1
JZ4                 1
SOL             10228
CL                  6

第五步：能量最小化

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

invoke mdrun to carry out the EM

gmx mdrun -v -deffnm em

Steepest Descents converged to Fmax < 1000 in 143 steps
Potential Energy  = -4.9014547e+05
Maximum force     =  8.7411469e+02 on atom 27
Norm of force     =  5.6676244e+01

Now that our system is at an energy minimum, we can begin real dynamics.

第六步：平衡

对配体施加约束

温度耦合组的处理

配体约束

为了约束配体，需要为其生成位置约束拓扑。首先为JZ4创建一个仅包含非氢原子的索引组：

gmx make_ndx -f jz4.gro -o index_jz4.ndx
...
 > 0 & ! a H*
 > q

然后，执行genrestr模块并选择这个新创建的索引组（该组在index_jz4.ndx文件中将被标记为组3）。

gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000

现在需要将这些信息添加到拓扑结构中。根据要使用的条件，可以通过以下几种方式实现：如果希望在蛋白质受约束时配体也同步受约束，请在拓扑文件的指定位置添加以下内容：

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Ligand position restraints
#ifdef POSRES
#include "posre_jz4.itp"
#endif

; Include water topology
#include "./charmm36-jul2022.ff/tip3p.itp"

位置至关重要！必须按照指示将posre_jz4.itp的调用语句放置在拓扑文件中。jz4.itp文件中的参数定义了配体的[ moleculetype ]指令，该分子类型的定义会持续生效直至包含水拓扑（tip3p.itp）之前。若将posre_jz4.itp的调用语句置于其他位置，会导致位置约束参数错误地应用于其他分子类型。

若需要在平衡化过程中进行更精细的控制（例如对蛋白质和配体分别施加约束），可以通过其他#ifdef条件块来控制配体位置约束文件的调用，具体方式如下：

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_jz4.itp"
#endif

; Include water topology
#include "./charmm36-jul2022.ff/tip3p.itp"

在后一种情况下，若需要同时约束蛋白质和配体，必须在.mdp文件中明确定义参数：define = -DPOSRES -DPOSRES_LIG。具体采用何种处理方式取决自身的需求。这些示例仅用于演示GROMACS所提供的灵活性。

温度调控

常规做法是设置 tc-grps = Protein Non-Protein。但需要注意的是，"Non-Protein"组实际上包含了JZ4配体。由于JZ4与蛋白质在物理结构上紧密连接，最合理的处理方式是将它们视为整体——即在温度耦合时把JZ4归入蛋白质组。同理，我们添加的少量氯离子（Cl-）也应归属于溶剂组。为实现这种分组，需要创建特殊的索引组来合并蛋白质和JZ4。这可以通过make_ndx模块来实现，操作时需提供完整体系的坐标文件（此处使用能量最小化后的输出结构em.gro）：

gmx make_ndx -f em.gro -o index.ndx
...
  0 System              : 33506 atoms
  1 Protein             :  2614 atoms
  2 Protein-H           :  1301 atoms
  3 C-alpha             :   163 atoms
  4 Backbone            :   489 atoms
  5 MainChain           :   651 atoms
  6 MainChain+Cb        :   803 atoms
  7 MainChain+H         :   813 atoms
  8 SideChain           :  1801 atoms
  9 SideChain-H         :   650 atoms
 10 Prot-Masses         :  2614 atoms
 11 non-Protein         : 30892 atoms
 12 Other               :    22 atoms
 13 JZ4                 :    22 atoms
 14 CL                  :     6 atoms
 15 Water               : 30864 atoms
 16 SOL                 : 30864 atoms
 17 non-Water           :  2642 atoms
 18 Ion                 :     6 atoms
 19 JZ4                 :    22 atoms
 20 CL                  :     6 atoms
 21 Water_and_ions      : 30870 atoms

通过以下操作合并"Protein"和"JZ4"组（其中">"符号表示make_ndx模块的命令提示符）：

> 1 | 13
> q

可以设置 tc-grps = Protein_JZ4 Water_and_ions 来实现预期的"蛋白质-非蛋白质"分组效果。

使用此.mdp文件继续进行NVT平衡化操作：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt

NPT平衡

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

gmx mdrun -deffnm npt

完成两个平衡化阶段后，系统现已在所需的温度和压力下达到充分平衡。接下来可以移除位置约束并开始正式的数据采集分子动力学模拟。该流程与之前操作类似——将利用包含压力耦合信息的检查点文件作为grompp的输入，并运行10纳秒的MD模拟（相关脚本可在此处获取）。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr

gmx mdrun -deffnm md_0_10

第七步：数据分析

坐标重居中与周期边界重包裹

在使用周期性边界条件的模拟过程中，分子可能出现"断裂"或穿越晶胞边界"跳跃"的现象。为使蛋白质重居中并将分子重新包裹在晶胞内（恢复理想的菱形十二面体结构），请调用trjconv命令：

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact

选择"Protein"作为居中基准组，选择"System"作为输出组。请注意：对于涉及多次周期边界跳跃的长时模拟，复合物（蛋白质-配体、蛋白质-蛋白质）的居中操作可能较为困难。此类情况（特别是蛋白质-蛋白质复合物）下，可能需要创建自定义索引组作为居中基准，例如选择某个蛋白质的活性位点或复合物中单体的界面残基。

要提取轨迹的初始帧（t = 0 ns），对重居中后的轨迹使用trjconv -dump命令：

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0

为实现更流畅的可视化效果，建议执行旋转和平移拟合。按以下方式运行trjconv命令：

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans

选择"Backbone"基于蛋白质主链进行最小二乘拟合，并选择"System"作为输出组。请注意：坐标的周期边界重包裹与拟合操作在数学上不兼容。如需执行拟合（该操作有助于可视化呈现，但大多数分析流程并不必需），请按此处所示分步进行坐标处理。

蛋白质-配体相互作用与配体动力学分析

计算轨迹过程中的距离变化

gmx distance -s md_0_10.tpr -f md_0_10_center.xtc -select 'resname "JZ4" and name OAB plus resid 102 and name OE1' -oall

平均距离为0.31±0.05 nm，符合氢键形成条件。

GROMACS的hbond模块将其定义为氢-供体-受体夹角，该值应≤30°。执行分析前需先创建索引组——供体原子组（需同时包含供体重原子和键合氢原子）及受体原子组：

gmx make_ndx -f em.gro -o index.ndx
...
 > 13 & a OAB | a H12
 (creates group 23)
 > 1 & r 102 & a OE1
 (creates group 24)
 > 23 | 24
 > q

使用angle模块计算三原子夹角：

gmx angle -f md_0_10_center.xtc -n index.ndx -ov angle.xvg

计算RMSD

gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg

计算得RMSD约为0.1 nm（1 Å），表明配体位置仅发生微小变化。

蛋白质-配体相互作用能

为量化JZ4与T4溶菌酶的结合强度，可计算两者间的非键相互作用能。GROMACS能分解任意定义组间的短程非键作用能（需注意：该值并非自由能或结合能）。CHARMM力场因针对显式水相互作用能进行参数化，故其相互作用能具参考价值。

通过.mdp文件中的energygrps关键词计算相互作用能（注意：此计算不应作为常规模拟部分。能量分解不支持GPU运行且会拖慢计算速度）。创建包含energygrps = Protein JZ4的新.tpr文件：

gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr

使用-rerun选项基于现有轨迹重新计算能量：

gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu

通过energy模块提取Coul-SR:Protein-JZ4和LJ-SR:Protein-JZ4能量项：

gmx energy -f ie.edr -o interaction_energy.xvg

测得短程库仑作用能为-20.5±7.4 kJ mol⁻¹，短程Lennard-Jones作用能为-99.1±7.2 kJ mol⁻¹。

生活很好，有你更好

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据分析

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

登录后参与评论

0 条评论

热度