GROMACS Tutorial 5-Protein-Ligand Complex

Introduction

本例将指导使用者建立一个包含与配体配合的蛋白质（T4溶菌酶L99A/M102Q，T4 lysozyme L99A/M102Q）的模拟系统的过程。本教程特别关注处理与配体相关的问题，假设你熟悉基本的GROMACS操作和拓扑的内容。如果你还不甚了解，请在尝试这个教程之前参考一个更基本的教程。

本教程需要2018年的GROMACS.x系列版本（或更新）。

pdb2gmx - The Protein Topology

我们必须下载将要处理的蛋白质结构文件。对于本教程，我们将使用T4溶菌酶L99A/M102Q (PDB：3HTB)。去RCSB网站下载PDB文件。

一旦您下载了这个结构，您就可以使用一个可视化程序(如VMD、Chimera、PyMOL等)查看这个结构。看完结构之后你也许会想去掉结晶水，PO4和BME。注意，这种方法并不是普遍适用的（例如，结合活性位点水分子的情况）。对于本教程而言，我们不需要结晶水或其他物质，它们只是结晶的共溶剂。我们将转而关注被称为JZ4（2-丙基苯酚）的配体。

如果您想要一个已经处理好的.pdb文件来检查自己做的对不对，您可以在这里下载它。我们现在面临的问题是，JZ4配体在GROMACS提供的任何力场中都不是一个可以被识别的实体，所以如果您试图通过pdb2gmx传递这个文件，它将给出一个致命报错。只有在力场的.rtp（residue topology，各残基的拓扑结构）文件中存在基本构成要素（building block）的条目时，拓扑才能自动组装。由于情况并非如此，我们将分两个步骤准备我们的系统拓扑:

1. 使用pdb2gmx准备蛋白质拓扑
2. 使用外部工具准备配体拓扑

因为我们将分别准备这两个拓扑，所以我们必须将蛋白质和JZ4配体保存到单独的坐标文件中。保存JZ4坐标如下:

grep JZ4 3HTB_clean.pdb > jz4.pdb

然后从3HTB_clean.pdb中删除JZ4行。这时，准备蛋白质拓扑结构是很简单的。我们在本教程中使用的力场是CHARMM36，从MacKerell实验室网站获得。在那里，下载最新的CHARMM36力场 tar 包和 “cgenff_charmm2gmx.py” 转换脚本，我们将在后面使用它。下载与您安装的Python版本相对应的转换脚本版本（python2.x 或 3. x）（译者在尝试使用2021年6月的力场版本时，生成.gro文件出现了严重的报错，显示为末端GLY缺少CB原子，但是一般力场中GLY没有CB；在换用2021年4月版本之后能够流畅运行，还请读者注意）。

在你的工作目录中解压力场 tar 包：

tar -zxvf charmm36-mar2019.ff.tgz

现在应该有一个 “charmm36-mar2019.ff” 在您的工作目录下的子目录。使用pdb2gmx编写T4溶菌酶的拓扑:

gmx pdb2gmx -f 3HTB_clean.pdb -o 3HTB_processed.gro

根据提示选择默认的水模型。这个结构将由pdb2gmx处理，并且会提示您选择一个力场：

Select the Force Field:
From current directory:
 1: CHARMM36 all-atom force field (July 2017)
From '/usr/local/gromacs/share/gromacs/top':
 2: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
 3: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
 4: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
 5: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
 6: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
 7: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
 8: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
 9: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
10: GROMOS96 43a1 force field
11: GROMOS96 43a2 force field (improved alkane dihedrals)
12: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
13: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
15: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
16: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)

对于本教程，选择CHARMM36力场（选项1），“From current directory” 下的第一个。

The Ligand Topology

现在我们必须处理配体。但是，对于力场无法自动识别的那些物质，人们是如何得出相应参数的呢？配体的正确处理是分子模拟中最具挑战性的任务之一。力场的作者花费了数年的时间来推导一个自洽的力场，而将一个新物质引入这个框架并不是一项简单的任务。任何新物质的力场参数必须以与原始力场一致的方式推导和验证。

