【GROMACS教程】水中溶菌酶

前言

教程原址

推荐版本

教程原作者推荐使用的GROMACS版本为2018.x。经验证，使用2023.2版本同样能够得到预期的结果。

获取帮助

在本教程中，我们将用到GROMACS的不同模块。如需获取任何模块的帮助信息，请使用以下命令：

gmx help <module>

或者：

gmx <module> -h

---

准备蛋白质的拓扑文件

获取鸡蛋清溶菌酶的PDB文件

鸡蛋清溶菌酶的PDB代码为1AKI。首先，我们需要从RCSB网站搜索并下载这个蛋白质的结构。你可以通过以下链接下载鸡蛋清溶菌酶PDB文件：
点击下载 1aki.pdb
使用PyMOL软件打开下载的1aki.pdb文件，对蛋白质结构进行可视化：

图1：蛋白质结构可视化

去除结晶水

蛋白质结构中红色的“+”代表结晶水，我们需要去除这些结晶水。使用grep命令可以轻松地去除所有包含"HOH"的行：

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

需要注意的是如果结构中存在紧密结合的或者处于其他功能性活性位点的水分子，你可能需要使用纯文本编辑器来手动编辑，以防去除这些对蛋白质功能有贡献的水分子。
使用PyMOL检查去除结晶水后的蛋白质结构，这时的PDB文件只包含组成蛋白质的原子：

图2：去除结晶水后的蛋白质结构

将PDB文件转化成GROMACS可识别的拓扑结构和坐标文件

接下来，我们将蛋白质结构输入到GROMACS的第一个模块——pdb2gmx。这一步将生成三个文件：

• 分子的拓扑结构 (topol.top)：包含了模拟中定义分子所需的所有信息，包括非键参数（原子类型和电荷）以及键参数（键、角度和二面角）
• 位置限制文件
• 后处理结构文件

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

图3：选择力场

这里需要选择力场。力场包含了将要写入拓扑文件的信息，本例中我们使用全原子OPLS力场，即$17$号力场：

17

这一步得到三个文件：

• topol.top：系统拓扑文件
• posre.itp：包含了用于限制重原子位置的信息
• 1AKI_processed.gro：GROMACS格式的结构文件，包含了定义在力场内的所有原子

pdb2gmx的主要目的是产生符合力场的拓扑，另外输出力场中的蛋白质结构。

检查输出的蛋白质结构

我们用PyMOL检查输出的蛋白质结构（1AKI_processed.gro）：

图4：输出的蛋白质结构

---

检查分子的拓扑文件

拓扑文件概览

现在，让我们深入了解拓扑文件（topol.top）的内容。使用纯文本编辑器打开并检查其内容。在几行注释之后，我们可以看到以下行：

图5：`topol.top`中的力场参数

这一行调用了力场OPLS-AA的参数。它位于文件的开头，表明接下来所有的参数都源自这个力场。

分子类型定义

下一个重要的部分是[ moleculetype ]：

图6：`topol.top`中的分子类型

这里，“Protein_chain_A”定义了分子的名字（在PDB文件中蛋白质被标记为A链）。

原子定义

接下来的部分[ atoms ]定义了蛋白质中的原子：

图7：`topol.top`中的原子定义

• nr：原子编号
• type：原子类型
• resnr：氨基酸残基序号
• residue：氨基酸残基名
• atom：原子名
• cgnr：电荷组序号
• charge：电荷
• mass：质量
• typeB，chargeB，massB：用于自由能微扰

后续部分

后续的部分包括[ bonds ]，[ pairs ]，[ angles ]和[ dihedrals ]等，这些定义了分子中的化学键和角度。

位置限制

该文件的其余部分定义了一些其他有用/必要的拓扑结构。首先是位置限制：

图8：`topol.top`中的位置限制

posre.itp定义了一个力常数，用于在平衡过程中将原子保持在原位。

分子类型“Protein_chain_A”的定义到此结束，拓扑文件的其余部分定义了其他粒子并提供了系统水平的描述。

溶剂定义

下一个分子类型是溶剂，在本例中是“SPC/E”（水），通过向pdb2gmx传递参数-water spce来进行选择。水的其他典型的选择包括SPC，TIP3P和TIP4P。

图9：`topol.top`中的溶剂定义

