【2.5.1.1】MM-PBSA

December 06, 2019 molecular_dynamics 阅读量：次

在本教程中，我们将使用MM-PBSA方法来计算两个蛋白质缔合的结合自由能。

MM-PBSA方法及其补充MM-GBSA方法的总体目标是计算两个状态之间的自由能差，这两个状态最常表示两个溶剂化分子的结合和未结合状态，或者比较两个不同溶剂化分子的自由能同一分子的构象。

理想情况下，我们希望直接计算该结合自由能，如下图所示：

但是，在这些溶剂化状态的模拟中，大部分能量贡献将来自溶剂与溶剂的相互作用，并且总能量的波动将比结合能大一个数量级。因此，该计算将花费非常多的时间来收敛。因此，一种更有效的方法是根据以下热力学循环对计算进行划分：

从该图显然可以看出，结合自由能delta-Gbind,solv 可以通过以下公式计算：

在MM-PBSA方法中，对上述结合自由能的不同贡献是以多种方式计算的：

通过为三个状态的每一个求解线性化的Poisson Boltzman方程或广义的Born方程来计算溶剂化自由能（这为溶剂化自由能提供了静电作用），并为疏水作用添加了一个经验项：

通过计算受体与配体之间的平均相互作用能，并在必要/期望的情况下考虑结合时的熵变来获得delta-Gvacuum：

可以通过对这三种物质进行正常模式分析来找到熵的贡献，但是实际上，如果只希望比较具有类似熵的状态，例如两个配体结合到同一蛋白质上，则可以忽略熵的贡献。其原因是，正常模式分析计算的计算量很大，并且往往具有较大的误差范围，从而在结果中引入了很大的不确定性。

受体和配体的平均相互作用能通常是通过对从平衡分子动力学（MD）模拟中收集的不相关快照的整体进行计算而获得的。

在本教程中，我们将演示如何使用Amber和AmberTools随附的MM / PB（GB）SA脚本自动执行所有必要步骤，以估算蛋白质-蛋白质复合物（RAS和RAF）与蛋白质的结合自由能。 -MM配体（MM-GBSA和MM-PBSA）以串联和并联方式结合-配体（雌激素受体和雷洛昔芬）。此外，我们将演示使用该脚本使用丙氨酸扫描和普通模式熵计算。原则上，上述结合自由能的计算将需要对复合物和两个单个蛋白质进行三个独立的MD模拟。但是，通常一个近似值是结合后不会发生显着的构象变化，因此可以从单个轨迹获得所有三个物种的快照。这是“单一轨迹方法”，也是本教程中将使用的方法。

第1节：构建起始结构并运行仿真以获得平衡的系统。
第2部分：运行生产模拟并获取快照集合。
第3节：计算结合自由能并分析结果（此处是MM / PBSA不同版本之间的教程叉）。

一、建立起始结构并运行模拟以获得平衡的系统

我们将在此模拟中建模的系统是人类H-Ras蛋白与C-Raf1的Ras结合域（Ras-Raf）之间的复合体，它是信号转导级联的核心。这是RAS-RAF复合系统的部分平衡的，预先准备的pdb文件。

ras-raf.pdb

该结构包含ras和raf蛋白，以及生理必需的GTP核苷酸，如下图所示：

出于本教程的目的，并且为了简化起见，我们将避免在计算中处理GTP分子，因为这将需要为此化合物设置新参数，并且不在本教程的讨论范围之内。因此，我们只需从pdb文件中删除它即可将其从计算中删除。从上图可以看出，尽管不是严格正确，但这种近似在一定程度上是合理的，因为GTP并不直接参与绑定接口。

该蛋白质中还存在一个镁离子，该镁离子本质上与GTP分子结合，因此我们也将其去除。因此，您应该从pdb文件中删除残基243和244。

下一步是将该pdb文件拆分为两个单独的结构，以便您拥有ras-raf.pdb，ras.pdb和raf.pdb。我们将使用这三种结构来创建用于MM-PBSA计算的三个气相prmtop和inpcrd文件对，以及用于溶剂化复合物的一对，用于运行MD模拟：