😇该教程使用的生成小分子参数的方法目前看来有些繁琐而且不能覆盖大多数情况，所以在参数生成上推荐使用sobtop，可以参考sobtop官网和这篇博客😇

对于OPLS、AMBER和CHARMM力场，这种推导通常采用各种量子力学计算的形式。这些力场的原始文献描述了所需的程序。对于GROMOS力场，参数化方法不太明确，依赖于凝聚相行为的经验拟合。也就是说，对每种原子类型计算出一些初始电荷和 Lennard-Jones 参数，评估其准确性，并加以改进。虽然最终结果非常令人满意，即流体模型与现实世界十分相似，但推导过程可能是费力和令人崩溃的。

因此，自动化工具是非常受欢迎的。对于每一个力场，都有一些方法或软件程序可以提供与各种力场兼容的参数。并不是所有的结果都一样准确。有几个例子，请参考下表：

力场	自动化工具	备注
AMBER	Antechamber	借助GAFF参数化分子
	acpype	连接Antechamber的Python接口，用于编写GROMACS拓扑
CHARMM	CGenFF	官方的CHARMM通用力场服务器（网站化服务）
GROMOS87/GROMOS96	PRODRG 2.5	用于拓扑生成的自动化服务器（网站化服务）
	ATB	用于生成适用于GROMOS96 54A7的拓扑的新服务器（网站化服务）
OPLS-AA	Topolbuild	将Tripos .mol2文件转换为拓扑
	TopolGen	一个Perl脚本，将一个全原子的.pdb文件转换为一个拓扑
	LigParGen	来自Jorgensen课题组的服务器，用于生成OPLS拓扑

Add Hydrogen Atoms to JZ4

在本教程中，我们将使用 CGenFF 生成JZ4的拓扑。点击上表中的链接访问该网站。注册一个（免费的）帐户并激活它。 CGenFF 需要一个 Sybyl .mol2 文件作为输入，以收集基本的原子类型信息和键连信息。 CHARMM 也是一个全原子力场，这意味着所有的H原子都是明确表示的。晶体结构通常不分配H坐标，所以它们必须内置。通过 Avogadro program 来生成一个.mol2文件并添加H原子。打开jz4.pdb，然后从“Build”菜单中，选择“Add Hydrogens”。 Avogadro 会在JZ4配体上建立所有的H原子。保存.mol2文件（File -> Save As…，从下拉菜单中选择Sybyl Mol2）并命名为“jz4.mol2”。

jz4. mol2必须作几处修正才能使用。在你最喜欢的纯文本编辑器中打开jz4. mol2，你会发现:

@<TRIPOS>MOLECULE
*****
 22 22 0 0 0
SMALL
GASTEIGER

