2.1极化力场概论
利用经典力学进行分子模拟,其结果的准确性关键取决于力场的质量,目前大部分常用力场,如OPLS[67~71],CHARMM[26,50~51],AMBER[22,24~25,49],MMFF[59~63]和GROMOS[72]等,都不同程度上存在着理论模型的局限,一个共同的不足就是,它们采用固定的原子中心固定点电荷来计算静电相互作用势能,把分子中的原子根据一定的化学环境进行分类,各类型电荷参数通过拟合气态模型分子的量子化学数据和一些实验结果得到。如果把溶质放入类似水的介电常数较大的溶剂中,或者当带有较大电荷的离子接近一个中性分子时,都会发生强烈的静电极化现象,这样的极化作用会影响体系的结构和能量,这时利用固定电荷静电势模型无法体现这种极化影响。Brooks和Forsman通过计算得出结果:氯离子与水分子簇的离解能不等于两体能的差值,从头计算的结果与力场有一定的偏差[73]。这样的弊端在凝聚态的模拟当中表现极为突出,以往为了减小误差,固定电荷力场在处理相关体系时,往往利用系数标度或者增加气态时原子点电荷的方法[74~76]得到比较平衡的溶剂-溶剂,溶质-溶剂和溶质-溶质静电相互作用势,尽管做了一定的校正,但是介电环境、极化现象广泛地存在于各种体系中,在气液界面和生物分子与溶剂水的界面处表现尤其突出,这一缺陷就不能以简单的方式来处理,许多文献已经报道和分析了在各种物理化学体系中准确处理极化的重要性,所以正确的处理静电极化是非常必要的[245~247]。
为了克服传统固定电荷力场的弊端,发展极化力场已成为近年来力场发展的焦点,近20年来得到了快速的发展,极化力场从方法上大体可分为四种类型:(1)诱导偶极(多极)模型(induciblepointdipolemodel);(2)浮动电荷模型(fluctuatingchargemodel,FQ);(3)浮动电荷与诱导偶极相结合模型(FQ-DP);(4)DrudeShell(dispersionoscillator)模型;DrudeShell(dispersionoscillator)模型近来引起了人们一定的重视,该模型已经应用到了液态水和离子水溶液体系。但是,无论何种方法都应该满足极化力场的目的和要求:(1)极化力场中的电荷必须随着化学环境的变化而随之改变,例如一些像生物大分子一样自由度很大的分子,分子间的极化强烈地依赖于其分子的构象和空间构型;(2)能够正确地计算多体相互作用,如气态分子团簇和溶液中分子间的极化作用;(3)极化力场的参数应具有良好的可转移性。(4)能准确计算静电长程相互作用能。今年来,各种极化力场在许多领域展开了初步的探讨,尤其在水和一些小分子体系进行了充分的研究,下面我们就对在各种体系中极化力场的发展作以简单的介绍。
由于液态水具有很好的氢键网络,在大多数的生物过程中都起着重要的作用,可以作为很好的体系来测试模型,所以,大多数极化力场都首先把水作为研究对象。
Caldwell等人的诱导偶极模型[77]采用各向同性的原子极化度,分子极化度等于各原子极化度的简单加和。Bernardo等人把极化区域限制在1,2或1,3键连的原子范围内[78];有的力场更是直接地把等于实验值1.4443的单极化度点放在氧原子或H-O-H的角平分线上[79~81]。这些诱导偶极力场大多都采用以实验为依据的可转移分子间作用势(TIP3P/TIP4P)的几何构型,或者是理想的简单四面体点电荷几何构型,H-O-H键角为109.5,但是,从目前的结果来看,还不清楚哪种方案更合理,通常极化力场通过拟合量子化学数据获取参数,TIP几何构型是较好的选择。Jedlovsky和Richardi通过比较三种水模型[79-81],结果发现,极化水模型获得的发散常数比TIP4P和SPC水模型更接近实验值[82];Sorensen等人通过X-ray衍射实验结果比较了两种力场模型下的水的径向分布函数,发现Chialvo-Cummings(CC)[79]和TIP4P-FQ模型[84]相对TIP3P和SPC三点模型更符合实验结果,但是不及TIP5P模型,这也表明,加入氧的孤对电子或非中心电荷会得到更好的结果[83]。浮动电荷极化力场(FQ)[84]
通过拟合水二聚体、三聚体等团簇从头计算方法的二体和三体相互作用能来获取极化参数[85],这种模型有一个最大的优点是高效性,缺点在于极化只能是限制在水分子的平面,虽然计算的介电常数很合理,但是它们的极化度是经过拟合得到的,而不是预先设定的。
极化力场在溶液和溶液表面也取得了一定的进展,Dang等人通过DC模型[81],发现水分子的平均偶极矩从液态接近并通过气液交界面的过程中,从2.75Debye下降到1.85Debye,有趣的是,在这个区域的几个埃范围内,这种过渡是平滑的,在气液交界面下几个埃的水分子的平均偶极矩有明显的减小;第二个研究发现[86],随着水分子进入有机相中,在水和CH2Cl2的界面处,水的平均偶极矩下降30%,类似的,CH2Cl2分子进入水中时平均偶极矩下降大约10%,同样的现象在水和CCl4界面也有类似现象发生[87]。第三个研究发现[88],(Dang等人)在碘离子和水组成的团簇中碘离子总是倾向于在团簇的表面,而不是在团簇的中心,如果使用非极化力场,这种现象就会消失。第四个研究发现[89],(Dang等人)与以往的非极化力场不同的是,在氯离子或者铯离子通过水和CCl4界面时没有自由能的最小值出现,表明它们是没有表面活性的;另外,他们发现[90],当苯分子通过气液界面时却有-4kcal/mol的自由能最小值出现,证明苯是具有表面活性的。
极化力场对小分子体系也做了开拓性的研究[91~93],例如CCl4、CH2Cl2、苯、胍离子、甲醛、胺、酰胺、N-甲基乙酰胺等分子,这些研究都表明,加入极化使电荷分布等结果与量子化学计算结果更加吻合,如果忽略了极化效应,在氢键的静电相互作用方面将导致一定的误差。值得注意的是,Levy等人计算的胺和酰胺的水合自由能,利用甲基连续的替换氨基上的氢,从氨水到三甲基氨,其自由能差为3kcal/mol,尽管仍旧比实验值1.1kcal/mol大,但是比固定电荷力场的4.4~6.6kcal/mol的结果准确性提高了很多。
金属蛋白质和金属酶的研究一直是力场研究的重要领域,所以极化力场也对其进行了重新的探讨。Gresh和他的合作者自1995年以来已经发表了十多篇关于SIBFA方法的论文[94~100],大部分都是用来建立和参数化应用于蛋白质金属主型配位化合物的力场,如Zn2+,对近距离的静电相互作用采用了屏蔽势能,他们把极化度和vanderWaals质心位点都放置在键和孤对电子上,而不是在原子上。SIBFA方法又加入了极化作用,来研究在配位键形成的过程中从蛋白质配体到金属电子转移对几何构型和能量的影响,虽然已经得出了一些重要的结果,但是它存在两个不足:第一,在分子动力学模拟过程中,不能计算解析梯度;第二,刚体近似,在能量优化的过程中键长和键角值是被固定的,势必会影响计算结果的合理性。
2.2应用于蛋白质模拟的极化力场的发展
在讨论应用于蛋白质体系的极化力场时,应该首先回顾一下传统固定电荷力场在该领域的发展历史,它们对此已经作出了突出贡献,为极化力场的发展奠定了稳固的基础。对蛋白质体系的分子动力学和MonteCarlo模拟已经有大约近30年的历史。1984年,Weiner等人建立AMBER力场[24],他们利用实验测得的肽片段的光谱数据来确定化学键伸缩和键角弯曲势能函数参数,开始时他们利用United-atom(以原子团为基本单元)模型,但是这一模型在处理芳香环、受力中心和计算光谱等方面存在着一定的弊端,后来随着计算机运算速度的加快,发展了All-atom(以原子为基本单元)模型,又利用从头计算6-31G基组下RESP静电势方法拟合净电荷的方法代替了HF/STO-3G小基组方法,结合小分子的液态模拟数据来确定Lennard-Jones势能参数,利用精密的从头计算方法来重新拟合二面角扭转势能参数,他们在几个重要的方面都作了突破性的改进,使AMBER力场能够更好地模拟蛋白质体系。几乎同AMBER力场同时发展的CHARMM力场[26,50~51],该力场更重视蛋白质在溶液中的模拟,他们利用RESP静电势方法拟合模型分子同水分子的二聚体的原子净电荷确定电荷参数,同时他们在拟合的过程中还考虑到了模型小分子二聚体的能量和偶极矩等物理量来最终优化电荷参数。他们为了减少非键长程相互作用的计算量,限制分子中基团呈电中性。该力场更重视溶液中性质的模拟,所以他们在二面角势能和其它参数确定的情况下大量地参照液态模拟的结果。CHARMM力场能更好地确定蛋白质在溶液中的结构。在AMBER和CHARMM力场之后,OPLS力场[67~71]也在AMBER力场的基础上发展起来,该力场能更好地计算液态的热力学物理量,后来他们完全拟合精密从头计算结果来拟合二面角扭转势能参数,建立了完全的蛋白质力场,能很好地重复从头计算的多肽构象结构和相对能量,在对实际蛋白质体系的模拟中也取得了不错的结果。除了以上三个力场之外,如GROMOS[72]、MMFF[59-63]和MM力场[7,17,27~28,44~48],在蛋白质体系的模拟中也取得了一定的进展。
以上蛋白质力场同大多数的力场一样,主要的弊端就在于没有充分考虑静电的极化影响,而这一因素在蛋白质和其溶液体系中是至关重要的,因为在蛋白质分子中存在着大量的极化基团和氢键网络,在溶液中溶剂水分子与溶质蛋白质之间的相互极化也普遍存在,因此对该体系必须考虑其极化效应,发展极化力场模拟蛋白质体系势在必行。极化力场已经对有机小分子和水进行了比较深入的研究,对多肽和蛋白质等生物大分子的研究也有了一定的进展,但仍旧还不够具体和全面。Kollman和其合作者创建的AMBER力场,根据Applequist模型加入了点极化度,通过抵消和减少一定的固定电荷,一般是减少到原来固定电荷的88%,该种方法应用到一般的有机分子时十分方便,但是在一定程度上,如果直接采用原来力场的其它参数,就会影响构象能的计算[101-102]。根据Banks和Stern提出的极化力场方法[103-104],Friesner、Berne和其合作者发展了应用于多肽和蛋白质的极化力场,对20种氨基酸二肽的构象进行了研究,很好地重复了从头计算的结果,提供了一个完全的蛋白质极化力场,并能够很好地描述多体影响和计算多体能量[105]。Banks首先应用了浮动电荷方法[104],通过拟合量子化学静电响应数据来拟合参数;Stern等人在原来的基础上,应用浮动电荷和诱导偶极相结合的方法来研究多肽构象[103],同时克服了单独使用浮动电荷方法的极化仅限于苯分子平面方向的弊端,后来他们用该方法完成了整个蛋白质极化力场[105]。
他们已经发展了一套完整的自动参数化方案拟合量子化学计算结果产生静电和极化参数。近来,他们采用Buckinghamexp-6的vanderWaals势能函数,通过高水平的从头计算,利用小分子氢键二聚体拟合力场的vanderWaals参数,然后通过计算与氨基酸片段相关的有机分子溶液来调整其势能函数中的色散吸引项参数。初步的工作表明,该方法计算得到的汽化热结果与实验值相当吻合,同时也很好地重复了高水平从头计算的构象能。Gresh等人也发展了蛋白质和多肽极化力场[106],利用SIBFA(sumofinteractionsbetweenfragmentsabinitiocomputed)方法对甲酰胺、氮甲基乙酰胺二聚体和丙氨酸、甘氨酸残基二肽的氢键能进行了研究,研究中发现,加入孤对电子或者更高级次的极矩比简单原子中心电荷能更好地描述和计算酰胺基之间的氢键取向和能量;同时,Dixon和Kollman的研究发现[107],水分子与吡啶形成氢键时,水分子更倾向于平行于吡啶平面,而不是垂直方向,两者相差3~4kcal/mol,但是用标准的固定电荷的AMBER力场却只能得到十分之几个千卡,从而更进一步地证明非中心电荷在处理氢键问题的重要性。Patel,etal.等人在2004年发展了第一代CHARMM浮动电荷力场[108~109],主要用来研究蛋白质和多肽体系,其静电模型参数,如原子的电负性、硬度,他们利用虚拟的偶极来模仿水分子放置在模型小分子的周围,通常是范德华半径处,从而利用密度泛函方法计算模型小分子的电荷响应,在不同方向和不同的距离放置虚拟偶极,以得到大量的计算数据,从而拟合以上参数。利用从头计算的气态的水分子和模型分子的二聚体和多聚体的结构和能量,以及对模型分子液态性质的模拟,来进一步优化其非键参数。