g_mmpbsa：Gromacs 插件，计算MM/PBSA

g_mmpbsa适用于 Gromacs 4.5.x-5.1.x 版本计算 MM/PBSA。原文链接及下载地址：http://rashmikumari.github.io/g_mmpbsa/。“😇”内的内容为译者添加。

如果你在模拟中使用了该脚本，请务必在发表时引用：

Kumari et al (2014) g_mmpbsa - A GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model. 54:1951-1962.
.
Baker et al. (2001) Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. USA 98:10037-10041.

😇有关于MMPBSA的基础原理可以参考这篇文章😇

Installation

在此处下载对应版本的预编译程序包。随后解压缩，并写入环境变量：

tar -zxvf g_mmpbsa.tar.gz

vim .bashrc
export PATH=$PATH:/path_to_your_file/g_mmpbsa/bin

随后就可以使用了。😇关于 g_mmpbsa 的具体使用方法以及需要的参数文件，建议跟着下面这个教程走一遍，你就可以理解了（耗时不会很长）😇

Tutorial 1: Binding Energy of Single Protein-Inhibitor Complex

点击此处下载教程所需文件。解压后进入文件夹：

tar -zxvf tutorial.tar.gz
cd tutorial/1EBZ

该文件夹包含HIV-1蛋白酶抑制剂复合物的拓扑参数（tpr）、索引（ndx）和轨迹（xtc）文件。结合能由三个能量项组成：真空中的势能（MM项，包括所选两组的 vdw and electrostatic）、极性溶剂化能（GB项）、非极性溶剂化能（SA项）。这些能量项可以用三步或一步来计算。

Three steps calculation

1. 计算真空中的势能

运行下述指令，并选择“1”和“13” ：

g_mmpbsa -f 1EBZ.xtc -s 1EBZ.tpr -n 1EBZ.ndx -pdie 2 -decomp

将会生成两个文件 energy_MM.xvg 和 contrib_MM.dat。这两个文件都可以由-mm filename1.xvg和-mmcon filename2.dat自主设定名称。energy_MM.xvg 文件包含范德华、静电相互作用以及蛋白质和抑制剂之间的净非键合势能。contrib_MM.dat 包含每个残基对计算的净非键相互作用能的贡献。

2. 计算极性溶剂化能

要计算极性溶剂化能，需要一个输入文件（例如 tutorial/polar.mdp）。该文件包含用于计算极性溶剂化能的输入参数。要获取输入参数的详细信息，请单击此链接（在本文后面也翻译出来了）。执行以下命令，选择1和13：

g_mmpbsa -f 1EBZ.xtc -s 1EBZ.tpr -n 1EBZ.ndx -i ../polar.mdp -nomme -pbsa -decomp

将会生成两个文件 polar.xvg 和 contrib_pol.dat。-pol filename1.xvg和-pcon filename2.dat可以用不同的名称生成这两个文件。polar.xvg 包含未结合蛋白质、未结合抑制剂和蛋白质抑制剂复合物的极性溶剂化能。contrib_pol.dat 包含每个残基对计算的净极性溶剂化能的贡献。

3. 计算非极性溶剂化能

为了计算非极性溶剂化能，需要一个输入文件（例如tutorial/apolar_sasa.mdp）。该文件包含用于计算非极性溶剂化能的参数。要获取输入参数的详细信息，请单击此链接。该文章中讨论了几种类型的非极性模型（数值见表1）。这里为只使用 SASA 和只使用 SAV 的模型提供了输入参数文件。

SASA-only model：

运行并选择“1”和“13” ：

g_mmpbsa -f 1EBZ.xtc -s 1EBZ.tpr -n 1EBZ.ndx -i ../apolar_sasa.mdp -nomme -pbsa -decomp -apol sasa.xvg -apcon sasa_contrib.dat

将会生成两个文件 sasa.xvg 和 sasa_contrib.dat。sasa.xvg 含有未结合蛋白、未结合抑制剂和蛋白抑制剂复合物的非极性溶剂化能。sasa-contrib.dat 包含每个残基对计算的净极性溶剂化能的贡献。

SAV-only model：

运行并选择“1”和“13” ：

g_mmpbsa -f 1EBZ.xtc -s 1EBZ.tpr -n 1EBZ.ndx -i ../apolar_sav.mdp -nomme -pbsa -decomp -apol sav.xvg -apcon sav_contrib.dat

输出与 SASA 类似。

One step calculation

运行并选择“1”和“13” ：

g_mmpbsa -f 1EBZ.xtc -s 1EBZ.tpr -n 1EBZ.ndx -i ../pbsa.mdp -pdie 2 -pbsa -decomp

pbsa.mdp 包含极性和仅SASA的非极性溶剂化能的输入参数。通过使用上述单个命令计算所有三个能量项，并生成所有输出文件。

Average Binding Energy Calculation

为了计算平均结合能，g_mmpbsa 包中提供了 python 脚本MmPbSaStat.py（点击此链接进入下载）。有关此脚本的详细信息，请单击此链接。要执行此脚本，需要以上获得的文件作为输入。

运行（所需三个文件均可由上一步获得）：

python MmPbSaStat.py -m energy_MM.xvg -p polar.xvg -a apolar.xvg (or sasa.xvg)

你将会获得两个输出文件 full_energy.dat 和 summary_energy_dat。summary_energy.dat 包含所有能量项的平均值和标准差，同样包含结合能：

#Complex Number:    1 
=============== 
   SUMMARY   