@<TRIPOS>ATOM
      1 C4         24.2940  -24.1240   -0.0710 C.3   167  JZ4167     -0.0650
      2 C7         21.5530  -27.2140   -4.1120 C.ar  167  JZ4167     -0.0613
      3 C8         22.0680  -26.7470   -5.3310 C.ar  167  JZ4167     -0.0583
      4 C9         22.6710  -25.5120   -5.4480 C.ar  167  JZ4167     -0.0199
      5 C10        22.7690  -24.7300   -4.2950 C.ar  167  JZ4167      0.1200
      6 C11        21.6930  -26.4590   -2.9540 C.ar  167  JZ4167     -0.0551
      7 C12        22.2940  -25.1870   -3.0750 C.ar  167  JZ4167     -0.0060
      8 C13        22.4630  -24.4140   -1.8080 C.3   167  JZ4167     -0.0245
      9 C14        23.9250  -24.7040   -1.3940 C.3   167  JZ4167     -0.0518
     10 OAB        23.4120  -23.5360   -4.3420 O.3   167  JZ4167     -0.5065
     11 H          25.3133  -24.3619    0.1509 H       1  UNL1        0.0230
     12 H          23.6591  -24.5327    0.6872 H       1  UNL1        0.0230
     13 H          24.1744  -23.0611   -0.1016 H       1  UNL1        0.0230
     14 H          21.0673  -28.1238   -4.0754 H       1  UNL1        0.0618
     15 H          21.9931  -27.3472   -6.1672 H       1  UNL1        0.0619
     16 H          23.0361  -25.1783   -6.3537 H       1  UNL1        0.0654
     17 H          21.3701  -26.8143   -2.0405 H       1  UNL1        0.0621
     18 H          21.7794  -24.7551   -1.0588 H       1  UNL1        0.0314
     19 H          22.2659  -23.3694   -1.9301 H       1  UNL1        0.0314
     20 H          24.5755  -24.2929   -2.1375 H       1  UNL1        0.0266
     21 H          24.0241  -25.7662   -1.3110 H       1  UNL1        0.0266
     22 H          23.7394  -23.2120   -5.1580 H       1  UNL1        0.2921
@<TRIPOS>BOND
     1     4     3   ar
     2     4     5   ar
     3     3     2   ar
     4    10     5    1
     5     5     7   ar
     6     2     6   ar
     7     7     6   ar
     8     7     8    1
     9     8     9    1
    10     9     1    1
    11     1    11    1
    12     1    12    1
    13     1    13    1
    14     2    14    1
    15     3    15    1
    16     4    16    1
    17     6    17    1
    18     8    18    1
    19     8    19    1
    20     9    20    1
    21     9    21    1
    22    10    22    1

首先需要改变的是分子标题。将 ***** 替换为 JZ4 ：

@<TRIPOS>MOLECULE
JZ4

接下来，修正残留的名称和编号（第6-7列），使它们都相同。此.mol2文件只包含一个分子，因此应该只有一个确切的残基名称和编号。修正后jz4. mol2的 ATMO 部分应该与下文相似：

@<TRIPOS>ATOM
      1 C4         24.2940  -24.1240   -0.0710 C.3     1  JZ4        -0.0650
      2 C7         21.5530  -27.2140   -4.1120 C.ar    1  JZ4        -0.0613
      3 C8         22.0680  -26.7470   -5.3310 C.ar    1  JZ4        -0.0583
      4 C9         22.6710  -25.5120   -5.4480 C.ar    1  JZ4        -0.0199
      5 C10        22.7690  -24.7300   -4.2950 C.ar    1  JZ4         0.1200
      6 C11        21.6930  -26.4590   -2.9540 C.ar    1  JZ4        -0.0551
      7 C12        22.2940  -25.1870   -3.0750 C.ar    1  JZ4        -0.0060
      8 C13        22.4630  -24.4140   -1.8080 C.3     1  JZ4        -0.0245
      9 C14        23.9250  -24.7040   -1.3940 C.3     1  JZ4        -0.0518
     10 OAB        23.4120  -23.5360   -4.3420 O.3     1  JZ4        -0.5065
     11 H          25.3133  -24.3619    0.1509 H       1  JZ4         0.0230
     12 H          23.6591  -24.5327    0.6872 H       1  JZ4         0.0230
     13 H          24.1744  -23.0611   -0.1016 H       1  JZ4         0.0230
     14 H          21.0673  -28.1238   -4.0754 H       1  JZ4         0.0618
     15 H          21.9931  -27.3472   -6.1672 H       1  JZ4         0.0619
     16 H          23.0361  -25.1783   -6.3537 H       1  JZ4         0.0654
     17 H          21.3701  -26.8143   -2.0405 H       1  JZ4         0.0621
     18 H          21.7794  -24.7551   -1.0588 H       1  JZ4         0.0314
     19 H          22.2659  -23.3694   -1.9301 H       1  JZ4         0.0314
     20 H          24.5755  -24.2929   -2.1375 H       1  JZ4         0.0266
     21 H          24.0241  -25.7662   -1.3110 H       1  JZ4         0.0266
     22 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

