谈谈分子模拟中的隐式溶剂模型与GB模型

本文主要讨论分子模拟中隐式溶剂模型和GB模型（=Generalized Born，广义Born模型）的各个方面。


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

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

溶质分子在溶剂环境下的总能量Etot=溶剂化能deltaG(solv)+相同构象但是在真空中的能量Evac。
deltaG(solv)=deltaG(elec)+deltaG(nonpolar)=deltaG(elec)+deltaG(vdw)+deltaG(cavity)

注：有时明确包含氢键项等近程强相互作用，但一般算进deltaG(elec)。另外，在真空和溶剂环境下，相同构象的溶质的振动模式改变（以低频为主）对配分函数的影响而导致的自由能变也应当加进去，但是由于难以计算，比较耗时，结果亦不准确，而且对自由能贡献很小，所以一般不考虑。至于在两种环境下同构象溶质的平动、转动变化对配分函数的影响以及P*delta(V)项引入的自由能变，则完全忽略。对于可极化力场，溶剂对溶质的极化作用造成的能变也算进在deltaG(solv)里。

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

因为deltaG(nonpolar)的两项都主要涉及第一溶剂化层，它正比于SASA（溶剂可及表面积），故最简单、常用的表达是a*SASA+b，比例系数a（有时也称为表面张力参数）由拟合实验值得到，一般为7.2cal/mol/Å^2，b一般为0，这样计算结果可以重现碳氢化合物的水和能。也有其它一些隐势溶剂模型使用的表达式是∑a(i)SASA(i)，即不同类型原子使用不同的表面张力参数，然后乘以每个原子的SASA，最后加和，结果比全部原子使用相同表面张力参数一般更为准确。另有其它更严格的方案计算，比如包括溶质体积、溶质直径、溶剂密度等信息在内的方程，还有比如deltaG(vdw)通过对SAS的积分得到，而deltaG(cavity)仍用正比于SASA的方案等等。在GB模型下的MD模拟中，对于原始状态的生物大分子，由于SASA变化很小，原子因非极性作用导致的受力往往可以忽略，若需计算的话则是求它对坐标的导数。由于a*SASA+b的a为正值，可见deltaG(nonpolar)的效果是减小溶质的SASA，此项体现了疏水效应。

deltaG(elec)在量子化学中使用自洽反应场方法(SCRF)精确计算，溶剂效应对溶质的影响体现在对哈密顿算符添加的微扰项，会使溶质电子云被极化，它对溶剂产生影响，又使微扰项改变，故需要反复迭代求解。在一般的分子模拟中不需要太高计算精度，也为了运算速度，并不直接表达出电子，而是将原子抽象为一个数值不变的点电荷，多数情况也不再用额外方案考虑极化效应，故溶剂化自由能可以一步解出。最常用的方法是GB和PB方法计算，GB速度快但一般准确度低于PB。

Born方程是指在某种介电常数ε的溶剂中，使其中半径为a的粒子带上q电荷所引起自由能的改变为0.5*q^2/ε/a。Born方程加上库仑静电能就得到了体系静电能，求处在真空与处在ε溶剂中的静电能差，就得到了广义Born方程：-0.5*(1-1/ε)*∑[i]q(i)^2/a(i)-(1-1/ε)∑[j>i]∑[i]q(i)q(j)/r(ij)，其中∑[i]代表对i下标的项加和。但是这样计算的结果高估了溶剂对处于分子内部原子之间相互作用的影响。目前常用的GB方程由Still于1990年提出（原文JACS,112,6127），依据是库仑定律和Born方程，但以另一种相似的形式出现。在这个方程中，当粒子相离较远时，方程等价于上述原始GB方程能量，相距为0时等价于Born方程，当距离较近时等价于Onsager反应场能量。方程中还可以再引入单价粒子的屏蔽校正项。GB方程得到的deltaG(elec)加上a*SASA得到的deltaG(nonpolar)来计算溶解自由能称为GB/SA方案，但往往直接简称为GB方案。GB方程使用的原子电荷q适合使用拟合静电势方法来得到，即ESP电荷。

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方程也被用于量子化学的溶剂模型中，称为SMx系列方法，它将GB方程作为微扰项，电荷以Mulliken布居、CMx等方法得到。

但GB方法也有缺点。与显式水模型仍有一定差距，将水分子进行了“连续化”的近似，故不能表现与溶剂的强相互作用（如氢键）、特殊相互作用（如水桥）。对于膜体系，溶剂、溶质、膜都有着不同的介电常数，这样复杂的环境下，只用单一介电常数的GB模型也难以表达。由于GB模型下的模拟并不抗衡离子，所以对于带电体系，不能表达抗衡粒子与它的相互作用，尤其是对于多价离子与带电荷较多的核酸间的作用。一些研究（如PNAS,vol.99,12777)表明用GB模拟多肽会产生与使用显式溶剂模型不同的构象分布，与实验结论有较大差异，主要是由于在GB下体系倾向生成盐桥使静电相互作用增强，同时削弱了疏水作用导致易破坏疏水核，这和缺乏抗衡离子的静电屏蔽有关，有人提出将电性相反残基间介电常数加大或加入惩罚函数解决。计算deltaG(elec)项时，相对于原理严格的PB也有差距（即便使用最佳的有效Born半径），但是好处是GB计算速度快得多。不过PB的结果亦不可作为绝对的金标准，因为一些误差在PB引入的近似中就已经出现了。在计算速度上，显式溶剂体系计算量是O(N^2)，但使用PME等方法后，可降至O(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是基于多层近似后的溶剂模型。