水的位置也可以被约束，力常数（$k_{pr}$）为$1000kJ\cdot mo{l}^{-1}\cdot n{m}^{-2}$。

离子参数

接下来是离子参数：

图10：`topol.top`中的离子参数

系统定义

最后是系统水平的定义。[ system ]直接给出了模拟过程中写入输出文件的系统的名称，[ molecules ]直接列出了系统中所有的分子：

图11：`topol.top`中的系统的定义

---

定义盒子并填充溶剂

在本例中，我们将模拟一个简单的水系统。定义盒子并填充溶剂分为两步：

1. 使用editconf模块定义盒子的尺寸
2. 使用solvate模块将盒子填充水

定义盒子

本例中我们使用简单的立方盒子作为晶胞：

gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

• -c：以蛋白质作为盒子的中心
• -d 1.0：将蛋白质放置于距离盒子边界至少$1.0nm$处
• -bt cubic：盒子的类型为立方形

在周期性边界条件下，蛋白质永远不应看到其周期性图像。溶质盒子的距离为$1.0nm$意味着两个蛋白质的周期性图像之间至少间隔$2nm$。

填充溶剂

接下来将盒子填充溶剂：

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

• -cp 1AKI_newbox.gro：蛋白质的配置信息
• -cs spc216.gro：溶剂的配置信息

填充溶剂将会在拓扑文件的[ molecules ]部分增加新内容：

图12：填充溶剂后`topol.top`中的新增部分

结果可视化

使用PyMOL或其他可视化工具查看填充溶剂后的蛋白质结构（1AKI_solv.gro）：

图13：填充溶剂后的蛋白质结构

---

添加离子

在分子动力学模拟中，添加离子是平衡系统电荷的重要步骤。

蛋白质的净电荷

拓扑文件中[ atoms ]部分的最后一行有着“qtot 8”的字样，这表示蛋白质具有净电荷$+8e$。

嵌入离子

genion模块将扫描拓扑文件，将其中的水分子替换成指定的离子。

生成.tpr文件

.tpr文件是这一步的输入文件，通过grompp模块（GROMACS pre-processor）产生。这个模块处理坐标文件和拓扑文件以产生原子级别的输入（.tpr），.tpr文件包含了系统中所有原子的所有参数。

获取参数文件

.mdp文件（molecular dynamics parameter file）通常用于运行能量最小化或进行分子动力学模拟。本例中仅用于产生对整个系统的原子描述。你可以从以下链接下载示例参数文件：
点击下载示例参数文件 ions.mdp

使用grompp模块

当我们有了参数文件后，我们就可以通过grompp模块产生.tpr文件：

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

• -p topol.top：过程中需要对拓扑文件进行处理，去除水分子，添加离子

使用genion模块

产生的二进制文件ions.tpr中包含了对系统的原子级描述，将这个文件传递给genion模块：

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

• -s ions.tpr：结构/状态文件，作为输入
• -pname：定义正离子
• -nname：定义负离子
• -neutral：系统呈电中性
• -conc：指定离子的浓度
  

图14：选择用离子替换的组分

选择替换溶剂（13 ( SOL)）来嵌入离子：

13

genion模块加入必要的离子来平衡蛋白质的电荷，使系统净电荷为$0$，本例中为$8$ Cl^-。

系统中的粒子组成

现在系统中的粒子组成如下：

图15：添加离子后`topol.top`中的新增部分

结果可视化

使用PyMOL或其他可视化工具查看添加离子后的蛋白质结构（1AKI_solv_ions.gro）：

图16：添加离子后的蛋白质结构

---

能量最小化

溶剂化的、电中性的系统已经组装好了，但是在开始模拟前，必须确保没有空间冲突或不恰当的几何结构。通过能量最小化（energy minimization，EM）来使结构松弛。
能量最小化（EM）过程与添加离子的过程相似，仍需使用grompp模块组装结构和拓扑，并将模拟参数写入二进制输入文件（.tpr）。但在这一步，我们将.tpr文件传递给mdrun模块进行能量最小化。

获取参数文件

你可以从以下链接下载示例参数文件：
点击下载输入参数文件 minim.mdp

使用grompp模块

将输入参数文件传递给grompp模块来组装二进制输入文件。

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

使用mdrun模块

然后使用mdrun来进行能量最小化：

gmx mdrun -v -deffnm em