接下来使用 jz4_fix.mol2（是正确的）。

Generate the JZ4 Topology with CGenFF

现在我们足以用 jz4_fix.mol2 来制作拓扑图了。访问 CGenFF 服务器，登录您的账户，点击页面顶部的“Upload molecule”。上传 jz4_fix.mol2 ，随后 CGenFF 服务器将以CHARMM“流（stream）”文件(扩展名.str)的形式快速返回一个拓扑。将其内容从web浏览器中保存到一个名为“jz4.str”的文件中。您也可以在这里下载该文件的副本（不过还是建议使用者亲自尝试一下转换过程）。

CHARMM流文件包含所有拓扑信息——原子类型、电荷和键连信息。它通过类比推理的方式，涵盖了所有未被力场覆盖的内部相互作用的键合参数。CGenFF还为每个参数提供惩罚分数（penalty scores），即对所分配参数的可靠性进行评估。任何低于10的都被认为是可以立即使用的。10-50意味着需要对拓扑进行一些验证，大于50通常需要手动重新参数化。这种惩罚得分是CGenFF服务器最重要的功能之一。许多其他服务器生成拓扑，然后就成了“黑盒”，用户只能隐忍地信任它们。在没有经过验证的情况下把你的整个研究项目都押在一个自动程序上是非常危险的。一个糟糕的配体拓扑结构会导致大量的时间浪费和不可靠的结果。一定要验证新参数化物质的拓扑！至少，要根据力场中存在的基团，检查配体的电荷大小和原子类型。

检查 jz4. str 的内容，看看电荷的得分和新的二面体参数。它们都是非常低的，表明这个拓扑非常好的，可以直接用于我们的模拟。

CHARMM格式的JZ4拓扑非常好，但是如果我们试图在GROMACS中运行我们的模拟，它就没有用处了。使用我们之前从 MacKerell 网站下载的 cgenff_charmm2gmx.py 脚本。脚本名称内含了版本信息，_py2或_py3，但为了简单起见，这里省略了该信息，但您需要使用Python版本下对应的脚本的实际名称。执行以下转换:

python cgenff_charmm2gmx.py JZ4 jz4_fix.mol2 jz4.str charmm36-mar2019.ff

请注意：这个脚本需要 NetworkX 包，版本1.11。它与NetworkX 2.x不兼容，运行只会退出并显示一个明显的错误。

请注意转换结束时下面的屏幕输出表示成功：

============ 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 ============

这个配体引入了不属于现有的力场的新的键合参数，这些参数将被写入一个名为“jz4. prm”的文件。该文件的格式为GROMACS .itp文件。我们待会就处理这个文件，务必注意它的存在。配体拓扑现在写入“jz4. itp”。它包含了配体[ moleculetype ]的定义。

Build the Complex

从pdb2gmx中，我们有一个名为“3HTB_processed.gro”的文件。它包含了我们蛋白质中经过处理的、符合力场的结构。我们还有从 cgenff_charmm2gmx.py 获得的 jz4_ini. pdb ，它具有所有必要的H原子，并匹配了JZ4拓扑中的原子名称。使用editconf将这个.pdb文件转换为.gro格式:

gmx editconf -f jz4_ini.pdb -o jz4.gro

将 3HTB_processed.gro 复制到一个新的文件中进行操作，可以将其命名为“complex.gro”（因为在蛋白质中加入JZ4将生成我们的蛋白质配体复合体）。接下来，复制 jz4.gro 的坐标部分，并将其粘贴到 complex.gro 中，蛋白质原子的最后一行下面，box vector的上面，如下所示:

