Materials Studio官方教程（Help-Tutorials）- 脂质双分子层的粗粒化分子动力学

背景：粗粒化分子动力学能够模拟比原子分子动力学更长的时间和更大的长度尺度。这是通过将几个原子的基团表示为单个单元或珠来实现的。为体系中每种类型的珠之间的相互作用参数化力场。一般来说，粗粒化力场必须为每个待研究的体系参数化。然而，Marrink等人已经生成了一个粗粒化力场MARTINI，该力场针对广泛的有机部分进行参数化，专注于模拟生物分子材料，如脂质。

二棕榈酰磷脂酰胆碱（DPPC）是一种磷脂，是肺表面活性剂的主要成分。它还用于研究脂质体、脂质双层和生物膜模型。

Materials Studio包括MARTINI力场的略微修改版本，MS Martini，可用于运行珠结构的粗粒化分子动力学模拟。

简介：在本教程中，您将使用 Mesocite 通过 MS Martini 力场运行脂质双层的粗粒化分子动力学模拟。您将使用mesostructure building工具来构建初始起始结构并执行双层的一些平衡。

目的：演示如何使用mesoscale modeling工具设置和模拟磷脂双分子层

本教程重要节点：

建立磷脂和溶剂分子模型-应用力场类型并修改力场-建立磷脂双分子层模型-优化和平衡结构-模拟磷脂双分子层并分析结果

1. 建立磷脂和溶剂分子模型

您可以通过Build菜单或使用Mesostructure工具栏访问中间结构构建工具。要使用工具栏，请确保它是可见的。

从菜单栏中选择View | Toolbars，并确保启用了Mesostructure工具栏。

第一步是定义用于构建磷脂和溶剂分子的珠子类型。将使用与Marrink等人相同的表示法，其中1个珠子表示大约四个重原子，平均质量为72 g/mol，有效半径为2.35 Å。

单击Mesostructure工具栏上的Bead Types按钮，打开Bead Types对话框。单击Defaults…按钮，打开Bead Type Defaults对话框。

将Mass设置为72，Radius设置为2.35，然后关闭Bead Type Defaults对话框。

在Bead Types对话框中，定义以下新的珠子类型：C、GL、NC、PO和W。关闭该对话框。

使用定义的珠子类型，可以构建磷脂和溶剂的初始珠子表示。

单击Mesostructure工具栏上的Mesomolecule按钮，打开Build Mesomolecule对话框。

为了构建具有其特征构象的磷脂分子，可以分三步构建该分子。从磷酸盐和胆碱的头部基团开始。

定义1单位PO和1单位NC的重复单元。单击Build按钮。

接下来将两个甘油珠子添加到磷酸基团中。

删除现有的重复单元定义，并定义一个1单位GL的新重复单元。在结构中选择PO珠子，并勾选Add to branch points复选框。点击More…按钮打开Mesomolecule Branches对话框，并确认勾选Branch from selected beads，且Number of branches to attach设置为2。

关闭对话框并单击Build按钮。

提示：使用Label对话框可在文档中添加BeadTypeName标签。

最后将碳尾部原子添加到两个甘油珠子中。将添加四个尾部原子并删除不需要的两个尾部原子，以便生成合适的初始几何构型结构。

删除现有的重复单元定义，并定义一个包含4单位C的新重复单元。选择结构中的两个GL珠子，然后单击Build按钮。

在生成的结构中，删除指向头部原子组方向的两个支链。

现在介观分子的结构应如下图所示。

磷脂的初始珠子结构为球棍模型显示样式，珠子由BeadTypeName标记

Save文档并close视图。在Project Explorer中，将文档rename为DPPC.xsd。

最后构建一个溶剂分子。这将是一个W类型单个珠子，代表四个水分子。

删除现有的重复单元定义，并定义一个新的重复单元，即1个W单元。单击Build按钮并关闭对话框。

Save并关闭新文档，然后将其rename为solvent.xsd。

2. 应用力场类型并修改力场

您将使用MS Martini力场来执行本教程中的所有计算。MS Martini基于Marrink等人开发的MARTINI力场，并已被参数化以研究磷脂双分子层和其他有机化合物。对于所有力场类型，基本力场使用4.7 Å的平衡键长和180°的平衡键角。对于非线性键（例如顺式不饱和键），Marrink等人使用了120°的修正角，并使用该非线性顺式键角来表示GL-PO-GL键角。默认情况下，MS Martini中不包含此非线性顺式键角度，为应用到该数据，需要对提供的力场进行细微修正。

当要修改现有的MS Martini力场时，需要将其导入到项目中。