> $AMBERHOME/bin/tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

注意：对于AMBER 14，请在使用 -f $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99 来加载ff99力场。

com = loadpdb ras-raf.pdb
ras = loadpdb ras.pdb
raf = loadpdb raf.pdb

确保为要使用的计算方法选择正确的半径。有关详细信息，请参见手册“ LEaP”部分中的PBRadii设置段落以及此处提供的建议。

set default PBRadii mbondi2

saveamberparm com ras-raf.prmtop ras-raf.inpcrd
saveamberparm ras ras.prmtop ras.inpcrd
saveamberparm raf raf.prmtop raf.inpcrd

现在，在退出tleap之前，您应该创建用于运行MD模拟的溶剂化复合物：

charge com
> Total unperturbed charge: -0.000000
> Total perturbed charge: -0.000000    (Hence there is no need to add counter ions)
solvatebox com TIP3PBOX 12.0
saveamberparm com ras-raf_solvated.prmtop ras-raf_solvated.inpcrd
quit

这是文件: ras-raf.prmtop, ras-raf.inpcrd, ras.prmtop, ras.inpcrd, raf.prmtop, raf.inpcrd, ras-raf_solvated.prmtop, ras-raf_solvated.inpcrd

具体下载地址参见：https://ambermd.org/tutorials/advanced/tutorial3/section1.htm

1.1 平衡溶剂化物

我们将通过短暂的最小化，50ps的加热和50ps的密度平衡（对复合物的约束较弱），然后在300K下进行500ps的恒压平衡来平衡溶剂化的复合物。所有模拟都将在氢原子摇动，2 fs的时间步长和用于温度控制的Langevin动力学下进行。输入文件如下：

min.in

minimise ras-raf
 &cntrl
  imin=1,maxcyc=1000,ncyc=500,
  cut=8.0,ntb=1,
  ntc=2,ntf=2,
  ntpr=100,
  ntr=1, restraintmask=':1-242',
  restraint_wt=2.0
 /

heat.in

heat ras-raf
 &cntrl
  imin=0,irest=0,ntx=1,
  nstlim=25000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=1,
  ntpr=500, ntwx=500,
  ntt=3, gamma_ln=2.0,
  tempi=0.0, temp0=300.0, ig=-1,
  ntr=1, restraintmask=':1-242',
  restraint_wt=2.0,
  nmropt=1
 /
 &wt TYPE='TEMP0', istep1=0, istep2=25000,
  value1=0.1, value2=300.0, /
 &wt TYPE='END' /

density.in

heat ras-raf
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=25000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=2, ntp=1, taup=1.0,
  ntpr=500, ntwx=500,
  ntt=3, gamma_ln=2.0,
  temp0=300.0, ig=-1,
  ntr=1, restraintmask=':1-242',
  restraint_wt=2.0,
 /

equil.in

heat ras-raf
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=250000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=2, ntp=1, taup=2.0,
  ntpr=1000, ntwx=1000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0, ig=-1,
 /

注意：在本教程的示例中，我们不会更改用于随机数生成器的随机种子的值。这由名称列表变量ig控制。这主要是针对教程设置中结果的可重复性问题。但是，在运行生产模拟时，尤其是在使用ntt = 2或3（Anderson或Langevin thermostats）时，必须在每次MD重新启动时将默认值更改为随机数种子。如果您使用的是AMBER 10（bugfix.26或更高版本）或AMBER 11或更高版本，则可以通过在cntrl名称列表中设置ig = -1来自动执行此操作。否则，您可以在每次重启计算时为ig指定一个正数。有关不执行此操作的陷阱的更多详细信息，请参考以下出版物：

Cerutti DS, Duke, B., et al., “A Vulnerability in Popular Molecular Dynamics Packages Concerning Langevin and Andersen Dynamics”, JCTC, 2008, 4, 1669-1680

您应该使用以下命令来运行所有4种模拟：

