【2.5.1.5】并行计算蛋白蛋白结合自由能
一、建立起始结构并运行模拟以获得平衡的系统。
我们将在此模拟中建模的系统是人类H-Ras蛋白与C-Raf1的Ras结合域（Ras-Raf）之间的复合体，它是信号转导级联的核心。 这是RAS-RAF复合系统的部分平衡的，预先准备的pdb文件。
该结构包含ras和raf蛋白，以及生理必需的GTP核苷酸，如下图所示：
有关如何构建起始结构和运行仿真以获得平衡系统的指导，请参阅第1节和第2节。
使用MMPBSA.py计算绑定自由能的重要文件是拓扑文件和mdcrd文件（ ras-raf_top_mdcrd.tgz ）
二、并行计算Ras-Raf的结合自由能。
现在，我们将计算复合物，受体和配体的相互作用能和溶剂化自由能，并对结果求平均值，以获得结合自由能的估计值。请注意，在本教程的这一部分中，我们将不会计算熵对绑定的贡献，因此严格来说，我们的结果将不是真正的自由能，但可用于与类似系统进行比较。有关使用正态分析（Nmode）来计算系统的熵贡献或取消注释常规部分中最后一行以使用AMBER中的ptraj模块执行准谐波熵计算的示例，请参见第5节。
我们将同时使用MM-GBSA和MM-PBSA方法进行结合能计算，以与从3.1节 获得的结果进行比较。 MMPBSA.py.MPI通过为每个线程（进程）分配相等数量的帧来并行化计算。因此，当处理的帧数是启动的线程数的倍数时，它最有效地运行。但是，这不是必需的。如果帧数不是线程数的倍数，则“剩余”帧将在启动的线程数的子集中平均分配。例如，在3个线程上运行50帧将导致2个线程计算17个帧，而最后一个线程仅计算16个帧。因此，第三个线程将必须等待前两个线程完成计算才能继续。因此，例如5个线程将是一个明智的选择（因为每个线程占用10帧）。但是，线程数不能超过已处理的帧数，否则MMPBSA.py.MPI将终止并显示一条错误消息。 MMPBSA.py.MPI的输入文件与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.MPI的输入文件被设计为类似于AMBER的sander模块中使用的mdin文件的设置。 每个名称列表的开头均以＆符号表示，后跟名称列表的名称。 此外，反斜杠（/）或'＆end’可用于结束名称列表。 有关所有变量的完整列表，请参见amber手册。 该输入文件分为三个名称列表：general，pb和gb。 通用名称列表旨在指定并非特定于计算的特定部分而是特定于所有部分的变量。 在此设置中，我们已将RAS定义为受体，将RAF定义为配体。 ‘endframe’变量设置mdcrd的停止帧。 “＆gb”和“＆pb”名称列表标记使脚本知道使用这些名称列表中定义的给定值执行MM-GBSA和MM-PBSA计算。
mpirun -np 4 $AMBERHOME/bin/MMPBSA.py.MPI -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 > progress.log 2>&1
或通过将脚本提交到排队系统（例如PBS）
qsub parallel.job
parallel.job
#!/bin/sh
#PBS -N rasraf_parallel
#PBS -o parallel.out
#PBS -e parallel.err
#PBS -m abe
#PBS -M email@address.com
#PBS -q brute
#PBS -l nodes=1:node:ppn=4
#PBS -l pmem=900mb
cd $PBS_O_WORKDIR
mpirun -np 4 MMPBSA.py.MPI -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 bigprod.mdcrd > progress.log 2>&1
这会将计算进度打印到文件 progress.log 中。 计算过程中的所有错误也将被打印到该文件中（这是2>＆1的目的）。 最后，一旦计算完成，将打印计时，显示脚本每个步骤所花费的时间。
progress.log
MMPBSA.py.MPI being run on 4 processors
ptraj found! Using /scr/arwen_3/swails/i686/amber11/exe/ptraj
sander found! Using /scr/arwen_3/swails/i686/amber11/exe/sander
Preparing trajectories with ptraj...
50 frames were read in and processed by ptraj for use in calculation.
Starting calculations
Starting gb calculation...
Starting pb calculation...
Calculations complete. Writing output file(s)...
Timing:
Processing Trajectories With Ptraj: 0.126 min.
Total GB Calculation Time (sander): 4.782 min.
Total PB Calculation Time (sander): 28.407 min.
Output File Writing Time: 0.053 min.
Total Time Taken: 33.379 min.
MMPBSA Finished. Thank you for using. Please send any bugs/suggestions/comments to mmpbsa.amber@gmail.com
命令行上的’-y * .mdcrd’告诉脚本读取工作目录中所有以'.mdcrd’结尾的文件，并将它们用作要分析的轨迹。
将’keep_files’设置为其默认值1，以下是所有输出文件：Parallel_output.tgz。
该脚本使用ptraj创建三个未溶剂化的mdcrd文件（复合物，受体和配体），这是在GB和PB计算期间分析的坐标。 * .mdout文件包含所有指定帧的能量。 如果使用ptraj执行准谐波熵计算，则会将平均pdb文件创建为平均结构。 除最终输出文件FINAL_RESULTS_MMPBSA.dat以外，由MMPBSA.py.MPI创建的所有文件都应以前缀'_MMPBSA_开头。
FINAL_RESULTS_MMPBSA.dat
| Run on Sun Feb 14 19:10:43 EST 2010
|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB in serial
|&general
| endframe=50, verbose=1,
| mpi_cmd='mpirun -np 3', nproc=3
|/
|&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): bigprod.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)
统计文件的开头包括日期/时间，输入文件的副本，配体和受体掩膜的最佳猜测，脚本使用的文件，分析的帧数(frames)以及哪个PB solver（如果有）被使用了。统计信息文件的其余部分包括GB的所有平均能量，标准偏差和均值的标准误差，其后是PB。在每个部分之后，结合的ΔG随误差值一起给出。该文件中不同术语的含义如下：
VDWAALS = MM的van der Waals贡献。
EEL =由MM力场计算出的静电能。
EPB / EGB =分别由PB或GB计算得出的对溶剂化自由能的静电贡献。
ECAVITY =通过经验模型计算得出的对溶剂化自由能的非极性贡献。
DELTA G binding=根据上述条件计算得出的最终估计结合自由能。 （千卡/摩尔）
请注意，尚未报告总气相能，因为使用单轨方法，受体和配体的键合项的值应恰好抵消了络合物的键合项的值。如果能量没有抵消到允许的公差（0.001 kcal / mol）之内，则会出现错误消息。
人们通常期望找到一种非常有利的静电能和一种不利的溶剂化自由能。这象征着人们必须使用的能量来使结合颗粒去溶剂化并对齐其结合界面。
从负的总结合自由能为-86.36 kcal / mol，我们可以清楚地看到，这是纯水中的一种有利的蛋白质-蛋白质复合物，但请记住，由于我们没有估算（不利于）熵对结合的贡献。请注意，GB方法给出的结合能略低，但仍表明这是一种有利的结合态。
参考资料