从菜单栏中选择Modules | Mesocite | Forcefield Manager以打开Mesocite Forcefield Manager。选择MS Martini并单击>>按钮打开力场。关闭对话框。

Save并关闭文档，将其rename为MS Martini CIS.off。

首先，将定义珠子的力场类型。为MARTINI文件中列出的DPPC磷脂参数化类型。

由于Mesocite中没有珠子类型规则，因此需要手动将力场类型分配给珠子。可以使用Properties Explorer或Mesocite对话框执行此操作。在下文中，将使用Mesocite对话框。

从菜单栏中选择Modules | Mesocite | Calculation。

将打开Mesocite Calculation对话框。Setup选项卡可用于定义计算任务，Energy选项卡定义力场和非键项，Job Control选项卡用于设置计算参数。在自定义forcefield时，应选择导入的Martini版本，而不是标准版本。

选择Energy选项卡，然后从Forcefield下拉列表中选择Browse…。在Choose Forcefield对话框中，选择MS Martini CIS.off。

单击Forcefield后面的More…按钮以打开Mesocite Preparation Options对话框。

Mesocite Preparation Options对话框显示当前力场中所有力场类型的列表。现在，可以将选定类型分配给磷脂和溶剂结构中的选定珠子。

打开DPPC.xsd文档。在3D Viewer中，按下ALT键并double-click任意C珠子。

将选择BeadTypeName为C的所有珠子。

在力场类型列表中选择C1，然后单击Assign按钮。对GL、PO和NC珠子类型重复此操作，按上表分配力场类型。

两种珠子类型，NC和PO，也分配了电荷。需要使用Properties Explorer或Charges对话框分配电荷。

选择PO珠子，并在Properties Explorer中将Charge设置为-1。对NC珠子重复此操作，将Charge设置为1。保存并关闭文档。

最后，需要将溶剂珠子的力场类型指定为高极性（P4）。

打开solvent.xsd文档。在Mesocite Preparation Options对话框中选择P4，单击Assign按钮并关闭对话框。保存并关闭文档。

现在，将通过为GL-PO-GL键角添加一个显式键角项来修改力场，该键角对应于力场类型Na-Qa-Na之间的角度。

打开力场文件MS Martini CIS.off，选择Interactions选项卡。从Show interaction下拉列表中选择Angle Bend。在空白一行中，将Fi和Fk设置为Na，将Fj设置为Qa。将Functional Form修改为Cosine Harmonic。

对于GL-PO-GL键角，Marrink等使用了一个120°的平衡键角，或T0，及一个10.8 kcal/mol (45 kJ/mol)的力常数。

将T0设置为120，K0设置为10.8。保存并关闭力场文件。

现在可以使用介观结构模板Mesostructure Template创建工具将构建介观分子体块结构。

3. 建立磷脂双分子层模型

在这一部分的教程中，将要使用Mesostructure builder工具构建包含bilayer和water体相结构。块体的盒子大小为64×64×100 Å。向盒子中添加DPPC lipid构建slab形式的结构。

使用实验中每个lipid的面积计算slab的厚度，值为65 Å2。事实上，在Martini中一个solvent代表了四个水分子 (Marrink et al.)。 首先利用水的密度和摩尔质量确定珠子的体积。

水的密度为1 g/cm3，摩尔质量为18 g/mol。因此一个水分子的体积为18/Nav cm3，或者为18/Nav 1024 Å3 = 29.9 Å3，这里Nav = 6.02×1023是阿佛加德罗常数。因为每个珠子中包含四个水分子，所以一个珠子的体积为 4×29.9 = 119.6 Å3。

下面利用实验中每个lipid的面积，确定lipid珠子的数目。

Slab的总面积为64×64 = 4096 Å2。因为每个lipid的面积为65 Å2 ，所以在双分子层的每一边均有 4096/65 = 63个lipids，共有 2×63 = 126个lipids。每一个lipid含有12个珠子，因此slab共含有126×12 = 1512个珠子。

利用珠子数和它们的体积，便可以计算slab的厚度.

共有1512个珠子，每个珠子的体积为119.6 Å，所以总的体积为1512×119.6 = 180835 Å3。因为面积为4096 Å2，所以计算的厚度为180835/4096 = 44 Å。

现在开始搭建介观模板，从之前体系开始。

在主菜单中选择 Build/Build Mesostructure/Mesostructure Template 打开Build Mesostructure Template对话框，设置X和Y为64，Z为100. 在Filler文本框中输入solvent，单击Build按钮。

生成了一个新的template文件，此文件包含一个蓝色体系。向这个template中添加厚度为44 Å的slab。

