分子动力学--蛋白配体模拟
作者，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>ATOM1 C4 24.2940 -24.1240 -0.0710 C.3 167 JZ4167 -0.06502 C7 21.5530 -27.2140 -4.1120 C.ar 167 JZ4167 -0.06133 C8 22.0680 -26.7470 -5.3310 C.ar 167 JZ4167 -0.05834 C9 22.6710 -25.5120 -5.4480 C.ar 167 JZ4167 -0.01995 C10 22.7690 -24.7300 -4.2950 C.ar 167 JZ4167 0.12006 C11 21.6930 -26.4590 -2.9540 C.ar 167 JZ4167 -0.05517 C12 22.2940 -25.1870 -3.0750 C.ar 167 JZ4167 -0.00608 C13 22.4630 -24.4140 -1.8080 C.3 167 JZ4167 -0.02459 C14 23.9250 -24.7040 -1.3940 C.3 167 JZ4167 -0.051810 OAB 23.4120 -23.5360 -4.3420 O.3 167 JZ4167 -0.506511 H 25.3133 -24.3619 0.1509 H 1 UNL1 0.023012 H 23.6591 -24.5327 0.6872 H 1 UNL1 0.023013 H 24.1744 -23.0611 -0.1016 H 1 UNL1 0.023014 H 21.0673 -28.1238 -4.0754 H 1 UNL1 0.061815 H 21.9931 -27.3472 -6.1672 H 1 UNL1 0.061916 H 23.0361 -25.1783 -6.3537 H 1 UNL1 0.065417 H 21.3701 -26.8143 -2.0405 H 1 UNL1 0.062118 H 21.7794 -24.7551 -1.0588 H 1 UNL1 0.031419 H 22.2659 -23.3694 -1.9301 H 1 UNL1 0.031420 H 24.5755 -24.2929 -2.1375 H 1 UNL1 0.026621 H 24.0241 -25.7662 -1.3110 H 1 UNL1 0.026622 H 23.7394 -23.2120 -5.1580 H 1 UNL1 0.2921
@<TRIPOS>BOND1 4 3 ar2 4 5 ar3 3 2 ar4 10 5 15 5 7 ar6 2 6 ar7 7 6 ar8 7 8 19 8 9 110 9 1 111 1 11 112 1 12 113 1 13 114 2 14 115 3 15 116 4 16 117 6 17 118 8 18 119 8 19 120 9 20 121 9 21 122 10 22 1
需要修改的第一个地方是MOLECULE标题行。将"*****"替换为"JZ4"
@<TRIPOS>MOLECULE
JZ4
接下来需要统一修正残基名称和编号。由于该.mol2文件仅包含一个分子，因此应只指定一个残基名称和编号。编辑完成后，jz4.mol2文件中修正后的ATOM部分应显示为：
@<TRIPOS>ATOM1 C4 24.2940 -24.1240 -0.0710 C.3 1 JZ4 -0.06502 C7 21.5530 -27.2140 -4.1120 C.ar 1 JZ4 -0.06133 C8 22.0680 -26.7470 -5.3310 C.ar 1 JZ4 -0.05834 C9 22.6710 -25.5120 -5.4480 C.ar 1 JZ4 -0.01995 C10 22.7690 -24.7300 -4.2950 C.ar 1 JZ4 0.12006 C11 21.6930 -26.4590 -2.9540 C.ar 1 JZ4 -0.05517 C12 22.2940 -25.1870 -3.0750 C.ar 1 JZ4 -0.00608 C13 22.4630 -24.4140 -1.8080 C.3 1 JZ4 -0.02459 C14 23.9250 -24.7040 -1.3940 C.3 1 JZ4 -0.051810 OAB 23.4120 -23.5360 -4.3420 O.3 1 JZ4 -0.506511 H 25.3133 -24.3619 0.1509 H 1 JZ4 0.023012 H 23.6591 -24.5327 0.6872 H 1 JZ4 0.023013 H 24.1744 -23.0611 -0.1016 H 1 JZ4 0.023014 H 21.0673 -28.1238 -4.0754 H 1 JZ4 0.061815 H 21.9931 -27.3472 -6.1672 H 1 JZ4 0.061916 H 23.0361 -25.1783 -6.3537 H 1 JZ4 0.065417 H 21.3701 -26.8143 -2.0405 H 1 JZ4 0.062118 H 21.7794 -24.7551 -1.0588 H 1 JZ4 0.031419 H 22.2659 -23.3694 -1.9301 H 1 JZ4 0.031420 H 24.5755 -24.2929 -2.1375 H 1 JZ4 0.026621 H 24.0241 -25.7662 -1.3110 H 1 JZ4 0.026622 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.126163ASN O1 1692 0.624 -0.616 -0.140163ASN O2 1693 0.683 -0.703 -0.0111JZ4 C4 1 2.429 -2.412 -0.0071JZ4 C7 2 2.155 -2.721 -0.4111JZ4 C8 3 2.207 -2.675 -0.5331JZ4 C9 4 2.267 -2.551 -0.5451JZ4 C10 5 2.277 -2.473 -0.4301JZ4 C11 6 2.169 -2.646 -0.2951JZ4 C12 7 2.229 -2.519 -0.3081JZ4 C13 8 2.246 -2.441 -0.1811JZ4 C14 9 2.392 -2.470 -0.1391JZ4 OAB 10 2.341 -2.354 -0.4341JZ4 H1 11 2.531 -2.436 0.0151JZ4 H2 12 2.366 -2.453 0.0691JZ4 H3 13 2.417 -2.306 -0.0101JZ4 H4 14 2.107 -2.812 -0.4071JZ4 H5 15 2.199 -2.735 -0.6171JZ4 H6 16 2.304 -2.518 -0.6351JZ4 H7 17 2.137 -2.681 -0.2041JZ4 H8 18 2.178 -2.476 -0.1061JZ4 H9 19 2.227 -2.337 -0.1931JZ4 H10 20 2.458 -2.429 -0.2141JZ4 H11 21 2.402 -2.577 -0.1311JZ4 H12 22 2.374 -2.321 -0.5165.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.0gmx 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 atoms1 Protein : 2614 atoms2 Protein-H : 1301 atoms3 C-alpha : 163 atoms4 Backbone : 489 atoms5 MainChain : 651 atoms6 MainChain+Cb : 803 atoms7 MainChain+H : 813 atoms8 SideChain : 1801 atoms9 SideChain-H : 650 atoms10 Prot-Masses : 2614 atoms11 non-Protein : 30892 atoms12 Other : 22 atoms13 JZ4 : 22 atoms14 CL : 6 atoms15 Water : 30864 atoms16 SOL : 30864 atoms17 non-Water : 2642 atoms18 Ion : 6 atoms19 JZ4 : 22 atoms20 CL : 6 atoms21 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.tprgmx 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.tprgmx 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.tprgmx 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⁻¹。