Gromacs：水中溶菌酶的分子动力学模拟

前言

看了几篇文献后，先对MD模拟进行了个简单的实操，对水中溶菌酶进行分子动力学模拟，pdb文件已经通过Pymol进行了处理，删除了水分子，对这个实操做了笔记。

操作过程

产生拓扑文件

使用pdb2gmx模块产生拓扑文件。在WSL中打开windows系统中的文件系统：

cd /mnt

输入命令：

gmx pdb2gmx -f 1aki_clean.pdb -o 1aki_processed.gro -water spce -ignh

出现立场选择提示：

选择15(现在不懂)。输出结果：

工作目录产生三个新的文件： 1aki_processed.gro:Gromacs的结构文件格式，包含力场中定义的所有原子 topol.top：体系的拓扑文件 posre.itp：包含了用于限制重原子位置的信息

定义盒子并溶剂化

模拟一个简单的蛋白水溶液体系，可在适当的情况下模拟其他不同溶剂中的蛋白质或其他分子，只要有合适的力场参数即可。首先使用editconf模块定义盒子的大小，输入以下命令：

gmx editconf -f 1aki_processed.gro -o 1aki_newbox.gro -c -d 1.0 -bt cubic

使用solvate模块给盒子中填充溶剂(水)，输入以下命令：

gmx solvate -cp 1aki_newbox.gro -cs spc216.gro -o 1aki_solv.gro -p topol.top

其中spc216.gro为溶剂构型文件，为Gromacs自带的。可视化软件Pymol查看溶剂化体系：

添加抗衡离子

由pdb2gmx模块的输出结果，可以看到蛋白质带有+8e的净电荷。分定动力学模拟通常是在电中性下进行的，因此需要加入抗衡离子。

接下来添加抗衡离子：首先建立参数文件ions.mdp，内容如下：(文件内容含义待学习)

; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm （steep = steepest descent minimization）
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of （minimization） steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme  = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list （simple, grid）
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

使用grompp模块生成.tpr文件:

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

产生了.tpr类型的未见，提供了体系的原子级别的描述，包含了模拟参数说明。

之后使用genion模块添加抗衡离子，将一些水分子替换为离子，输入命令：

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

选择溶液13

能量最小化

分子动力学模拟之前，要保证体系结构正常，几何构型合理，原子无冲撞，此时需要先将溶剂化体系弛豫至能量较低的状态，即能量最小化。需要一个参数文件minim.mdp，内容如下：

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm （steep = steepest descent minimization）
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of （minimization） steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list （simple, grid）
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

使用grompp模块将结构、拓扑和模拟参数整合到.tpr文件中，输入命令：

gmx grompp -f minim.mdp -c 1aki_solv_ions.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

得到四个文件： em.log：EM过程log文件 em.edr：二进制能量文件 em.trr：二进制全精度坐标文件 em.gro：能量最小化后的解构有两个因素评价EM是否成功：第一个是势能：Epot为负，丙炔对于水中的蛋白质在10^5~10^6量级，具体取决于系统大小和水分子数量第二个是Fmax，其目标在minim.mdp中设置，emtol=1000.0，意思是目标Famx的值要低于1000KJ/(mol nm)。分析em.edr文件中包含了Gromacs在EM期间收集的所有能量项。采用energy模块：

gmx energy -f em.edr -o potential.xvg

生成图如下：

这是系统处于能量最小值，开始模拟。

平衡

EM确保我们在构象和溶剂取向方面具有合理的起始结构。在开始正式的动力学模拟之前，必须平衡蛋白质周围的溶剂和离子。如果这时候直接进行都力学模拟，体系可能会崩溃。原因是溶剂一般都在自身内部优化结构，而不一定需要跟溶质一起优化。需要将溶剂调整到模拟所需的温度然后在溶质周围形成温度的构象。在达到准确温度之后（基于动能），我们会给体系施加压力指定体系密度达到合理值。

平衡过程锋味两个阶段，第一阶段实在NVT系综(正则系综，等温等容系综)下进行。其中体系的粒子数(N)，体积(V)和平均温度(T)保持不变。如果温度尚未稳定，则需要额外的时间，一般50-100ps.进行100ps的NVT平衡。需要一个nvt.mdp参数文件。(Index of /gmx/lysozyme/Files (mdtutorials.com)) 输入以下命令执行NVT平衡：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt

第一条命令运行结果如下：

第二条命令执行需要时间，结果如下：

分析温度变化过程，使用energy模块：

gmx energy -f nvt.edr -o temperature.xvg

输入16 0已选择体系温度并退出。

结构作图：

体系温度迅速到达目标值(300K)，然后保持稳定。上一步稳定体系的温度，接下来稳定体系的压力。压力的平衡是在NPT系综中进行的，此时粒子数目、压力和温度保持不变。NPT系综也叫等温等压系综，是最解决实验条件的系综。使用100ps的NPT平衡的npt.mdp参数文件，其中增加了压力耦合选项，使用的是Parrinello-Rahman压浴，其他一些更改：

continuation=yes：我们正在从NVT平衡阶段继续仿真
gen_vel=no：不生成速度，而是从轨迹文件中读取跟NVT平衡模拟一样，还是使用grompp和mdrun。

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tp
gmx mdrun -deffnm npt

分析压力变化过程：

gmx energy -f npt.edr -o pressure.xvg

选择18 0选择体系压力并退出，作图如下：

分析密度，使用energy模块输入24 0：

gmx energy -f npt.edr -o density.xvg

作图如下：

上述两图说明压力和体系达到了很多的平衡。

MD模拟

完成两步平衡后，体系处于温度和压力下的平衡状态，现在我们将位置限制取消，进行MD模拟。运行一个1ns的MD模拟，参数文件为md.mdp：

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr

执行mdrun：

gmx mdrun -deffnm md_0_1

分析

对体系进行分析。介绍一些工具：第一个是trjconv，该后处理工具可以提取坐标、纠正周期性或者手动调整轨迹(时间单位、帧频等)。在本例中，将使用trjconv消除体系的周期性。在模拟中，蛋白质在扩散过程中可能会穿过盒子边界，这可能会导致蛋白质的一部分出现在盒子的另一侧，考虑到这种行为，使用下面的命令：

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center

居中的组为1(Protein)，输出的组为0(System)。对修正的轨迹进行所有的分析。

先看结构稳定性，Gromacs有计算rmsd计算的内置使用程序，成为rms。输入命令：

gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns

计算结晶结构的RMSD，输入以下命令：

gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns

两个结果放在一起：

表明结构非常稳定，图之间的细微差异表明，t=0ns时的结构与该晶体结构略有不同。蛋白质的均方回转半径(RG)可以描述其紧密程度。如果蛋白质处于稳固的折叠状态，Rg会处于有一个相对稳定值。如果蛋白没有折叠，Rg会一直变化。分析Rg：

gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg

选择1(Protein)进行分析：

绘图如下：

蛋白质在300K下以1ns的速度保持非常稳定，以这得形式存在。这一结果并不出乎意料。

总结

经过实操，对MD模拟有了基本的认识，也深感自己的不足，需要继续看文献以及阅读Gromacs手册，学习软件的使用方法以及丰富自己的理论知识。其次，想要利用Pymol做一个模拟后的轨迹动画，还在学习，对软件的熟练使用仍需加强。
PS:Pymol是世界上最烂的软件

参考

[1] GROMACS Tutorial

前言