最近一直在忙着我的论文修改，现在总算要告一段落了，马上要开启新的篇章，其中一个部分就是膜的体验，先准备学习以下gromacs李继存老师翻译的教程，这里记录一下学习的过程：

DPPC膜中的KALP15：

教程使用的文章是06年的一篇3.9分的文章，准备学习完以后参考一下他的分析步骤：文献查看

第一步：准备拓扑

多肽如何准备的参考资料并未详细给出，由于我也是第一次学，自己也先跳过，你可以在文章底部中文教程下载也可以KALP-15_princ.pdb

1.pdb进行处理

1
gmx pdb2gmx -f KALP-15_princ.pdb -o KALP-15_processed.gro -ignh -ter -water spc

ignh选项指示gmx pdb2gmx忽略输出中的H原子. 默认情况下, xLeap会给出全原子结构(由于AMBER力场使用显式的氢). 在AMBER的命名约定中, 这些H原子的命名可能与GROMOS96力场中的不同. 如果让gmx pdb2gmx忽略输入的所有H原子, 它只会添加需要的那些.所以如果在蛋白配体复合物模拟中选用amber的话那么也不应该加-ignh选项，否则可能会干掉一些氢

选择13：GROMOS96 53A6力场
两端不需要封端，所以均选择：None

2.修改拓扑

下载膜：连接地址

dppc128.pdb - 含128个DPPC脂分子的双脂膜
dppc.itp - DPPC的分子类型定义
lipid.itp - Berger脂分子参数

你也可以在我这里dppc下载好了的

将gromos53a6.ff目录中文件复制出来，也可以在我这里下载已经打包好了的，gromos53a6_lipid.ff

接下来, 将lipid.itp中的[ atomtypes ], [ nonbond_params ]和[ pairtypes ]部分复制粘贴到ffnonbonded.itp中相应的行. 你可以发现lipid.itp中的[ atomtypes ]节缺少原子编号(at.num列), 因此添加上编号. 修改后的行应该是:

1
2
3
4
5
6
7
8
9
10
11
12
LO    8    15.9994      0.000     A  2.36400e-03 1.59000e-06 ;carbonyl O, OPLS
LOM    8    15.9994      0.000     A  2.36400e-03 1.59000e-06 ;carboxyl O, OPLS
LNL    7    14.0067      0.000     A  3.35300e-03 3.95100e-06 ;Nitrogen, OPLS
 LC    6    12.0110      0.000     A  4.88800e-03 1.35900e-05 ;Carbonyl C, OPLS
LH1    6    13.0190      0.000     A  4.03100e-03 1.21400e-05 ;CH1, OPLS
LH2    6    14.0270      0.000     A  7.00200e-03 2.48300e-05 ;CH2, OPLS
 LP   15    30.9738      0.000     A  9.16000e-03 2.50700e-05 ;phosphor, OPLS
LOS    8    15.9994      0.000     A  2.56300e-03 1.86800e-06 ;ester oxygen, OPLS
LP2    6    14.0270      0.000     A  5.87400e-03 2.26500e-05 ;RB CH2, Bergers LJ
LP3    6    15.0350      0.000     A  8.77700e-03 3.38500e-05 ;RB CH3, Bergers LJ
LC3    6    15.0350      0.000     A  9.35700e-03 3.60900e-05 ;CH3, OPLS
LC2    6    14.0270      0.000     A  5.94700e-03 1.79000e-05 ;CH2, OPLS

在[ nonbond_params ]节中, 你可以看到一行;; parameters for lipid-GROMOS interactions., 删除这一行以及所有后续行, 因为这些非键组合特定地用于已经废弃的ffgmx力场. 移除这些行后, 会根据GROMOS96 53A6力场的标准组合规则生成蛋白质和脂分子之间的相互作用参数. 文件中也有涉及HW原子类型的非键相互作用, 由于它们都是零, 你也可以删除这些行, 否则请将HW重命名为H, 以便与GROMOS96 53A6的命名约定兼容. 如果你不进行重命名或移除这些行, 后面执行gmx grompp时会导致致命错误.

追加[ dihedraltypes ]的内容到ffbonded.itp中相应的节. 这些行看起来差异很大, 但无关紧要. 它们是Ryckaert-Bellemans二面角, 与默认GROMOS96 53A6力场中使用的标准周期性二面角形式不同. 最后, 将topol.top文件中的#include语句由

1
#include "gromos53a6.ff/forcefield.itp"

修改为

1
#include "gromos53a6_lipid.ff/forcefield.itp"

最后, 我们需要包含DPPC分子的特定参数. 做法非常简单. 在你的topol.top文件中添加一行#include “dppc.itp”即可, 放在蛋白质位置限制部分的后面, 位置限制部分是蛋白质[ moleculetype ]定义的结束. 这些做法与在拓扑文件中添加任何其他小分子或溶剂类似. 在本节及整个教程中, 绿色文本表示你要添加的行, 其他文本(黑色)表示修改前拓扑中已有的内容.

1
2
3
4
5
6
7
8
9
10
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include DPPC chain topology
#include "dppc.itp"
; Include water topology
#include "gromos53a6_lipid.ff/spc.itp"

如果你正确地遵循了上面的所有步骤, 你就会得到一个全功能的力场, 可与gmx pdb2gmx一起用于处理其他膜蛋白. 这样做不必手动修改topol.top. 将新的gromos53a6_lipid.ff目录放于$GMXLIB下就可以在整个系统中使用这个力场了.

3.定义单位晶胞并添加溶剂

比起水中的蛋白质来, 为膜蛋白定义单位晶胞更复杂一些. 构建单位晶胞时有几个关键步骤:

