上期四川大学李建宗博士为大家分享了》。
。
。
。

今天，，，将继续为大家介绍如何在豆听云计算平台上利用Amber完成蛋白-小分子体系的动力学模拟，，以及利用MMGBSA计算小分子与蛋白质（Abl和伊马替尼）之间的结合自由能。。。。

Amber是美国加州大学Peter Kollman等开发的一款著名的分子动力学模拟软件包。。
。Amber主要适用于蛋白质
，，，小分子和多糖等生物分子体系的模拟。。。

建议现在自己电脑上安装一个UCSF Chimera小程序，
，用于处理分子模拟所需的结构文件。。。。该软件学术免费，，而且具备和PyMOL类似的图形显示功能。。获取地址：

https://www.cgl.ucsf.edu/chimera/download.html

• 获取晶体结构结构复合物

本教程旨在利用Amber在豆听云平台上完成蛋白质Abl和靶向药物伊马替尼体系的动力学模拟，
，并且计算他们之间的结合自由能
。。。。幸运的是PDB数据库中包含了Abl和伊马替尼的复合物的晶体结构，，，，编号为6E4F。。。通过UCSF Chimera的fetch工具可以方便的获取晶体结构。。。在File—Fetch by ID中的PDB后输入编号6E4F即可下载复合物晶体结构，
，并将其保存为6E4F.pdb到自定义工作路径中
。。。

 Abl和伊马替尼的相互作用细节

• 处理蛋白质

 下载复合物晶体结构以后，，先对蛋白质进行处理，，，删除复合物中所有的除开蛋白质外的原子
。。。在Residue中选择所有非标准残基，
，
，然后通过Actions—Atoms/bonds—delete进行删除。。。。将蛋白质保存为protein.pdb文件。。。

• 处理配体

复合物中伊马替尼的名称为HRA，，，，重新打开6E4F.pdb文件，，，，然后选中HRA
，
，，在Select菜单中用Invert选中其余分子，，
，并且删除
，
，只保留伊马替尼的结构。。。如下所示。。。。

然后在Tools中找到Structure Editing分别对小分子进行加氢（Add H）和加电荷（Add charge）处理,电荷选择AM1-BCC。。。
。并将处理后的小分子保存为ligand.mol2文件
。
。

登录豆听云计算平台，
，，，将处理好的结构文件上传至豆听云计算平台，，，然后在提交作业中选择图形界面提交或者命令界面提交。
。。豆听云可快速帮你安装诸如Amber在内的模拟软件，，，免去了安装和环境配置等繁琐步骤，
，而且上传输入文件和下载计算结果文件都是图形操作，，，直观且方便。。

Amber支持GPU加速，，，在NVIDIA GPU 上的运行速度比仅使用 CPU 的系统快 15 倍，，，，因此，，
，在提交作业的时候注意选择GPU计算区
，，这里的GPU计算资源十分丰富
。
。。。对于不需要GPU的作业，，，，可以选择通用算区节省成本。。。
。

通过图形界面提交任务，，
，
，然后启动计算节点后
，，
，进入到我们刚才上传的文件夹下，，可以看到如下所示的操作系统界面
，
，，豆听云计算平台目前使用的是Centos 7系统，
，鼠标右键选择Open in Terminal然后就可以开始进行动力学模拟和结合自由能计算了。。。

本教程展示在豆听云计算平台上进行动力学模拟的具体细节，，计算量大的任务请参考本教程第04小结进行任务提交。。。。

利用Antechamber处理配体，，，并且利用parmchk2检查转换为Amber识别的参数结构文是否正确，
，命令如下：

antechamber -i ligand.mol2 -fi mol2 -o ligand.ante.mol2 -fo mol2 
parmchk2 -i ligand.ante.mol2 -f mol2 -o ligand.frcmod

其中-i指定输入文件，
，-fo或则-f指定输入文件格式，，
，，-o指定输出文件
，，-fo指定输出文件格式。
。。。检查是否有ligand.ante.mol2和ligand.frcmod文件生成。。。。利用Chimera打开ligand.ante.mol2，，，，可观测到如下参数化后的结构。
。。。

然后利用tleap模块生成配体、、、蛋白以及复合物的拓扑文件和模拟初始文件
。。。首先调用tleap

tleap

然后在tleap窗口中分别加载小分子、、蛋白质和溶剂所需要的力场参数