将Former type的类型改为Slab. 将Depth值设44，Orientation设置为Along Z。

接下来需要利用surface packing向两个表面中堆积 lipid珠子。

勾选Enable surface packing复选框，在Filler对话框中输入lipid，单击Add按钮，然后关闭对话框，Save文件。

现在您已经构建了模板，您需要用介观分子填充它。

在主菜单中选择Build | Build Mesostructure | Mesostructure打开Build Mesostructure对话框。

此时需要指定与fillers相关联的components。

在与solvent filler相对应的Mesoscale Molecule列中选择solvent.xsd. 在lipid filler中选择DPPC结构。

设置目标 Density为1 g/cm3。

在Options 选项中，不勾选Randomize conformations复选框。

在构建介观结构之前，您需要在脂质上定义头部和尾部珠粒，这样脂质头部就会粘在slab former的表面上。

单击Build Mesostructure对话框中的More…选项，打开Bead Packing Options对话框。此时需要让3D文件显示的结构为DPPC.xsd。

选择NC头部珠子，然后单击Create bead Head set from selection按钮。选择两个终止的尾珠，将Bead tag标记更改为Tail，然后单击Create bead Tail set from selection按钮。关闭对话框。

标记的DPPC珠结构如下所示：

您现在可以构建介观结构。

激活template文件，在Build Mesostructure 对话框中，单击Build按钮，关闭对话框。

此时生成了一个新的3D文件，结构中含有solvent和lipid分子且堆积在template中不同的区域。

在Project Explorer中，将此结构rename为 bilayer.xsd。

在主菜单中选择 File | Save Project，而后在主菜单中中选择 Window | Close All。

教程进行到这一阶段，已经构建的初始结构。接下来将要优化平衡结构和进行生产运行。

4. 优化和平衡结构

由Mesostructure Builder构建的结构是非常粗糙的，会包含重叠珠子。因此，在对这一结构进行动力学模拟前，需要首先对构建的结构进行优化，以移除过密堆积的现象和结构的初始应力。

进行结构优化最好分三个步骤进行：

优化分子间相互作用，保持内部构象和晶胞刚性
优化固定晶胞中的所有分子相互作用
优化单元和所有其他自由度。
这减少了通过引入不切实际的粘结长度或角度来解决珠子重叠。

打开Mesocite Calculation对话框，在Setup选项中，将Task改为Geometry Optimization并单击More… 按钮打开Mesocite Geometry Optimization对话框，勾选Keep motion groups rigid复选框.

通过 keeping motion groups rigid，使优化过程中所有运动组中的原子之间的距离是固定的。为使所有的分子为刚性分子，需要为每一个分子定义motion group。

激活bilayer.xsd结构文件,单击More…选项，打开Motion Groups对话框，单击Assign automatically按钮，关闭对话框。

现在可以开始进行初始结构优化。

在Mesocite Calculation对话框中，点击Run按钮。

优化这一结构可能需要几分钟的时间。优化过程中能量文件是实时更新的，可以观察能量的变化。

打开 bilayer Energies.xcd文件。

从能量的变化曲线中可以发现，初始结构的能量很高，并呈现快速下降的趋势。这是因为初始构建的结构中重叠的珠子在优化的过程中被逐渐地移除。

结构优化完成后需要进行的下一步工作是进行另一个最小化过程。这一过程没有约束 motion groups，并且对晶胞参数进行优化。

打开bilayer Mesocite GeomOpt\bilayer.xsd。

打开Mesocite Geometry Optimization对话框，不勾选Keep motion groups rigid 复选框。

点击Run按钮。

最后还要优化单元参数。

打开bilayer Mesocite GeomOpt\bilayer Mesocite GeomOpt\bilayer.xsd.在Mesocite Geometry Optimization对话框中选中Optimize cell复选框。关闭对话框。

单击Run按钮。

要减少文件夹嵌套，请将生成的结构保存到项目的根文件夹中。

打开bilayer Mesocite GeomOpt\bilayer Mesocite GeomOpt\bilayer Mesocite GeomOpt\bilayer.xsd。在File菜单上选择Save As….输入bilayer_opt作为File name，导航到Documents文件夹，然后单击Save。

您可以从结构中删除不再需要的motion groups，然后在继续下一步之前清除工作区。

按住ALT键并双击bilayer_opt.xsd中的运动组以选择所有运动组。按DELETE。

从菜单栏中选择File | Save Project，然后选择Window | Close All。重新打开bilayer_opt.xsd.

现在继续运行动力学。

