AMBER高级教程A3:MM-PBSA

类别:    标签: amber   阅读次数:   版权: (CC) BY-NC-SA

本教程使用mm_pbsa脚本计算RAS-RAF蛋白复合物的结合能, 并对计算中的每一步进行解释说明. 同时教程还包括了对如何使用mmpbsa_py脚本进行这种计算的说明.

AMBER高级教程A3: MM-PBSA

MM-PBSA方法可用于计算蛋白-配体复合物间的结合自由能, 其中配体可以是蛋白, 也可以是小分子, 多肽等等. 在本教程中, 我们将使用MM-PBSA方法计算两个蛋白结合的结合能.

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

\[[A]_{aq} +[B]_{aq} \iff [A^* B^*]_{aq^*}\]

理想情况下, 我们想直接计算如下图所示的结合自由能:

然而, 在这些溶剂化状态的模拟中, 大部分能量贡献将来自溶剂-溶剂相互作用, 并且总能量的波动将比结合能大一个数量级, 因此结合能的计算需要非常长的时间才能收敛. 因而更有效的方法是按照下面的热力学循环将计算划分开来:

显然, 从上图可以看出结合自由能 $\D G_\text{bind,solv}$ 可以通过下式计算:

\[\D G_\text{bind,solv}^0=\D G_\text{bind,vacuum}^0+\D G_\text{solv,complex}^0-(\D G_\text{solv,ligand}^0+\D G_\text{solv,receptor}^0)\]

在MM-PBSA方法中, 对上述结合自由能的不同贡献以不同方式进行计算:

\[\D G_\text{solv}^0=\D G_{electrostatic,\e=80}^0-\D G_{electrostatic,\e=1}^0+\D G_{hydrophobic}^0\] \[\D G_\text{vacuum}^0=\D E_{molecula\ mechanics}^0-T \D S_{norma\ mode\ analysis}^0\]

在本教程中, 我们将演示使用Amber和AmberTools自带的MM/PB(GB)SA脚本自动执行所有必要步骤, 对蛋白-蛋白复合物(RAS和RAF)和蛋白-配体复合物(雌激素受体和雷洛昔芬)的结合自由能进行计算, 计算时MM-GBSA和MM-PBSA采用串行和并行两种模式. 此外, 我们还将演示使用脚本进行丙氨酸扫描和简正模式熵计算. 原则上, 上述结合自由能的计算需要对复合物以及两种单独蛋白质进行三次独立的MD模拟. 不过, 通常我们近似地认为, 在结合时不会发生显著的构象变化, 从而可以从单一轨迹获得所有三个物种的快照, 这称为”单轨迹方法”. 本教程就采用这种方法.

1. 构建初始结构, 运行模拟获得平衡体系

在模拟中我们将要建模的体系是人类H-Ras蛋白与C-Raf1的Ras结合结构域之间的复合物(Ras-Raf), 它是信号转导级联的中心. 这里是一个部分平衡的, 预先准备好的RAS-RAF复合物的pdb文件ras-raf.pdb.

结构中包含了ras和raf蛋白, 另外还有一个生理上必需的GTP核苷酸, 如下图所示:

出于本教程的目的, 为简单起见, 我们在计算中不处理GTP分子, 因为这需要为该化合物设置新的参数, 内容超出了本教程的范围. 因此, 我们从pdb文件中删除它, 这样就简单地将它从计算中删除了. 虽然并非严格正确, 但这种近似是合理的, 因为从上图可以看出, GTP并不直接参与结合面.

另外, 这个蛋白质中还含有一个镁离子, 它主要与GTP分子结合, 因此我们也将其删除. 因此, 你应该从pdb文件中删除残基243和244.

下一步是将res-raf.pdb文件分成两个独立的结构, 这样就有了ras-raf.pdb, ras.pdbraf.pdb. 我们将使用这三个结构创建三个气相prmtopinpcrd文件组用于MM-PBSA计算, 以及溶剂化复合物的一个文件组用于运行MD模拟:

tleap -s -f oldff/leaprc.ff14SB

注意: 对AMBER 14, tleap调用中请使用-f $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99来加载ff99力场

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

确保为要使用的计算方法选择了正确的半径. 详细信息见手册中LEaP节的PBRadii部分. 如果需要设定特殊的原子半径, 可以参考推荐recommended_settings.pdf.

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    (因此无需添加抗衡离子)

solvatebox com TIP3PBOX 12.0
saveamberparm com ras-raf_solvated.prmtop ras-raf_solvated.inpcrd
quit

上述命令的输入可以使用脚本进行简化:

tleap -s -f tleap.in > tleap.out

其中tleap.in文件内容为:

tleap.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
source leaprc.protein.ff14SB
source leaprc.water.spce
source leaprc.gaff
com = loadpdb ras-raf.pdb
ras = loadpdb ras.pdb
raf = loadpdb raf.pdb
set default PBRadii mbondi2
saveamberparm com ras-raf.prmtop ras-raf.inpcrd
saveamberparm ras ras.prmtop ras.inpcrd
saveamberparm raf raf.prmtop raf.inpcrd
charge com
solvatebox com TIP3PBOX 12.0
saveamberparm com ras-raf_solvated.prmtop ras-raf_solvated.inpcrd
quit

下载相关文件:

1.1 平衡溶剂化复合体

我们将对溶剂化复合物按以下步骤进行平衡: 短的能量最小化, 50 ps的升温, 对复合物进行弱限制条件下50 ps的密度平衡, 300 K下500 ps的等压平衡. 所有的模拟在运行时都对氢原子使用SHAKE约束, 并使用2 fs的时间步长和Langevin动力学控制温度. 输入文件如下:

min.in
1
2
3
4
5
6
7
8
9
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
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
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
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
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=23(Anderson或Langevin恒温器)时, 每次 重新启动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, 16核的IBM P690机器上, 完成上面的模拟大于需要5小时.

可以在这里下载输出文件: equil.tar.gz

在我们进行MM-PBSA的成品模拟之前, 我们需要验证体系是否已经平衡. 为此, 我们从温度, 密度, 总能量和RMSD四个方面进行分析. 我们可以使用perl脚本(process_mdout.pl)从输出文件中提取有用的信息. 命令如下:

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

此外, 我们还要检查相对于最小化结构的蛋白质骨架算RMSD, 以确定在平衡期间构象是否已经稳定. 这可以使用measure_equil_rmsd.ptraj脚本借助ptrajcpptraj完成:

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

由于第一次对复合物体系的升温模拟是在等容条件下进行的, 并没有记录体系的密度数据. 因此需要编辑summary.DENSITY文件删除前50行(之所以这么做是因为xmgrace太愚蠢无法处理这种情况).

对复合物体系的温度, 密度, 总能量和RMSD进行绘图:

xmgrace summary.DENSITY
xmgrace summary.TEMP
xmgrace summary.ETOT
xmgrace equil.rmsd

结果如下:

在平衡期结束时, 密度, 温度和总能量曲线都明显地收敛. 看起来开始稳定的RMSD似乎没有完全收敛, 但鉴于本教程的目的也是可以接受的. 在实际计算中, 根据体系的大小, 我们需要运行更长的平衡时间. 接下来我们将运行成品模拟.

2. 运行成品模拟获得轨迹快照集

运行成品模拟时应该使用与平衡的最后阶段相同的条件, 以防止由于模拟条件变化而引起的势能突变.

我们将运行总共2 ns的成品模拟, 每10 ps记录一次坐标. 10 ps的时间间隔足够大能保证结构彼此是不相关的. 取决于你研究的体系, 使用更小时间间隔的快照也可能获得好的结果. 只要你获得的所有结构彼此不相关, 快照数目越多, 结果的统计误差应该越低. 对于RAS-RAF复合物这样的体系, 2 ns的模拟时间可能太短, 不足以获得足够多的不相关快照以对平衡系综进行充分采样. 20 ns左右的模拟时长可能更合适. 但是, 2 ns对于本教程来说已经足够了.

输入文件prod.in如下:

prod.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
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来进行模拟. PMEMD通常性能更好, 并且可以并行扩展. 下面是我在San Diego超算中心的Teragrid集群上使用96个核来运行上面的作业所用的PBS脚本run.x:

run.x
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#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

绘图结果如下:

