恒定PH值分子动力学模拟

自然界中往往存在一些对PH值敏感的蛋白，当我们想要探究这些蛋白在不同PH条件下构象变化时自然离不开MD的辅助。遗憾的是Gromacs暂不支持恒定PH条件下的模拟，因为这涉及到了蛋白表面残基质子化状态需要时刻改变。一种折衷的方法是：用propka等预测工具预测蛋白极性/可质子化残基的pka值（因为实际蛋白表面的pka值可能受到邻近残基的影响而发生偏移），然后根据想要模拟的PH环境，自主确定各个可质子化残基质子化状态（比如某个残基pKa=4，pH环境为7，那么这个残基自然是处于离去质子的解离状态）。然后修改对应pdb文件中的残基名，比如脱去质子的CYS改名为CYM（这个问题文末有提及）。这种方法存在的最大缺点在于：质子化状态一但设定好就会一成不变，这和实际情况显然是不符的。因此我们今天不介绍这种方法，而是介绍在Amber中的实现方法。

1. 准备结构

有这样一张表(下图所示)，罗列了常见质子化残基，包括ASP、GLU、HIS、CYS、LYS、TYR. 第三列数据对应它们各自的pKa值。第二列对应它们在恒定PH环境模拟时Amber识别的名字，如果想让Amber在模拟时对某个残基上面的质子化状态"多多关照",总得披上"马甲"好让Amber认识它吧。比如我想要模拟PH=7的环境，查看下表可以看出ASP、GLU、HIS在此条件下质子会发生解离，那么就把pdb文件中它们的名称改成对应的"马甲"。ASP→AS4、GLU→GL4、HIS→HIP

本次就以11个氨基酸长度的小肽test.pdb在pH为7的环境下模拟为例：

第一步，修改pdb文件

首先修改需要更换“马甲”的残基名，这里我们做了ASP→AS4、GLU→GL4修改，由于没有HIS残基，所以就不做考虑。此外该分子种有两个距离较靠近的CYS，我们想让他们之间成二硫键，所以修改做出修改：CYS→CYX。修改后的文件名为test_fix.pdb

#用pdb4amber不加氢，补全缺失原子，连接二硫键
pdb4amber -i test_fix.pdb --nohyd -o test_fix2.pdb
#不知道什么原因pdb4amber除氢总是失败，只能用下条命令自行除去了
grep -v '.............H' test_fix2.pdb > test_fix3.pdb

下图所示，我们想要连接的二硫键并没有形成，什么原因呢？这是因为只有硫原子间距小于2.5Å时pdb4amber才会自动连接。不过也没关系，后面我们还能在构建拓扑时强制将它们连上。

第二步，构建拓扑

#唤醒tleap程序
tleap
#加载constph力场（底层用的amber_ff10力场，同时于蛋白而言等价于ff99SB）
source leaprc.constph
source leaprc.water.tip3p  #为后面的抗衡离子添加力场
#加载修改好的pdb
mol = loadPDB test_fix3.pdb
#连接二硫键
#遵循语法bond <原子1> <原子2> [bondtype]
#bondtype分三类：“-”单键、“=”双键、“#”三键、“:”芳香键,若不指定则默认是单键
bond mol.4.SG mol.8.SG
#维持体系电中性
addions mol Cl- 0
addions mol Na+ 0
#保存拓扑、坐标
saveAmberParm mol test.prmtop test.inpcrd
#退出tleap程序
quit

注意，这里我们构建拓扑时没有添加水分子，这是因为本例要使用的是’隐式水‘模型。’隐式水‘不是真正的水分子，你可以把它理解为一个个数据点。启用’隐式水‘模型的参数为icnstph=1,需要将它添加到后文用到的以*.in为后缀的文件中。

第三步，准备恒定pH输入文件(cpin file)

#因为本例中我们只考虑了GLU、ASP的质子化状态改变，所以resname后面的参数改为它们对应的'马甲'——GL4、AS4
cpinutil.py -resnames GL4 AS4 -p test.prmtop -o test.cpin

2.模拟

第一步，EM

#em.in文件在分享文件中的em文件夹内
pmemd -O -i em/em.in -p test.prmtop -c test.inpcrd -o em/test_min.mdout -r em/test_min.rst -ref test.inpcrd -cpin test.cpin -inf em/em_mdinfo

第二步，体系加热

#heat.in文件在分享文件中的heat文件夹内
pmemd -O -i heat/heat.in -p test.prmtop -c em/test_min.rst -o heat/test_heat.mdout -r heat/test_heat.rst -ref em/test_min.rst -cpin test.cpin -inf heat/heat_mdinfo