• -v参数会将模拟过程打印在屏幕上
• -deffnm参数指定了输入和输出文件的文件名
• 如果输入文件不是em.tpr，那么可以通过-s参数来指定。

使用GPU进行运算

如果在这一步使用GPU进行运算，那么需要通过-ntmpi指定使用GPU的数量以及-ntomp指定每个GPU处理几个线程，例如以下命令表示使用$4$个GPU：

gmx mdrun -v -deffnm em -ntmpi 4

结果文件

在本例中，我们将得到以下文件：

• em.log：ASCII编码形式的文本文档，记录了EM的过程
• em.edr：二进制能量文件
• em.trr：二进制完整的轨迹文件
• em.gro：能量最小化后的结构

结果可视化

使用PyMOL或其他可视化工具查看能量最小化后的结构（em.gro）：

图17：能量最小化后的结构

---

分析能量最小化

分析能量最小化的结果时，有两个关键点需要关注：

• 势能分析
• 最大力分析

势能分析

首先是势能（$E_{pot}$），它应该为负，并且对于水中的简单蛋白质，其数量级应为（$10^5\sim 10^6$），这与系统的大小和水分子的数量有关。

最大力分析

另一点是最大力（$F_{max}$），其目标在minim.mdp中设定为emtol = 1000.0，表明$F_{max}$不大于$1000kJ\cdot mo{l}^{-1}\cdot n{m}^{-1}$。
当$F_{max}>emtol$时，仍有可能得到一个合理的$E_{pot}$，但这样系统可能不够稳定。

分析能量文件

em.edr 文件包含了能量最小化过程中所有与能量有关的内容。使用 energy 模块对其进行分析：

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

图18：选择要对`em.edr`分析的数据

选择分析势能（10 Potential），0 没有特殊含义，表示输入终止：

10 0

这将会显示$E_{pot}$的平均值，同时得到文件 potential.xvg。

图19：势能分析结果

势能图像绘制

使用Xmgrace软件绘图，图像展示了$E_{pot}$良好、稳定的收敛：

图20：势能收敛图

结论

现在系统达到了能量最小状态，可以进行下一步分子动力学模拟。

---

平衡

平衡分为两步进行：

1. NVT平衡
2. NPT平衡

NVT平衡

能量最小化（EM）确保了我们在几何和溶剂方向上有了正确的结构。为了开始模拟，我们必须平衡蛋白质周围的溶剂和离子，否则如果直接在这里尝试无限制动力学，系统将崩溃。原因是溶剂主要进行了自身内部的优化，而不一定与溶质一起优化，因此必须被调整到我们希望进行模拟的温度并建立溶质（蛋白质）的正确方向。当达到正确的温度（基于动能）后我们将对系统施加压力直到达到合适的密度。

使用位置限制

这里我们将用到pdb2gmx产生的文件posre.itp（position restrain），这个文件对蛋白质中的重原子（所有非氢原子）施加位置限制力。只有在克服相当大的能量惩罚之后这些原子才允许移动。位置限制允许我们在蛋白质周围平衡溶剂，而不会因为蛋白质的结构变化而增加变量。位置限制的起点（即限制势能为零的坐标）是通过传递给grompp的-r参数的坐标文件提供的。

平衡过程

平衡通常分两个阶段进行。第一个阶段在NVT系综（恒定粒子数量、体积和温度，也称“等温等容”或“正则”系综）下进行。这一过程所用的时间长度取决于系统的内容。但是在NVT系综中，系统的温度应该稳定在期望值，如果温度还没有稳定，就需要额外的时间。通常$50\sim 100ps$足够了，本例中我们将进行$100ps$ NVT平衡。

获取参数文件

从以下链接下载参数文件：
点击下载参数文件 nvt.mdp

使用grompp和mdrun模块进行平衡

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

• gen_vel = yes：启动速度生成，使用不同的随机种子（gen_seed）会给出不同的初始速度，这样就能从相同的起始结构进行多次不同的模拟
• tcoupl = V-rescale：速度重置恒温器，是对Berendsen弱耦合方法的改进，后者没有正确地再现正确的动力学系综

可能会花点时间。这一步我使用了$1$个GPU和$11$个OpenMP线程，花费$14.623s$的墙钟时间来完成。

分析温度变化

现在可以利用energy模块分析温度的变化过程：

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

图21：选择要对`nvt.edr`分析的数据

选择分析系统的温度（16 Temperature）：

