通过Gromacs做Lysozyme在水中的模拟

1. 问题描述

是个课程实验，要求分析Lysozyme(溶菌酶)在水中的状态

2. 环境

安装过程参考 我的另一篇文章 在Windows的Subsystem Linux（WSL）下安装gromacs，普通Linux也一样的。

3. 模拟过程

模拟过程大致如下

3.1 PDB文件下载

在 PDB的资源网站中找到1AKI的下载地址，记录。

打开服务器的shell，到你喜欢的路径。
用wger命令下载刚刚找到的 下载地址，-O是大写的O，后面的参数表示保存的文件名，最后的参数是下载网址。

wget -O 1aki.pdb https://files.rcsb.org/download/1AKI.pdb

然后就会如下所示

--2021-05-17 10:46:12--  https://files.rcsb.org/download/1AKI.pdb
正在解析主机 files.rcsb.org (files.rcsb.org)... 132.249.210.142
正在连接 files.rcsb.org (files.rcsb.org)|132.249.210.142|:443... 已连接。
已发出 HTTP 请求，正在等待回应... 200 OK
长度：未指定 [application/octet-stream]
正在保存至: “1aki.pdb”

1aki.pdb                                      [    <=>                                                                            ] 113.59K   180KB/s  用时 0.6s    

2021-05-17 10:46:19 (180 KB/s) - “1aki.pdb” 已保存 [116316]

输入 ls 即可看到刚刚下载到的pdb文件

$ ls
1aki.pdb

3.2 准备拓扑结构

3.2.1 清除结构中的水分子

这边简单介绍一下VMD，全称Visual Molecular Dynamics，是一种pdb文件的可视化软件，功能很多，可以用作模拟的可视化中。
需要的可以在官网下载界面找到适合自己的版本下载，总有一款适合你，使用方式自行查阅官网说明。

下载完之后通过VMD看起来是这样的，里面的蓝点是文件中的水分子。

第一步先去除pdb文件中的水分子，这个步骤不是通用的，可能水分子会和一些活性位点紧密结合或者是有重要作用时，是不可以的。
直接通过文本编辑器删除文件中的HOH即可。

grep -v HOH 1aki.pdb > 1AKI_clean.pdb

去除水后的可视化结果如下。

3.2.2 pdb to gmx

gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce

得到的图跟前面的一样，换了个展示方法。

这边的输出需要关注一下

Making bonds...
Number of bonds was 1985, now 1984
Generating angles, dihedrals and pairs...
Before cleaning: 5142 pairs
Before cleaning: 5187 dihedrals
Keeping all generated dihedrals
Making cmap torsions...
There are 5187 dihedrals,  426 impropers, 3547 angles
          5106 pairs,     1984 bonds and     0 virtual sites
Total mass 14313.332 a.m.u.
Total charge 8.000 e
Writing topology

Writing coordinate file...
		--------- PLEASE NOTE ------------
You have successfully generated a topology from: 1AKI_clean.pdb.
The Amber99 force field and the spce water model are used.

里面的total charge表示蛋白质基团总的带电荷量。

3.3 Define box and Solvate

在构建好拓扑结构后，需要构建一个空间来做模拟的空间，这里构建一个十二面体的空间，因为这样可以在给蛋白质分子更大的自由空间的情况下，容纳更少的水分子。

gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 0.8 -bt dodecahedron

现在还看不出来它是方的，这里用solvate填充一下水分子，看着还是一个正方形。

gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

3.4 离子平衡

接下来要做一下离子平衡，因为蛋白质里面的电荷是不平衡的，总体来讲带了一点负电荷，在topol.top里面可以看到，这边就不详细展开了。
在3.2中也有看到蛋白质带有8个电子的负电。
因此，接下来需要对电荷进行一个平衡。

在这儿下载ions.mdp这个参数文件，然后利用grompp命令来生成添加离子时需要用到的坐标文件和拓扑文件，以及包含了所有原子级别参数的ions.tpr文件。

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

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

之后会问你，选择什么换成离子，直接看上面的东西就可以看出来，这里要选溶液对吧，