$AMBERHOME/bin/sander -O -i min.in -o min.out -p ras-raf_solvated.prmtop -c ras-raf_solvated.inpcrd \
-r min.rst -ref ras-raf_solvated.inpcrd
$AMBERHOME/bin/sander -O -i heat.in -o heat.out -p ras-raf_solvated.prmtop -c min.rst \
-r heat.rst -x heat.mdcrd -ref min.rst
gzip -9 heat.mdcrd
$AMBERHOME/bin/sander -O -i density.in -o density.out -p ras-raf_solvated.prmtop -c heat.rst \
-r density.rst -x density.mdcrd -ref heat.rst
gzip -9 density.mdcrd
$AMBERHOME/bin/sander -O -i equil.in -o equil.out -p ras-raf_solvated.prmtop -c density.rst \
-r equil.rst -x equil.mdcrd
gzip -9 equil.mdcrd

在1.7GHz IBM P690的16个处理器上，这大约需要5个小时。

以下是输出文件：equil.tar.gz

在继续进行MM-PBSA生产MD运行之前，我们需要验证系统是否平衡。为此，我们将研究温度，密度，总能量和RMSD。我们可以使用以下perl脚本（ process_mdout.pl ）开始，该脚本将从输出文件中提取有用的信息。

./process_mdout.pl heat.out density.out equil.out

由于第一次加热是在恒定体积条件下进行的，因此不会记录密度数据。因此，您将需要编辑summary.DENSITY文件并删除前50行。（因为xmgrace是愚蠢的，否则会感到困惑）。

xmgrace summary.DENSITY
xmgrace summary.TEMP
xmgrace summary.ETOT

另外，我们将针对最小化结构检查蛋白质骨架RMSD，以查看平衡过程中是否已实现构象稳定性。可以使用ptraj或cpptraj通过以下脚本完成此操作：

measure_equil_rmsd.ptraj

trajin equil.mdcrd.gz 1 250 1
reference ras-raf_solvated.inpcrd
rms reference out equil.rmsd @CA,C,N

xmgrace equil.rmsd

平衡期结束时，密度，温度和总能量图都已清晰收敛。虽然RMSD似乎开始趋于水平，但似乎尚未完全收敛，但就本教程而言，是可以接受的。在实际计算中，您可能希望运行更长的平衡时间，具体取决于您的系统。现在，我们准备执行生产运行。

二、运行生产模拟并获得完整的快照。

模拟的生产阶段应在与平衡的最终阶段相同的条件下运行，以防止由于模拟条件的变化而导致势能突然跳跃。

我们将总共运行2 ns或每10 ps记录一次生产。这应该足够远，以使结构不相关。取决于您的系统，将快照靠近在一起可以获得良好的效果。只要获得的所有结构都不相关，快照数量越多，结果的统计误差就应该越低。请注意，对于诸如RAS-RAF复杂系统之类的系统，此处的2 ns模拟时间很可能太短而无法获得一组足以对平衡集合进行充分采样的不相关快照。大约20 ns的值可能更合适。但是，这足以满足本教程的目的。

这是输入文件：

prod.in

prod ras-raf
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=250000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=2, ntp=1, taup=2.0,
  ntpr=5000, ntwx=5000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0, ig=-1,
 /

然后应运行4次以获得2 ns的模拟时间。由于这是简单的周期性边界PME模拟，因此可以根据需要使用PMEMD进行模拟。通常，这将提供更好的性能和并行扩展。以下是我在圣地亚哥超级计算机中心的Teragrid群集上使用的示例脚本，以在96个处理器上运行此作业。计算总共花了10个小时。

run.x

#SDSC Teragrid PBS Script
#PBS -j oe
#PBS -l nodes=48:ppn=2
#PBS -l walltime=12:00:00
#PBS -q dque
#PBS -V
#PBS -M name@email.com
#PBS -A account_no
#PBS -N run_pmemd_96

