MD基本教程-溶菌酶的水溶液的MD模拟

溶菌酶水溶液的简单分子动力学模拟教程GROMACS是一个使用经典分子动力学理论研究蛋白质动力学的工具[1]。这个软件包遵守GNU协议，因此可以免费从其主页（）上下载。GROMACS可以在linux、unix和Windows上使用。第一步建立拓扑结构本教程使用鸡蛋清溶菌酶（PDBcode：1AK1）作为目标蛋白质。首先需要从RCSB（）网站下载其三维结构文件。然后通过诸如SYBYL、DiscoveryStudio、VMD、Chimera和PyMOL等图形显示软件分析其三维结构。根据模拟的需要去除或保留该文件中的结晶水分子。例如，在研究蛋白质和配基之间的亲和作用时，结合在活性位点附近的水分子就必须保留。但如果利用MD模拟研究蛋白质的构象转换时，蛋白质中的结晶水就可以删除。删除PDB文件中分子可以利用简单的文本编辑软件如vi(Linux)、emacs(Linux/Mac)和Notepad(Windows)删除这些分子所对应的条目。例如，如果要删除文件中高的结晶水就可以将PDB文件中的所有“HOH”残基删除。切记不要用文字处理软件Word处理该文件。在运行pdb2gmx程序之前，一定要检查pdb文件中的MISSING条目。这些条目会标注出该晶体结构中缺失哪些原子或残基。这一步是必须做的，因为不完全氨基酸残基或分子都会导致pdb2gmx命令的失败。这些遗失的原子或残基必须利用其他的软件补充完整。但末端的遗失可能不会对动力学模拟的顺利进行造成影响，但有可能会对结果有影响。此外，还要注意pdb2gmx程序并不是万能的，它并不能为所有的分子建立拓扑结构，仅能对一些力场中定义好的分子才有效，例如大部分蛋白质，核酸，还有一些辅因子如NAD、NADH、ATP和ADP等。此外，一些糖类（海藻糖等）的力场也都定义好，可以从力场文件中获得。因此，如果目标蛋白质含有小分子配基的，在运行该软件之前就需要先用简单的文本编辑软件将配基删除。最后除去结晶水和小分子配基，同时确保所有必需原子存在，此时的PDB文件中只含蛋白质原子，现在可以运行pdb2gmx命令了。常用pdb2gmx命令的详细参数：pdb2gmx-f1AKI.pdb-o1conf.gro–pconf.top–ignh此时程序会提示你选择力场（力场的选择类型与你安装的力场密切相关）：SelecttheForceField:From'/usr/local/gromacs/share/gromacs/top':1:AMBER03forcefield(Duanetal.,J.Comp.Chem.24,1999-2012,2003)2:AMBER94forcefield(Cornelletal.,JACS117,5179-5197,1995)3:AMBER96forcefield(Kollmanetal.,Acc.Chem.Res.29,461-469,1996)4:AMBER99forcefield(Wangetal.,J.Comp.Chem.21,1049-1074,2000)5:AMBER99SBforcefield(Hornaketal.,Proteins65,712-725,2006)6:AMBER99SB-ILDNforcefield(Lindorff-Larsenetal.,Proteins78,1950-58,2010)7:AMBERGSforcefield(Garcia&Sanbonmatsu,PNAS99,2782-2787,2002)8:CHARMM27all-atomforcefield(withCMAP)-version2.09:GROMOS9643a1forcefield10:GROMOS9643a2forcefield(improvedalkanedihedrals)11:GROMOS9645a3forcefield(SchulerJCC2001221205)12:GROMOS9653a5forcefield(JCC2004vol25pag1656)13:GROMOS9653a6forcefield(JCC2004vol25pag1656)14:OPLS-AA/Lall-atomforcefield(2001aminoaciddihedrals)15:[DEPRECATED]Encadall-atomforcefield,usingfullsolventcharges16:[DEPRECATED]Encadall-atomforcefield,usingscaled-downvacuumcharges17:[DEPRECATED]Gromacsforcefield(seemanual)18:[DEPRECATED]GromacsforcefieldwithhydrogensforNMR力场中包含将写进拓扑文件中的所有信息，所以此选择非常重要，应仔细衡量每个力场且选择与你的模拟情况最相符的力场。在此教程中，我们选择全原子的OPLS力场，因此，选择14并回车。此外，在pdb2gmx软件中还有很多其他操作参数可以选择：如-ignh：在PDB文件中忽略非极性氢原子，尤其在处理利用核磁共振获得的三维结构时尤其需要利用该参数；-ter：交互指定N末端和C末端的带电情况；-inter：交互指定带电氨基酸（精氨酸，赖氨酸，天冬氨酸，谷氨酸和组氨酸）的电荷，指定半胱氨酸的二硫键。运行pdb2gmx命令会获得GROMACS运行所需要的三个文件力：1．conf.top分子拓扑文件2.posre.itp位置限制文件3.conf.gro分子结构文件conf.top为分子拓扑文件（后缀名为top）包含所模拟分子的所有信息，例如原子类型、电荷数、键长、键角和二面角等。posre.itp包含限制重原子位置的信息。conf.gro为一个GROMACS格式的结构文件，它包含在已定义力场中上所有原子的结构文件（氢原子已经加入其中）。conf.top是整个系统文件的拓扑文件。下面首先来看看分子的拓扑信息。利用文本编辑器打开conf.top，几行之后你会发现如下句子：#includeoplsaa.ff/forcefield.itp，这一行指定了本次MD模拟选用OPLS-AA力场。即，本次MD模拟的所有拓扑参数皆来自于此力场。第二个重要行是[moleculetype]，在它之后可以看到;NamenrexclProtein_A3Protein_A是分子的名字，即蛋白质在PDB文件中被标为A链。大部分内容为蛋白质分子中的原子信息，以下所示：[atoms];nrtyperesnrresidueatomcgnrchargemasstypeBchargeBmassB1opls_2871LYSHN1-0.314.0067;qtot-0.32opls_2901LYSHH110.331.008;qtot0.033opls_2901LYSHH210.331.008;qtot0.364opls_2901LYSHH310.331.008;qtot0.695opls_293B1LYSHCA10.2512.011;qtot0.946opls_1401LYSHHA10.061.008;qtot1注释如下：-nr：原子数；-type：原子类型；-resnr：氨基酸残基数；-residue：氨基酸残基名；注意，在pdb文件中为”LYS”；而在.rtp文件中改为LYSH表明赖氨酸质子化。-atom：原子名；-cgnr：带电基团数，带电基团定义整个基团所带电荷为整数，这样有利于加速计算；-charge：原子所带电荷；-mass：原子质量。bonds：成键信息；pairs：成对信息；angles：键角信息；dihedrals：二面角信息TypeB，chargeB和massB：在自由能微扰计算过程中需要指定，在本章中不予讨论。posre.itp：建立位置限制，它定义了在平衡中用于保持原子位置的一种平衡常数;IncludePositionrestraintfile#ifdefPOSRES#includeposre.itp#endif至此对蛋白质A的分子类型设定基本结束。拓扑文件中剩余的内容都是用来定义其他分子和描述系统的等级。1.定义水分子类型：;Includewatertopology#includeoplsaa.ff/spce.itp#ifdefPOSRES_WATER;Positionrestraintforeachwateroxygen[position_restraints];ifunctfcxfcyfcz11100010001000#endif水分子也可以被限定位置，一般使用力常数为1000kJmol-1nm-22.离子参数描述：Includegenerictopologyforions#includeoplsaa.ff/ions.itp3.最终来到系统等级设定[system]指令给出系统名称并将在模拟中写入模拟文件，[molecules]指令列出系统中所有分子。[system];NameLYSOZYME[molecules];Compound#molsProtein_A1一些关于[molecules]指令的提示：1.molecules中所列分子顺序必须与轨迹文件（.gro）中分子的顺序一一对应，否则会出错；2.名字列表必须与[moleculetype]名字一致，不只是残基名或其它。第二步定义模拟盒子大小本章模拟蛋白质在简单的水溶液中构象转换情况。当然我们也可以模拟蛋白质在任意溶液中的构象转换情况，如果能够获得这些溶剂分子的拓扑参数的话，如何获得分子的力场参数我们会在下面的内容中提到。在GROMACS软件中，定义盒子尺寸是用editconf命令来完成的。现有多种晶胞类型供选择：立方体，长方体和菱形十二面体。本章中我们选择最简单的立方体盒子作为晶胞。当你对周期性边界条件和盒子类型更加熟悉以后，可以使用菱形十二面体盒子，因为其体积是同周长立方体的~71%，这样就可以降低模拟体系内的水分子个数，从容大大节省计算量。Editconf的命令范例：editconf-f1AKI.gro-o1AKI_box.gro-c-d1.0-btcubic–box6.56.56.5参数如下所示：-c：使蛋白质位于盒子中央；-d：使蛋白质与盒子距离最近的原子与盒子边缘的距离大于等于1纳米。-bt：盒子类型。蛋白质距盒子边缘的距离是一个重要的参数。因为MD模拟中使用周期性边界条件，因此原子之间必须满足原子之间的距离最小规则，即蛋白质不能与其镜像发生作用，否则MD模拟过程中在计算力的过程中会发生错误。将蛋白质与盒子边缘距离设置为1纳米就意味着任何两个蛋白质镜像分子间的距离至少为2纳米。第三步加入溶剂分子现在可以向模拟盒子中添加溶剂水分子了。genbox命令可以完成以上操作genbox-cp1AKI_box.gro-csspc216.gro-o1AKI_solv.gro-pconf.top此处溶剂水的结构是用GROMACS中的默认值（spc216.gro）。它是一个预平衡好的的SPC水模型。输出结果为1AKI_solv.gro，此外需要指定genbox中拓扑文件的名字，以便对其进行修改，注意比较运算前后条目[molecules]的改变。[molecules];Compound#molsProtein_A1SOL10832Genbox仅需要记录加入的水分子数，稍后将其写入拓扑文件，从而改变模拟环境，因为水的信息更新是事先编程好的。注意：如果模拟非水溶剂，genbox不会修改拓扑文件，此时就需要自己去修改添加。第四步加入离子从第一步pdb2gmx的输出结果可以看出，所模拟的蛋白质溶液系统的净电荷为-22，这主要基于其氨基酸成分。如果你在pdb2gmx输出结果中没注意看到此信息，请转向拓扑文件中[atoms]指令的最后一行，应显示qtot8.。因为生命体是不带净电