談談分子模擬中的隱式溶劑模型與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是基于多層近似后的溶劑模型。