16 0

图22：温度分析结果

温度变化图像

利用Xmgrace软件绘图。从图中可以看出，系统的温度迅速达到目标值（$300K$），并在之后的平衡过程中保持稳定：

图23：温度变化图

NVT平衡到此完成，接下来我们要进行NPT平衡。

NPT平衡

NVT平衡稳定了系统的温度，接下来我们要在NPT系综（恒定粒子数量、压力和温度，也称“等温等压”系综，最接近实验条件）下稳定系统的压力（密度）。

获取参数文件

在这里下载用于$100ps$ NPT系综的参数文件：
点击下载 npt.mdp

参数文件差异

它与用于NVT平衡的参数文件并没有很大的不同。注意增加了使用Parrinello-Rahman压力控制器的压力耦合部分。一些其他的改变：

• continuation = yes：从NVT平衡阶段继续进行模拟
• gen_vel = no：速度从轨迹文件中读取

使用grompp和mdrun模块进行平衡

同样地，使用grompp和mdrun模块进行平衡，但是注意我们要使用-t参数来包含NVT平衡的检查点文件（包含了用于继续模拟所需的所有状态变量），这是为了保持NVT期间产生的速度。-c参数指定的坐标文件是NVT模拟的最终输出。

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

这一步同样要花点时间。我使用了$1$个GPU和$11$个OpenMP线程，花费$15.668s$的墙钟时间完成了模拟。

分析压力变化

再次通过energy模块分析压力的变化过程：

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

图24：选择要对`npt.edr`分析的数据

分析压力的变化过程：

18 0

图25：压力分析结果

压力变化图像

使用Xmgrace作图。整个$100ps$平衡过程中压力值波动很大。运行平均值以红线标出：

图26：压力变化图

平衡过程中，压力的平均值是$14.7\pm 205.2bar$。尽管参考压力为$1bar$，但是压力是分子动力学模拟过程中波动非常大的一个量，因此从统计学上将，这并不表明我们获得的压力（$14.7\pm 205.2bar$）和参考压力值（$1bar$）之间存在显著差异。

分析密度变化

现在再来看一下密度的变化过程：

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

图27：选择要对`npt.edr`分析的数据

选择分析密度：

24 0

图28：密度分析结果

密度变化图像

利用Xmgrace进行绘图分析：

图29：密度变化图

$100ps$模拟过程中，密度的平均值为$1019\pm 5kg\cdot m^{-3}$，接近实验值（$1000kg\cdot m^{-3}$）和SPC/E模型的期望密度为（$1008kg\cdot m^{-3}$），并且随着时间的推移，密度值非常稳定，表明系统在压力和密度上已经达到了良好的平衡。

---

MD模拟

经过两个阶段的平衡，系统达到了预期的温度和压力并保持良好的平衡，现在我们可以解除位置限制并进行分子动力学（MD）模拟以收集数据。

获取参数文件

在这里下载参数文件：
点击下载 md.mdp

使用grompp准备模拟

使用检查点文件（包含保持压力耦合信息）来进行grompp，运行$1ns$分子动力学模拟：

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

grompp将会给出PME（粒子-网格-Ewald）负载的估计值，这个信息指出多少处理器将会被用于PME计算，以及多少处理器用于PP（粒子-粒子）计算：

图30：PME负载估计值

对于立方盒子，最佳设置的PME负载为$0.25$（$PP:PME=3:1$），非常接近我们所得到的值。在执行mdrun时，程序会自动决定用于执行PME和PP运算的的最佳处理器数量，因此为了获得最佳的性能，要确保为计算指定了合适的线程/核心数（通过-nt参数指定）。

执行mdrun

现在可以执行mdrun了：

gmx mdrun -deffnm md_0_1

模拟过程将需要较长的时间，其中非键相互作用和PME计算可以在GPU上进行来加速模拟过程（键力计算在CPU上运行）。如果你有可用的GPU，可以通过以下指令来使用：

gmx mdrun -deffnm md_0_1 -nb gpu

这一步我使用$1$个GPU和$11$个线程，经过$125.108s$的墙钟时间来完成模拟过程。

---

分析

trjconv模块是一个后处理工具，用于剔除坐标、校正周期性或者手动改变轨迹（时间单位、帧频等）。在本例中，我们将用trjconv来解释系统中的任何周期性。蛋白质会通过晶胞扩散，可能会看起来“断裂”或者可能会“跳”到盒子的另一边。为了解释这种行为，使用以下指令：

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

