专利名称：一种用于建立蛋白质体系分子模拟力场的方法
技术领域：
本发明涉及一种用于蛋白质体系的分子模拟力场。
背景技术：
自从分子力学诞生以来，不同的计算机模拟方法已经发展起来。如研究生物分子动力学的MD方法，研究大分子和小分子之间的相互结合的docking方法，从头预测蛋白质结构的Rosetta方法。对于所有这些计算模拟方法，正确描述研究体系中原子之间的相互作用是基础和核心，它的精确性直接关系到模拟结果的可靠性。原则上我们可以 通过使用量子化学计算方法(如分子轨道，密度泛函)求解薛定谔方程来得到原子之间相互作用，但对于一个生物分子体系(包括大分子本身和周围环境)求解薛定谔方程是不现实的，基于分子力学的经验力场成为行之有效的手段。当今普遍应用于生物大分子模拟的力场有AMBER、CHARMM, OPLS和GROMOS力场。在所有这些力场中，描述原子间的相互作用的物理模型基本一致，主要包括键的伸缩、键角弯曲、二面角扭转、范德瓦尔斯和静电相互作用能(式I)。不同力场的差别在于它们应用不同的方法所得到的不同参数组。E = [Kr(r — req f +YjKei0 ~ d,qf +S yt1 + cos(^° 一 /)](式I)
W A" B" \ - hQ j+X(7T-+
丨〈丨Klj,<j ^1J取样的完整性和力场的精度是分子力学模拟中的两个瓶颈问题。在发展的初始阶段，由于极为有限的计算机速度，发展快速有效的取样方法成为当时的研究重点。一些人甚至认为，只要给予足够快的计算机和足够长的时间，我们就可以解决复杂的蛋白质折叠问题。当初IBM公司的Bluegene工程的主要任务之一就是用来研究蛋白质折叠机理。然而随着取样完整性的改善(如并行计算，Replica exchange MD)和计算机速度的大大提高，越来越多的计算机模拟实例揭示力场的精度还不够。现在普遍应用的力场都具有构象偏向的(Bias)缺点。如AMBER力场往往倾向于右螺旋构象，而OPLS和GROMOS偏爱P -折叠构象。当前这些力场研究小组试图用non-additive力场或重新拟合主链参数,虽然对所检测的体系，这些方法看起来有了一些改善，但对真正的蛋白质折叠研究，效果并不清楚。直到现在还没有一个真正的从头折叠a/0蛋白质的成功实例。发展具有足够精度的力场已成为这个领域的共识和当务之急，它直接影响该方法应用的广泛性和预测性能力。
发明内容
一套经验力场的建立包括参数拟合和精度检验。一般地，参数通过拟合由量子化学计算得到的氨基酸二肽模型的构象能(ACE-XXX-NME，XXX = ALA, GLY，或PRO)及静电势得到；力场精确性的检验采用个别的氨基酸二肽模型的构象能作为评价标准；这种评估方法存在不足之处，主要包括(I) 二肽模型中的ACE和NME端基在蛋白质肽链中并不出现，蛋白质计算仅使用模型中心氨基酸的参数，由于端基带来的差别，能较好描述二肽模型的参数并不一定适合描述蛋白质中氨基酸单元。(II)蛋白质中氢键相互作至关重要，但二肽模型不存在氢键，该模型所拟合的力场参数不能描述由氢键作用导致的构象能差异。(III)精度检验通常局限于几个二肽或多肽模型，力场是否适合其它氨基酸并不清楚。首先我们对AMBER极化模型进行了参数化，虽然极化力场较传统的力场精确，但它比传统力场大约慢三倍，这大大地限制了该力场在大分子模拟中的应用。考虑到生物过程(主要是蛋白质折叠)发生所需时间(通常在微秒级别)和现阶段计算机的能力，基于近年来发展起来的快速广义玻恩(Generalized Born)溶剂模型，我们对Amber蛋白质力场进行参数化，具体工作如下重新拟和点电荷，静电场的量子化学计算选用的是M052x/6-31G**方法；在原有的AMBER力场中，使用傅里叶级数表示扭转二面角的能量，对于主链其系数是通过数值拟合量子化学计算得到的(0/V) 二维势能面，在这里我们直接使用量子化学M052X/6-311++G**//M052X/6-31G**方法计算的(￠/ /)势能面。利用数值插值的方法得到能量和力，这个方法克服了傅里叶级数拟合时所带来的误差。
我们要建立的力场是针对连续介质溶剂模型的，Gaussian程序中有众多的连续介质溶剂模型可供选择。力场参数根据考虑了溶剂效应的量子化学计算结果进行拟合。在保证计算工作量合理的同时，选取最优的溶剂模型。为了解决这个问题，我们选取类似于氨基酸侧链并有溶剂化自由能实验数据的分子作为模型，通过与实验数据比较，确定最优的量子化学计算方法和溶剂模型。在Gaussian 03中，考察的量子化学计算方法包括MP2，M05-2X, B3LYP，基组包括6-31G，6-31G#，6-311++G#，cc-pVTZ，aug-cc-pVTZ。在气相条件下，使用相应的理论方法和6-31G**基组优化模型分子的结构，然后在考察的连续介质溶剂中做单点计算。考察的溶剂模型包括 IEFPCM，CPCM, DPCM, IPCM, SCIPCM 和 C0SM0。在 IEFPCM，CPCM 和 DPCM 的计算中，进一步考察了原子半径的选择(包括UAO，UAHF, UAKS,UFF, PAULING和BONDI半径)。在Gaussian 09中，考察了同样的量化方法和基组，溶剂模型只考虑与上述六种原子半径结合的IEFPCM模型和Gaussian 09中新的SMD模型。和Gaussian 03计算不同的是,Gaussian09中分子的构型是采用6-31G**基组在液相中进行优化得到的。我们将针对某一特定的GB溶剂模型建立力场，无疑针对的GB模型应该是最优的。GB模型的精度依赖于三个因素模型本身、GB原子半径以及原子电荷。在未得到新力场之前，使用Amber力场并用上节提到的相似于氨基酸侧链的分子作为模型，考察这些因素对精度的影响。考察的GB方法及相关的GB半径包括IGB1 (Bondi, Mbondi，Mbondi2)，IGB2(Bondi，Mbondi，Mbondi2)，IGB5(Bondi，Mbondi，Mbondi2)，IGB7(Bondi，Mbondi2)和IGB8 (Bondi，Mbondi2)。考察原子电荷的影响时，计算静电势方法包括HF，M05_2X，B3LYP和MP2，基组包括 6-31G，6-31G#，6-311++G#，cc-pVTZ，aug-cc-pVTZ。介质包括气相和液相(共80套)，用RESP方法拟合模型的原子电荷，计算这些分子的溶剂化自由能。根据以上研究，在MAP的方法中，我们以氨基酸的二肽链结构(ACE-X-NME，X = 20个氨基酸之一)为模型，使用30 X 30的网格密度，格点上的能量通过计算得到，而其它任意点的能量通过插值的方法计算。在计算二维势能面时，对于ALA和GLY选择的的量子化学方法是M052x/6-311++G#//M052x/6-31G#，ASP和GLU选择的量子化学方法是M052x/6-311++G#//M052x/6-31+G*。模型分子首先在低级别(M052x/6_31G** 或M052x/6-31+G*)气相下进行构型优化，然后用更高级别的方法(M052x/6_311++G#)在液相下计算单点能量。量化计算选用的水溶剂模型是IEFPCM，原子半径用的是UAKS联合原子拓扑模型，力场计算用的水溶剂模型是A. Onufriev, D. Bashford和D. A. Case发展的GB模型,对应AMBER9中的关键词igb = 5,原子半径采用的mbondi2模型。在新的力场方法中，ALA、GLY、ASP和GLU用的是各自的势能面，除PRO外，其它的氨基酸用的是ALA的势能面。电荷拟和用的分子模型是氨基酸的二肽结构(ACE-XXX-NME，X = 20个氨基酸之一)，对紧凑和扩展构象选取了两种结构右螺旋(a K，= -57. 0, ￥ = -47. 0)和P -折叠(P，4) =-119. 0, V = 113. 0),侧链的二面角从Dunbrack等人统计出来的数据库中获得，当主链的值确定后，侧链值选用在时布局数最大构型的二面角值，分子构型优化以及静电势的计算在M052x/6-31G**的级别、气相下进行。在用RESP方法拟和电荷时，每个氨基酸的电荷总量限制为0或±1，ACE和NME电荷总量限制为I，模型中具有相同化学环境的原子设定具有相同的电荷。
具体实施方式
我们以氨基酸的四肽链结构(ACE-ALA-XXX-ALA-NME，XXX = 20个氨基酸之一)为模型，选取了五种具有代表性的结构(右螺旋(a K，0 = -57. 0, ￥ = -47. 0)，左螺旋(a L,
= 57. 0, ￥ = 47. 0), ￠-折叠，小=-119.0，￥ = 113.0)，反平行 3 -叠片(3a，小=-140.0，￥ = 135.0), PPII (PPII, =-79.0，V = 150. 0)，侧链二面角的选取和二肽模型侧链二面角的选取相同)，以MP2/cc-pVTZ计算的结果为标准来检验力场的准确性。为了对比同时用不同的量子化学和其它版本的AMBER力场方法计算了这20个氨基酸5种构型的相对能量。在密度泛函方法M052x和B3LYP的计算中，模型分子的构型首先在气态用6-31G**基组进行优化，其中ASP和GLU模型优化用的是6-31+G*基组，优化时把主链和侧链的二面角固定。然后在液态用6-31G**，6-311++G**和cc-pVTZ基组做单点校正，MP2/cc-pVTZ方法下用M052x/6-31G**优化的构型做了单点能量计算。液态计算采用的水溶剂模型是IEFPCM，构造空穴时选择的是UAKS联合原子拓扑模型。所有的量子化学计算用Gaussian 03软件完成。力场计算用的是AMBER9软件包，选择的力场有AMBER94，AMBER96，AMBER99，AMBER99SB和AMBER03，溶剂模型用的是GB模型(igb = 5)，原子半径采用的是mbondi2模型。优化时把主链和侧链的二面角固定。在数据处理时计算了两种均方根差RMS，每个氨基酸的五种构型的相对能量分别计算然后在20个氨基酸中平均；RMS-C，对20个氨基酸的同一种构型进行进行计算，然后对5种构型平均。RMS的计算公式如下
权利要求
1.一种用于建立蛋白质体系分子模拟力场的方法，该方法中力场的建立以氨基酸的四肽链结构为模型，选取了五种具有代表性的结构，其特征在于首先将模型分子的构型在气态用6-31G**基组进行优化，然后在液态用6-31G**，6-311++G料和cc-pVTZ基组做单点校正，MP2/cc-pVTZ方法下用M052x/6-31G**优化的构型做了单点能量计算。液态计算采用的水溶剂模型是IEFPCM，构造空穴时选择的是UAKS联合原子拓扑模型。
全文摘要
本发明涉及一种用于建立蛋白质体系分子模拟力场的方法。该力场的建立以氨基酸的四肽链结构为模型，选取了五种具有代表性的结构，模型分子的构型首先在气态用6-31G**基组进行优化，然后在液态用6-31G**，6-311++G**和cc-pVTZ基组做单点校正，MP2/cc-pVTZ方法下用M052x/6-31G**优化的构型做了单点能量计算。液态计算采用的水溶剂模型是IEFPCM，构造空穴时选择的是UAKS联合原子拓扑模型。用这种的方法得到的力场，准确性大大提高，计算的RMS值明显小于其它的力场，接近QM方法的M052x的结果，计算的RMS-C值表明新的力场方法克服了原有力场中构象偏向的缺点。
文档编号G06F19/12GK102779239SQ201110117428
公开日2012年11月14日 申请日期2011年5月9日 优先权日2011年5月9日
发明者吴春, 姜金良, 汪志祥 申请人:中国科学院研究生院