=============== 

 van der Waal energy      =        -334.587   +/-   15.514 kJ/mol 

 Electrostattic energy    =        -159.380   +/-   15.810 kJ/mol 

 Polar solvation energy   =         313.698   +/-   10.174 kJ/mol 

 SASA energy              =         -30.431   +/-    0.996 kJ/mol 

 SAV energy               =           0.000   +/-    0.000 kJ/mol 

 WCA energy               =           0.000   +/-    0.000 kJ/mol 

 Binding energy           =        -210.699   +/-   19.745 kJ/mol 

=============== 
    END     
=============== 

full_energy.dat 包含作为时间函数的能量项的值。最后四列包含作为时间函数的Δ_E_MM、Δ_G_polar、Δ_G_nonpolar 和 Δ_G_binding。这些数量可以用 xmgrace/matplotlib/gnuplot 绘制（你也可以直接复制到 GraphPad 或者用 Excel 处理）。xmgrace格式（.agr）的四个文件在tutorial/1EBZ/output中提供。

要使用bootstrap分析计算平均结合能，请执行以下命令：

python MmPbSaStat.py -bs -nbs 2000 -m energy_MM.xvg -p polar.xvg -a apolar.xvg (or sasa.xvg)

同样，生成两个输出文件 full_energy.dat 和 summary_energy_dat。full_energy.dat 与上面的类似。然而，summary_energy.dat 包含包括结合能在内的所有含能分量的平均误差和标准误差。summary_energy.dat 中的平均值与上述值略有不同。有关此方法的更多详细信息，请参阅此文献。

Tutorial 2: Energy Contribution of Residues to Binding

😇暂时咕咕😇

Usage

😇与一般 gromacs 的指令选项相似（该部分较为简单，不再翻译）😇：

Option     Filename  Type         Description
------------------------------------------------------------
  -f       traj.xtc  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -s      topol.tpr  Input        Run input file: tpr tpb tpa
  -i     grompp.mdp  Input, Opt.  grompp input file with MD parameters
  -n      index.ndx  Input, Opt.  Index file
 -mm  energy_MM.xvg  Output, Opt. xvgr/xmgr file
-pol      polar.xvg  Output, Opt. xvgr/xmgr file
-apol    apolar.xvg  Output, Opt. xvgr/xmgr file
-mmcon contrib_MM.dat  Output, Opt. Generic data file
-pcon contrib_pol.dat  Output, Opt. Generic data file
-apcon contrib_apol.dat  Output, Opt. Generic data file


Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   yes     Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-tu          enum   ps      Time unit: fs, ps, ns, us, ms or s
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   xmgrace  xvg plot formatting: xmgrace, xmgr or none
-[no]silent  bool   no      Display messages, output and errors from external APBS program. Only works with external APBS program
-rad         enum   bondi   van der Waal radius type: bondi, mbondi, mbondi2 or amber
-rvdw        real   1       Default van der Waal radius (in nm) if not found
-[no]mme     bool   yes     To calculate vacuum molecular mechanics energy
-pdie        real   1       Dielectric constant of solute. Should be same as of polar solvation
-[no]incl_14 bool   no      Include 1-4 atom-pairs, exclude 1-2 and 1-3 atom pairs during MM calculation. Should be "yes" when groups are bonded with each other.
-[no]focus   bool   no      To enable focusing on the specfic region of molecule, group of atoms must be provided in index file
-[no]pbsa    bool   no      To calculate polar and/or non-polar solvation energy
-ndots       int    24      Number of dots per sphere in the calculation of SASA, more dots means more accuracy
-[no]diff    bool   yes     Calculate the energy difference between two group otherwise only calculates for one group
-[no]decomp  bool   no      Decomposition of energy for each residue

😇译者个人关于MMPBSA的一些感想😇

😇MMPBSA是一种端点（End point）式的结合自由能计算方法，在译者看来，传统的MMPBSA有一些局限。它并不能全面描述结合过程中动态的能量变化，以至于一昧的对所有轨迹进行运算后得到的平均结合能很可能是不可靠的。

通常，MMPBSA更适合于计算静态结合位点的结合自由能，这需要在模拟过程中底物稳定的结合在结合位点上。但是，对于大部分蛋白及其底物（或者说配体），结合位点未必是单一的，我们可以设想这样一种情况：有一个底物L，它与蛋白P的结合发生在P活性口袋的底部，但实际上L也可以结合在活性口袋侧壁上的一个亚结合位点上，尽管后者结合不如前者稳定；而在模拟过程中，很有可能发生L在结合位点与亚结合位点之间来回交换的现象（这里插一嘴，理论上可以统计结合在两个不同位点的帧所占总模拟时间的比例来比较结合的强弱），如果此时我们不对轨迹加以聚类，然后分别计算MMPBSA，那么所得到的值会是两位点结合能的加权均值；这个平均值不能代表任何一个结合位点的结合自由能，因此一昧的对所有轨迹进行运算后得到的平均结合能很可能是不可靠的。

此外，考虑到MM/PBSA方法使用单一的静态结构来表示蛋白质-配体复合物，它不能完全捕捉结合过程固有的动态变化。或者换句话说，就是我们无法通过MMPBSA知道底物是如何结合到蛋白相应的结合位点上的，这一点也倒是显而易见，毕竟MMPBSA计算的轨迹大都来自于常规MD，采样空间有限，结合过程中是否有能垒存在也无法得知。

从这两点来看，一些增强采样的方法如SMD、TMD、伞形采样或者MtMD，就能更好地反应结合过程的能量变化，由此得到的自由能形貌图或者PMF曲线能提供延结合途径和反应坐标的更多信息。

但是这并不意味着MMPBSA鸡肋，它简单快速，适合用在大规模筛选中；对于结合口袋较浅或者结合非常稳定的情况一般是可靠的；而且一般稍加处理，对轨迹进行个聚类，也可以消除结合不稳定带来的影响。译者的一些拙见，如有不对还请指出😇