分子动力学模拟

一、准备和程序运行

1.1预处理受体蛋白、配体分子

Pymol打开分子对接结果，依次进行以下结果

//将受体蛋白另存为pdb，AutoDock中去H，再另存为pdb

protein.pdb

//将配体分子另存为pdb，AutoDock去除氢离子，然后openBabel转化为mol2

mol.mol2

1.2为小分子生成top文件和gro文件

打开sobtop，把mol.mol2拖入sobtpo

//选择1，生成top文件

//选择3，尽可能使用GAFF力场

//选择0，进入下一步

//选择4，如果可能，预先构建成键参数，任意猜测缺少的选项

//回车，生成的top文件生成在sobtop软件根目录下

//回车，生成的itp位置限制文件在sobtop软件根目录下

//选择2，生成gro文件

//回车，生成的gro文件在sobtop软件根目录下

//回车，退出sobtop软件

//将sobtop软件中的

mol.gro

mol.itp

mol.top

//三个文件剪切到实验工作环境目录下

1.3产生蛋白质连带一个离子的拓扑文件

gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top

//选择Amber99sb-ildn力场和TIP3P水(注意，如果静电荷不为零: Total charge in system x.000 e 则需要添加抗衡离子)

//得到的topol.top，protein.gro，posre.itp

//Protein.gro对应蛋白，topol.itp对应离子,posre.itp开头的文件对应限制势

1.4合并gro文件，另存为complex.gro文件

1.4.1//将mol.gro加入到pdb2gmx产生的protein.gro的末尾，并将第二行的原子数改为（蛋白质原子数+小分子原子数）,建议作图确认结构合理性，注意文件中的空格和回车问题

//注意，protein.gro文件末尾的三个小数是晶格的坐标，不要删除或大幅修改

1.4.2//蛋白质的限制势itp文件在pdb2gmx的时候已经产生，但小分子的没有，genrestr是对输入的结构产生坐标或距离限制势itp文件的工具，接下来运行命令，进行限制势的产生

gmx genrestr -f mol.gro -o posre_mol.itp

//选择组的时候选择0，system默认的位置限制势常数是1000kJ/mol/nm2，已经足够大

1.4.3将下列语句插入到mol.itp文件的末尾，注意，复制时连“#”井号一同复制，最好在末尾添加之前空一行，方便检查文件错误

#ifdef POSRES

#include “posre_mol.itp”

#endif

//这样当mdp中使用define = -DPOSRES的时候配体的位置也会被限制了

1.4.4把配体的itp文件引入整体的拓扑文件topol.top，把分子数也设置好

//注意，在引入的时候需要将小分子的mol.itp文件引入到蛋白质链之前，因为mol.itp最开头定义了[atomtypes]因此，这个itp要最优先被引入

//将下列语句插入到引入蛋白质的itp文件引入之前；

#include “mol.itp”

//并在末尾的[molecules]中引入mol 1（注意空格），将topol.top的格式与complex.gro中分子出现的顺序对应

//即topol.top中应该是类似这样的顺序：

; Include forcefield parameters

#include “amber99sb-ildn.ff/forcefield.itp”

; Include chain topologies

#include “mol.itp”

#include “topol_Protein_chain_A.itp”

; Include water topology

#include “amber99sb-ildn.ff/tip3p.itp”

#ifdef POSRES_WATER

; Position restraint for each water oxygen

[ position_restraints ]

;  i funct       fcx        fcy        fcz

   1    1       1000       1000       1000

#endif

; Include topology for ions

#include “amber99sb-ildn.ff/ions.itp”

[ system ]

; Name

CATIONIC TRYPSIN in water

[ molecules ]

; Compound        #mols

Protein_chain_A     1

mol 1

//注意，末尾的mol和1之间有一个空格，且下一步添加溶剂后有可能会缺失回车，需要纠正

1.5设定盒子，加水，加离子，能量极小化

1.5.1//设定盒子

gmx editconf -f complex.gro -o complex_box.gro -d 0.8 -bt cubic

//设定的盒子是立方盒子，可能会增加计算量，但是不容易产生边界相互作用的问题

1.5.2//加水

gmx solvate -cp complex_box.gro -o complex_SOL.gro -p topol.top

//注意这一步加水后有可能topol.top文件最后一行的SOL可能会串行，需要手动添加回车，避免其与mol 1连在同一行，容易在后续处理中报错

1.5.3//将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录

//产生临时tpr文件

gmx grompp -f em.mdp -c complex_SOL.gro -p topol.top -o em.tpr -maxwarn 1

//这里如果警告较多可以将maxwarn的数值改大一些

//出现以下报错：

将

和

; Include ligand topology

#include “mol.itp”

一起加到topol.top的立场文件“; Include forcefield parameters

#include “amber99sb-ildn.ff/forcefield.itp””之后，即

1.5.4//加离子，使得整个体系变为电中性（与第三步中的电荷数相对应）

gmx genion -s em.tpr -p topol.top -o system.gro -neutral

//这里选择分组时选择SOL对应的分组

//产生的带有离子且电中性的体系为system.gro

1.5.5//能量极小化

gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

1.6限制性动力学100ps

gmx grompp -f pr.mdp -c em.gro -p topol.top -r em.gro -o pr.tpr

如果出现以下报错：

则在代码最后加上-maxwarn 10，忽略警告

gmx mdrun -v -deffnm pr

//注意，对于复杂体系如果用常规2fs步长可能模拟刚开始就会崩溃，因此此处做限制性动力学期间，步长用较小的1fs以求稳妥。

1.7产生索引文件