在Setup选项卡上，将Task更改为Dynamics，然后单击More…按钮以打开Mesocite Dynamics对话框。将Time step设置为 40 fs，并将Ensemble更改为 NPT。

时间步长（time step）40 fs 是Marrink等人参数化MARTINI力场时使用的time step。 选择NPT系综来平衡体系压力和温度。首先利用 Velocity scaling thermostat 来平衡体系温度，而后利用Andersen barostat来优化体系的残余应力。 Velocity scaling是非常合适的控温thermostat，可迅速使整个体系更接近平衡状态

选择Thermostat选项，设置Thermostat为Velocity Scale.在 Barostat选项中设置 Barostat为Andersen，关闭对话框。

在Mesocite Calculation对话框中点击Run按钮

进行动力学计算得到的结构适于进行产出运行并分析结果。

5. 模拟磷脂双分子层并分析结果

最终的产出运行是一个长时间的NPT 动力学，使用Nose恒温和 Andersen 恒压.

激活bilayer_opt Mesocite Dynamics文件夹中的 bilayer_opt.xtd.

打开Mesocite Dynamics 对话框设置Thermostat选项为Nose.在 Dynamics选项中设置Frame output every 为250 STEPS，关闭对话框。

提示：实际的计算中，需要进行一个更长时间的production run, 通过密度的变化曲线判断是否收敛。

可以重新开始进行计算，以使先前模拟得到的原子速度作为此次模拟的初始速度。当已经进行了速度的平衡计算，重新开始计算可以保证重新利用得到的速度。

在 Mesocite Calculation 对话框中, 勾选Restart复选框，然后单击Run按钮。在出现的警告对话框中，单击Yes按钮。关闭Mesocite Calculation对话框。

从开始到计算完成可能需要几分钟的时间。计算完成后，将关注的目光转向轨迹文件并进行一些分析工作。

从主菜单上选择 File | Save Project,而后选择 Window | Close All。

下面对结构进行分析.

从分析的结果，会得到GL- PO – GL键角和C – PO- C键角的角分布的差异。首先，需要编辑DPPC分子，添加GL-PO-GL和C-PO-C 角度；需要进行针对轨迹的参数匹配寻找所有的GL-PO-GL和C-PO-C 角，并分析角度分布。

在bilayer_opt Mesocite Restart文件夹中打开轨迹文件和初始的DPPC结构。

激活DPPC.xsd结构文件，按如下图所示定义两个键角，选择 GL-PO-GL键角。

可以在分析中使用的两个键角

下面使用参数匹配工具为轨迹文件中的键角创建两个设定：

从主菜单中选择Edit | Find Patterns打开 Find Patterns对话框。在Pattern document中选择DPPC.xsd，同时结构中选择了GL-PO-GL 键角.将Match property 设置为BeadTypeName.激活bilayer_opt.xtd结构文件并点击Find按钮。

体系中所有的GL-PO-GL 键角均创建了Angle monitors，并且都被选定；可以将它们设置为一个set以便于选择和之后的分析。

单击New Sets…按钮，打开 Define New Set对话框。输入GLPOGLAngles，单击OK按钮。在bilayer.xtd结构文件中，取消键角的选择。

现在，可以利用Find Patterns工具寻找 C-PO-C键角并对这些键角创建set。

回到 DPPC.xsd并选择C-PO-C键角；激活bilayer_opt.xtd结构文件，在Find Patterns对话框中点击Find按钮。

单击New Sets…按钮，打开Define New Set对话框，输入CPOCAngles点击OK按钮，然后关闭 Find Patterns对话框。在bilayer_opt.xtd结构文件中任意位置单击以取消键角的选择。

通过以上步骤已经创建了包含两个不同键角的两个sets。接下来可以对键角的分布进行分析。

选择主菜单中的Modules | Mesocite | Analysis打开Mesocite Analysis对话框。

在分析对话框中可以对轨迹文件进行多种类型的分析工作。在这个例子中，将要分析键角分布。

在Analysis列中选择Angle distribution。在Sets下拉菜单中选择GLPOGLAngles ，然后单击 Analyze按钮。

下一步，在Sets下拉菜单中选择 CPOCAngles，然后单击 Analyze按钮。

进行两次Analyze之后将会产生两个不同的图表文件。GL-PO-GL键角分布的峰值出现在100°左右，对于C-PO-C键角，在10°到80°之间出现一个展宽的峰。这表明尾部角度相比约束的头部来说具有更大的灵活性。

参考文献

S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de Vries., “The MARTINI forcefield: coarse grained model for biomolecular simulations.”, J. Phys. Chem. B, 111:7812-7824, 2007.

本入门教程到此结束。