cd /gpfs/projects/prod/
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod1.out \
-p ras-raf_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod2.out \
-p ras-raf_solvated.prmtop -c prod1.rst -r prod2.rst -x prod2.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod3.out \
-p ras-raf_solvated.prmtop -c prod2.rst -r prod3.rst -x prod3.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod4.out \
-p ras-raf_solvated.prmtop -c prod3.rst -r prod4.rst -x prod4.mdcrd
gzip -9 prod*.mdcrd

以下是输出文件：prod.tar.gz（84.8 MB）

为了获得良好的结果，至关重要的是，我们的系统在生产阶段仍要探索平衡相空间。我们将通过绘制密度，温度，总能量和主链RMSD，以与最后一个平衡步骤相同的方式进行检查。

请注意，生产RMSD看起来并没有真正达到平衡，而其他属性基本上是恒定的（请注意小比例尺）。理想情况下，我们应该运行更长的生产时间（约20 ns）。但是，出于本教程的目的，我们将继续介绍现有内容。

现在，我们可以进行第3节，计算结合自由能。

第一个链接将带您进入使用（和安装）Python脚本MMPBSA.py的说明。
第二个链接将带您进入使用Perl脚本mm_pbsa.pl的说明。

三、计算结合自由能并分析结果

3.1 Python Script MMPBSA.py

在本教程中，我们将演示如何使用AmberTools中发布的MM-PBSA方法来计算结合自由能，运行丙氨酸扫描并计算熵校准的正常模式。本教程分为以下几类：

计算蛋白质-蛋白质复合物（Ras-Raf）的结合自由能。
计算蛋白质-配体复合物（雌激素受体和雷洛昔芬）的结合自由能。
计算Ras-Raf的结合自由能，并使用丙氨酸扫描与残基突变为丙氨酸的突变型Ras-Raf复合物的结合能进行比较，并分析结果。
使用三个处理器并行计算Ras-Raf的结合自由能。
使用正态模式分析（Nmode）计算雌激素受体和雷洛昔芬复合物的熵。
以每个残基或成对的每个残基为基础，将自由能的贡献分解为Ras-Raf的结合自由能。

除了6，其他几个的详情见后续的博文

3.2 计算结合自由能并分析结果。（此处使用mm_pbsa.pl执行的所有计算）

现在，我们需要从生产运行中提取快照（无水）以用于MM-PBSA计算。 mm_pbsa.pl脚本（在$ AMBERHOME / bin目录中）为我们自动执行此提取过程。请注意，如果未安装gzcat，则需要在运行mm_pbsa之前解压缩轨迹文件。我们必须提供以下输入文件（请参阅链接的文件以获取各个术语的解释）：

extract_coords.mmpbsa

@GENERAL

PREFIX                  snapshot
PATH                    ./
COMPLEX                 1
RECEPTOR                1
LIGAND                  1
COMPT                   ./ras-raf.prmtop
RECPT                   ./ras.prmtop
LIGPT                   ./raf.prmtop
GC                      1
AS                      0
DC                      0
MM                      0
GB                      0
PB                      0
MS                      0
NM                      0
@MAKECRD
BOX                     YES
NTOTAL                  42193
NSTART                  1
NSTOP                   200
NFREQ                   1
NUMBER_LIG_GROUPS       1
LSTART                  2622
LSTOP                   3862
NUMBER_REC_GROUPS       1
RSTART                  1
RSTOP                   2621
@TRAJECTORY
TRAJECTORY              ./prod1.mdcrd
TRAJECTORY              ./prod2.mdcrd
TRAJECTORY              ./prod3.mdcrd
TRAJECTORY              ./prod4.mdcrd
@PROGRAMS

这仅指定哪些原子是受体，配体和复合物的一部分，并指定与未溶剂化结构相对应的prmtop文件，弹道中快照的总数，步幅长度和弹道文件的名称。我们已指定每个输出文件都应以前缀“快照”编写。在此设置中，我们已将RAS定义为受体，将RAF定义为配体。这纯粹是一个命名约定。

$AMBERHOME/bin/mm_pbsa.pl extract_coords.mmpbsa > extract_coords.log

这将需要几分钟才能运行。以下是相关的输出文件：extract_coords.tar.gz（14 MB）