Reading file ions.tpr, VERSION 2019.4 (single precision)
Reading file ions.tpr, VERSION 2019.4 (single precision)
Will try to add 0 NA ions and 8 CL ions.
Select a continuous group of solvent molecules
Group     0 (         System) has 15037 elements
Group     1 (        Protein) has  1960 elements
Group     2 (      Protein-H) has  1001 elements
Group     3 (        C-alpha) has   129 elements
Group     4 (       Backbone) has   387 elements
Group     5 (      MainChain) has   517 elements
Group     6 (   MainChain+Cb) has   634 elements
Group     7 (    MainChain+H) has   646 elements
Group     8 (      SideChain) has  1314 elements
Group     9 (    SideChain-H) has   484 elements
Group    10 (    Prot-Masses) has  1960 elements
Group    11 (    non-Protein) has 13077 elements
Group    12 (          Water) has 13077 elements
Group    13 (            SOL) has 13077 elements
Group    14 (      non-Water) has  1960 elements
Select a group: 

上面说要添加0个钠离子和8个氯离子，选择跟什么替换，所以输入13就行了，或者12也一样，因为可以看到，SOL和water是一样多的，这里可以认为他们两个指的是同一个内容，然后就整完了。

效果图如下，其中，水分子被隐藏，蛋白质展示为条状，白珠为添加后的氯离子：

3.5 能量最小化

在这里下载参数文件，里面有注释，需要的话要打开来改一改。

vim minim.mdp

vim打开可以看到，其中的emtol就是最小能量，改这个参数可以调整模拟停止是能量的阈值。

; 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       = 500.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

然后用mdrun来执行em

gmx mdrun -v -deffnm em

接下来选择一些参数来输出

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

在这个界面里能选择参数来输出，输出的结果存在了potential.xvg中，可以看到，势能减小的同时，压力也减小了。

能量最小化后的如图所示，可以看到，水分子已经均匀的分布开，整个水分子的排布呈现十二面体的样子。

3.6 运动平衡（温度平衡）

最开始的水溶液的动能不确定，需要让它有和蛋白质相同的动能才可以进行仿真。

从这里下载参数文件。

其中的

gen_vel = yes

指的是生成随机的初始速度。

最开始生成蛋白质时候的posre.itp文件的作用是锁定蛋白质的位置，这里用它来作为约束，然后用mdrun命令来进行模拟。

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

gmx mdrun -v -deffnm nvt

最后，输入命令，选择需要输出的内容

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

16 0

得到如图所示的温度曲线

3.7 压力平衡

从这里下载参数文件，比上一个多出了两个参数：

continuation = yes：从NVT平衡的结构开始继续进行模拟
gen_vel = no：不生成速度，而是从轨迹文件中读取

在进行完温度平衡后，需要对体系的压力进行平衡。
与温度平衡类似，执行下列命令进行平衡：

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

gmx mdrun -v -deffnm npt

而后，通过如下命令进行数据输出压力和密度

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

18 24 0

此时得到的分子的速度特性符合了温度平衡和压力平衡。
下图是平衡过程中的密度和压力。

3.8 进行分子在水中的模拟

在这里下载参数文件，其中，参数文件的steps可以根据所需要模拟的时间来调整下。
文件长这样，我需要调整时间到2ns，因此采用了这个。

title                   = OPLS Lysozyme NPT equilibration
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 1000000    ; 2 * 1000000 = 100000 ps (2 ns)
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 0         ; suppress bulky .trr file by specifying
nstvout                 = 0         ; 0 for output frequency of nstxout,
nstfout                 = 0         ; nstvout, and nstfout
nstenergy               = 5000      ; save energies every 10.0 ps
nstlog                  = 5000      ; update log file every 10.0 ps
nstxout-compressed      = 5000      ; save compressed coordinates every 10.0 ps
compressed-x-grps       = System    ; save the whole system
; Bond parameters
continuation            = yes       ; Restarting after NPT
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighborsearching
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K

然后，仍然用grompp将不同数据整合到tpr文件中。

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

然后进行模拟

gmx mdrun -v -deffnm md_0_1

接下来使用trjconv命令将轨迹文件输出

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

可以看到RMSD先升高再稳定住，在0.1ns左右的位置稳定住了，说明结构比较稳定。