从成品模拟RMSD的波动趋势上看, 体系还没有达到平衡状态, 但其他性质基本上是恒定的(注意纵轴标尺的大小). 理想情况下, 我们应该延长成品模拟的时长(约20 ns). 然而, 鉴于本教程的目的, 我们就使用2 ns的轨迹进行下面的步骤.

在下一节中我们将计算结合自由能, 可采用的方式有两种: 第一种使用(需要先安装)Python脚本MMPBSA.py, 第二种使用Perl脚本mm_pbsa.pl.

3. 计算结合自由能并分析结果

使用Perl脚本mm_pbsa.pl计算结合自由能

我们需要从成品模拟轨迹中抽取快照(不含溶剂水)用于MM-PBSA计算. mm_pbsa.pl脚本(位于$AMBERHOME/bin目录)可以自动完成这个提取过程. 请注意, 如果没有安装gzcat, 则需要在运行mm_pbsa之前将轨迹文件解压缩. 我们必须提供如下的输入文件extract_coords.mmpbsa(下载的文件中包含对每项的解释说明):

extract_coords.mmpbsa
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
@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文件, 轨迹中的快照总数, 提取步长和轨迹文件的名称. 我们还指定了每个输出文件都以snapshot作为前缀. 在本教程中, 我们定义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实现的(下载的文件中有注释):

binding_energy.mmpbsa
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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

该输入文件的各个部分指定了需要运行的计算, 计算时需要的文件以及计算结合自由能的不同贡献时所需要的所有特殊参数. 如果你打开下载的文件, 可以看到其中对不同项的解释说明. 不同参数的数值是根据经验数据确定的, 有待于当前研究的验证. 对于非极性溶剂化自由能的计算, Recommendation_setting.pdf文件给出当前的推荐设置. 计算结合自由能的示例输入文件以及常用的参数设置可以在$AMBERHOME/AmberTools/src/mm_pbsa/Examples/TEMPLATE_INPUT_SCRIPTS目录中找到. 更多信息请参考相关文献. 请注意, 早期版本的AMBER需要额外的Poisson-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

binding_energy.log文件只显示了计算是否成功完成. all.out文件给出了每个物种每一快照的单独的能量贡献, 而statistics.out文件包含了平均结合自由能的最终结果:

snapshot_statistics.out
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#                  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, PBCAL/GBCALVDW部分, 前两项大致相互抵消, 这可以通过查看PBELE/GBELE的值是否远小于ELEPBCAL/GBCAL对它的贡献进行检查.

通常我们会期望能发现非常有利的静电能和不利的溶剂化自由能, 这意味着结合粒子去溶剂化并对齐到结合界面的能量.

总结合自由能-89.13 kcal/mol为负值, 这表明在纯水中, 该蛋白-蛋白复合体的结合是有利的. 但是, 要记住, 该结果并不等于真正的结合自由能, 因为我们没有考虑(不利的)熵对结合自由能的贡献. 注意到, GB方法给出的结合能稍低, 但仍然表明这是一个有利的结合状态.

若要扩展本教程, 可以研究更改溶剂的盐含量, 修改特定残基或选择不同的起始结构对于结合自由能的影响. 另外, 也可以考虑利用nmode来计算熵变, 但要注意, 对于这种规模的复合物来说, 计算熵将会耗费大量的内存和时间.

利用Python脚本MMPBSA.py计算结合自由能

在本教程中, 我们将演示使用AmberTools中发布的MM-PBSA方法计算结合自由能, 运行丙氨酸扫描, 进行简正模式分析以计算熵变. 教程分为如下几部分:

3.1 计算蛋白-蛋白复合物(Ras-Raf)的结合自由能

本小节中, 我们将模拟人类H-Ras蛋白与结合于Ras结合结构域的C-Raf1形成的复合物(Ras-Raf), 该复合物是信号转导级联的中心. 部分平衡的经过预处理的RAS-RAF复合体的pdb结构如下图所示. 该结构中含有ras和raf蛋白, 另外还有一个生理必需的GTP核酸.

关于如何构建初始结构和运行动力学模拟以获得平衡体系, 请参考前面几节.

使用MMPBSA.py计算结合自由能的重要文件是拓扑文件和mdcrd文件(ras-raf_top_mdcrd.tgz).

我们现在将计算复合物, 受体和配体的相互作用能和溶剂化自由能, 并对结果进行平均以获得结合自由能的估计值. 请注意, 在教程的这一部分我们不计算熵对结合能的贡献, 因此严格来说, 所得的结果并不是真正的自由能, 但可以用来对相似的体系进行比较. 有关使用简正模式分析(Nmode)计算体系的熵贡献, 可参考3.5节. 也可以取消下面输入文件中&general命名列表中最后一行的注释, 以便使用AMBER的ptraj模块进行准简谐熵计算.

我们将分别使用MM-GBSA方法和MM-PBSA方法进行结合自由能的计算并进行比较. 以下为MMPBSA.py的输入文件mmpbsa.in:

mmpbsa.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
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. mmpbsa.in被分成三个命名列表: general, pbgb. general命名列表旨在指定并非仅适用于计算的特定部分, 而是针对所有部分的变量. 在本教程中, 我们将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

这将交互式地运行脚本, 并将计算进度输出至标准输出终端, 将任何错误或警告输出至标注错误终端. 最后, 计算完成后会在终端上输出总耗时以及计算过程中每一步骤的耗时.

另外, 命令行参数可以用shell识别的通配符(即bash的*?). 例如, 命令行中的-y *.mdcrd会告知程序读入工作目录中所有以.mdcrd结尾的文件, 并将其作为待分析的轨迹.