gmx make_ndx -f pr.gro

依次输入：

x1|x2|x3（定义“蛋白-离子-配体”组，新的组号接在原有的组号之前，组号为R）

（注释：|为右大括号右侧的按键，用“shift”+“、”便可输入该符号）

!R(把其他部分定义为一个组，新的组号为R+1)

name R protein_lig(改名为蛋白_配体)

name R+1 envir（改名为环境）

q

//得到的index。ndx里的组名就和模板文件里面的md.mdp中的组对应了

（单纯蛋白的话运行命令后直接输入q退出即可）

1.8常规动力学1ns

gmx grompp -f md.mdp -c pr.gro -p topol.top -o md.tpr -n index.ndx -maxwarn 10

//如果这里警告多，可以将maxwarn后的数值改大一些

1.9正式运行

//如果需要跑较长时间的动力学，需要修改md.mdp文件，默认是100 ns

gmx mdrun -v -deffnm md

//由于配体、离子与蛋白的相互作用较为明显，所以这里任务之中将他们一起作为一个控温组和同一个消除屏东转动的组

//-v 为可视化

暂停：control+c

续跑：gmx mdrun -s md.tpr -cpi md.cpt -deffnm md -v   

二、结果分析

2.1 蛋白定位至中心

在计算RMSD以及回旋半径过程，大分子（蛋白）以及小分子会在盒子中不断变换位置和构象，甚至会出现跨盒子的现象；所以在整个计算开始之前，我们需要将整个体系定位到盒子中心；使用gmx trjconv命令完成操作：

gmx trjconv -s md.tpr -f md.xtc -o md_0_1_noPBC.xtc -pbc mol -center

选择以蛋白质为中心（Group 1）；输出蛋白质（Group 1）

gmx trjconv -s md.tpr -f md_0_1_noPBC.xtc -o center.pdb -dt 1000

选择蛋白质（Group 1）

运行后得到md_0_1_noPBC.xtc文件，这是居中后的体系文件，用于后续分析，这里我们可以可视化一下刚刚生成的center.pdb文件（pymol），这里记录蛋白运动1000帧的动画，时长00:10

2.2 去除平衡转动

上面的动画我们发现整个蛋白复合物运动幅度都比较大，所以下一步需要去除体系中的平衡转动；

gmx trjconv -s md.tpr -f md_0_1_noPBC.xtc -o fit.xtc -fit rot+trans

选择蛋白骨架，（Group 4），选择蛋白（Group 1）

gmx trjconv -s md.tpr -f fit.xtc -o fit.pdb -dt 1000

选择蛋白（Group 1）

去除平衡转动后，我们再继续可视化一下fit.pdb，并使用fit.xtc进行后续的计算时长00:11

2.3 计算RMSD、RMSF、回旋半径

使用fit.xtc作为输入文件，进行下一步的计算：

2.3.1计算RMSD

gmx rms -s md.tpr -f fit.xtc -o rmsd.xvg -tu ns

选择蛋白骨架（Group 4），随后选择蛋白（Group 1）

dit xvg_show -f rmsd.xvg

2.3.2接下来计算RMSF，使用gmx rmsf命令完成：

计算RMSF

gmx rmsf -s md.tpr -f fit.xtc -o rmsf.xvg -res

选择Group 1，protein

dit xvg_show -f rmsf.xvg

(在其后加上-res)便可以以氨基酸残基为横坐标

2.3.3计算回旋半径，使用gmx gyrate进行计算：

gmx gyrate -s md.tpr -f fit.xtc -o gyrate.xvg

选择Group 1，protein

dit xvg_show -f gyrate.xvg

2.3.4 计算氢键

gmx hbond -s md.tpr -f fit.xtc -n index.ndx -num hbond_num.xvg -dt 100 -life hbond_life.xvg -ac hbond_ac.xvg

选择Group 1，protein，然后选择Group 19，Mol

gmx hbond -f fit.xtc -s md.tpr -num hb_num.xvg -ac hbond_ac.xvg -dist hbdist.xvg -hx hbhelix.xvg

三、结果可视化

用pymol处理统一轨道中蛋白和小分子。在计算RMSD以及回旋半径过程，大分子（蛋白）以及小分子会在盒子中不断变换位置和构象，甚至会出现跨盒子的现象；所以在整个计算开始之前，我们需要将整个体系定位到盒子中心；使用gmx trjconv命令完成操作：

3.1对于蛋白和小分子：

gmx trjconv -s md.tpr -f md.xtc -o P_L_noPBC.xtc -pbc mol -center

选择蛋白质骨架为中心（Group 4）；输出除了水之外的原子组（Group 17/14）

运行后得到P_L_noPBC.xtc文件，这是居中后的除了水之外的体系文件

gmx trjconv -s md.tpr -f P_L_noPBC.xtc -o P_L_protein.pdb -dt 1000

选择蛋白质（Group 1），得到蛋白质的pdb文件

gmx trjconv -s md.tpr -f P_L_noPBC.xtc -o P_L_ligend.pdb -dt 1000

选择小分子Group 13），得到小分子的pdb文件

先后将两个文件导入到PyMol，随后便可作图。

3.2对于单纯蛋白

gmx trjconv -s md.tpr -f md.xtc -o P_noPBC.xtc -pbc mol -center

选择蛋白质骨架为中心（Group 4）；输出除了水之外的原子组（Group 17/14）

gmx trjconv -s md.tpr -f P_noPBC.xtc -o P_protein.pdb -dt 1000

选择蛋白质（Group 1）1，得到蛋白质的pdb文件

先后将文件导入到PyMol，随后便可作图