【2.5.1.2】计算蛋白-蛋白的结合自由能
使用MMPBSA.py计算绑定自由能的重要文件是拓扑文件和mdcrd文件（ ras-raf_top_mdcrd.tgz ）
现在，我们将计算复合物，受体和配体的相互作用能和溶剂化自由能，并对结果求平均值，以获得结合自由能的估计值。请注意，在本教程的这一部分中，我们将不会计算熵对绑定的贡献，因此严格来说，我们的结果将不是真正的自由能，但可用于与类似系统进行比较。有关使用正态分析（Nmode）来计算系统的熵贡献或取消注释以下输入文件的＆general名称列表中的最后一行以使用的ptraj模块执行准谐波熵计算的示例，请参见第3.5节。
为了进行比较，我们将同时使用MM-GBSA方法和MM-PBSA方法进行结合能计算。这是通过以下MMPBSA.py输入文件完成的：
mmpbsa.in
Input file for running PB and GB
&general
endframe=50, verbose=1,
# entropy=1,
/
&gb
igb=2, saltcon=0.100
/
&pb
istrng=0.100,
/
MMPBSA.py的输入文件被设计为类似于AMBER的sander模块中使用的mdin文件的设置。每个名称列表的开头均以＆符号表示，后跟名称列表的名称。此外，反斜杠（/）或'＆end’可用于结束名称列表。有关所有变量的完整列表，请参见Amber手册或键入 $ $AMBERHOME/bin/MMPBSA.py --input-file-help $。该输入文件分为三个名称列表：general，pb和gb。通用名称列表旨在指定并非特定于计算的特定部分而是特定于所有部分的变量。在此设置中，我们已将RAS定义为受体，将RAF定义为配体。 ‘endframe’变量设置mdcrd的停止帧。 “＆gb”和“＆pb”名称列表标记使脚本知道使用这些名称列表中定义的给定值执行MM-GBSA和MM-PBSA计算。 “ verbose”变量允许用户指定将多少输出写入到输出文件。
$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd
这将以交互方式运行脚本，并将计算进度打印到STDOUT，并将任何错误或警告打印到STDERR。最后，一旦计算完成，将打印计时，显示计算的每个步骤所花费的时间。
命令行参数可以使用shell识别的通配符（即*和？用于bash）。例如，命令行上的’-y * .mdcrd’告诉脚本读取工作目录中所有以'.mdcrd’结尾的文件，并将它们用作要分析的轨迹。
这是此脚本创建的所有输出文件：pb_gb_output1.tgz。
该脚本使用ptraj创建三个未溶剂化的mdcrd文件（复合物，受体和配体），这是在GB和PB计算期间分析的坐标。 * .mdout文件包含所有指定帧的能量。创建平均结构的PDB文件（通过RMS）对齐所有快照，以准备使用ptraj进行准谐波熵(quasi-harmonic entropy)计算（如果需要）。除最终输出文件FINAL_RESULTS_MMPBSA.dat外，由MMPBSA.py创建的所有文件均应以前缀“ MMPBSA”开头。
FINAL_RESULTS_MMPBSA.dat
| Run on Thu Feb 11 12:18:37 EST 2010
|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB
|&general
| endframe=50, verbose=1,
|# entropy=1,
|/
|&gb
| igb=2, saltcon=0.100
|/
|&pb
| istrng=0.100,
|/
|--------------------------------------------------------------
|Solvated complex topology file: ras-raf_solvated.prmtop
|Complex topology file: ras-raf.prmtop
|Receptor topology file: ras.prmtop
|Ligand topology file: raf.prmtop
|Initial mdcrd(s): prod.mdcrd
|
|Best guess for receptor mask: ":1-166"
|Best guess for ligand mask: ":167-242"
|Calculations performed using 50 frames.
|Poisson Boltzmann calculations performed using internal PBSA solver in sander.
|
|All units are reported in kcal/mole.
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
GENERALIZED BORN:
Complex:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1863.7944 17.1704 2.4283
EEL -17200.7297 75.9366 10.7391
EGB -3249.6511 65.2075 9.2217
ESURF 91.3565 1.3938 0.1971
G gas -19064.5240 77.8536 11.0102
G solv -3158.2946 65.2224 9.2238
TOTAL -22222.8186 51.0216 7.2155
Receptor:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1268.1888 14.2342 2.0130
EEL -11557.0773 71.7127 10.1417
EGB -2532.0669 57.7003 8.1600
ESURF 64.2843 1.1143 0.1576
G gas -12825.2661 73.1118 10.3396
G solv -2467.7826 57.7110 8.1616
TOTAL -15293.0487 35.3527 4.9996
Ligand:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -529.3090 9.4198 1.3322
EEL -4684.4720 36.1449 5.1117
EGB -1688.9631 26.5353 3.7527
ESURF 37.0493 0.6185 0.0875
G gas -5213.7811 37.3522 5.2824
G solv -1651.9138 26.5425 3.7537
TOTAL -6865.6949 25.8878 3.6611
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -66.2966 4.2751 0.6046
EEL -959.1803 34.9190 4.9383
EGB 971.3789 33.0497 4.6739
ESURF -9.9770 0.3759 0.0532
DELTA G gas -1025.4769 35.1797 4.9752
DELTA G solv 961.4018 33.0518 4.6742
DELTA G binding = -64.0750 +/- 6.3729 0.9013
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
POISSON BOLTZMANN:
Complex:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1863.7944 17.1704 2.4283
EEL -17200.7297 75.9366 10.7391
EPB -3207.7160 66.4023 9.3907
ECAVITY 67.8762 0.7818 0.1106
G gas -19064.5240 6061.1875 857.1813
G solv -3139.8399 66.4069 9.3914
TOTAL -7686.8660 52.5400 7.4303
Receptor:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1268.1888 14.2342 2.0130
EEL -11557.0773 71.7127 10.1417
EPB -2483.7242 56.4551 7.9840
ECAVITY 47.1495 0.4737 0.0670
G gas -12825.2661 5345.3320 755.9441
G solv -2436.5747 56.4571 7.9842
TOTAL -5250.2060 38.5188 5.4474
Ligand:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -529.3090 9.4198 1.3322
EEL -4684.4720 36.1449 5.1117
EPB -1670.4169 27.6694 3.9131
ECAVITY 28.0328 0.4133 0.0584
G gas -5213.7811 1395.1865 197.3092
G solv -1642.3841 27.6725 3.9135
TOTAL -2350.3020 25.1197 3.5525
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -66.2966 4.2751 0.6046
EEL -959.1803 34.9190 4.9383
EPB 946.4251 34.5128 4.8808
ECAVITY -7.3062 0.3004 0.0425
DELTA G gas -1025.4769 1237.6138 175.0250
DELTA G solv 939.1189 34.5141 4.8810
DELTA G binding = -86.3579 +/- 8.3264 1.1775
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
WARNINGS:
igb=2 should be used with mbondi2 pbradii set. Yours are modified Bondi radii (mbondi)
统计文件的开头包括日期/时间，基于给定值和文件的任何警告，mmpbsa.in文本，脚本使用的文件，分析的帧数以及哪个PB solver（如果有）用过的。统计信息文件的其余部分包括GB的所有平均能量，标准偏差和均值的标准误差，其后是PB。在每个部分之后，结合的ΔG随误差值一起给出。该文件中不同术语的含义如下：
VDWAALS = MM的van der Waals贡献。
EEL =由MM力场计算出的静电能。
EPB / EGB =分别由PB或GB计算得出的对溶剂化自由能的静电贡献。
ECAVITY =通过经验模型计算得出的对溶剂化自由能的非极性贡献。
DELTA G binding =根据上述条件计算得出的最终估计结合自由能。 （千卡/摩尔）
请注意，尚未报告总气相能，因为使用单轨方法，受体和配体的键合势项的值应恰好抵消了配合物的键合势项的值。如果能量未在允许的公差范围内消除，则会出现错误消息。
人们通常期望找到一种非常有利的静电能和一种不利的溶剂化自由能。这象征着人们必须使用的能量来使结合颗粒去溶剂化并对齐其结合界面。
从负的总结合自由能为-86.36 kcal / mol，我们可以清楚地看到，这是纯水中的一种有利的蛋白质-蛋白质复合物，但请记住，由于我们没有估算（不利于）熵对结合的贡献。请注意，GB方法给出的结合能略低，但仍表明这是一种有利的结合态。
参考资料