下载脚本生成的输出文件: `pb_gb_output1.tgz.

脚本使用ptraj生成了三个非溶剂化的mdcrd文件(复合物, 受体, 配体), 它们是在GB和PB计算过程中分析过的坐标. *.mdout文件包含所有指定帧的能量. 还会生成一个平均结构的PDB文件, 其中的结构(通过RMS)对齐到所有快照. 如果需要, 可以使用ptraj对这个结构进行准简谐熵计算. MMPBSA.py生成的所有文件的名称均以前缀_MMPBSA_开始, 除了最终的输出文件FINAL_RESULTS_MMPBSA.dat之外:

FINAL_RESULTS_MMPBSA.dat
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
| 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求解器(如果使用的话). 统计文件的其余内容包括所有的平均能量, 标准偏差和平均值的标准误差, 先列出GB对应的值, 再列出PB对应的值. 每个部分之后, 会给出结合自由能ΔG及误差. 文件中不同项的含义如下:

值得注意的是, 结果中并未报告总的气相能量, 因为使用单一轨迹计算方法时, 受体和配体的成键势能项的值会精确地与复合体中的对应值互相抵消. 如果这两项能量在允许的误差范围内不能抵消, 将会给出错误信息.

通常我们会期望能发现非常有利的静电能和不利的溶剂化自由能, 这意味着结合粒子去溶剂化并对齐到结合界面的能量.

总结合自由能-86.36 kcal/mol为负值, 这表明在纯水中, 该蛋白-蛋白复合体的结合是有利的. 但是, 要记住, 该结果并不等于真正的结合自由能, 因为我们没有考虑(不利的)熵对结合自由能的贡献. 注意到, GB方法给出的结合能稍低, 但仍然表明这是一个有利的结合状态.

3.2 使用三个处理器并行计算Ras-Raf的结合自由能

在本节中, 我们将并行地计算结合自由能. MMPBSA.py.MPI通过为每个线程(进程)分配相同数目的帧进行并行化计算. 因此, 当处理的帧数是所启动的线程数的倍数时, 它的运行效率最高. 但是, 这并不是必须的. 如果帧数不是线程数的倍数, 则”剩余”帧将在启动线程数的一个子集中均匀分配. 例如, 在3个线程上计算50帧将导致2个线程计算17帧, 最后一个线程只计算16帧. 因此, 第三个线程将不得不等待前两个线程完成计算后才能继续进行下面的计算. 出于这个原因, 使用5个线程将是更明智的选择(每个线程处理10帧). 但是, 线程数不能超过要处理的帧数, 否则MMPBSA.py.MPI将会终止, 给出错误信息. MMPBSA.py.MPI的输入文件与MMPBSA.py的输入文件完全相同:

mmpbsa.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
Input file for running PB and GB
&general
   endframe=50, verbose=1,
#  entropy=1,
/
&gb
  igb=2, saltcon=0.100
/
&pb
  istrng=0.100,
/

运行命令如下:

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内容如下:

parallel.job
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#!/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
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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

keep_files设置为默认值1, 输出文件为Parallel_output.tgz

最终结果为FINAL_RESULTS_MMPBSA.dat:

FINAL_RESULTS_MMPBSA.dat
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
| 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)

3.3: 计算蛋白-配体复合物(雌激素受体和雷洛昔芬)的结合自由能

1. 建立起始结构, 通过模拟获得平衡体系

本小节中, 我们将模拟雌激素受体蛋白和雷洛昔芬配体形成的蛋白-配体复合物. 其预处理过的pdb文件如下Estrogen_Receptor-Raloxifene.pdb.

该结构包含雌激素受体蛋白以及配体雷洛昔芬, 雷洛昔芬已经锚定在了受体蛋白上, 如下图所示:

这个体系的构建方式与第1节中的Ras-Raf体系类似. 关于如何构建起始结构和运行动力学模拟以获得平衡体系, 请参考第1节和第2节. 另外, 必须利用antechamber获得雷洛昔芬的正确参数. 详细说明请参考AMBER基础教程B4.

MD模拟中, 使用MMPBSA.py计算结合自由能所需的重要文件是MD模拟所用的拓扑文件和mdcrd文件, 下载Est_Rec_top_mdcrd.tgz.

2. 计算雌激素受体和雷洛昔芬的结合自由能

此节说明文件与前几节几乎相同, 故此省略, 只给出命令和结果.

MMPBSA.py计算的输入文件mmpbsa.in:

mmpbsa.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
Input file for running PB and GB
&general
   endframe=50, keep_files=2,
/
&gb
  igb=2, saltcon=0.100,
/
&pb
  istrng=0.100,
/

运行命令:

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y *.mdcrd

使用keep_files=2, 得到的输出文件为pb_gb_output2.tgz

最终结果

FINAL_RESULTS_MMPBSA.dat
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
| Run on Thu Feb 11 12:44:26 EST 2010

|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB
|&general
|   endframe=50, keep_files=2,
|/
|&gb
|  igb=2, saltcon=0.100,
|/
|&pb
|  istrng=0.100,
|/
|--------------------------------------------------------------
|Solvated complex topology file:  1err.solvated.prmtop
|Complex topology file:           complex.prmtop
|Receptor topology file:          receptor.prmtop
|Ligand topology file:            ligand.prmtop
|Initial mdcrd(s):                1err_prod.mdcrd
|
|Best guess for receptor mask:   ":1-240"
|Best guess for  ligand  mask:   ":241"
|Ligand residue name is "RAL"
|
|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                  -2013.3801               20.3021              2.8712
EEL                     -16938.6450               85.7631             12.1287
EGB                      -3507.0086               67.7839              9.5861
ESURF                       97.5448                1.3301              0.1881

G gas                   -18952.0251               88.1333             12.4639
G solv                   -3409.4639               67.7969              9.5879

TOTAL                   -22361.4889               47.1982              6.6748

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.2272               19.2311              2.7197
EEL                     -16895.0354               85.5797             12.1028
EGB                      -3528.7276               68.3585              9.6673
ESURF                      101.2613                1.3071              0.1849

G gas                   -18850.2626               87.7138             12.4046
G solv                   -3427.4663               68.3710              9.6691

TOTAL                   -22277.7288               48.1057              6.8032

 
Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.8595                2.0516              0.2901
EEL                         -5.5796                2.0333              0.2876
EGB                        -28.4863                0.6040              0.0854
ESURF                        4.4326                0.0462              0.0065

 
G gas                       -7.4391                2.8885              0.4085
G solv                     -24.0538                0.6058              0.0857

TOTAL                      -31.4929                5.0748              0.7177

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.2934                2.9265              0.4139
EEL                        -38.0300                3.2114              0.4542
EGB                         50.2053                2.5869              0.3658
ESURF                       -8.1491                0.2589              0.0366

 
DELTA G gas                -94.3234                4.3449              0.6145
DELTA G solv                42.0562                2.5999              0.3677

 
 DELTA G binding =        -52.2672     +/-      2.4568                 0.3475
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

POISSON BOLTZMANN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -2013.3801               20.3021              2.8712
EEL                     -16938.6450               85.7631             12.1287
EPB                      -3329.1708               67.0354              9.4802
ECAVITY                     68.2656                0.5195              0.0735

G gas                   -18952.0251             7767.4837           1098.4881
G solv                   -3260.9052               67.0374              9.4805

TOTAL                    -5265.0831               49.0426              6.9357

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.2272               19.2311              2.7197
EEL                     -16895.0354               85.5797             12.1028
EPB                      -3355.4746               67.3299              9.5219
ECAVITY                     70.1184                0.5285              0.0747

G gas                   -18850.2626             7693.7163           1088.0558
G solv                   -3285.3562               67.3320              9.5222

TOTAL                    -5279.4509               50.4067              7.1286

 
Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.8595                2.0516              0.2901
EEL                         -5.5796                2.0333              0.2876
EPB                        -31.3364                0.6953              0.0983
ECAVITY                      3.1896                0.0288              0.0041

G gas                       -7.4391                8.3434              1.1799
G solv                     -28.1468                0.6959              0.0984

TOTAL                       56.0934                5.0476              0.7138

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.2934                2.9265              0.4139
EEL                        -38.0300                3.2114              0.4542
EPB                         57.6402                3.0642              0.4333
ECAVITY                     -5.0423                0.1683              0.0238

DELTA G gas                -94.3234               18.8778              2.6697
DELTA G solv                52.5978                3.0688              0.4340

 
 
 DELTA G binding =        -41.7256     +/-      2.9618                 0.4189
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

3.4 计算Ras-Raf的结合自由能, 并利用丙氨酸扫描(Alanine Scanning)比较分析Ras-Raf复合物中某个残基突变为丙氨酸后的结果

1. 设置pdb文件, 为后续的tleap做准备

部分平衡后预处理过的RAS-RAF复合物的结构文件ras-raf.pdb

接下来我们必须准备突变体结构的pdb文件供tleap使用. 为了保证prmtop文件的一致性, 我们强烈建议在运行任何模拟之前准备该突变体的pdb及其拓扑文件, 并与第1节中创建的初始拓扑文件一起. 本教程中, 我们选择将残基21(异亮氨酸, I21)突变为丙氨酸, 因为I21位于受体和配体结合的界面处, 应该对结合能有明显的影响. 请注意, 本教程涉及的代码目前只适用于丙氨酸突变.

由于I21仅处于受体中, 所以我们不需要准备突变配体的pdb文件. 因此, 我们只需要修改ras-raf.pdbras.pdb. 为此, 你需要了解一些所涉及到的氨基酸结构的知识. 异亮氨酸的侧链是-CH(CH3)CH2CH3, 而丙氨酸的侧链是-CH3. 由于异亮氨酸侧链的原子数比丙氨酸侧链的多, 因此我们必须从pdb文件中去除原子及其相应的信息(名称, 编号, 坐标等). 这种突变涉及除β-C(CB)之外的所有侧链原子. 在I21中, 这意味着我们必须删除Ras-Raf和Ras的pdb文件中与原子294到305相对应的行. 我们不需要添加β-H(HB)原子, 因为基于为体系选择的特定库文件, tleap会将这些原子添加到合适的位置. 最后, 将残基21所有剩余的原子对应的残基名称从ILE更改为ALA. 此过程将生成两个突变的pdb文件: ras-raf_mutant.pdbras_mutant.pdb. RAS-RAF及其I21A突变结构如下图所示:

其他残基的突变也可以使用类似的方法, 其中位于CB之后羰基碳(C)之前的原子组可以从pdb文件中去除. 值得注意的是, 单次计算中只能进行一次突变.

2. 构建起始拓扑和坐标文件, 模拟获得平衡体系

基于已经生成的突变pdb文件, 可用tleap为这些结构构建相应的拓扑和坐标文件. 首先, 我们将生成对应非突变复合体的文件:

$AMBERHOME/exe/tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99SB

com = loadpdb ras-raf.pdb
ras = loadpdb ras.pdb
raf = loadpdb raf.pdb
saveamberparm com ras-raf.prmtop ras-raf.inpcrd
saveamberparm ras ras.prmtop ras.inpcrd
saveamberparm raf raf.prmtop raf.inpcrd

构建用于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

在退出tleap之前, 构建与突变pdb文件相应的拓扑和坐标文件:

com_mut = loadpdb ras-raf_mutant.pdb
ras_mut = loadpdb ras_mutant.pdb
saveamberparm com_mut rasraf_mutant.prmtop rasraf_mutant.inpcrd
saveamberparm ras_mut ras_mutant.prmtop ras_mutant.inpcrd
quit

上述步骤共生成了12个文件(6个.prmtop文件和6个.inpcrd文件). 非突变体.prmtop.inpcrd文件用于运行MD模拟以获得平衡体系, 具体方法可参考第1节和第2节.

利用MMPBSA.py计算结合自由能的重要文件是(非突变体和突变体的)拓扑文件, 以及使用非突变体的拓扑文件和坐标文件运行得到的mdcrd文件ras-raf_alascan.tgz

3. 对Ras-Raf的结合自由能进行丙氨酸扫描计算

我们现在将计算复合物, 受体和配体的相互作用能和溶剂化自由能, 并对结果进行平均以获得结合自由能的估计值. 然后, 为了与”野生型”进行比较, 我们会对突变后的结构进行相同的计算, 计算前需要先将mdcrd中的坐标突变. 请注意, 在教程的这一部分我们不计算熵对结合能的贡献, 因此严格来说, 所得的结果并不是真正的自由能, 但可以用来对相似的体系进行比较. 有关使用简正模式分析(Nmode)计算体系的熵贡献, 可参考3.5节. 也可以取消下面输入文件中&general命名列表中最后一行的注释, 以便使用AMBER的ptraj模块进行准简谐熵计算.

计算方法与前面的相同, MMPBSA.py的输入文件mmpbsa.in如下:

mmpbsa.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
sample input file for running alanine scanning
 &general
   startframe=1, endframe=50, interval=1,
   verbose=1,
/
&gb
  saltcon=0.1
/
&pb
  istrng=0.100
/
&alanine_scanning
/

文件中的&alanine_scanning命名列表表示初始化脚本中的丙氨酸扫描, 其中唯一可使用的变量是mutant_only, 其详细说明见手册.

运行命令如下:

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -sp ras-raf_solvated.prmtop -cp rasraf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd -mc rasraf_mutant.prmtop -mr ras_mutant.prmtop

输出文件ALASCAN_output.tgz

最终结果FINAL_RESULTS_MMPBSA.dat

FINAL_RESULTS_MMPBSA.dat
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
| Run on Thu Feb 11 13:11:48 EST 2010

|Input file:
|--------------------------------------------------------------
|sample input file for running alanine scanning
| &general
|   startframe=1, endframe=50, interval=1,
|   verbose=1,
|/
|&gb
|  saltcon=0.1
|/
|&pb
|  istrng=0.100
|/
|&alanine_scanning
|/
|--------------------------------------------------------------
|Solvated complex topology file:  ras-raf_solvated.prmtop
|Complex topology file:           rasraf.prmtop
|Receptor topology file:          ras.prmtop
|Ligand topology file:            raf.prmtop
|Initial mdcrd(s):                bigprod.mdcrd
|Mutant complex topology file:    rasraf_mutant.prmtop
|Mutant receptor topology file:   ras_mutant.prmtop
|Mutant ligand topology file:     raf.prmtop
|
|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                      -3142.2247               63.1977              8.9375
ESURF                       91.3565                1.3938              0.1971

G gas                   -19064.5240               77.8536             11.0102
G solv                   -3050.8682               63.2131              8.9397

TOTAL                   -22115.3922               51.5332              7.2879

 
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                      -2444.8629               54.9156              7.7662
ESURF                       64.2843                1.1143              0.1576

 
G gas                   -12825.2661               73.1118             10.3396
G solv                   -2380.5786               54.9269              7.7678

TOTAL                   -15205.8447               36.8422              5.2103

 
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                      -1661.8286               26.5442              3.7539
ESURF                       37.0493                0.6185              0.0875

 
G gas                    -5213.7811               37.3522              5.2824
G solv                   -1624.7794               26.5514              3.7549

TOTAL                    -6838.5604               25.6515              3.6277

 
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                        964.4668               32.9201              4.6556
ESURF                       -9.9770                0.3759              0.0532

 
DELTA G gas              -1025.4769               35.1797              4.9752
DELTA G solv               954.4898               32.9223              4.6559

 
 DELTA G binding =        -70.9871     +/-      6.6875                 0.9457
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
I21A MUTANT:
GENERALIZED BORN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1855.4226               17.0765              2.4150
EEL                     -17210.2882               75.8866             10.7320
EGB                      -3145.1010               63.2477              8.9446
ESURF                       91.8639                1.3913              0.1968

G gas                   -19065.7108               77.7842             11.0003
G solv                   -3053.2370               63.2630              8.9467

TOTAL                   -22118.9478               50.9582              7.2066

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1261.9126               14.1817              2.0056
EEL                     -11566.4419               71.5475             10.1183
EGB                      -2447.4831               55.0008              7.7783
ESURF                       64.5090                1.1105              0.1570

 
G gas                   -12828.3545               72.9394             10.3152
G solv                   -2382.9741               55.0120              7.7799

TOTAL                   -15211.3287               36.2055              5.1202

 
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                      -1661.8286               26.5442              3.7539
ESURF                       37.0493                0.6185              0.0875

 
G gas                    -5213.7811               37.3522              5.2824
G solv                   -1624.7794               26.5514              3.7549

TOTAL                    -6838.5604               25.6515              3.6277

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -64.2010                4.0841              0.5776
EEL                       -959.3742               34.9114              4.9372
EGB                        964.2108               32.9092              4.6541
ESURF                       -9.6943                0.3800              0.0537

 
DELTA G gas              -1023.5752               35.1495              4.9709
DELTA G solv               954.5164               32.9114              4.6544

 
 DELTA G binding =        -69.0588     +/-      6.5302                 0.9235
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

RESULT OF ALANINE SCANNING: (I21A MUTANT:) DELTA DELTA G binding =    1.9283  +/-    9.3470
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

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                      -3227.2145               64.4523              9.1149
ECAVITY                     68.4754                0.7567              0.1070

G gas                   -19064.5240             6061.1875            857.1813
G solv                   -3158.7391               64.4568              9.1156

TOTAL                    -7522.2032               51.2973              7.2545

 
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                      -2485.3559               54.5638              7.7165
ECAVITY                     47.5088                0.4610              0.0652

G gas                   -12825.2661             5345.3320            755.9441
G solv                   -2437.8471               54.5658              7.7168

TOTAL                    -5118.8075               38.9610              5.5099

 
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                      -1684.5802               28.2572              3.9962
ECAVITY                     28.1687                0.3939              0.0557

G gas                    -5213.7811             1395.1865            197.3092
G solv                   -1656.4114               28.2599              3.9966

TOTAL                    -2313.4381               24.9082              3.5225

 
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                        942.7215               33.8861              4.7922
ECAVITY                     -7.2022                0.3069              0.0434

DELTA G gas              -1025.4769             1237.6138            175.0250
DELTA G solv               935.5194               33.8875              4.7924

 
 
 DELTA G binding =        -89.9575     +/-      8.2480                 1.1664
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
I21A MUTANT:
POISSON BOLTZMANN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1855.4226               17.0765              2.4150
EEL                     -17210.2882               75.8866             10.7320
EPB                      -3229.1405               64.8100              9.1655
ECAVITY                     68.5521                0.7596              0.1074

G gas                   -19065.7108             6050.3755            855.6523
G solv                   -3160.5884               64.8144              9.1661

TOTAL                    -7520.1586               50.6710              7.1660

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1261.9126               14.1817              2.0056
EEL                     -11566.4419               71.5475             10.1183
EPB                      -2487.5603               54.6289              7.7257
ECAVITY                     47.6466                0.4632              0.0655

G gas                   -12828.3545             5320.1609            752.3844
G solv                   -2439.9137               54.6309              7.7260

TOTAL                    -5118.8820               38.4370              5.4358

 
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                      -1684.5802               28.2572              3.9962
ECAVITY                     28.1687                0.3939              0.0557

G gas                    -5213.7811             1395.1865            197.3092
G solv                   -1656.4114               28.2599              3.9966

TOTAL                    -2313.4381               24.9082              3.5225

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -64.2010                4.0841              0.5776
EEL                       -959.3742               34.9114              4.9372
EPB                        942.9999               34.0350              4.8133
ECAVITY                     -7.2632                0.3107              0.0439

DELTA G gas              -1023.5752             1235.4872            174.7243
DELTA G solv               935.7367               34.0364              4.8135

 
 
 DELTA G binding =        -87.8385     +/-      8.0665                 1.1408
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

RESULT OF ALANINE SCANNING: (I21A MUTANT:) DELTA DELTA G binding =    2.1190  +/-   11.5368
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

该结果文件与前面的同类文件的不同在于, 在每种方法之后会报告结合的ΔΔG, 表明突变对于复合物的结合自由能ΔG的相对影响. 特定的突变也在文件末尾输出. 在本例中, 我们将残基21从异亮氨酸突变为丙氨酸(即I21A).

3.5 分解单个残基或每对残基对Ras-Raf结合自由能的贡献(适用于amber11)

MMPBSA.py计算结合自由能所需的拓扑文件和mdcrd: ras-raf_top_mdcrd.tgz

我们将对3.1节获得的Ras-Raf体系进行结合自由能分解. AMBER支持两种类型的分解方法: 分解到每对残基和分解到单个残基. 单个残基分解方法通过计算单个残基与体系中其他所有残基的相互作用总和来计算单个残基的能量贡献. 成对分解计算体系中残基对之间的相互作用能. 下面将分别展示这两种类型的能量分解. 需要注意的是, 只有 MMPBSA.py正确识别出输入的残基后, 输出的每个残基的DELTA贡献才是有效的.

a. 单个残基的自由能分解

要进行分解, 必须在MMPBSA.py的输入文件中指定&decomp命名列表. 此外, 必须指定idecomp变量(没有默认值). 未能为此变量赋值将导致程序终止, 并给出错误信息. idecomp看指定4个允许的值, 其中两个用于单个残基的分解, 另外两个用于残基对分解. 值12决定了单个残基分解的方式: 值1将把1-4非键相互作用能(1-4 EEL1-4 VDW)加到内部势能项上; 值2将把1-4 EEL相互作用能加到静电势项上, 将1-4 VDW相互作用能加到范德华势能项上.

下面的MMPBSA.py输入文件mmpbsa_per_res_decomp.in使用PB和GB隐式溶剂模型进行单个残基分解:(注意:PB的非极性溶剂化能目前不可分解)

mmpbsa_per_res_decomp.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
Per-residue GB and PB decomposition
&general
   endframe=50, verbose=1,
/
&gb
  igb=5, saltcon=0.100,
/
&pb
  istrng=0.100,
/
&decomp
  idecomp=1, print_res="5; 30-40; 170-200"
  dec_verbose=1,
/

文件中dec_verbose变量用于指定能量分解输出文件的详细程度.

运行命令如下:

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -do FINAL_DECOMP_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd

注意, 上述命令可以用MMPBSA.py.MPI并行计算, 详情可参考3.3节.

输出文件per_res_output.tgz

最终结果文件FINAL_RESULTS_MMPBSA.dat, FINAL_DECOMP_MMPBSA.dat

FINAL_RESULTS_MMPBSA.dat
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
| Run on Thu May 20 14:55:43 EDT 2010

|Input file:
|--------------------------------------------------------------
|Per-residue GB and PB decomposition
|&general
|   endframe=50, verbose=1,
|/
|&gb
|  igb=5, saltcon=0.100,
|/
|&pb
|  istrng=0.100,
|/
|&decomp
|  idecomp=1, print_res="5; 30-40; 170-200"
|  dec_verbose=1,
|/
|--------------------------------------------------------------
|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               16.9979              2.4039
EEL                     -17200.7297               75.1734             10.6311
EGB                      -2918.9628               65.1000              9.2065
ESURF                       92.2138                0.9782              0.1383

G gas                   -19064.5240               77.0712             10.8995
G solv                   -2826.7490               65.1073              9.2076

TOTAL                   -21891.2730               52.3724              7.4066

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.1888               14.0912              1.9928
EEL                     -11557.0773               70.9920             10.0398
EGB                      -2314.8693               56.2410              7.9537
ESURF                       64.4513                0.6128              0.0867

G gas                   -12825.2661               72.3770             10.2356
G solv                   -2250.4181               56.2443              7.9542

TOTAL                   -15075.6842               36.8322              5.2089

 
Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.3090                9.3251              1.3188
EEL                      -4684.4720               35.7816              5.0603
EGB                      -1587.3051               26.8494              3.7971
ESURF                       38.5992                0.5158              0.0730

G gas                    -5213.7811               36.9768              5.2293
G solv                   -1548.7058               26.8544              3.7978

TOTAL                    -6762.4869               26.1943              3.7044

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.2966                4.2321              0.5985
EEL                       -959.1803               34.5681              4.8887
EGB                        983.2116               33.0175              4.6694
ESURF                      -10.8367                0.3832              0.0542

DELTA G gas              -1025.4769               34.8262              4.9252
DELTA G solv               972.3749               33.0197              4.6697

 DELTA G binding =        -53.1020     +/-      6.8437                 0.9678
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

POISSON BOLTZMANN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1863.7944               16.9979              2.4039
EEL                     -17200.7297               75.1734             10.6311
EPB                      -3216.4587               65.8638              9.3146
ECAVITY                     67.8762                0.7739              0.1094

G gas                   -19064.5240               77.0712             10.8995
G solv                   -3148.5825               65.8684              9.3152

TOTAL                   -22213.1066               51.7402              7.3172

 
Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.1888               14.0912              1.9928
EEL                     -11557.0773               70.9920             10.0398
EPB                      -2489.5955               55.9343              7.9103
ECAVITY                     47.1495                0.4689              0.0663

G gas                   -12825.2661               72.3770             10.2356
G solv                   -2442.4460               55.9363              7.9106

TOTAL                   -15267.7121               38.0243              5.3774

 
Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.3090                9.3251              1.3188
EEL                      -4684.4720               35.7816              5.0603
EPB                      -1673.2574               27.4055              3.8757
ECAVITY                     28.0328                0.4091              0.0579

G gas                    -5213.7811               36.9768              5.2293
G solv                   -1645.2246               27.4085              3.8761

TOTAL                    -6859.0057               24.7882              3.5056

 
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.2966                4.2321              0.5985
EEL                       -959.1803               34.5681              4.8887
EPB                        946.3942               34.1674              4.8320
ECAVITY                     -7.3062                0.2973              0.0420

DELTA G gas              -1025.4769               34.8262              4.9252
DELTA G solv               939.0881               34.1687              4.8322

 
 
 DELTA G binding =        -86.3888     +/-      8.1817                 1.1571
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

WARNINGS:
igb=5 should be used with either mbondi2 or bondi pbradii set. Yours are modified Bondi radii (mbondi)
FINAL_DECOMP_MMPBSA.dat
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
| Run on Thu May 20 14:55:43 EDT 2010
idecomp = 1: Decomposition per-residue adding 1-4 interactions added to Internal.
Energy Decomposition Analysis (All units kcal/mol): Generalized Born solvent

 
DELTAS:
Total Energy Decomposition:
Residue |  Location |       Internal      |    van der Waals    |    Electrostatic    |   Polar Solvation   |   Non-Polar Solv.   |       TOTAL
-------------------------------------------------------------------------------------------------------------------------------------------------------
LYS   5 | R LYS   5 |    0.000 +/-  4.870 |   -0.156 +/-  1.465 |   69.267 +/-  9.154 |  -67.061 +/-  9.601 |   -0.009 +/-  0.156 |    2.040 +/- 14.208
ASP  30 | R ASP  30 |    0.000 +/-  5.623 |   -0.065 +/-  0.961 |  -52.559 +/- 11.072 |   52.622 +/-  9.530 |    0.000 +/-  0.084 |   -0.003 +/- 15.684
GLU  31 | R GLU  31 |    0.000 +/-  5.174 |   -0.247 +/-  1.099 |  -79.946 +/- 10.550 |   80.630 +/-  9.693 |   -0.227 +/-  0.125 |    0.210 +/- 15.272
TYR  32 | R TYR  32 |    0.000 +/-  4.615 |   -0.290 +/-  1.515 |    0.639 +/-  4.431 |   -0.076 +/-  3.229 |   -0.012 +/-  0.175 |    0.261 +/-  7.327
ASP  33 | R ASP  33 |    0.000 +/-  4.464 |   -0.556 +/-  1.073 | -103.116 +/-  5.820 |  103.788 +/-  5.821 |   -0.459 +/-  0.094 |   -0.343 +/-  9.426
PRO  34 | R PRO  34 |    0.000 +/-  3.388 |   -1.829 +/-  0.869 |   -3.383 +/-  2.647 |    3.854 +/-  1.944 |   -0.308 +/-  0.130 |   -1.666 +/-  4.800
THR  35 | R THR  35 |    0.000 +/-  5.702 |   -1.829 +/-  1.049 |    0.376 +/-  4.365 |    0.947 +/-  2.028 |   -0.204 +/-  0.070 |   -0.709 +/-  7.536
ILE  36 | R ILE  36 |    0.000 +/-  4.650 |   -2.987 +/-  1.918 |    0.092 +/-  2.149 |    0.991 +/-  0.697 |   -0.377 +/-  0.072 |   -2.282 +/-  5.515
GLU  37 | R GLU  37 |    0.000 +/-  5.221 |   -1.627 +/-  1.388 | -126.728 +/-  6.441 |  120.528 +/-  4.686 |   -0.745 +/-  0.048 |   -8.573 +/-  9.624
ASP  38 | R ASP  38 |    0.000 +/-  3.750 |   -1.583 +/-  1.560 | -104.899 +/-  6.925 |   99.370 +/-  5.710 |   -0.254 +/-  0.037 |   -7.367 +/-  9.852
SER  39 | R SER  39 |    0.000 +/-  3.447 |   -2.184 +/-  1.086 |  -13.696 +/-  3.959 |    8.918 +/-  1.800 |   -0.504 +/-  0.035 |   -7.466 +/-  5.655
TYR  40 | R TYR  40 |    0.000 +/-  4.687 |   -4.403 +/-  1.682 |   -3.076 +/-  2.884 |    1.652 +/-  1.092 |   -0.366 +/-  0.042 |   -6.193 +/-  5.858
ARG 170 | L ARG   4 |    0.000 +/-  4.987 |   -0.094 +/-  1.646 |  -86.951 +/- 10.352 |   82.074 +/-  6.005 |   -0.147 +/-  0.073 |   -5.118 +/- 13.069
VAL 171 | L VAL   5 |    0.000 +/-  3.812 |   -0.183 +/-  1.390 |    2.128 +/-  2.460 |   -2.010 +/-  0.555 |    0.000 +/-  0.008 |   -0.065 +/-  4.778
PHE 172 | L PHE   6 |    0.000 +/-  4.289 |   -0.217 +/-  0.944 |    0.037 +/-  1.743 |    0.132 +/-  0.939 |    0.000 +/-  0.064 |   -0.048 +/-  4.818
LEU 173 | L LEU   7 |    0.000 +/-  4.907 |   -0.398 +/-  1.241 |   -0.940 +/-  3.050 |    1.683 +/-  1.446 |    0.000 +/-  0.022 |    0.345 +/-  6.084
PRO 174 | L PRO   8 |    0.000 +/-  3.433 |   -0.188 +/-  1.422 |    2.303 +/-  3.219 |   -2.589 +/-  1.289 |    0.000 +/-  0.051 |   -0.474 +/-  5.083
ASN 175 | L ASN   9 |    0.000 +/-  4.796 |   -1.671 +/-  1.017 |   -1.833 +/-  4.740 |    4.535 +/-  2.409 |   -0.354 +/-  0.119 |    0.678 +/-  7.233
LYS 176 | L LYS  10 |    0.000 +/-  4.403 |   -1.848 +/-  0.810 |  -33.879 +/-  7.269 |   36.798 +/-  6.704 |   -0.315 +/-  0.107 |    0.756 +/- 10.856
GLN 177 | L GLN  11 |    0.000 +/-  4.261 |   -3.791 +/-  1.560 |   -1.910 +/-  3.338 |    4.530 +/-  2.016 |   -0.359 +/-  0.050 |   -1.530 +/-  5.983
ARG 178 | L ARG  12 |    0.000 +/-  6.180 |   -2.462 +/-  1.321 |  -77.671 +/-  6.496 |   73.669 +/-  4.608 |   -0.386 +/-  0.076 |   -6.850 +/- 10.167
THR 179 | L THR  13 |    0.000 +/-  4.716 |   -1.277 +/-  1.200 |  -10.976 +/-  3.020 |    9.344 +/-  0.977 |   -0.158 +/-  0.031 |   -3.068 +/-  5.810
VAL 180 | L VAL  14 |    0.000 +/-  4.196 |   -3.837 +/-  1.389 |   -3.014 +/-  2.541 |    2.972 +/-  0.804 |   -0.501 +/-  0.041 |   -4.379 +/-  5.161
VAL 181 | L VAL  15 |    0.000 +/-  4.333 |   -1.791 +/-  1.119 |   -3.565 +/-  2.809 |    3.472 +/-  0.656 |   -0.155 +/-  0.055 |   -2.039 +/-  5.324
ASN 182 | L ASN  16 |    0.000 +/-  4.282 |   -1.978 +/-  0.859 |   -3.199 +/-  5.507 |    3.645 +/-  2.886 |   -0.369 +/-  0.085 |   -1.900 +/-  7.598
VAL 183 | L VAL  17 |    0.000 +/-  4.088 |   -0.187 +/-  1.149 |    1.057 +/-  4.557 |   -0.672 +/-  2.388 |    0.000 +/-  0.073 |    0.199 +/-  6.671
ARG 184 | L ARG  18 |    0.000 +/-  4.797 |   -0.183 +/-  1.450 |  -90.812 +/-  7.977 |   87.306 +/-  6.336 |   -0.335 +/-  0.109 |   -4.023 +/- 11.353
ASN 185 | L ASN  19 |    0.000 +/-  5.744 |   -0.018 +/-  0.966 |   -0.268 +/-  7.498 |    0.303 +/-  4.029 |    0.000 +/-  0.099 |    0.017 +/- 10.315
GLY 186 | L GLY  20 |    0.000 +/-  2.371 |   -0.008 +/-  0.701 |   -0.334 +/-  2.324 |    0.379 +/-  1.810 |    0.000 +/-  0.057 |    0.037 +/-  3.846
MET 187 | L MET  21 |    0.000 +/-  3.770 |   -0.156 +/-  1.254 |   -1.692 +/-  3.999 |    1.697 +/-  2.588 |   -0.031 +/-  0.089 |   -0.181 +/-  6.204
SER 188 | L SER  22 |    0.000 +/-  5.828 |   -0.013 +/-  1.008 |    2.808 +/-  4.910 |   -2.793 +/-  1.893 |    0.000 +/-  0.061 |    0.002 +/-  7.917
LEU 189 | L LEU  23 |    0.000 +/-  4.943 |   -0.021 +/-  1.312 |    1.683 +/-  2.195 |   -1.464 +/-  0.671 |    0.000 +/-  0.013 |    0.197 +/-  5.606
HIP 190 | L HIP  24 |    0.000 +/-  5.252 |   -0.024 +/-  1.131 |  -43.617 +/-  5.567 |   43.652 +/-  4.925 |    0.000 +/-  0.083 |    0.011 +/-  9.172
ASP 191 | L ASP  25 |    0.000 +/-  3.724 |   -0.058 +/-  0.723 |   62.413 +/-  8.165 |  -61.719 +/-  8.199 |    0.000 +/-  0.107 |    0.636 +/- 12.178
CYS 192 | L CYS  26 |    0.000 +/-  5.318 |   -0.098 +/-  1.398 |    1.937 +/-  3.894 |   -1.552 +/-  1.485 |    0.000 +/-  0.042 |    0.287 +/-  6.900
LEU 193 | L LEU  27 |    0.000 +/-  4.324 |   -0.108 +/-  1.390 |    0.884 +/-  2.119 |   -0.740 +/-  0.648 |    0.000 +/-  0.006 |    0.036 +/-  5.054

... 以下250行省略

FINAL_DECOMP_MMPBSA.dat输出文件包含有关每个残基与体系其余部分相互作用的信息, 这些信息分为: 内能(idecomp=1时, 势能项包括键, 角, 二面角和1-4相互作用), 范德华能(idecomp=2时, 包括VDW1-4 VDW), 静电能(idecomp=2时, 包括EEL1-4 EEL), 极性溶剂化能和非极性溶剂化能. 这个文件被分成如下几个部分:

复合物, 受体, 配体和DELTA(定义为复合物减去受体减去配体)中每个残基的分解能量会分别输出到各自的部分中. 此外, 每个残基的能量贡献会进一步分解为对骨架, 侧链及总能量贡献. Backbone是骨架原子与体系中所有其他原子之间的相互作用能, Sidechain是侧链原子与体系中所有其他原子之间的相互作用能. Total是残基中每个原子与系统中所有其他原子之间的相互作用能(因此是该残基的BackboneSidechain值的总和). 每个残基项会分解成上述的组成部分, 由相互作用的平均值+/-该项的标准偏差组成. DELTA部分额外包含一个名为Location的列, 列出复合物中特定残基的位置(R代表受体, L代表配体). 变量dec_verbose控制分解输出文件的详细程度(详细信息请参阅手册).

b. 残基对能量分解

注意: 根据我们的经验, 用PB隐式溶剂模型进行残基对能量分解分析需要非常长的时间才能完成. 下面的分析, 同时使用了GB和PB对50帧进行分析, 在9个处理器(9个独立的32位单核2.8 GHz Xeon处理器)上花费了61小时. 而GB分析只花了3分钟. 所以如果你选择进行PB残基对能量分解, 可能需要很长的模拟时间.

在本节中, 我们将修改输入文件以执行残基对能量分解. 与上面的单个残基能量分解的输入文件大部分相同, 残基对能量分解输入文件mmpbsa_pairwise_decomp.in如下所示:

mmpbsa_pairwise_decomp.in
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
Pairwise GB and PB decomposition
&general
   endframe=50, verbose=1,
/
&gb
  igb=5, saltcon=0.100,
/
&pb
  istrng=0.100,
/
&decomp
  idecomp=1, print_res="5; 30-40; 170-200"
  dec_verbose=0,
/

用于启动MMPBSA.py的命令与单个残基能量分解的命令相同. 但是, 定义用于残基对能量分解的&decomp命名列表中的print_res时需格外小心. 残基对能量分解需要计算的项数为n2, 其中n为print_res指定的残基数. 默认情况下, print_res对应于复合物中的每个残基, 对Ras-Raf会创建大约65 MB(超过450,000行)的分解输出文件. 而且, 由sander创建的mdout文件也将非常大(可能有几个GB, 具体取决于分析的帧数和残基对数), 并且解析所需的内存/时间需求变得相当巨大(即, 解析输出可能需要几分钟时间). 计算的残基对是print_res中指定的残基与print_res中指定的每个其他残基之间的残基对. 有关print_res变量的语法说明, 请参阅手册.

FINAL_DECOMP_MMPBSA.dat的部分输出结果如下图所示:

FINAL_DECOMP_MMPBSA.dat
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
| Run on Sun May 23 05:36:28 EDT 2010
idecomp = 3: Pairwise decomposition adding 1-4 interactions added to Internal.
Pairwise Energy Decomposition Analysis (All units kcal/mol): Generalized Born solvent

 
DELTAS:
Total Energy Decomposition:
Resid 1 | Resid 2 |       Internal      |    van der Waals    |    Electrostatic    |   Polar Solvation   |   Non-Polar Solv.   |       TOTAL
-----------------------------------------------------------------------------------------------------------------------------------------------------
LYS   5 | LYS   5 |    0.000 +/-  0.000 |    0.000 +/-  1.075 |    0.000 +/-  3.408 |    1.601 +/-  7.229 |    0.000 +/-  0.051 |    1.601 +/-  8.064
LYS   5 | ASP  30 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.341 |   -0.000 +/-  0.339 |    0.000 +/-  0.000 |   -0.000 +/-  0.480
LYS   5 | GLU  31 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.350 |    0.000 +/-  0.348 |    0.000 +/-  0.000 |    0.000 +/-  0.494
LYS   5 | TYR  32 |    0.000 +/-  0.000 |    0.000 +/-  0.001 |    0.000 +/-  0.071 |    0.000 +/-  0.070 |    0.000 +/-  0.000 |    0.000 +/-  0.099
LYS   5 | ASP  33 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.434 |    0.000 +/-  0.431 |    0.000 +/-  0.000 |    0.000 +/-  0.612
LYS   5 | PRO  34 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.064 |    0.000 +/-  0.064 |    0.000 +/-  0.000 |    0.000 +/-  0.090
LYS   5 | THR  35 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.189 |    0.008 +/-  0.183 |    0.000 +/-  0.000 |    0.008 +/-  0.262
LYS   5 | ILE  36 |    0.000 +/-  0.000 |    0.000 +/-  0.001 |    0.000 +/-  0.120 |    0.021 +/-  0.130 |    0.000 +/-  0.000 |    0.021 +/-  0.177

... 以下1800行省略

ARG 200 | LEU 197 |    0.000 +/-  0.000 |    0.000 +/-  0.423 |    0.000 +/-  0.779 |   -0.226 +/-  0.442 |    0.000 +/-  0.022 |   -0.226 +/-  0.991
ARG 200 | LYS 198 |    0.000 +/-  0.000 |    0.000 +/-  0.284 |    0.000 +/-  0.793 |    0.237 +/-  0.572 |    0.000 +/-  0.011 |    0.237 +/-  1.018
ARG 200 | VAL 199 |    0.000 +/-  0.000 |    0.000 +/-  0.291 |    0.000 +/-  0.832 |    1.460 +/-  0.587 |    0.000 +/-  0.019 |    1.460 +/-  1.059
ARG 200 | ARG 200 |    0.000 +/-  0.000 |    0.000 +/-  0.562 |    0.000 +/-  3.888 |   14.394 +/-  2.388 |   -0.000 +/-  0.035 |   14.394 +/-  4.598
idecomp = 3: Pairwise decomposition adding 1-4 interactions added to Internal.
Pairwise Energy Decomposition Analysis (All units kcal/mol): Poisson Boltzmann solvent

 
DELTAS:
Total Energy Decomposition:
Resid 1 | Resid 2 |       Internal      |    van der Waals    |    Electrostatic    |   Polar Solvation   |   Non-Polar Solv.   |       TOTAL
-----------------------------------------------------------------------------------------------------------------------------------------------------
LYS   5 | LYS   5 |    0.000 +/-  0.000 |    0.000 +/-  1.075 |    0.000 +/-  3.408 |    0.670 +/-  7.128 |    0.000 +/-  0.000 |    0.670 +/-  7.974
LYS   5 | ASP  30 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.341 |    0.007 +/-  0.334 |    0.000 +/-  0.000 |    0.007 +/-  0.477
LYS   5 | GLU  31 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.350 |    0.009 +/-  0.344 |    0.000 +/-  0.000 |    0.009 +/-  0.491
LYS   5 | TYR  32 |    0.000 +/-  0.000 |    0.000 +/-  0.001 |    0.000 +/-  0.071 |    0.002 +/-  0.069 |    0.000 +/-  0.000 |    0.002 +/-  0.098
LYS   5 | ASP  33 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.434 |    0.007 +/-  0.423 |    0.000 +/-  0.000 |    0.007 +/-  0.606
LYS   5 | PRO  34 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.064 |   -0.000 +/-  0.063 |    0.000 +/-  0.000 |   -0.000 +/-  0.090
LYS   5 | THR  35 |    0.000 +/-  0.000 |    0.000 +/-  0.000 |    0.000 +/-  0.189 |    0.024 +/-  0.175 |    0.000 +/-  0.000 |    0.024 +/-  0.257
LYS   5 | ILE  36 |    0.000 +/-  0.000 |    0.000 +/-  0.001 |    0.000 +/-  0.120 |    0.009 +/-  0.102 |    0.000 +/-  0.000 |    0.009 +/-  0.158
LYS   5 | GLU  37 |    0.000 +/-  0.000 |    0.000 +/-  0.008 |    0.000 +/-  1.649 |   -0.223 +/-  1.624 |    0.000 +/-  0.000 |   -0.223 +/-  2.315
LYS   5 | ASP  38 |    0.000 +/-  0.000 |    0.000 +/-  0.009 |    0.000 +/-  1.802 |   -0.191 +/-  1.383 |    0.000 +/-  0.000 |   -0.191 +/-  2.272

... 以下1800行省略

值得注意的是, FINAL_RESULTS_MMPBSA.dat输出的结果与单个残基能量分解的结果完全相同, 因为能量分解的方式不会影响总值. 因此, 这里忽略了文件以避免重重.

3.6 利用简正模式分析(Nmode)计算雌激素受体和雷洛昔芬复合物的熵

0. 简正模式计算的参数设置

请注意, 截至2010年5月, 简正计算是通过nab编译的nab程序mmpbsa_py_nabnmode完成的, 正常的安装过程可参考这里. 该程序必须存在于PATH$AMBERHOME/exe中才能运行计算. mmpbsa_py_nabnmodenmode之间的主要区别在于nab程序可以计算广义波恩(GB)溶剂中的简正模式, 因此另外又增加了两个输入变量(参见手册中的nmode_igbnmode_istrng). 这个实现还应该改进了早期版本脚本中出现的最小化收敛问题.

1. 构建起始拓扑和坐标文件, 模拟获得平衡体系

方法见3.3节.

2. 简正模式分析计算结合熵(normal mode)

我们将计算复合物, 受体和配体的简正模式, 并对结果进行平均以估计结合熵. 请注意, 复合物体系的熵贡献可以通过取消&general命名列表中最后一行的注释进, 利用AmberTools中的ptraj模块进行准简谐熵计算.

我们将利用MMPBSA.py的输入文件mmpbsa.in进行简正模式计算:

mmpbsa.in
1
2
3
4
5
6
7
8
Input file for running entropy calculations using NMode
&general
   endframe=50, keep_files=2,
/
&nmode
   nmstartframe=1, nmendframe=50,
   nminterval=5, nmode_igb=1, nmode_istrng=0.1,
/

文件中的endframe变量设置当用ptraj构建无溶剂体系的mdcrd文件时, 读入轨迹时在哪一帧停止. keep_files变量允许用户指定计算结束后删除哪些文件. &nmode命名列表用于定义用于简正模式计算的变量. 对要分析的无溶剂体系, nmstartframe变量定义简正模式分析开始时对应的帧. nmendframenminterval分别定义了简正模式分析的结束帧和帧间隔. 值得注意的是, nmstartframe/endframe/interval这些变量对应的”轨迹”是由&general命名列表中定义的startframe, endframeinterval变量提取出来的快照. 因此, 只有这些帧的一部分才能用于简正模式计算(最多是_MMPBSA_complex.mdcrd中的每帧).

运行命令如下:

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y *.mdcrd > progress.log

或者使用下面的命令将作业脚本提交到排队系统, 如PBS系统:

qsub nmode.job

作业脚本nmode.job看上去和bash类似, 内容如下:

nmode.job
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#!/bin/sh
#PBS -N nmode
#PBS -o nmode.out
#PBS -e nmode.err
#PBS -m abe
#PBS -M email@domain.edu
#PBS -q brute
#PBS -l nodes=1:surg:ppn=3
#PBS -l pmem=1450mb

cd $PBS_O_WORKDIR

mpirun -np 3 MMPBSA.py.MPI -O -i mmpbsa_nm.in -o FINAL_RESULTS_MMPBSA.dat -sp 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y 1err_prod.mdcrd > progress.log

计算进度会输出到文件progress_nm.log. 计算过程中的所有错误也会输出到这个文件. 最后, 计算完成后会显示每个步骤所用的时间. 如果串行计算10帧大约需要40个小时. MMPBSA.py.MPI可以用来并行(运行细节请参考3.3节). 计算时会根据线程数量平均分配帧数. 体系的大小对计算时间和所需内存(RAM)有显著影响. 段错误(segfaults)通常是由于系统中RAM不足造成的.

progress_nm.log
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
ptraj found! Using /share/local/lib/amber10/i686/bin/ptraj
sander found! Using /share/local/lib/amber10/i686/bin/sander (serial only!)
nmode found! Using /share/local/lib/amber10/i686/bin/nmode

Preparing trajectories with ptraj...
50 frames were read in and processed by ptraj for use in calculation.

Starting sander calls

 
Starting nmode calculations...
Timing:
Processing Trajectories With Ptraj:       0.240 min.
Total Harmonic nmode Calculation Time: 2363.906 min.
Output File Writing Time:                 0.018 min.

Total Time Taken:                      2364.165 min.

 
MMPBSA Finished. Thank you for using. Please send any bugs/suggestions/comments to mmpbsa.amber@gmail.com

输出文件Nmode_output.tgz

上述命令利用ptraj生成了三个非溶剂化mdcrd文件(复合物, 受体和配体), 其中包含计算熵变时分析过的坐标. 简正模式计算时利用每个快照的PDB文件, 因此从复合物, 受体和配体的无溶剂mdcrd文件生成了10个PDB文件. 如果entropy=1, 还会生成一个平均结构的PDB文件, 其中的结构(通过RMS)对齐到所有快照. 如果需要, 可以使用ptraj对这个结构进行准简谐熵计算.

最终结果文件FINAL_RESULTS_MMPBSA.dat:

FINAL_RESULTS_MMPBSA.dat
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
| Run on Thu Feb 11 12:44:26 EST 2010

|Input file:
|--------------------------------------------------------------
|Input file for running entropy calculations using NMode
|&general
|   endframe=50, keep_files=2,
|/
|&nmode
|   nmstartframe=1, nmendframe=50,
|   nminterval=5,
|/
|--------------------------------------------------------------
|Solvated complex topology file:  1err.solvated.prmtop
|Complex topology file:           complex.prmtop
|Receptor topology file:          receptor.prmtop
|Ligand topology file:            ligand.prmtop
|Initial mdcrd(s):                1err_prod.mdcrd
|
|Best guess for receptor mask:   ":1-240"
|Best guess for  ligand  mask:   ":241"
|Ligand residue name is "RAL"
|
|Calculations performed using 50 frames.
|Poisson Boltzmann calculations performed using internal PBSA solver in sander.
|
|All units are reported in kcal/mole.
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ENTROPY RESULTS (HARMONIC APPROXIMATION) CALCULATED WITH NMODE:

Complex:
Entropy Term          Average        Std. Dev.
-----------------------------------------------------------
Translational:        16.9389           0.0000
Rotational:           17.3953           0.0038
Vibrational:        2784.7967           2.3982
Total:              2819.1307           2.3964

 
Receptor:
Entropy Term          Average        Std. Dev.
-----------------------------------------------------------
Translational:        16.9233           0.0000
Rotational:           17.3911           0.0045
Vibrational:        2755.5693           3.0352
Total:              2789.8840           3.0342

 
Ligand:
Entropy Term          Average        Std. Dev.
-----------------------------------------------------------
Translational:        13.2972           0.0000
Rotational:           11.4991           0.0496
Vibrational:          33.7549           0.0442
Total:                58.5511           0.0058

 
DELTA S total=       -30.0192 +/-       3.4064

NOTE: All entropy results have units kcal/mol. (Temperature has already been multiplied in as 300. K)
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

该输出文件的开始部分包含有关计算的各种细节. 输出文件的其余部分包括平动, 转动, 振动和总熵贡献的平均值, 标准偏差和均值的标准误差. 在这些部分之后, ΔS以均值以及标准差和标准误差的形式给出.

通常我们会期望生物复合物的ΔS值为负. 这意味着蛋白质和配体结合形成复合物时可用微观状态数减少. 可用微观状态数的减少主要来自配体被限制在结合位点附近, 并且当配体与蛋白结合之后, 配体的运动受到了限制.

我们得到的ΔS值等于-30.02 kcal/mol, 这意味着在纯水中该蛋白-配体复合体的生成是一个熵不利的过程. 但是, 该结果不等于真实的结合自由能, 因为我们没有估计结合自由能中的(有利的) 焓贡献.

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 分子中原子编号的匹配
后一篇: AMBER高级教程A1(旧版):模拟含有非标准残基的溶剂化蛋白(简单版本)

访问人次(2017年1月27日起): | 最后更新: 2021-04-11 10:19:28 CST | 版权所有 © 2008 - 2021 Jerkwin