layout: post title: Martini实例教程：新分子的参数化 categories:
在本教程中, 我们将讨论如何为一个含有已知片段, 但尚未被描述过的新分子构建Martini拓扑. 随后我们将对该Martini模型进行自组装模拟以评估其行为, 但不会进行模型的精细化. 如果你觉得有必要对其进行精细化, 请参见本教程第二部分(不过它是基于另一个分子来说明的), 或者参照脂质和聚合物教程中类似的精细化步骤.
假设现在你开始对某类呈非对称"流星锤"形状的两亲分子开始感兴趣, 图1所示即为此类分子中的一个实例, 它来自论文Masuda and Shimizu, Langmuir 20, 5969 (2004). 这种"流星锤"分子链的两端均为亲水基团, 但两个基团大小明显不一: 一端是相对较大的葡萄糖环, 另一端则是相对较小是羧酸. 根据Masuda和Shimizu在论文里的描述, 这种"流星锤"分子可自组装成管状结构, 管径大小取决于连接臂的长度. 从图1中你可以观察到分子所含的化学基团与标准脂质/蛋白质的基本组成单元相差不大, 因此, 这些化学基团很可能可以使用Martini力场进行描述. 事实上, Martini力场可以处理蛋白和糖类, 也确有文献报道过糖脂分子(与"流星锤"分子的结构类似)的模拟, 相关参数也已经发表, 具体可见Lopez et al., Journal of Chemical Theory and Computation, 9, 1694 (2013).
【付彦凯 注】原文此处没给出论文出处, 但在下文又反复提到了这篇论文, 因此在这里加上了出处.
图1 非对称"流星锤"脂分子族, 大亲水端葡萄糖通过烷烃链与小亲水端连接
【付彦凯 注】原文插图分辨率较低, 译文插图系从原始文献重新截图. 下同.
不幸的是, 你所感兴趣的"流星锤"分子与已报道的糖脂的构成方式并不完全一样, 你需要自己动手构建它的Martini拓扑结构! 不过幸运的是, 你可以最大限度地利用Martini力场中已有的组成单元. 接下来, 我们将为你展示构建过程.
首先, 将图1中的分子拆分为若干个合理的基本组成单元, 该过程有时会被称为(从原子描述到粗粒级描述的)映射(mapping). 如果你熟悉Martini力场中已有的脂质和蛋白质模型, 并考虑到各个分子片段的化学本性, 那么会很自然地把含有十二烷连接臂的"流星锤"分子粗粒化为8个珠子(图1, n=12时). 你可以通过画图来辅助实现这个想法, 把需要粗粒化成一个珠子的所有原子圈在一起, 并写下珠子类型, 正如上面那篇糖脂模拟论文中所做的一样. 图2所示为论文中对糖脂分子GCER的一种合理映射方式. GCER同时代表葡糖苷(脂)酰鞘氨醇(glucosylceramide)和半乳糖苷(脂)酰鞘氨醇(galactosylceramide), 关于两者结构的更详细信息可参见脂库网站. 尽管GCER中含有酰胺基团, 却是通过醚键与糖基连接的, 而在我们的脂分子中, 酰胺基团是直接与糖基相连的. 因此, 对我们的"流星锤"分子, 合理的映射方案可参考GCER, 但键和键角参数将会有所不同.
图2 糖脂模拟论文中GCER分子的映射, 以及非对称"流星锤"脂分子族可能的(部分)映射
让我们给将要Martini化的"流星锤"分子起个新名字GDAL(glucose-dodecane-acid-lipid), 当然你可以给它起个更好的名字, 不过从现在开始我们还是叫它GDAL好了...
你需要给GDAL写一个拓扑文件, 可以自己从头开始写, 也可以在其他分子拓扑文件的基础上进行改造, GCER分子的拓扑文件就是一个很好的选择. GCER自身只是连接到脂酰基的头基的名字, 如果去检索Martini力场官网提供的脂分子的拓扑文件, 我们会发现一个名为DPGS的带有两条脂肪链尾巴的糖脂分子. 接下来我们将以DPGS的拓扑文件为模板, 其.itp文件部分内容如下:
[ moleculetype ]
DPGS 1
[ atoms ]
1 P4 1 DPGS C1 1 0 72.0000
2 P1 1 DPGS C2 2 0 72.0000
3 P1 1 DPGS C3 3 0 72.0000
4 P1 1 DPGS AM1 4 0 72.0000
5 P5 1 DPGS AM2 5 0 72.0000
6 C3 1 DPGS D1A 6 0 72.0000 ; corrected (C3 instead of C1), oct 2013
7 C1 1 DPGS C2A 7 0 72.0000
8 C1 1 DPGS C3A 8 0 72.0000
9 C1 1 DPGS C4A 9 0 72.0000
10 C1 1 DPGS C1B 10 0 72.0000
11 C1 1 DPGS C2B 11 0 72.0000
12 C1 1 DPGS C3B 12 0 72.0000
13 C1 1 DPGS C4B 13 0 72.0000
根据 1) 其他脂类分子的一般映射规则 2) DPGS分子.itp文件中的成键规则 3) 图2所展示信息, 不难猜测, DPGS分子.itp文件中, 1-3号珠子代表葡萄糖环, 4号珠子代表甘油部分, 5号珠子代表酰胺基团, 其他珠子代表脂肪链尾端. 如果据此编写GDAL分子的.itp文件, 5号珠子应该和3号珠子相连, 而不是与DPGS中一样, 和4号珠子相连, 然后脂肪链尾端(10-13号珠子) 应该与5号珠子相连. 由于GDAL是十二烷连接臂(连接臂占用12/4=3个珠子), 原来的13号珠子应该修改为羧酸末端. 按照上述方式对原DPGS分子的.itp文件进行修改, 我们会发现, 除了13号珠子外, 其他珠子的设定均不用改动. 根据2007年发表的Martini力场论文, 13号珠子应该从C1改为P3类型(论文中P3代表乙酸). 综上所述, GDAL分子.itp文件编写过程如下:
.itp文件, 将DPGS替换为GDAL;C1修改为P3[ bonds ], [ angles ]和[ dihedrals ]中所有涉及4, 6, 7, 8, 9号珠子的条目都删去如此修改完成后, 仍有一个地方需要修改, 即涉及4号酰胺基团珠子(原5号珠子)的成键参数, 包括环-酰胺的键长, 环-酰胺-烷烃的键角. 这些参数要设为合理值. 酰胺基团整体非常刚性, 因此以其为顶点的键角应接近180°, 且具有高的力常数; 参与成环的珠子与酰胺基团珠子成键的键长可参考蛋白质中骨架粒子-骨架粒子的键长. 你可以从已发表的Martini力场相关论文中查找这些信息.
有了拓扑文件之后, 你还可以直接下载DPGS的坐标文件, 删去4, 6, 7, 8, 9号珠子, 使用它做真空中的能量最小化, 然后利用得到的坐标随机地部分填充空的模拟盒子. 可以再向盒子中添加水分子, 10%的水分子可修改成抗冻水, 然后就可按照脂质教程的步骤进行自组装模拟了.
【付彦凯 注】普通水名称W, 粒子类型P4; 抗冻水名称WF, 粒子类型BP4. 两者均已在Martini力场文件中定义.
如果你完成了上述所有步骤, 可点击这里下载我们做的自组装模拟, 与你的进行比较. 200个分子足够形成双层结构, 但是双层结构弯曲度没有达到预期效果, 其原因可能是, 两种亲水头端并未像Masuda和Shimizu所预测的那样各自分布在两亲性双层结构一侧. 而导致两种亲水头端未能分离开来的原因又可能是: 1) 插入分子个数不够; 2) 周期性边界条件和/或 参数设置; 3) GDAL分子从双层结构中实现180°翻转需要较长时间. 这些问题也正是研究最有趣的地方! 你可以随意改动珠子类型设定从而比较不同相互作用强度下的模拟结果. 也许, 你需要在全原子力场下模拟单个GDAL分子在水中的行为, 确认其柔性正确后再进行映射. 在接下来的教程中, 我们将以甲苯为例来讲述类似的过程, 在脂质和聚合物教程中也有类似步骤. 好好享受吧!
在本部分教程中, 我将向你展示如何为一个未参数化的新分子构建Martini拓扑结构. 构建过程高度依赖于分子的化学结构以及与其相关的(实验)数据. 我们将使用一个小分子甲苯作为示例分子来讲述该过程. 甲苯的全原子拓扑是已知的, 其水-辛醇分配系数(logP)为2.73(数据来自http://dx.doi.org/10.1021/ci050282s). 教程假设你对Martini力场有基本了解(若否, 请先学习Martini官网上的其他教程), 对GROMACS有基本了解(若否, 请学习GROMACS官网提供的教程), 同时对文本编辑器有基本了解.
对新分子进行可靠的Martini化是基于其在溶液中的全原子力场模拟数据来开展的. 我们自己做Martini化通常基于Gromos力场, 但你也可以自由选择任何力场. 如果没有拓扑文件, 你可以基于ATB或GAFF力场自己来构建一个(这一步的工作量可能也很大), 或者其他人已经发表了你要研究的分子的全原子模拟结果, 你可以将论文中的数据作为参考材料.
溶剂的选择很重要: 要Martini化的分子在极性或无极性疏水溶剂中的表现可能会有所不同, 构建的粗粒化模型可能无法两者兼顾. 选择与你做粗粒化模拟时的环境匹配的溶剂, 如果极性和非极性环境都很重要, 那么只能妥协, 选择表现最好的一个.
你可以输入下面的命令进行甲苯的模拟:
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">grompp -f gromos.mdp -c gromos.gro -p gromos.top -o gromos mdrun -v -deffnm gromos </pre></div>点击下载相应的文件: gromos.mdp, gromos.gro, gromos.top
你刚才运行的是单个甲苯分子在癸烷中的模拟, 模拟时长1 ns. 需要指出的是, 由于甲苯只是个简单的分子, 1 ns的模拟时长足够了.
此步是粗粒化的核心步骤, 需要依靠你的经验, 化学知识以及不断试错. 我最喜欢的做法是: 把分子结构画在纸上, 复制10次, 然后像小学生画画那样尝试不同的映射方法. 该过程没有规律可循, 但我可以提供一些建议:
依据你刚才建立的映射规则, 把在步骤1)中开展的全原子模拟转换为粗粒的模拟. 方法有很多种, 你可以使用转换工具或initram.sh脚本. 不过我认为最简单的方法是, 创建一个GROMACS索引文件, 每一索引组代表一个珠子, 包含该珠子映射原子的编号. 写完索引文件后运行如下命令:
点击下载相应的文件: mapping.ndx
命令中, -ng 3表示索引分组(即珠子个数)有3个, 如果命令前端不加seq 0 2, 你需要挨个输入3个索引分组的名称. 你可以打开我写的索引文件mapping.ndx, 看看是如何对甲苯分子进行映射的.
再去看看你在步骤2)中绘制的映射图, 确定这些珠子之间的成键关系(键, 约束, 键角, 正常和异常二面角). 键长/键角等信息可以从步骤3)中得到的粗粒化轨迹文件(mapped.xtc)中提取. 你可以再制作一个名为bonded.ndx的索引文件, [bonds]下两个珠子一行, [angles]下三个珠子一行, [dihedrals]下四个珠子一行. 使用g_bond来获取键长信息, g_angle来获取键角和二面角信息.
点击下载相应的文件: bonded.ndx
根据g_angle命令生成的angles.xvg文件可以绘制键角的分布直方图(angles.xvg文件第一列y值是三个键角的平均值, 可舍弃), 你可以使用g_analyze命令:
某些情况下, 仅仅计算键角分布的平均值和标准差就足够了, 然而键角分布可能存在双峰, 最明智的选择是绘制分布图观察下, 可使用xmgrace命令来完成:
一些注意事项: 如果你的分子被周期性边界条件分成两半, 你可能会得到诡异的键长分布结果. 此种情况下, 可首先用trjconv -pbc whole处理一下全原子力场的模拟轨迹, 然后再进行映射. 另一点需要注意的是, 在执行上文中的g_bond命令时, 每条键长的平均值和标准差会存储在bonds.log文件里, 而distance.xvg存储的是所有键长平均后的分布情况. 若想得到每条键长的分布情况, 要么为每一条键单独运行一次g_bond命令, 要么参见后文提供的方法. 一旦你有了CG.tpr文件, 即可在g_bond命令中使用, 从而可以直接比较二者的键长分布.
.itp文件尽管你可以对已有的.itp文件进行复制粘贴(通常是为了保持文件格式正确), 大部分情况下, 编写.itp文件的工作还是要手动完成的. 你的.itp文件中(可能)将包括以下几项:
[moleculetype]: 一行即可, 包括分子的名字和键合近邻的排除数. Martini力场下, 键合近邻的排除数为1.[atoms]: 每个珠子一行, 信息包括珠子编号, 珠子类型, 残基编号, 残基名称, 珠子名称, 电荷组编号以及电荷. 对小分子, 残基编号和残基名称和粗粒化之前一样. 珠子名称可自定义. 在Martini力场中, 每个珠子的电荷都有自己的电荷组, 带电量可以为1, 0和-1. 珠子质量已经预定义好了, 在.itp文件中不予注明.[bonds]: 每条键一行, 信息包括珠子1, 珠子2, 键型、键长和力常数. 在Martini力场中, 键型为1. 键长可采用步骤4)中所获得的平均值. 力常数应该根据键长分布的宽度来调整: 窄峰用大的力常数, 宽峰用小的力常数.[constraints] [angles] [dihedrals] 正确格式参见GROMACS手册.从mapped.xtc文件中提取出一帧, 使用与全原子模拟相同的溶剂(粗粒化的)癸烷重新溶剂化. 创建一个.top文件(里面应包含Martini力场的.itp文件, 新创建的粗粒化分子的.itp文件, 并添加正确的分子和溶剂个数). 创建一个Martini力场适用的.mdp文件, 可参考官网上的示例. 在正式开始粗粒化模拟之前, 最好先进行能量最小化, 然后用小的时间步长(如5 fs)跑几千步, 让分子弛豫一下. 除此之外, 你还需要.itp文件, .itp文件以及.itp文件. 最好新建一个文件夹来运行粗粒化模拟.
点击下载相应的文件: em.mdp, martini.gro, martini.top, relax.mdp, md.mdp
模拟结束后, 重新提取键长和键角. 你仍然可以使用步骤4)中那个索引文件.
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">g_bond -f md.xtc -s md.tpr -n bonded.ndx -noaverdist -d bonds_martini.xvg g_angle -f md.xtc -n bonded.ndx -ov angles_martini.xvg -type angle -all </pre></div>注意和在步骤4)中执行g_bond命令时有什么不同. 这里你使用了-s md.tpr和-noaverdist.
步骤6)可能会有两种结果: 模拟顺利进行, 你得到了键长和键角分布; 模拟崩溃. 如果模拟顺利, 你可以对从步骤6)和步骤4)中获得的分布进行比较, 有必要的话, 回到步骤5)调整.itp中的参数以得到更好的分布. 如果步骤6)中得到的分布峰比较宽, 意味着你需要增大力常数; 如果平均值整体有所偏移, 意味着你需要做相应的回调; 诸如此类. 如果模拟崩溃, 你同样需要返回步骤5), 此行目的是让分子更稳定. 可行的方法有: 减小力常数: 去除不必要的二面角, 将约束(constraint)改变为键或反之. 导致系统不稳定的常见因素有:
反复调整直到你对结果满意后, 进入最后一步.
两相分配自由能固然是一个好的评估指标(我们目前正在编写一个教程, 主要讲解针对Martini分子如何有效地计算此参数), 但对于纯的液态体系, 密度也可以用来评估. 此外, 任何实验数据都可以用来评估.