您应该检查日志文件中是否有任何错误。还要确保盒子的尺寸看起来合理。如果有任何奇怪的盒子尺寸，则可能是您选择的原子序号有问题或其中一个轨迹文件损坏了。

从我们刚刚提取的快照开始，我们现在将计算复合物，受体和配体的相互作用能和溶剂化自由能，并对结果取平均值以获得对结合自由能的估计。请注意，在本教程中，我们将不会计算熵对绑定的贡献，因此严格来说，我们的结果将不是真正的自由能，但可用于与类似系统进行比较。例如。可以对沿结合界面的氨基酸点突变的影响进行分析。常见的解决方法称为丙氨酸扫描。

为了演示的目的，我们将同时使用MM_PBSA方法和MM_GBSA方法进行结合能计算。这可以通过mm_pbsa.pl的以下输入文件来完成（请参阅链接的文件以获取注释）：

binding_energy.mmpbsa

VERBOSE               0
PARALLEL              0
PREFIX                snapshot
PATH                  ./
START                 1
STOP                  200
OFFSET                1
COMPLEX               1
RECEPTOR              1
LIGAND                1
COMPT                 ./ras-raf.prmtop
RECPT                 ./ras.prmtop
LIGPT                 ./raf.prmtop
GC                    0
AS                    0
DC                    0
MM                    1
GB                    1
PB                    1
MS                    1
NM                    0
@PB
PROC                  2
REFE                  0
INDI                  1.0
EXDI                  80.0
SCALE                 2
LINIT                 1000
ISTRNG                0.0
RADIOPT               0
ARCRES                0.0625
INP                   1
SURFTEN               0.005
SURFOFF               0.00
IVCAP                 0
CUTCAP                -1.0
XCAP                  0.0
YCAP                  0.0
ZCAP                  0.0
@MM
DIELC                 1.0
@GB
IGB                   2
GBSA                  1
SALTCON               0.00
EXTDIEL               80.0
INTDIEL               1.0
SURFTEN               0.005
SURFOFF               0.00
@MS
PROBE                 0.0
@PROGRAMS

此输入文件的各个部分指定要运行的计算。在哪些文件上运行它们以及计算对结合自由能的不同贡献所必需的任何特殊参数。如果您在上方打开链接的文件，则会发现不同术语的解释。不同参数的值是根据经验数据确定的，并有待不断研究。请在此处找到用于计算非极性溶剂化自由能贡献值的当前推荐设置的列表。可以在 $AMBERHOME/AmberTools/src/mm_pbsa/Examples/TEMPLATE_INPUT_SCRIPTS 目录中找到将自由能计算与常用参数设置绑定在一起的示例输入文件。有关更多信息，请参考该领域的出版物。请注意，AMBER的早期版本需要外部泊松Boltzmann解算器（例如Delphi），但从AMBER 8开始，可以使用内部PBSA程序。该程序可显着提高速度，并更易于集成到AMBER框架中。您可以使用以下方法运行计算：

$AMBERHOME/bin/mm_pbsa.pl binding_energy.mmpbsa > binding_energy.log

您可以通过以下方式查看计算进度：

tail -f binding_energy.log

该计算将花费大约2个小时来运行（P4 3.2 GHz）。在600个快照中，每个快照的PBSA部分都需要花费大部分时间。计算的GBSA部分以秒为单位。为了加快计算速度，您还可以通过在PARALLEL下指定可用处理器的数量，来请求一次并行执行多个快照的mm_pbsa分析。

完成后，您应该找到以下输出文件：binding_energy.log，snapshot_statistics.out，snapshot_com.all.out，snapshot_rec.all.out，snapshot_lig.all.out

生成文件见： https://ambermd.org/tutorials/advanced/tutorial3/section3.htm

日志文件仅显示计算是否成功完成。 all.out文件为每种物种的每个快照提供了各自的能量贡献，而statistics.out文件包含最终的平均结合自由能结果：

snapshot_statistics.out