source leaprc.protein.ff14SBsource leaprc.gaffsource leaprc.water.tip3p

其中protein.14SB指定蛋白质力场，，gaff指定小分子力场，，，tip3p指定水分子模型。。。。

然后生成配体初始文件ligand.prmtop和ligand.inpcrd，
，，并且检测配体结构完整性。。

LIG = loadmol2 ligand.ante.mol2check LIGloadamberparams ligand.frcmodsaveoff LIG ligand.libsaveamberparm LIG ligand.prmtop ligand.inpcrd

生成蛋白初始文件receptor.prmtop和receptor.inpcrd

REC = loadpdb rec.pdbcheck RECsaveamberparm REC receptor.prmtopreceptor.inpcrd

生成蛋白-配体复合物初始文件com_solvated.prmtop和com_solvated.inpcrd

COM = combine {REC, LIG}saveamberparm COM com.prmtop com.inpcrdcharge COMsolvatebox COM TIP3PBOX 12.0saveamberparm COM com_solvated.prmtopcom_solvated.inpcrd

tleap会有类似提示

Building topology.Building atom parameters.Building bond parameters.Building angle parameters.Building proper torsion parameters.Building improper torsion parameters.total943 improper torsions appliedBuilding H-Bond parameters.Incorporating Non-Bonded adjustments.Not Marking per-residue atom chain types.Marking per-residue atom chain types.(Residues lacking connect0/connect1 -these don't have chain types marked:res total affected

完成上述操作后
，，，没有错误提示
，，，，则可退出tleap

quit

如果成功生成下面这些这些文件，
，则可以开始下一步的能量优化了。。。

• 能量最小化

动力学模拟起始需要对体系进行初始能量最小化，，，，以消除不合理原子间的碰撞和不规范的作用力，，命令如下
：

pmemd.cuda -O -i min.in -o min.out -p com_solvated.prmtop -c com_solvated.inpcrd -r min.rst -ref com_solvated.inpcrd

pmemd.cuda,调用pmemd的GPU版本,这里-i指定参数文件min.in, 对体系的平衡进行了限定，，参考官网指南具体min.in内容

minimise com&cntrlimin=1,maxcyc=2000,ncyc=1000,cut=8.0,ntb=1,ntc=2,ntf=2,ntpr=100,ntr=1, restraintmask=':1-288',restraint_wt=2.0/

imin=1 指定模拟步骤为能量最小化

maxcyc表明采用共轭梯度算法，，，且指定最大循环数，，

ncyc=1000表明采用最陡下降算法

ntpr=100表明保存模拟信息间隔

restraintmask指定限制氨基酸范围

restraint_wt=2.0表示限制力的大小

• 升温模拟

对体系进行升温，
，
，，温度从0K开始到300K结束，，
，命令如下
，，，并且检测是否有heat.rst和 heat.mdcrd文件生成。
。
。。

pmemd.cuda -O -i heat.in -o heat.out -p com_solvated.prmtop -c min.rst -r heat.rst -x heat.mdcrd -ref min.rst

heat.in参数内容

heat ras-raf&cntrlimin=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.rst和density.mdcrd文件生成 ,命令如下
：

pmemd.cuda -O -i density.in -o density.out -p com_solvated.prmtop -c heat.rst -r density.rst -x density.mdcrd -ref heat.rst

所用参数文件 density.in,

heat ras-raf&cntrlimin=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.rst和equil.mdcrd生成

pmemd.cuda -O -i equil.in -o equil.out -p com_solvated.prmtop -c density.rst -r equil.rst -x equil.mdcrd

equil.in参数文件内容如下
：

heat ras-raf&cntrlimin=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,/

利用process_mdout.pl脚本检查上述平衡是否达到合理水平，，，
，命令如下

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

利用xmgrace作图：

xmgrace summary.DENSITYxmgrace summary.TEMPxmgrace summary.ETOT

温度平衡结果

密度平衡结果

体系平衡结果

• 成品模拟
  

从上面分析可知体系基本上达到平衡，，
，因此可以进行成品模拟步骤，，，，命令如下

pmemd.cuda -O -i prod.in -o prod.out -p com_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrd

prod.in参数为：

heat ras-raf&cntrlimin=0,irest=1,ntx=5,nstlim=500000,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, /

同样的利用Chimera的MD movie工具打开模拟结果文件的prod1.mdcrd和com_solvated.prmtop，，，
，观看动力学模拟动画效果。。

Amber动力学模拟动态展示效果

通过MD movie--Analysis--Plot--RMSD工具绘制体系的RMSD图,分析结果表明整个体系在模拟的时间段内处于2.0 Å范围内波动，，
，是合理范围，，，，可用于后续结合自由能计算。。

动力学模拟体系RMSD结果

Amber集成了很多种结合自由能计算，，例如MMGBSA或者MMPBSA,这里我们用MMGBSA来计算Abl和伊马替尼的结合自由能，，Amber中计算结合自由能的模块为MMPBSA.py命令如下：

python $AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o result.dat -sp com_solvated.prmtop -cp com.prmtop -rp receptor.prmtop -lp ligand.prmtop -y prod*.mdcrd

mmpbsa.in参数文件指定了计算结合自由能类型，，具体参数内容如下：

Input file for running PB and GB&general   startframe=5, verbose=1,# entropy=1,/&gb  igb=2, saltcon=0.100/&pb  istrng=0.100,/

运行完毕会生成一系列记载能量信息的文件，，，最终结果包含在命令指定的dat后缀的文件中，，，，可提取数值进行作图分析。。。。通过计算可知Abl蛋白质与imatinib的结合自由能为-62.73 kcal/mol,具体能量项如下图所示
，，，，主要包括了范德华相互作用、、、
、静电相互作用、、、溶剂自由能等。。。。

Abl和imatinib结合自由能计算结果

上述教程是在豆听云计算平台上运行的细节，，在计算平台上提交作业，，请不要在管理节点逐条需要消耗计算资源的命令，，，管理节点只有2核4G内存的计算资源，
，进行动力学模拟必须将作业提交到运算节点进行。。

因此需要将上述命令写进脚本中
，，然后利用如下命令提交即可：

sbatch -N 1 -p g-v100-1 -c 12 amber.sh

N为节点的数量，
，
，这里输入的是1。。
。。p为选择的PARTITION，，，，这里使用的是V100卡（g-v100-1）
。
。查看作业运行情况及参数详细介绍可使用slurm命令进行查看
。。。。

amber.sh内容是上述教程中所有命令的脚本形式，
，内容为：

#!/bin/bashmodule add Anaconda3/2020.02source /public/software/.local/easybuild/software/amber/aber20/amber.shulimit -s unlimitedulimit -l unlimitedantechamber -i ligand.mol2 -fi mol2 -o ligand.ante.mol2 -fo mol2parmchk2 -i ligand.ante.mol2 -f mol2 -o ligand.frcmodpdb4amber -i protein.pdb -o rec.pdb --reduce --dry -y --nohydtleap -s -f tleap.inpmemd.cuda -O -i min.in -o min.out -p com_solvated.prmtop -c com_solvated.inpcrd -r min.rst -ref com_solvated.inpcrdpmemd.cuda -O -i heat.in -o heat.out -p com_solvated.prmtop -c min.rst -r heat.rst -x heat.mdcrd -ref min.rstpmemd.cuda -O -i density.in -o density.out -p com_solvated.prmtop -c heat.rst -r density.rst -x density.mdcrd -ref heat.rstpmemd.cuda -O -i equil.in -o equil.out -p com_solvated.prmtop -c density.rst -r equil.rst -x equil.mdcrdpmemd.cuda -O -i prod.in -o prod1.out -p com_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrdpython $AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o result.dat -sp com_solvated.prmtop -cp com.prmtop -rp receptor.prmtop -lp ligand.prmtop -y prod*.mdcrd

上述脚本中前面5行命令是告诉豆听云计算平台需要调用Amber程序，，后面的命令就是本教程所使用的具体命令了。。。

因此只需要将protein.pdb和ligand.mol2文件、、、所有in参数文件上传到豆听云计算平台，，，
，然后提交amber.sh作业脚本即可完成本教程所有内容。。。本教程参数文件获取地址
：

 https://pan.baidu.com/s/1n3HuK9dp9o_uTdBjLuPpdQ   提取码: 49x8

• 高性能计算资源，，极大的节省计算成本

• 支持海量的计算工具，，且开箱即用

• 计算任务进度实时追踪提示的贴心服务

•  详细耐心的技术咨询