【3.3.5】gromacs的隐式溶剂模型和GB模型

本文主要讨论隐式溶剂模型和GB模型(广义Born模型)的各个方面。

隐式水模型优势主要有下:

  • 减少体系的粒子数从而降低计算开销。
  • 没有溶剂摩擦效应,可以加快构象变化。
  • 不显著、具体地表现溶剂运动而直接表现平均效果,便于计算自由能,不需要动力学很长时间采样。
  • 不需要显式水的预平衡过程。
  • 适用于常PH模拟、REMD等方法。
  • 适合计算能量面,结果比较光滑,真实的水会带来噪音,溶剂的运动、碰撞使主体分子产生许多能量局部极小点。等等……

总体来说,隐式水模型相对于显式水模型,是一个离散化到连续化的过程。溶剂模型的近似关系如下(越左边越真实,越右边近似程度越大): 显式溶剂模型–>隐式溶剂模型–>明确拆分为静电、非极性贡献–>PB–>GB

分子在溶剂环境下的总能量Etot=溶剂化能deltaG(solv)+溶质在真空条件下的能量Evac。(对于极化力场,溶剂对溶质的极化作用造成的能变也算进在deltaG(solv)里。) deltaG(solv)=deltaG(elec)+deltaG(nonpolar)=deltaG(elec)+deltaG(vdw)+deltaG(cavity)

注:有时明确包含氢键项。

*其中deltaG(elec)就是溶剂环境下的静电作用能与真空下的静电作用能之差,主要包括溶剂与溶质间的静电作用,以及溶剂对溶质原子间静电作用的屏蔽效果。 其中deltaG(nonpolar)就是除此以外,即体系各原子皆无电荷时,溶剂环境与真空环境下能量之差。它又分为deltaG(vdw)和 deltaG(cavity)。deltaG(vdw)是溶质与溶剂的范德华作用,由于范德华作用随距离衰减快,所以只考虑第一溶剂化层的溶剂,此项造成能量降低。deltaG(cavity)表现的是溶质反抗溶剂压力形成洞穴的做功以及溶质周围溶剂(第一溶剂化层为主)的重构造成的熵罚,此项造成能量升高。

因为deltaG(nonpolar)的两项都主要涉及第一溶剂化层,它正比于SASA,故最简单、常用的表达是a*SASA+b,而比例系数a由拟合实验值得到,b一般为0。另有其它方案计算,比如包括溶质体积、溶质直径等信息在内的方程。也有比如deltaG(vdw)通过SAS的积分得到,而 deltaG(cavity)仍用正比于SASA的方案。在GB模型下的MD模拟中,对于原始状态的生物大分子,由于SASA变化很小,原子因非极性作用导致的受力往往可以忽略,若需计算的话则是求它对坐标的导数。由于a*SASA+b的a为正值(一般为7.2cal/mol/Å^2),可见 deltaG(nonpolar)的效果是减小溶质的SASA,此项体现了疏水效应。