从：

  163ASN      C 1691   0.621  -0.740  -0.126
  163ASN     O1 1692   0.624  -0.616  -0.140
  163ASN     O2 1693   0.683  -0.703  -0.011
   5.99500   5.19182   9.66100   0.00000   0.00000  -2.99750   0.00000   0.00000   0.00000

变为（新添加段落为绿色，翻译后没弄上）：

  163ASN      C 1691   0.621  -0.740  -0.126
  163ASN     O1 1692   0.624  -0.616  -0.140
  163ASN     O2 1693   0.683  -0.703  -0.011
    1JZ4     C4    1   2.429  -2.412  -0.007
    1JZ4     C7    2   2.155  -2.721  -0.411
    1JZ4     C8    3   2.207  -2.675  -0.533
    1JZ4     C9    4   2.267  -2.551  -0.545
    1JZ4    C10    5   2.277  -2.473  -0.430
    1JZ4    C11    6   2.169  -2.646  -0.295
    1JZ4    C12    7   2.229  -2.519  -0.308
    1JZ4    C13    8   2.246  -2.441  -0.181
    1JZ4    C14    9   2.392  -2.470  -0.139
    1JZ4    OAB   10   2.341  -2.354  -0.434
    1JZ4     H1   11   2.531  -2.436   0.015
    1JZ4     H2   12   2.366  -2.453   0.069
    1JZ4     H3   13   2.417  -2.306  -0.010
    1JZ4     H4   14   2.107  -2.812  -0.407
    1JZ4     H5   15   2.199  -2.735  -0.617
    1JZ4     H6   16   2.304  -2.518  -0.635
    1JZ4     H7   17   2.137  -2.681  -0.204
    1JZ4     H8   18   2.178  -2.476  -0.106
    1JZ4     H9   19   2.227  -2.337  -0.193
    1JZ4    H10   20   2.458  -2.429  -0.214
    1JZ4    H11   21   2.402  -2.577  -0.131
    1JZ4    H12   22   2.374  -2.321  -0.516
   5.99500   5.19182   9.66100   0.00000   0.00000  -2.99750   0.00000   0.00000   0.00000

因为我们已经在.gro文件中添加了22个原子，所以增加 complex. gro 的第二行来反映这一变化。现在坐标文件中应该有2636个原子。

Build the Topology

在系统拓扑结构中包含JZ4配体的参数是非常容易的。只需在topol.top中插入一行 #include“jz4.itp” 来说明位置约束文件中已完成涵盖。位置约束声明的出现表明了 “Protein” moleculetype 部分的结束。

从：

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"

变成：

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"

配体引入了新的二面体参数，该参数由 cgenff_charmm2gmx.py 脚本写入“jz4.prm”。在topol.top的顶部，插入 #include 语句来添加以下参数:

; Include forcefield parameters
#include "./charmm36-mar2019.ff/forcefield.itp"

; Include ligand parameters
#include "jz4.prm"

[ moleculetype ]
; Name            nrexcl
Protein_chain_A     3

#include语句的位置很关键——它必须出现在任何 [moleculetype] 条目之前，因为所有的参数都必须在构建分子之前被定义。它还必须出现在父力场的 #include 语句之后，因为必须知道所有原子类型，才能引入使用它们的绑定参数。

最后要做的调整是在[ molecules ]。为了说明 complex.gro 中有一个新分子的事实，我们得把配体的名称加在这里，像这样：

[ molecules ]
; Compound        #mols
Protein_chain_A     1
JZ4                 1

现在拓扑和坐标文件与系统的内容一致。

Solvation

在这一点上，操作流程就像任何其他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

在可视化 solve.gro 时，你可能会想知道为什么 editconf 没有生成所要求的十二面体单元格形状，因为系统看起来会像这样:

GROMACS程序总是使用最有效的数字坐标表示，在此它已经将所有东西重新包装成一个三斜的单元格。mdrun 执行的物理计算可以用不同的坐标封装等价地执行，因此它会首选最高效的。在生成 .tpr 文件之后，可以恢复所需的单元格形状。

