gromacs manual_GROMACS蛋白配体分子动力学模拟结果分析简要笔记

0. 引言

本文以前文(https://zhuanlan.zhihu.com/p/149862369)为基础，对蛋白配体复合物分子模拟体系的结果进行一系列的粗浅分析，本文记述了简要的分析方法。

1 MD 体系参数与周期性校正

主体模拟结束之后，需要检查其温度、体系密度等体系参数符合预设要求。

通过 gmx energy命令提取模拟过程中随时间变化的体系参数数据，然后利用此命令给出的各参数平均值等数据，以及根据数据绘制得到的参数随时间变化的折线图，即可以判断体系是否在预设的参数条件下运行。

gmx energy -f md_results_1.edr -o output.xvg
# then choose the items that u wanna output

同样还可以利用上述命令，提取体系potential、total-energy等参数信息进行分析。

轨迹文件中分子可能存在跨过周期性边界的情况，需要校正模拟体系的周期性。

首先将蛋白质和配体组合为一组，保存为 prolig_center.ndx，在其末尾添加如下内容：

[ center ]
500

上述内容设定了周期性校正的中心原子，之后的校正会将此原子始终置于盒子中心，编号 500 的原子为蛋白质上某一个原子。

之后利用下述命令校正周期性：

gmx trjconv -s md_results_1.tpr -f md_results_1.xtc -n prolig_center.ndx -o prolig_fit.xtc -pbc atom -center  
# 如果需要消除热运动：
# gmx trjconv -f md.xtc -s  md.tpr -fit rot+trans -o mdfit.xtc

运行命令之后先选择校正中心为设置的center，然后选择对整个体系进行校正即可。

2 配体径向分布

径向分布函数可以表征两种类型的分子之间的距离关系。分析配体到蛋白质表面的径向分布函数，可以获得稳定时期配体到蛋白质表面的距离等信息。

GROMACS中A类型粒子和B类型粒子之间的径向分布函数(RDF)被定义为：

其中<rho B(r)> 表示距离 A 粒子 r 距离，厚度为 bin 的壳层内B粒子的密度， <rho B(r)> 表示以A粒子为中心，半径的球内B粒子的平均密度。

在本研究中 r_max 为4 nm，也即盒子边长的一半，bin 设置为 0.01 nm。

通常径向分布函数用于分析连续相(比如说模拟体系中的水分子在蛋白质周围的径向分布)，但是也可以利用其对蛋白质周围配体的存在进行分析。根据上述公式可得，将其应用到配体分析时，函数会呈现出与连续相的径向分布函数不一样的特征。对于半径为 r 的球体，配体的平均密度是一个固定值，但在某些壳层内，因为没有配体分布，故而会使得函数值为零，而在某些配体密集分布的地方，会使得函数值较高；而对于连续相，则函数最终应该会收敛到1附近(也即在半径很大的时候，壳层平均密度与总体平均密度接近)，且函数峰值不会太高。

生成包含1ZIN的组ZIN.ndx，利用如下命令分析模拟稳定时期(模拟的最后15ns)配体1ZIN到蛋白质表面的径向分布函数：

gmx rdf -s prolig.tpr -f prolig_fit.xtc -o rdf_1_surf.xvg -n ZIN.ndx -ref protein -sel 1ZIN -bin 0.01 -surf mol -b 85000 -e 100000 -pbc yes

结果图举例如下：