MMPB/GBSA结合自由能计算以残基贡献度分析

1.MMPB/GBSA结合自由能计算简介

在运行完MD之后,毫无疑问需要对其进行分析.其中结合自由能计算就是必不可少的一部分:

配体和受体的游离态到结合态之间存在一个能量差,用这个能量差可以衡量两者结合的亲密度.下图是结合MMGBSA的原理,其将结合自由能分为了分子力学能Emm,(范德华力Evdw,静电相互作用Eele),溶剂化自由能(Gsolvation)

2.输入文件准备

为了进行结合自由能计算 我们分别需要

1.复合物拓扑文件 - ->ligand_1_complex.prmtop
2.动力学模拟结果mdcrd文件 - ->07_prod.nc
3.输入参数文件–>mmpbsa.in

2.1 拓扑文件处理

在命令行执行如下命令:

ante-MMPBSA.py -p ligand_1_complex.prmtop -c ligand_1_complex-now.prmtop -r 4l6s_docked.prmtop -l ligand_1.prmtop -s ":WAT,Na+,Cl-" -n ":352"
#-p 输入符合物拓扑文件
#-c 处理后的复合物拓扑文件输出
#-r 处理后的靶点拓扑文件输出
#-l 处理后的配体拓扑文件输出
#-s 需要提出的元素 本例剔除水,钠离子,氯离子
#-n 配体所处的残基序列位置

2.2 mmpbsa.in的编译

&general             
startframe=1901, endframe=2000, interval=1,   
verbose = 1
ligand_mask = ':352'
receptor_mask = ':1-351'
/              
&gb              
igb=8, saltcon=0.100,  
/              
&pb              
inp=0, radiopt=0,           
istrng=0.100, 
/              
&decomp             
idecomp=1, print_res="44,46,51,99,102-113,200-205,208-212,215-222,224,228-249,258-260,290,324,326-331,334,337,339,352",
dec_verbose=1,
/
    '''
# start/end frame 起始结束帧 interval 步长(间隔10帧取一次) 一般帧数/步长为1000就够
# verbose =1, 输出能量信息。（=0：不输出，=2输出每步能量信息)
# &gb 采用 Generalized Born 方法，盐浓度 单位是摩尔浓度（Morality）
# &pb 离子强度,单位为摩尔浓度（Morality）
# idecomp=1 & 2 时输出结果为单个残基的结合自由能，其中=1将1-4非键相互作用能（1-4 EEL & 1-4 VDW）归为内部势能，而=2则将1-4非键相互作用各自归为静电势能和范德华势能上
# idecomp=3 & 4 时输出结果为残基对的结合自由能。
# print_res 为需要计算残基贡献度的残基,一般选择ligand 5A范围之内的残基以及其本身就行.

3.运行程序

执行下列代码

mpirun -np {mpi_num} MMPBSA.py.MPI -O -i mmgbsa.in -o 11_FINAL_RESULTS_MMPBSA.dat -do 11_FINAL_DECOMP_MMPBSA.csv -sp ligand_1_complex.prmtop -cp ligand_1_complex-now.prmtop -rp 4l6s.docked.prmtop -lp ligand_1.prmtop -y ./07_prod.nc > 11_mmgbsa.log'
# {mpi_num} 并行cpu核数
# -i 输入参数文件 -o 输出总能量变化文件 -do 残基贡献度输出文件 
# cp/rp/lp 各类输入拓扑文件
# y mdcrd 输入文件
# > 日志文件输出

4.绘图

运行完毕之后 取11_FINAL_DECOMP_MMPBSA.dat结果选择能量贡献度靠前的残基绘制能量贡献图(如下):