Adding Ions

我们现在有一个包含带电蛋白质的溶剂化系统。pdb2gmx的输出结果告诉我们，该蛋白质的净电荷为+6e(基于其氨基酸组成)。如果您在pdb2gmx输出中错过了这个信息，请查看topol.top中的[atoms]指令的最后一行;它应该（部分）写作“qtot 6”。由于生物体内不存在净电荷，我们必须向我们的系统中加入离子。

用grompp来生成一个.tpr文件，使用任何.mdp文件。作者习惯使用需要最少的参数的 .mdp，因为最容易维护。比如这个。

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

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

用 -pname 和 -nname 指定的离子名称在以前的GROMACS版本中是特定于力场的，但在4.5版本中进行了标准化。指定的原子名称总是元素符号，所有大写字母，连同[ moleculetype ]也是如此。残基名称可以加上或不加上电荷的符号(+/-)。如果遇到困难，请参考 ions.itp 以获得适当的名称。。

现在你的[ moleculetype ]应该长这样：

[ molecules ]
; Compound        #mols
Protein_chain_A     1
JZ4                 1
SOL             10228
CL                  6

Energy Minimization

现在系统已经组装好了，使用这个输入参数文件，使用grompp创建二进制输入:

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

当运行genbox和genion时，确保你一直在更新 topol.top 文件，否则你会得到很多糟糕的错误消息(“坐标文件中的坐标数不匹配拓扑”等)。

接下来调用mdrun：

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

正如在溶菌酶教程中，可以使用 energy 模块来监测势能的各种组成部分。我不会在这里说明这些原则。玩得开心。

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

Equilibration, Phase 1

平衡我们的蛋白质配体复合体就像平衡任何含有蛋白质的水系统一样。在这种情况下，有一些特殊的考虑：

1. 对配体施加约束
2. 温度偶联基团的处理

Restraining the Ligand

为了约束配体，我们需要为它生成一个位置约束拓扑。首先，为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

现在，我们需要在拓扑中包含这些信息。根据我们希望使用的条件，我们可以通过几种方式来实现这一点。如果我们只是想在蛋白质也被抑制的时候限制配体，在你的拓扑结构中指定的位置添加以下行（; Ligand position restraints 部分）：

; 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-mar2019.ff/tip3p.itp"

位置很重要!您必须按照指示将对 posre_jz4.itp 的调用放到拓扑中。jz4.itp 中的参数为我们的配体定义了一个[ molecular type ]定义。moleculetype 以包含水的拓扑结构(tip3p.itp)结束。如果将对posre_jz4.itp的调用行放到别处，系统可能会尝试将位置约束应用到错误的 moleculetype 上。

如果你想在平衡过程中进行更多的控制，即独立抑制蛋白质和配体，你可以在不同的 #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-mar2019.ff/tip3p.itp"

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

Thermostats

对温度耦合的适当控制是一个需要谨慎对待的问题。将每一种 moleculetype 与它自己的恒温基团（thermostatting group）耦合在一起是个坏主意。例如，如果你这样做：

tc-grps = Protein JZ4 SOL CL

你的系统可能会爆炸，因为温度耦合算法不够稳定，无法控制与少量原子(即JZ4和CL)组合的动能的波动。不要将系统中的每一个物质单独耦合。

通用的方法是（在 .mdp 文件中）设置 tc-grps = Protein Non-Protein，然后继续。不幸的是，“非蛋白质”组也包括JZ4。由于JZ4和蛋白质的物理联系非常紧密，最好将它们视为一个单一的实体。也就是说，为了温度耦合的目的，JZ4要与蛋白质绑定在一起。以同样的方式，我们插入的几个氯离子被认为是溶剂的一部分。要做到这一点，我们需要一个特殊的索引组来合并蛋白质和JZ4。我们使用 make_ndx 模块来完成这个任务，它提供了整个系统的任何坐标文件。在这里，我使用的是 em.gro，我们系统输出(最小化)的结构：