在相同的坐标系中对蛋白和膜进行取向
在蛋白周围堆积脂分子
使用水进行溶剂化

1.蛋白质和膜的取向
(1) 使用gmx grompp对仅含有DPPC的体系生成.tpr文件. 你可以使用任何有效的.mdp文件, 相应于纯DPPC的拓扑文件. 这里是一个示例.mdp文件, 以及一个拓扑文件. 注意拓扑文件非常简单, 仅包含了dppc.itp和spc.itp, 用以读入DPPC和水的参数. 就这么简单! 运行gmx grompp:

这里教程没有写命令，我暂时这么写：

1
editconf_mpi -f dppc128-q.pdb -o dppc128.gro -princ

选择2 DPP 正确与否待考

1
gmx grompp -f minim.mdp -c dppc128.gro -p topol_dppc.top -o em.tpr

你可能会得到一个致命错误, 像这样, 但在这种情况下对上面的命令可以安全地使用-maxwarn 1选项来忽略错误. 可以这样做的原因见这里. 请注意 这一步是使用topol_dppc.top的唯一一步, 此文件不能用于任何其他目的, 你不应该将它用于本教程中的 任何 剩余步骤

对复杂的分子, 去除周期性的子程序可能会崩溃, 在这种情况下你可以使用gmx trjconv.但是后面教程里做了

1
gmx grompp -f minim.mdp -c dppc128.gro -p topol_dppc.top -o em.tpr

(2) 使用gmx trjconv移除周期性

1
gmx trjconv -s em.tpr -f dppc128.gro -o dppc128_whole.gro -pbc mol -ur compact

这一步产生的结构有一定的飞线，原因待考

vmdscene
vmdscene

现在查看下这个.gro文件的最后一行, 它对应于DPPC单位晶胞的x/y/z盒向量. 我们需要在相同的坐标系内调整KALP多肽的取向, 并将多肽的质心置于盒子的中心:

1
gmx editconf -f KALP-15_processed.gro -o KALP_newbox.gro -c -box 6.41840 6.44350 6.59650

2. 在蛋白四周堆积脂分子

目前我发现, 围绕嵌入蛋白质堆积脂分子的最简单方法是InflateGRO方法(参考文献), 你可以在这里inflategro. 请注意, 我发布的代码是我自己保存的InflateGRO原始版本的副本, 而 不是 来自作者的InflateGRO2. 下载上面链接中的文件, 将其重命名为inflategro.pl再继续. 首先, 整合蛋白质和双脂层的结构文件:

1
cat KALP_newbox.gro dppc128_whole.gro > system.gro

移除不需要的行(来自KALP结构的盒向量, 来自DPPC结构的标题信息), 并相应地更新坐标文件的第二行(总原子数).

InflateGRO脚本的作者建议对蛋白质的重原子使用非常强的位置限制力, 以保证在EM过程中蛋白质的位置不发生改变. 在拓扑中添加新的#ifdef语句, 调用特殊的位置限制, 这样你的拓扑现在就包含类似下面的部分:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Strong position restraints for InflateGRO
#ifdef STRONG_POSRES
#include "strong_posre.itp"
#endif
; Include DPPC chain topology
#include "dppc.itp"
; Include water topology
#include "gromos53a6_lipid.ff/spc.itp"

现在就可以使用gmx genrestr生成新的位置限制文件:

1
gmx genrestr -f KALP_newbox.gro -o strong_posre.itp -fc 100000 100000 100000

在用于能量最小化的.mdp文件中, 添加一行define = -DSTRONG_POSRES以保证使用这些新的位置限制. 然后, 简单地按照InflateGRO的指示(包含在脚本自身中)做, 过程很简单. 对脂分子的位置使用因子4进行缩放:

perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat
由于InflateGRO脚本对命令行变量顺序的要求非常严格, 这里需要给出一个简略的说明:

system.gro - 将要施加缩放的输入坐标文件的名称
4 - 施加的缩放因子, 值>1表示扩增, 值<1表示收缩/压缩
DPPC - 施加缩放的脂分子的残基名称
14 - 搜索脂分子的截断半径(单位为Å), 处于此半径内的脂分子会被删除
system_inflated.gro - 输出文件的名称
5 - 计算每个脂分子面积时所用的格点间距(单位为Å)
area.dat - 输出文件, 包含每个脂分子的面积信息, 对确定结构是否合适很有帮助

记下删除了多少脂分子, 并相应地更新拓扑文件中[ molecules ]部分的内容.

There are 1 lipids within cut-off range…
1 will be removed from the upper leaflet…
0 will be removed from the lower leaflet…

貌似我的是1个
运行能量最小化. 然后, 使用因子0.95收缩脂分子(假定你使用了默认的名称, 能量最小化的结果文件为confout.gro):

perl inflategro.pl confout.gro 0.95 DPPC 0 system_shrink1.gro 5 area_shrink1.dat

在网上找到一个脚本，脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
perl inflategro.pl minimizedZ.gro 0.95 DPPC 0 shrinked1.gro 5 area_shrink1.dat
grompp -f minim.mdp -c shrinked1.gro -p topol.top -o minimizedZ2.tpr
mdrun -v -deffnm minimizedZ2
for i in $(seq 1  20)
do
perl inflategro.pl minimizedZ2.gro 0.95 DPPC 0 shrinked2.gro 5 area_shrink1.dat $i
grompp -f minim.mdp -c shrinked2.gro -p topol.top -o minimizedZ2.tpr $i
mdrun -v -deffnm minimizedZ2 $i
done