上一节我们对蛋白配体复合物体系进行一个能量最小化与动力学模拟的操作，这一节我们将对其结果进行分析。

重新定位和重新包装坐标

与任何使用周期性边界条件进行的模拟一样，分子可能会出现 "断裂 "或在单元格中来回 "跳跃"。要重新调整蛋白质的位置并将分子重新包覆在单元格内以恢复所需的菱形十二面体形状，请调用 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 "为输出。需要注意的是，如果模拟时间较长，涉及许多跨周期边界的跃迁，则可能难以对复合物（protein-ligand, protein-protein）进行居中。在这种情况下（尤其是在protein-protein复合物中），可能需要创建一个自定义索引组来进行定 心，该索引组与复合物中一个蛋白质的活性位点或一个单体的界面残基相对应。

要提取轨迹的第一帧（t = 0 毫微秒），请使用 trjconv -dump 和重新定心的轨迹：

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

为了获得更流畅的视觉效果，可能需要进行旋转和平移拟合。执行 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（系统）"进行输出。请注意，同时进行 PBC 重包和坐标拟合在数学上是不相容的。如果您希望进行拟合（这对可视化很有用，但对大多数分析程序来说并非必要），请按照此处的说明分别进行这些坐标操作。

分析蛋白质配体相互作用和配体动力学

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

2-propylphenol 配体可以与 Gln102 侧链形成氢键。我们将计算 JZ4 的羟基与 Gln102 的羰基 O 原子间的距离。形成氢键的典型标准是供体原子和受体原子之间的距离≤ 3.5 Å（0.35 nm）。使用命令行选择语法（示例和更多语法请参见 gmx 帮助选择），使用距离模块计算轨迹上的距离。

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 模块中，角度定义为氢-供体-受体，该角度应≤ 30°。要进行此分析，首先要为供体原子（必须包括供体重 原子和键合氢原子）和受体原子创建索引组：

gmx make_ndx -f em.gro -o index.ndx

使用角度模块分析这三个原子形成的角度：

gmx angle -f md_0_10_center.xtc -n index.ndx -ov angle.xvg

请注意，与其他大多数 GROMACS 分析模块不同，角度计算不使用 -s 参数。角度计算不需要质量或周期信息、原子名称等，因此在处理轨迹时不需要 .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

执行 rms 模块，选择 "Backbone "进行最小二乘拟合，选择 "JZ4_Heavy "进行 RMSD 计算。这样做可以通过拟合去除蛋白质的整体旋转和平移，所报告的 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 Å）左右，这表明配体的位置变化非常小。

蛋白质配体相互作用能量

为了量化 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.tpr，并将所有输出文件写入 ie.* 作为文件名。-rerun 选项使用了要重新计算能量的轨迹名称，而 -nb cpu 则告诉 mdrun 只尝试在 CPU 硬件上运行，而忽略任何可能可用的 GPU。如上所述，这种计算无法在 GPU 上执行。重新运行的速度应该很快，几分钟就能完成。

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

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

平均短程库仑相互作用能为 -20.5 ± 7.4 kJ mol-1，短程伦纳德-琼斯能为 -99.1 ± 7.2 kJ mol-1。从这些量的相对大小得出结论可能很有诱惑力，但尽管 CHARMM 是根据明确的水相互作用能进行参数化的，将相互作用能进一步分解为这些分量并不一定是真实的。没有办法通过实验验证这些量，因此无法知道它们是否有意义。不过，在这种情况下，总相互作用能是有用的。该值（根据两个量相加的标准公式传播误差后）为 -119.6 ± 10.3 kJ mol-1。

至此，分析部分结束。