gmx make_ndx -f em.gro -o index.ndx
...
  0 System              : 33506 atoms
  1 Protein             :  2614 atoms
  2 Protein-H           :  1301 atoms
  3 C-alpha             :   163 atoms
  4 Backbone            :   489 atoms
  5 MainChain           :   651 atoms
  6 MainChain+Cb        :   803 atoms
  7 MainChain+H         :   813 atoms
  8 SideChain           :  1801 atoms
  9 SideChain-H         :   650 atoms
 10 Prot-Masses         :  2614 atoms
 11 non-Protein         : 30892 atoms
 12 Other               :    22 atoms
 13 JZ4                 :    22 atoms
 14 CL                  :     6 atoms
 15 Water               : 30864 atoms
 16 SOL                 : 30864 atoms
 17 non-Water           :  2642 atoms
 18 Ion                 :     6 atoms
 19 JZ4                 :    22 atoms
 20 CL                  :     6 atoms
 21 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.tpr

gmx mdrun -deffnm nvt

Equilibration, Phase 2

一旦NVT模拟完成，使用此.mdp文件继续NPT:

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

gmx mdrun -deffnm npt

Production MD

在两个平衡阶段完成后，系统现在就在期望的温度和压力下很好地平衡了。我们现在已经准备好释放位置限制，并运行生产MD进行数据收集。这个过程就像我们之前看到的一样，我们将利用检查点文件（在本例中，它现在包含保存压力耦合信息）来进行 grompp。我们将运行一个10-ns MD模拟，脚本可以在这里找到。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr

gmx mdrun -deffnm md_0_10

Analysis

Recentering and Rewrapping Coordinates

在任何用周期性边界条件进行的模拟中，分子可能出现“断裂”或可能在盒子中来回“跳跃”。为了重新定位蛋白质并重新包裹单元细胞内的分子，以恢复所需的菱形十二面体形状，调用 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 = 0ns)，对重定位的轨迹使用trjconv -dump：

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0

为了更平滑的可视化，进行适当的旋转和平移拟合可能是有益的（请注意，同步PBC重新包装和拟合坐标在数学上是不兼容的。如果您希望执行拟合（这对可视化很有用，但对大多数分析例程不是必需的），请按此处所示，分别执行这些坐标操作）。按如下方式执行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”进行输出。

Analyzing Protein-Ligand Interactions and Ligand Dynamics

本教程不可能涵盖您希望执行的所有分析方法。这里将说明一些基本的操作。

2-丙基酚配体能与Gln102侧链形成氢键。GROMACS hbond 模块可以很容易地用来计算任何原子团之间的氢键数，但在我们的例子中，唯一的值将是1或0。为了更详细的了解配体与Gln102的相互作用，我们将计算JZ4的羟基与Gln102的羰基O原子之间的距离。对于氢键的形成，典型的标准是给体和受体原子之间的距离≤3.5 Å (0.35 nm)。使用 distance 模块计算轨迹过程中的距离，指令如下（有关示例和更多语法，请参阅 gmx help）：

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 模块中，这个角度被定义为 hydrogen-donor-acceptor ，这个角度应该≤30°。要进行这一分析，首先为给体原子(必须包括给体重原子和键合的氢)和受体原子创建索引组：

gmx make_ndx -f em.gro -n 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

请注意，angle 不接受 -s 参数，这与大多数其他GROMACS分析模块不同。角度计算不需要质量或周期信息、原子名称等，所以轨迹处理时不需要.tpr或坐标文件。命令返回的值应该大约为23±9°。这个结果可能有些出乎意料，因为我们构建的指标组顺序为OAB, H12, OE1，这对应于供体-氢-受体的距离，我们预计其接近线性(~180°)。检查索引文件的内容，你会发现:

[ JZ4_&_OAB_H12_Protein_&_r_102_&_OE1 ]
1616 2624 2636

