注册 |登录

Jerkwin分享 http://blog.sciencenet.cn/u/Jerkwin

博文

AMBER高级教程A3：MM-PBSA

已有 30323 次阅读 2018-1-28 21:03 |系统分类:科研笔记

  • 原始文档: AMBER ADVANCED TUTORIAL 3: MM-PBSA

  • Perl Version By Ross Walker & Thomas Steinbrecher

  • Python Version By Dwight McGee, Bill Miller III, & Jason Swails

  • 2018-01-28 06:39:37 翻译: 袁媛; 修订: 李继存

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

AMBER高级教程A3: MM-PBSA

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

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

[A]aq+[B]aq⟺[A∗B∗]aq∗

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

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

显然, 从上图可以看出结合自由能 ΔGbind,solv 可以通过下式计算:

ΔGbind,solv0=ΔGbind,vacuum0+ΔGsolv,complex0−(ΔGsolv,ligand0+ΔGsolv,receptor0)

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

  • 对三种状态的溶剂化自由能是通过求解线性化的泊松-波尔兹曼(Poisson-Boltzman)方程或广义波恩(Generalized Born)方程(给出了静电相互作用对溶剂化自由能的贡献), 并添加表示疏水贡献的经验项来计算的.

ΔGsolv0=ΔGelectrostatic,ϵ=800−ΔGelectrostatic,ϵ=10+ΔGhydrophobic0

  • ΔGvacuum 是通过计算受体和配体之间的平均相互作用能, 并且在必要或需要时考虑受体与配体结合时熵的变化来得到的:

ΔGvacuum0=ΔEmolecula mechanics0−TΔSnorma mode analysis0

  • 熵变对于 ΔGvacuum 的贡献可以通过对三个物种(即两个蛋白和蛋白-蛋白复合物)进行简正模式分析获得. 但实际上如果只是需要比较具有相似熵的状态, 比如两个配体结合到同一个蛋白, 那么熵对于 ΔGvacuum 的贡献可以忽略. 这样做的原因在于简正模式分析计算非常费时, 而且一般存在较大的误差, 会在结果中引入显著的不确定性.

  • 受体与配体间的平均相互作用能通常是通过对不相关的快照(snapshots)系综进行计算得到的, 这些快照来自一段平衡的分子动力学模拟.