图31：选择蛋白质作为要居中的组

这里选择蛋白质（1 ( Protein)）作为要居中的组：

1

图32：选择系统作为输出

选择系统（0 ( System)）作为输出：

0

我们要在这条“校正”轨迹上进行我们所有的分析。
首先来看一下结构稳定性。使用rms模块进行RMSD（Root-Mean-Square Deviation，均方根偏差）计算：

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

-tu ns参数将以ns为单位输出结果，这是为了提高图像清晰度（比如对于较长时间的模拟——在图像上显示$1e+05ps$并不如$100ns$好看）。

图33：选择主链作为最小二乘拟合的组

选择主链（4 ( Backbone)）作为最小二乘拟合和RMSD计算的组：

4

图34：选择主链作为RMSD计算的组

4

输出的图表将显示相对于在最小化平衡系统中存在的结构的RMSD。如果希望计算相对于晶体结构的RMSD，可以使用以下指令：

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

然后同样选择主链作为最小二乘拟合和RMSD计算的组。
利用Xmgrace软件绘图：

图35：RMSD分析图

两个时间序列都显示RMSD水平下降至$\sim 0.1nm$（$1\mathring{A}$），表明结构非常稳定。两图之间的细微差别表明$t=0ns$处的结构与该晶体结构略有不同。这是预料之中内的，因为系统中的结构被能量最小化了，而且位置限制并不是$100\%$完美的。

---

其他分析

回转半径分析

蛋白质的回转半径（radius of gyration，$R_g$）是衡量其紧密度的一个指标。如果蛋白质紧密折叠，其将可能保持一个相对稳定的$R_g$值。如果蛋白质展开，其$R_g$将会随着时间变化。现在我们计算模拟中溶菌酶的回转半径：

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

选择蛋白质作为待分析组（1 ( Protein)）：

1

同样地，使用Xmgrace进行绘图：

图36：回转半径分析图

我们从相对不变的$R_g$值可以看出蛋白质在$300K$下$1ns$的模拟过程中保持非常稳定，处于紧凑（折叠）状态。

制作轨迹动画

使用以下命令对MD模拟过程中粒子的轨迹进行抽帧：

gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o traj.pdb -dt 10

• -dt 10：抽帧间隔为$10ps$，对于本例中的$1ns$模拟，帧数为$101$
  

图37：选择一组作为抽帧的输出

这里需要选择作为输出的组，为了可视化整个系统中粒子的运动，我选择系统作为输出（0 ( System)）：

0

现在我们得到了包含MD模拟过程中系统内粒子运动的轨迹信息的文件——traj.pdb，可以使用PyMOL制作运动轨迹动画：

图38：MD模拟过程系统中粒子的运动轨迹

---

GROMACS分子动力学模拟流程概述：

1. 使用pdb2gmx处理.pdb文件：
   • 将PDB文件转换为GROMACS可识别的格式，并生成拓扑文件。
2. 使用editconf定义盒子：
   • 定义模拟盒子的尺寸和类型，以容纳蛋白质。
3. 使用solvate填充溶剂：
   • 在模拟盒子中添加溶剂分子。
4. 使用grompp产生.tpr文件：
   • 准备运行文件，将坐标文件、拓扑文件和模拟参数结合生成.tpr文件。
5. 使用genion添加离子：
   • 向系统中添加离子以平衡电荷。
6. 能量最小化（EM）：
   • 先用grompp产生.tpr文件，再使用mdrun进行能量最小化，确保系统没有几何冲突。
7. NVT平衡：
   • 先用grompp产生.tpr文件，再用mdrun进行热浴平衡，稳定系统的温度。
8. NPT平衡：
   • 与NVT平衡步骤基本相同，但还包括压力的平衡。
9. 分析：
   • 使用trjconv校正周期性：消除由于周期性边界条件引起的原子穿越盒子界面的假象。
   • 使用rms分析RMSD：计算均方根偏差，评估蛋白质结构的稳定性。
   • 使用gyrate分析回转半径：计算蛋白质的回转半径，评估其紧凑度。

如果在操作过程中遇到任何问题，你可以访问GROMACS用户论坛来获取帮助。
你也可以通过微信公众号 AnzThoughts 私信我：