Gromacs: 蛋白-配体复合物实例 (2)

2023-08-13 18:47 作者:纯粹之狐 0人读过 | 我要投稿

上回我们将拓扑文件都准备完全，并修改了top文件。接下来这一节是定义盒子添加溶剂。

此时，工作流程与其他 MD 模拟一样。我们将定义单元格并在其中注入水。

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

在查看 solv.gro 时，您可能会想，为什么 editconf 没有生成所需的十二面体单胞形状，因为系统看起来是这样的：

GROMACS 程序总是使用数值效率最高的坐标表示法，即把所有坐标重新包扎成三菱单元。mdrun 所执行的物理计算可以通过不同的坐标包裹方式等效进行，因此首选效率最高的方式。在生成 .tpr 文件后，可以恢复所需的单元格形状。

添加离子

现在我们有了一个包含带电蛋白质的溶解系统。pdb2gmx 的输出结果告诉我们，该蛋白质的净电荷为 +6e（基于其氨基酸组成）。如果您在 pdb2gmx 的输出中忽略了这一信息，请查看 topol.top 中 [ atoms ] 指令的最后一行；它应该是（部分）"qtot 6"。由于生命不是以净电荷形式存在的，因此我们必须在系统中添加离子。

使用任何 .mdp 文件，使用 grompp 生成 .tpr 文件。我使用 .mdp 文件来运行能量最小化，因为它们需要的参数最少，因此最容易维护。

该mdp文件长这样：

; LINES STARTING WITH ';' ARE COMMENTS

title = Minimization ; Title of run

; 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 < 10.0 kJ/mol

emstep = 0.01 ; Energy 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

ns_type = grid ; Method to determine neighbor list (simple, grid)

rlist = 1.0 ; Cut-off for making neighbor list (short range forces)

coulombtype = cutoff ; Treatment of long range electrostatic interactions

rcoulomb = 1.0 ; long range electrostatic cut-off

rvdw = 1.0 ; long range Van der Waals cut-off

pbc = xyz ; Periodic Boundary Conditions

接下来运行命令：

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

我们现在将 .tpr 文件传递给 genion

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

选择替换溶剂分子的group选择SOL，也就是输入15

在以前的 GROMACS 版本中，用 -pname 和 -nname 指定的离子名称是针对特定力场的，但在 4.5 版中已经标准化。指定的原子名称始终是大写的元素符号以及[ moleculetype ]。残基名称可以附加或不附加电荷符号（+/-）。如果遇到困难，请参阅 ions.itp 了解正确的命名方法。

[ molecular ] 现在应该如下所示：

[ molecules ]

; Compound #mols

Protein_chain_A 1

JZ4 1

SOL 10228

CL 6

能量最小化

现在系统已组装完毕，使用 grompp进行操作，em.mdp文件如下：

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

确保在运行 genbox 和 genion 时更新了 topol.top 文件，否则你会收到很多讨厌的错误消息（“坐标文件中的坐标数与拓扑不匹配”等）

现在我们可以调用 mdrun 来执行 EM

gmx mdrun -v -deffnm em

现在，我们的系统处于能量最小值，可以开始真正的动力学研究了。

平衡

平衡我们的蛋白质配体复合物与平衡任何其他含有水中蛋白质的系统非常相似。在这种情况下，有一些特殊的考虑因素：

1. 对配体施加限制

2. 温度耦合基团的处理

约束配体

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

gmx make_ndx -f jz4.gro -o index_jz4.ndx

然后，执行 genrestr 模块，选择新创建的索引组（即 index_jz4.ndx 文件中的第 3 组）：

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

现在，我们需要在拓扑中加入这些信息。根据我们希望使用的条件，有几种方法可以做到这一点。如果我们只想在限制蛋白质的同时限制配体，那么请在拓扑中的指定位置添加以下几行：

位置很重要！您必须在拓扑结构中调用 posre_jz4.itp。jz4.itp 中的参数定义了配体的 [ moleculetype ] 指令。分子类型以包含水拓扑（tip3p.itp）结束。在其他地方调用 posre_jz4.itp 会试图将位置限制参数应用于错误的分子类型

如果您希望在平衡过程中进行更多控制，即独立限制蛋白质和配体，您可以在另一个 #ifdef 块中控制配体位置限制文件，如图所示：

在后一种情况下，为了同时限制蛋白质和配体，我们需要在 .mdp 文件中指定 define = -DPOSRES -DPOSRES_LIG。如何处理您的系统取决于您。这些例子只是为了说明 GROMACS 提供的灵活性。对于标准平衡程序，同时限制蛋白质和配体可能就足够了。您自己的需求可能会有所不同。

温度调节

适当控制温度耦合是一个敏感问题。将每个分子类型与自己的恒温组耦合是个坏主意。例如，如果这样做

tc-grps = Protein JZ4 SOL CL

您的系统很可能会爆炸，因为温度耦合算法不够稳定，无法控制只有几个原子的原子团（即 JZ4 和 CL）所产生的动能波动。请勿单独耦合系统中的每一个种类。

典型的方法是设置 tc-grps = Protein Non-Protein，然后继续。但是，"非蛋白质 "组也包括 JZ4。由于 JZ4 和蛋白质的物理联系非常紧密，最好将它们视为一个整体。也就是说，为了温度耦合的目的，将 JZ4 与蛋白质归为一组。同样，我们插入的少量 Cl- 离子也被视为溶剂的一部分。为此，我们需要一个特殊的索引组来合并蛋白质和 JZ4。我们可以使用 make_ndx 模块，提供完整系统的任意坐标文件来实现这一目的。这里我使用的是 em.gro，即我们系统的输出（最小化）结构：