make_ndx 自动将原子数从低到高排序，因此计算的结果是受体-供体-氢的角度，与 hbond 模块计算的角度相同。这与氢键的形成是一致的，因为氢键的角度≤30°。为了得到理想的给体-氢-受体的角度，我们必须在文本文件中手动编辑索引组来重新排序原子号(2624 2636 1616)。用该指标组重新进行角度计算，得到的平均值为147±11°。

最后，我们可能对量化在模拟过程中配体结合姿势的变化感兴趣。为了计算一个仅为JZ4的重原子RMSD，为它创建一个新的索引组：

gmx make_ndx -f em.gro -n index.ndx
...
 > 13 & ! a H*
 > name 26 JZ4_Heavy
 > q

执行 rms 模块，最小二乘拟合选择“Backbone”，计算RMSD选择“JZ4_Heavy”。通过拟合去除蛋白质的整体旋转和平移，RMSD报告出JZ4位置相对于蛋白质变化的程度，这是一个说明在模拟过程中结合姿势保持得如何的很好的指标。

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 Å)，表明配体的位置变化很小。

Protein-Ligand Interaction Energy

为了量化JZ4和T4溶菌酶相互作用的强度，计算这两种物质之间的非键相互作用能可能是有用的。GROMACS具有分解任意数量的定义基团之间的短期非键能的能力。值得注意的是，这个量不是自由能或结合能。事实上，大多数力场都不是这样参数化的，这个量实际上是有物理意义的。CHARMM专门针对与水的量子力学相互作用能进行参数化，因此它在本质上与那个有意义的量相平衡，这种相互作用能是有用的。

相互作用能是通过.mdp文件中的 energygrps 关键字计算的。尽管它是一个 .mdp 关键字，交互作用能计算不应该被认为是正常模拟的一部分。短距离能量的分解在GPU上运行是不兼容的，也会不必要地降低计算速度。mdrun 模块不需要做这些额外的工作来执行有效的模拟。因此，只计算相互作用能作为分析的一部分，而不是动力学。从定义了 energygrps = Protein JZ4 的.mdp文件中创建一个新的.tpr文件，如下所示:

gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr

接下来，使用 -rerun 选项调用 mdrun ，从现有的模拟轨迹中重新计算能量：

gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu

注意使用 -deffnm 来读取 ie.pr 并将所有输出文件写入 ie.* 作为他们的文件名。-rerun 选项取要重新计算能量的轨迹的名称，-nb cpu告诉 mdrun 只尝试在cpu硬件上运行，忽略任何可能可用的GPU。如上所述，这种类型的计算不能在GPU上执行。重新运行应该非常快，只需几分钟即可完成。

通过 energy 模块提取感兴趣的能量项。我们感兴趣的术语有 Coul-SR:Protein-JZ4 和 LJ-SR:Protein-JZ4。

gmx energy -f ie.edr -o interaction_energy.xvg

近程库仑相互作用能平均为-20.5±7.4 kJ mol^-1，近程Lennard-Jones能平均为-99.1±7.2 kJ mol^-1。从这些量的相对大小中得出结论可能是很诱人的，但是即使CHARMM是根据明确的水相互作用能进行参数化的，将相互作用能进一步分解成这些分量并不一定是真实的。没有办法通过实验来验证这些量，所以不可能知道它们是否有意义。然而，总相互作用能在这种情况下是有用的。这个值(根据两个量相加的标准公式传播误差后)是-119.6±10.3 kJ mol^-1。

Summary

你现在已经进行了蛋白质配体复合体与GROMACS的分子动力学模拟。本教程不应被视为全面的。我们可以用GROMACS进行更多类型的模拟(自由能计算、非平衡MD和正态模态分析，仅举几个例子)。为了提高效率和准确性，您还应该查阅文献和GROMACS手册，以便对这里提供的.mdp文件进行调整。

略。

Happy simulating!