第三步，预平衡

这一步*in文件里需要添加或调整参数solvph=, ntcnstph=。

solvph=设定体系pH值，ntcnstph= 设置每隔多少步监测一次所选残基质子化状态。要注意的是，amber每次只能监测一个残基的状态，也就是说所选“滴定残基”(即需要质子化状态改变的残基）数越多，于单个残基而言监测间隔越长。这就需要我们进行适当取舍。amber官方举的例子：当选取滴定残基数为10时，监测步长设为5得到的结果比较好。可见即便amber可以模拟恒定pH条件，但这也是有限制的。本例中我们选取的滴定残基数总共就两个，所以ntcnstph=可以设得比较大，我设定为25。（ps:后来才发现前面设置了三个“滴定残基”，不过影响应该不大）

#pr.in文件在分享文件中的pr文件夹内
pmemd -O -i pr/pr.in -p test.prmtop -c heat/test_heat.rst -o pr/test_pr.mdout -cpin test.cpin -r pr/test_pr.rst -cpout pr/test_pr.cpout -cprestrt pr/test_pr.cpin -inf pr/pr_mdinfo

本步输出的test_pr.cpin是test.cpin文件更新版。在预平衡中，“滴定残基”质子化状态发生了改变，所以需要将更新后的信息保存到新结果中，以供下一步成品模拟调用。

无关紧要的内容

恒定PH值分子动力学模拟相关推荐

深度学习DL蒙特卡洛法平衡态分子动力学模拟并计算苯酚键值
接上文<用反向传导进行分子动力学模拟并比较NN二甲基苯胺,N甲基苯胺,苯胺,硝基苯的定位效应>继续用神经网络模拟分子,这次计算苯酚苯酚的网络结构算出来的数值画成图 . 对比前面算出来的 ...
vasp 模拟退火_【转】vasp的分子动力学模拟 - 第一原理 - 小木虫 - 学术科研互动社区...
vasp做分子动力学的好处,由于vasp是近些年开发的比较成熟的软件,在做电子scf速度方面有较好的优势. 缺点:可选系综太少. 尽管如此,对于大多数有关分子动力学的任务还是可以胜任的. 主要使用的系 ...
分子动力学模拟手把手教你
如果你是AMBER的新用户或对一般的分子动力学模拟毫无了解可通过此教程入门. 教程简介这个教程旨在介绍如何使用Amber进行分子动力学模拟,并假设您以前没有使用过Amber. 它专门为想要了解如 ...
分子动力学模拟软件_分子模拟软件Discovery Studio教程(十三)：构建PLS模型（3D-QSAR）...
Discovery Studio™ (简称DS)是专业的生命科学分子模拟软件,DS目前的主要功能包括:蛋白质的表征(包括蛋白-蛋白相互作用).同源建模.分子力学计算和分子动力学模拟.基于结构药物设计工 ...
gromacs manual_GROMACS蛋白配体分子动力学模拟结果分析简要笔记
0. 引言本文以前文(https://zhuanlan.zhihu.com/p/149862369)为基础,对蛋白配体复合物分子模拟体系的结果进行一系列的粗浅分析,本文记述了简要的分析方法. 1 M ...
分子动力学模拟之SETTLE约束算法
Python微信订餐小程序课程视频 https://edu.csdn.net/course/detail/36074 Python实战量化交易理财系统 https://edu.csdn.net/cou ...
去除em斜体的方法_鱼缸水体pH值对观赏鱼的影响，以及偏高或偏低的调节方法...
大家好,我是小罗,专注于观赏鱼饲养技术,欢迎关注. 最近有一位饲养七彩神仙鱼的鱼友经常跟小罗聊天,他被一个容易出现误差的酸度计和降不下去的亚硝酸盐搞得焦头烂额.交流良久,期间我建议他把珊瑚砂慢慢换走, ...
Amber进行分子动力学模拟以及计算mmpbsa
使用amber计算mmpbsa记录 1.文件处理 2.蛋白与分子处理 (1) 前处理 (2) 生成crd与prm文件 3.分子动力学模拟 (1)能量最小化 (2)体系加热 (3)均匀密度 (4)全局平 ...
Vasp进行分子动力学模拟关键词解析及计算示例1
针对周期性体系的分子动力学模拟(包含计算所需输入文件) 通过vasp进行分子动力学的计算需要进行以下具体步骤: 1 选择一个足够大的晶胞制作成POSCAR,原子数越多越好(100个以上),一般k点较小 ...

恒定PH值分子动力学模拟

1. 准备结构

2.模拟

恒定PH值分子动力学模拟相关推荐

最新文章

热门文章