#                  COMPLEX                RECEPTOR                  LIGAND        
#          ----------------------- ----------------------- -----------------------
#                  MEAN        STD         MEAN        STD         MEAN        STD
#          ======================= ======================= =======================
ELE            -8656.78      70.18     -5602.09      63.10     -2102.25      52.57
VDW             -984.99      24.34      -661.18      20.33      -256.02      12.93
INT             5085.33      50.22      3449.57      38.65      1635.76      29.42
GAS            -4556.44      75.96     -2813.70      65.21      -722.52      53.50
PBSUR             65.09       1.05        45.25       0.64        27.24       0.46
PBCAL          -3223.64      58.68     -2490.86      48.73     -1671.27      47.46
PBSOL          -3158.55      58.26     -2445.62      48.45     -1644.03      47.31
PBELE         -11880.42      34.25     -8092.96      29.34     -3773.52      17.30
PBTOT          -7714.99      48.25     -5259.32      36.97     -2366.55      26.61
GBSUR             65.09       1.05        45.25       0.64        27.24       0.46
GB             -3407.82      58.49     -2631.83      50.08     -1731.06      47.68
GBSOL          -3342.74      58.15     -2586.58      49.83     -1703.82      47.55
GBELE         -12064.60      26.94     -8233.92      23.57     -3833.31      13.40
GBTOT          -7899.17      47.07     -5400.28      35.65     -2426.34      26.80

#                    DELTA        
#          -----------------------
#                  MEAN        STD
#          =======================
ELE             -952.43      44.10
VDW              -67.79       5.18
INT               -0.00       0.00
GAS            -1020.22      44.58
PBSUR             -7.40       0.41
PBCAL            938.50      42.51
PBSOL            931.09      42.31
PBELE            -13.94       9.43
PBTOT            -89.13       7.94
GBSUR             -7.40       0.41
GB               955.07      41.30
GBSOL            947.66      41.10
GBELE              2.63       7.41
GBTOT            -72.56       6.40

该文件中不同术语的含义如下：

ELE =由MM力场计算出的静电能。 VDW = MM的van der Waals贡献。 INT =由MM力场中的键，角和二面角项产生的内部能量。（在单轨迹方法中，该项始终等于零）。 GAS =气相总能量（ELE，VDW和INT之和）。 PBSUR / GBSUR =通过经验模型计算得出的对溶剂化自由能的非极性贡献。 PBCAL / GB =分别由PB或GB计算得出的对溶剂化自由能的静电贡献。 PBSOL / GBSOL =对溶剂化的非极性和极性贡献之和。 PBELE / GBELE =静电溶剂化自由能和MM静电能之和。 PBTOT / GBTOT =根据上述条件计算得出的最终估计结合自由能。（千卡/摩尔）

这里有几件事要注意。通常，ELE，PBCAL/GBCAL和VDW部分是对结合自由能的主要贡献，并且这两个术语中的前两个趋于彼此抵消，可以通过PBELE/GBELE的值来检查比它的贡献小。

人们通常会期望找到一种非常有利的静电能和一种不利的溶剂化自由能。这象征着人们必须使用的能量来使结合颗粒去溶剂化并对齐其结合界面。

从负的总结合自由能为-89.13 kcal / mol，我们可以清楚地看到，这是纯水中的一种有利的蛋白质-蛋白质复合物，但请记住，由于我们没有估算（dis （有利）熵对结合的贡献。请注意，GB方法给出的结合能略低，但仍表明这是一种有利的结合态。

要扩展本教程，您可以研究如果更改溶剂的盐含量，修改特定的残基或选择不同的起始几何形状会发生什么情况。您也可以考虑运行nmode计算来寻找熵，但是要知道，对于这种大小的复合体，这将是非常大的内存，并且计算量很大。

参考资料

https://ambermd.org/tutorials/advanced/tutorial3/index.htm

这里是一个广告位，，感兴趣的都可以发邮件聊聊：tiehan@sina.cn

个人公众号，比较懒，很少更新，可以在上面提问题，如果回复不及时，可发邮件给我： tiehan@sina.cn