在本教程中, 我们将演示使用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.pdb和raf.pdb. 我们将使用这三个结构创建三个气相prmtop和inpcrd文件组用于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
setdefault 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=2或3(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脚本借助ptraj或cpptraj完成:

measure_equil_rmsd.ptraj
1
2
3
trajin equil.mdcrd.gz 1250 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

在平衡期结束时, 密度, 温度和总能量曲线都明显地收敛. 看起来开始稳定的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_96cd /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

从成品模拟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.7870.18-5602.0963.10-2102.2552.57
VDW             -984.9924.34-661.1820.33-256.0212.93
INT             5085.3350.223449.5738.651635.7629.42
GAS            -4556.4475.96-2813.7065.21-722.5253.50
PBSUR             65.091.0545.250.6427.240.46
PBCAL          -3223.6458.68-2490.8648.73-1671.2747.46
PBSOL          -3158.5558.26-2445.6248.45-1644.0347.31
PBELE         -11880.4234.25-8092.9629.34-3773.5217.30
PBTOT          -7714.9948.25-5259.3236.97-2366.5526.61
GBSUR             65.091.0545.250.6427.240.46
GB             -3407.8258.49-2631.8350.08-1731.0647.68
GBSOL          -3342.7458.15-2586.5849.83-1703.8247.55
GBELE         -12064.6026.94-8233.9223.57-3833.3113.40
GBTOT          -7899.1747.07-5400.2835.65-2426.3426.80

#                    DELTA
#          -----------------------
#                  MEAN        STD
#          =======================
ELE             -952.4344.10
VDW              -67.795.18
INT               -0.000.00
GAS            -1020.2244.58
PBSUR             -7.400.41
PBCAL            938.5042.51
PBSOL            931.0942.31
PBELE            -13.949.43
PBTOT            -89.137.94
GBSUR             -7.400.41
GB               955.0741.30
GBSOL            947.6641.10
GBELE              2.637.41
GBTOT            -72.566.40

该文件中不同项的含义如下:

  • ELE: 由分子力场(MM)计算的静电能

  • VDW: 来自MM的范德华贡献

  • INT: MM力场中由键, 角和二面角项引起的内能(在单轨迹方法中该项总是等于零)

  • GAS: 总的气相能量(ELE, VDW和INT的总和)

  • PBSUR/GBSUR: 由经验模型计算的溶剂化自由能中的非极性贡献

  • PBCAL/GB: 分别由PB或GB计算的溶剂化自由能中的静电贡献

  • PBSOL/GBSOL: 溶剂化的非极性和极性贡献之和

  • PBELE/GBELE: 静电溶剂化自由能和MM静电能的总和

  • PBTOT/GBTOT: 根据上述项得出的最终结合自由能的估计(kcal/mol)

值得注意的是, 通常情况下, 结合自由能的主要贡献来自ELE, PBCAL/GBCAL和VDW部分, 前两项大致相互抵消, 这可以通过查看PBELE/GBELE的值是否远小于ELE或PBCAL/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, pb和gb. 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 1112: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.794417.17042.4283
EEL                     -17200.729775.936610.7391
EGB                      -3249.651165.20759.2217
ESURF                       91.35651.39380.1971

G gas                   -19064.524077.853611.0102
G solv                   -3158.294665.22249.2238

TOTAL                   -22222.818651.02167.2155　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.188814.23422.0130
EEL                     -11557.077371.712710.1417
EGB                      -2532.066957.70038.1600
ESURF                       64.28431.11430.1576　
G gas                   -12825.266173.111810.3396
G solv                   -2467.782657.71108.1616

TOTAL                   -15293.048735.35274.9996　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.30909.41981.3322
EEL                      -4684.472036.14495.1117
EGB                      -1688.963126.53533.7527
ESURF                       37.04930.61850.0875　
G gas                    -5213.781137.35225.2824
G solv                   -1651.913826.54253.7537

TOTAL                    -6865.694925.88783.6611　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.29664.27510.6046
EEL                       -959.180334.91904.9383
EGB                        971.378933.04974.6739
ESURF                       -9.97700.37590.0532　
DELTA G gas              -1025.476935.17974.9752
DELTA G solv               961.401833.05184.6742　
 DELTA G binding =-64.0750+/-6.37290.9013--------------------------------------------------------------------------------------------------------------------------------------------------------------

POISSON BOLTZMANN:Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1863.794417.17042.4283
EEL                     -17200.729775.936610.7391
EPB                      -3207.716066.40239.3907
ECAVITY                     67.87620.78180.1106

G gas                   -19064.52406061.1875857.1813
G solv                   -3139.839966.40699.3914

TOTAL                    -7686.866052.54007.4303　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.188814.23422.0130
EEL                     -11557.077371.712710.1417
EPB                      -2483.724256.45517.9840
ECAVITY                     47.14950.47370.0670

G gas                   -12825.26615345.3320755.9441
G solv                   -2436.574756.45717.9842

TOTAL                    -5250.206038.51885.4474　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.30909.41981.3322
EEL                      -4684.472036.14495.1117
EPB                      -1670.416927.66943.9131
ECAVITY                     28.03280.41330.0584

G gas                    -5213.78111395.1865197.3092
G solv                   -1642.384127.67253.9135

TOTAL                    -2350.302025.11973.5525　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.29664.27510.6046
EEL                       -959.180334.91904.9383
EPB                        946.425134.51284.8808
ECAVITY                     -7.30620.30040.0425

DELTA G gas              -1025.47691237.6138175.0250
DELTA G solv               939.118934.51414.8810　　
 DELTA G binding =-86.3579+/-8.32641.1775--------------------------------------------------------------------------------------------------------------------------------------------------------------WARNINGS:
igb=2 should be used with mbondi2 pbradii set. Yours are modified Bondi radii (mbondi)

统计文件的开始部分包括日期/时间, 基于给定值和文件的任何警告, mmpbsa.in输入文件内容, 脚本所用的文件的名称, 分析的帧数以及使用了哪个PB求解器(如果使用的话). 统计文件的其余内容包括所有的平均能量, 标准偏差和平均值的标准误差, 先列出GB对应的值, 再列出PB对应的值. 每个部分之后, 会给出结合自由能ΔG及误差. 文件中不同项的含义如下:

  • VDWAALS: 来自分子力场(MM)的范德华贡献

  • EEL: 由MM力场计算的静电能

  • EPB/EGB: 分别由PB或GB计算的静电对溶剂化自由能的贡献

  • ECAVITY: 通过经验模型计算的非极性对溶剂化自由能的贡献

  • DELTA G binding: 根据上面各项计算出的最终结合自由能的估计值(kcal/mol)

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

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

总结合自由能-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=900mbcd$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 readinand 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 1419: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.794417.17042.4283
EEL                     -17200.729775.936610.7391
EGB                      -3249.651165.20759.2217
ESURF                       91.35651.39380.1971

G gas                   -19064.524077.853611.0102
G solv                   -3158.294665.22249.2238

TOTAL                   -22222.818651.02167.2155　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.188814.23422.0130
EEL                     -11557.077371.712710.1417
EGB                      -2532.066957.70038.1600
ESURF                       64.28431.11430.1576　
G gas                   -12825.266173.111810.3396
G solv                   -2467.782657.71108.1616

TOTAL                   -15293.048735.35274.9996　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.30909.41981.3322
EEL                      -4684.472036.14495.1117
EGB                      -1688.963126.53533.7527
ESURF                       37.04930.61850.0875　
G gas                    -5213.781137.35225.2824
G solv                   -1651.913826.54253.7537

TOTAL                    -6865.694925.88783.6611　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.29664.27510.6046
EEL                       -959.180334.91904.9383
EGB                        971.378933.04974.6739
ESURF                       -9.97700.37590.0532　
DELTA G gas              -1025.476935.17974.9752
DELTA G solv               961.401833.05184.6742　
 DELTA G binding =-64.0750+/-6.37290.9013--------------------------------------------------------------------------------------------------------------------------------------------------------------

POISSON BOLTZMANN:Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1863.794417.17042.4283
EEL                     -17200.729775.936610.7391
EPB                      -3207.716066.40239.3907
ECAVITY                     67.87620.78180.1106

G gas                   -19064.52406061.1875857.1813
G solv                   -3139.839966.40699.3914

TOTAL                    -7686.866052.54007.4303　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.188814.23422.0130
EEL                     -11557.077371.712710.1417
EPB                      -2483.724256.45517.9840
ECAVITY                     47.14950.47370.0670

G gas                   -12825.26615345.3320755.9441
G solv                   -2436.574756.45717.9842

TOTAL                    -5250.206038.51885.4474　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.30909.41981.3322
EEL                      -4684.472036.14495.1117
EPB                      -1670.416927.66943.9131
ECAVITY                     28.03280.41330.0584

G gas                    -5213.78111395.1865197.3092
G solv                   -1642.384127.67253.9135

TOTAL                    -2350.302025.11973.5525　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.29664.27510.6046
EEL                       -959.180334.91904.9383
EPB                        946.425134.51284.8808
ECAVITY                     -7.30620.30040.0425

DELTA G gas              -1025.47691237.6138175.0250
DELTA G solv               939.118934.51414.8810　　
 DELTA G binding =-86.3579+/-8.32641.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 1112: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.380120.30212.8712
EEL                     -16938.645085.763112.1287
EGB                      -3507.008667.78399.5861
ESURF                       97.54481.33010.1881

G gas                   -18952.025188.133312.4639
G solv                   -3409.463967.79699.5879

TOTAL                   -22361.488947.19826.6748　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.227219.23112.7197
EEL                     -16895.035485.579712.1028
EGB                      -3528.727668.35859.6673
ESURF                      101.26131.30710.1849

G gas                   -18850.262687.713812.4046
G solv                   -3427.466368.37109.6691

TOTAL                   -22277.728848.10576.8032　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.85952.05160.2901
EEL                         -5.57962.03330.2876
EGB                        -28.48630.60400.0854
ESURF                        4.43260.04620.0065　
G gas                       -7.43912.88850.4085
G solv                     -24.05380.60580.0857

TOTAL                      -31.49295.07480.7177　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.29342.92650.4139
EEL                        -38.03003.21140.4542
EGB                         50.20532.58690.3658
ESURF                       -8.14910.25890.0366　
DELTA G gas                -94.32344.34490.6145
DELTA G solv                42.05622.59990.3677　
 DELTA G binding =-52.2672+/-2.45680.3475--------------------------------------------------------------------------------------------------------------------------------------------------------------

POISSON BOLTZMANN:Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -2013.380120.30212.8712
EEL                     -16938.645085.763112.1287
EPB                      -3329.170867.03549.4802
ECAVITY                     68.26560.51950.0735

G gas                   -18952.02517767.48371098.4881
G solv                   -3260.905267.03749.4805

TOTAL                    -5265.083149.04266.9357　Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.227219.23112.7197
EEL                     -16895.035485.579712.1028
EPB                      -3355.474667.32999.5219
ECAVITY                     70.11840.52850.0747

G gas                   -18850.26267693.71631088.0558
G solv                   -3285.356267.33209.5222

TOTAL                    -5279.450950.40677.1286　Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.85952.05160.2901
EEL                         -5.57962.03330.2876
EPB                        -31.33640.69530.0983
ECAVITY                      3.18960.02880.0041

G gas                       -7.43918.34341.1799
G solv                     -28.14680.69590.0984

TOTAL                       56.09345.04760.7138　
Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.29342.92650.4139
EEL                        -38.03003.21140.4542
EPB                         57.64023.06420.4333
ECAVITY                     -5.04230.16830.0238

DELTA G gas                -94.323418.87782.6697
DELTA G solv                52.59783.06880.4340　　
 DELTA G binding =-41.7256+/-2.96180.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.pdb和ras.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.pdb和ras_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

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

b. 残基对能量分解

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


◆本文地址: https://jerkwin.github.io/2018/01/28/AMBER高级教程A3-MM-PBSA/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆



https://blog.sciencenet.cn/blog-548663-1097321.html

上一篇：分子中原子编号的匹配
下一篇：AMBER高级教程A1：模拟含有非标准残基的溶剂化蛋白(简单版本)
收藏 IP: 72.221.39.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-8 15:29

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部