deltaG(elec)项一般通过GB或PB方法计算,GB速度快但一般准确度低于PB。目前常用的GB模型计算方程由Still于1990年提出(原文 http://www.brsbox.com/filebox/down/fc/90f984299e8583a61e0f33942bd13260 ),依据是库仑定律和born方程,推导出另一种形式。当粒子相离较远时,方程等价于库仑定律+born方程能量,而相距为0时等价于born方程,而当距离较近时等价于Onsager反应场能量。方程中还可以再引入单价粒子的屏蔽校正项。

GB方程中重要参数是有效born半径。某个原子电荷与溶剂作用对deltaG(elec)的贡献实际上是把分子中其它原子电荷去掉后,溶剂化状态与真空状态能量差。如果用的有效born半径以GB方程算的结果与之一样,就是正确的born半径。或者粗略来说:某原子用有效born半径算得的GB式中的“自身项”的能量应当与真实的分子中此原子与溶剂的静电相互作用能一致。这样,有了正确的born半径,就可以用GB式得到合理的 deltaG(elec)。有效born半径对计算deltaG(elec)有重要影响,模拟过程中分子构象不断改变,依赖于它的有效born半径也不断改变,需要每一步重新计算,若构象变化不明显可每隔一定步数重新计算,这是GB计算量中的主要部分之一。计算有效born半径在不同程序中有不同方法。一般使用CFA(库仑场近似),有效born半径就可以对溶质分子表面内的空间进行积分获得。而分子内空间的范围不便于描述,简单的方法是直接将原子VDW球内区域叠加作为分子内空间,称GB-HCT,但这样就忽略了原子VDW空间之间的缝隙,可以乘以4/3来近似校正。也有人根据原子被埋程度的概念提出了基于经验的简单、快速的函数,其中引入可调参数,在amber里称为GB-OBC方法,结果明显优于GB-HCT,应注意可调参数是通过预先优化得到的,面向不同体系模拟应使用不同参数。amber中还支持更新的GBn方法,结果优于GB-OBC,适合蛋白而不适于核酸体系。

GB方程的优点是形式简单,计算快速,结果较准确,而且有简单的解析导数,可直接计算GB环境下原子在MD过程中的受力。GB可以结合分子力场及量子力学模拟各种分子体系的溶剂化效果。可以用在类似MM/PBSA的结合自由能计算的静电项中,即MM/GBSA。一些模拟方法则只能使用GB,比如 amber中可以实现的结合MC的常PH模拟,原理是根据Metropolis判据在MD过程中动态决定是否某点被质子化或去质子化。对于REMD,因为随着体系粒子数的增加,需要的副本数目飙升(否则能量重叠较差,交换概率低而起不到效果),所以总是结合GB方法来显著降低粒子数量。对于不大的体系,GB比显式溶剂模型有更快的速度。

但GB方法也有缺点。与显式水模型仍有一定差距,将水分子进行了“连续化”的近似,故不能表现与溶剂的强相互作用、特殊相互作用(如水桥)。对于膜体系,溶剂、溶质、膜都有着不同的介电常数,这样复杂的情况介电常数环境下GB也难以表达。对于电荷较多的核酸体系,也不能较好表达多价抗衡粒子与它的相互作用。一些研究也表明用GB模拟多肽会产生与显式溶剂模型不同的构象,与实验结果有一定偏差。计算deltaG(elec)项时,相对于原理严格的PB 也有差距(即便使用最佳的有效born半径),但是好处是GB计算速度快得多。不过PB的结果亦不可作为绝对的金标准,因为一些误差在PB引入的近似中就已经出现了。在计算速度上,显式溶剂体系计算量是O(N^2),但使用Ewald等方法后,可降至Nlog(N)以下,而GB模型如果在cutoff上不做近似(往往取无限大),对于大体系反倒更慢。

除GB、PB以外也有其它一些方案,最简单的隐式模型是介电常数依赖于作用距离的方法,其中库仑作用的介电常数等于粒子间距离(即静电作用能1/r 变为了1/(r^2)),以体现溶剂的屏蔽效果,此方法精度显然低于GB。近来还有一些新方法,如AGBNP(Analytical Generalized Born plus Nonpolar),在GB基础上引入另外形式的deltaG(nonpolar)项。ALPB(analytical linearized Poisson-Boltzmann),其方程类似GB(在无限介电常数下回归为GB方程),故有着基本一致的计算速度,但引入了有效分子静电尺寸参数,模拟水比GB有着更好的精度,结果更接近于PB。这些新模型理论上有着更好的结果,但仍需广泛地应用于模拟来检验。还有隐式与显式溶剂模型相结合的方法,也就是在主体分子外面包一层溶剂分子,往往以弹簧势限制住溶剂避免跑走,而更外面则用隐式溶剂模型,结合了两种方法的一些优点。

目前Amber对隐式溶剂模型支持较好,支持GB、ALPB,也有内建的PB模块pbsa(老版本挂的是第三方的delphi)。而Gromacs在最新版本4.0.5中尚不支持GB/PB,需要外接第三方软件如APBS,但在接下来的版本中即将加入GB。 对于GB模型,适用领域与不适用领域的界限尚为模糊,一般来说,对于必须用GB的地方,或者已证明有明显优势的地方可以使用,但如果不确定,尽量还是用显式溶剂模型,毕竟GB是基于多层近似后的溶剂模型。

个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn

Sam avatar
About Sam
专注生物信息 专注转化医学