文章_可極化力場（PPC）
摘要
對溶液中的蛋白質進行了從頭量子力學計算，以產生極化的特定蛋白質電荷（PPC），用於蛋白質的分子動力學（MD）模擬。通過開發基於片段的量子化學方法和隱式的連續溶劑模型，使蛋白質的量子計算成為可能。蛋白質的計算電子密度被用來推導PPCs，這些PPCs代表了蛋白質在原始結構附近的極化靜電狀態。這些PPCs與標準力場中的PPCs一樣是以原子為中心的，因此在計算上對蛋白質的分子動力學模擬具有吸引力。為了研究電子極化對硫氧還蛋白的結構和動力學的影響，已經進行了廣泛的MD模擬。我們的研究表明，通過明確比較使用PPC和AMBER指控的MD結果，硫氧還蛋白的動力學被電子極化所穩定。特別是，使用PPCs的MD自由能計算準確地再現了埋在硫氧還蛋白內的可電離殘基Asp26的pKa移位的實驗值，而以前使用標準力場的計算則高估了pKa移位的兩倍之多。通過嚴格的MD自由能模擬準確預測埋藏在蛋白質內部的可電離殘基的pKa移位，幾十年來一直是計算生物學的一個重大挑戰。這項研究提出了強有力的證據，證明蛋白質的電子極化在蛋白質動力學中起著重要作用。對溶液中的蛋白質進行了無源量子力學計算，以產生極化的蛋白質特定電荷（PPC），用於蛋白質的分子動力學（MD）刺激。通過開發基於片段的量子化學方法和隱含的連續溶劑模型，使蛋白質的量子計算成為可能。蛋白質的計算電子密度被用來推導PPCs，這些PPCs代表了蛋白質在原始結構附近的極化靜電狀態。這些PPCs與標準力場中的PPCs一樣是以原子為中心的。為了研究電子極化對硫氧還蛋白的結構和動力學的影響，已經進行了廣泛的MD模擬。我們的研究表明，通過明確比較使用PPC和AMBER電荷的MD結果，硫氧還蛋白的動力學被電子極化所穩定。特別是，使用PPCs的MD自由能計算準確地再現了埋在硫氧還蛋白內的可電離殘基Asp26的pKa移位的實驗值，而以前使用標準力場的計算則高估了pKa移位的兩倍之多。通過嚴格的MD自由能模擬準確預測埋藏在蛋白質內部的可電離殘基的pKa移位，幾十年來一直是計算生物學的一個重大挑戰。這項研究提出了強有力的證據，證明蛋白質的電子極化在蛋白質動力學中起著重要作用。
引言
蛋白質中相互作用的一個重要組成部分是靜電相互作用。蛋白質的摺疊、蛋白質與配體的結合、蛋白質與蛋白質的相互作用、電子轉移、質子的結合和釋放、酶的反應等過程主要是由靜電相互作用驅動的。特別是，質子結合是一個依賴於pH值的過程，並受到與蛋白質和區域性環境的靜電作用的強烈影響。蛋白質側鏈之間的靜電作用不僅取決於它們的距離，而且還取決於它們在蛋白質中的位置和它們的區域性溶劑環境。在質子結合中，蛋白質和溶劑經歷了涉及電子極化和原子團位移的介電鬆弛。蛋白質內部的區域性靜電環境是不均勻的和疏水的，這通常對電離不太有利。因此，與蛋白質表面或溶液中的孤立氨基酸相比，蛋白質內部的可電離（或帶電）殘基可能有很大的pKa偏移。
目前的標準力場，如CHARMM和AMBER，在應用上取得了巨大的成功，但存在著根本的侷限性。具體來說，這些FFs是氨基酸特異性的或類似於均值場的，因此不能準確表示高度不均勻和蛋白質特異性的特定蛋白質環境的靜電作用。例如，在相同或不同的蛋白質中的兩個相同型別的氨基酸，由於其不同的靜電環境，應該有相當不同的電荷狀態。而目前基於氨基酸的FFs無法描述特定蛋白質結構的極化狀態，例如，原生蛋白質結構。為了克服標準FFs的這一基本缺陷，人們努力開發肽和蛋白質的可極化力場模型。
在這篇文章中，我們提出了一個包含極化的特定蛋白電荷(The polarized protein-specific charge) 的替代力場，用於蛋白質的動力學研究。極化特定蛋白質電荷（PPC）在原子電荷中建立了蛋白質的極化。PPC本身並不是 "可極化的"，但它是由第一原理的量子溶解計算得出的原生（或特定）結構中的蛋白質。因此，PPCs正確地代表了蛋白質的電子極化狀態，因此提供了接近原生結構的準確的靜電作用。在我們的方法中，蛋白質溶解的量子化學計算是通過將最近開發的基於片段的方案，即molecular fractionation with conjugate caps （MFCC）
為了研究PPCs在蛋白質的MD模擬中可能產生的影響，我們進行了一系列的MD研究，研究了硫氧還蛋白的各種結構和動力學特性。使用PPCs的新的MD結果與相應的使用了AMBER電荷MD研究中得到的結果進行了明確的比較。特別是，進行了MD自由能模擬以預測埋藏Asp26的硫氧還蛋白的pKa移動。以前使用AMBER和CHARMM進行的MD自由能計算給出的pKa偏移是實驗值的兩倍。為了與這些結果進行直接比較，我們按照Simonson等人的計算方法和程式來計算pKa移位，只是用我們計算的硫氧還蛋白PPC取代標準AMBER電荷。這將消除結果比較中任何可能的動態不確定性。
理論方法
在連續體-溶劑模型中，溶質(蛋白質)由電荷分佈ρ(r)表示，它嵌入到一個空腔中，周圍是介電常數為ϵ的可極化介質。溶質電荷分佈ρ(r)使介電介質極化，併產生一個反應場，該反應場反作用於溶質，直到達到平衡。根據經典的靜電理論，作用於溶質的反應場可以有效地由腔體表面的感應電荷來表示。通過離散化空腔表面的誘導電荷並迭代求解表面電荷產生的外部反應場中的溶質的量子化學方程，可以得到流行的PCM方法，最近被推廣到蛋白質的溶解。然而，對於大的蛋白質，PCM模型需要許多離散的表面電荷，從而使線性方程的解決在計算上變得困難。在目前的方法中，我們通過數值求解泊松-波爾茲曼（PB）方程來獲得反應場。這避免了大型線性方程的求解。
在我們的方法中，擬合蛋白質原子電荷的基本程式可以描述如下。首先，用MFCC方法進行蛋白質的氣相計算，以獲得給定結構的蛋白質的初始電子密度。使用RESP程式將計算出的電子密度用於擬合原子電荷。這裡使用的電荷擬合理念與AMBER力場中使用的相同，這就保證了PPC與AMBER力場的其他引數一致。然後進行PB方程的求解，得到反應場，從中產生腔體表面的離散表面電荷。
蛋白質的MFCC/RESP電荷擬合方法與Poison-Boltzmann連續溶劑的耦合應提供蛋白質在水中的溶解自由能的準確估計。由MFCC/RESP程式產生的蛋白質每個原子上的Partial charges被傳遞給PB求解器DELPHI以確定自洽反應場。在電介質邊界上得出一組誘導表面電荷qind，代表溶劑分子的反應場效應。介質溶質/溶劑邊界由AMBER範德瓦爾斯半徑定義，溶質分子原子的探測半徑為1.4 Å。內部介電常數，表示為ϵsolute，被設定為unity，因為分子的極化性明確地包含在量子力學計算中。溶劑介電常數，表示為ϵsolvent，被設定為80。網格密度被設定為4.0網格/Å。然後，在接下來的量子力學計算中，每個被蓋住的片段（CF）的表面電荷被新增為背景電荷。如上所述，蛋白質的部分電荷和篩選的表面電荷相互極化，直到達到收斂。只有當蛋白質的偶極子和表面電荷都收斂到一定的數值精度內時，這個迴圈才會停止。
計算公式見原文。
從圖1的流程圖中可以很容易理解生成PPCs的過程。對於本文報道的結果，單個蛋白質片段的電子密度的量子化學計算是在B3LYP/6-31G*水平進行的。
結果與討論
蛋白質的溶化和極化
上一節所述的MFCC-PB計算協議需要進行廣泛的數值測試以驗證其功效。特別是，我們需要研究由PPC代表的蛋白質極化對蛋白質溶解的影響。我們進行了基準研究，考察了一些蛋白質系統的電子極化對溶解能的貢獻。一些蛋白質系統的靜電溶解能的計算貢獻列於表1，為了比較，之前用MFCC-CPCM方法計算的結果也列於表中。與經典的PB或廣義Born（GB）方法不同，溶質極化能在計算的靜電溶存自由能中是相當重要的。我們的計算表明，如表1所示，這個項對總的靜電溶解能的貢獻是7∼19%。極化能的這一重要部分表明，在目前的力場中需要包括極化效應。
為了研究蛋白質極化對各種蛋白質特性的影響，我們還將目前PPC計算出的溶液中蛋白質的偶極矩與AMBER03力場中的偶極矩進行了比較。由於大多數蛋白質具有非零淨電荷，而計算出的偶極矩取決於帶電實體在座標系中的位置，因此，蛋白質首先以其電荷中心為中心。電荷中心定義為rc=（∑ri∥qi∥）/∑∥qi∥，rc=（∑ri∥qi∥）/∑∥qi∥，其中ri是蛋白質中第i個原子的原始座標，qi是其部分電荷。然後，整個蛋白質的偶極子由D=∑（ri-rc）qi.D=∑（ri-rc）qi.為了研究更詳細的區域性靜電環境，我們計算了蛋白質單個殘基的偶極矩。圖2顯示了一些蛋白質系統中由PPC計算的單個殘基的偶極矩與AMBER03計算的偶極矩的對比圖。比較結果表明，從PPCs計算出的大多數殘基的偶極矩一般都比從AMBER電荷的大。這可能表明，在MFCC-PB電荷中，殘基的極化程度更高。
Protein-specific的原子電荷
由於MFCC-PB電荷正確地描述了蛋白質的靜電極化，所以在蛋白質不同位置的同一型別的氨基酸上partial charge通常是不同的。一個特定殘基上的partial charge是由它的特定構象和由於蛋白質的其他殘基而產生的化學環境決定的。為了說明這一特點，我們在表2中列出了溶菌酶的一些區域性電荷。如同在溶菌酶中，我們發現同一型別的氨基酸上的電荷相互之間有很大的不同。例如，ASP中骨幹原子N的原子電荷值從-0.083到-0.642，Ser中骨幹原子O從-0.727到-0.556，Ile中原子CG1從0.025到0.304。這是PPC最重要的特點，它與標準的基於氨基酸的對一個給定的殘基有相同的固定原子電荷的力場不同。
硫氧還蛋白中Asp26/Asp20的pKa移位計算（過）
儘管標準的基於氨基酸的部分電荷模型在MD模擬中對蛋白質的許多巨集觀特性建模時效果很好，但在模擬對區域性靜電環境更敏感的特性時預計會有困難。這是因為標準力場電荷是類似於均值場的，它們不包含蛋白質的極化和上述其他特定的蛋白質靜電資訊。例如，準確預測蛋白質中埋藏殘基的pKa位移是目前MD模擬中的一個難題。準確計算pKa位移對於我們理解酶的酸鹼催化機制至關重要。這種酶的活性需要催化殘基以適當的質子化狀態存在。pKa計算中常用的方法是基於求解PB方程，其中溶劑被視為連續介質。PB方法是一種平均場理論，對於一些簡單的情況可以給出有用的見解。然而，PB方法在更復雜的情況下失敗了，因為它不是一個微觀模型，因此不能說明詳細的分子過程。為了準確地預測pKa移動，人們需要正確地說明影響質子結合過程的分子因素。因此，需要有明確包括水分子的微觀方法來從第一原理上正確預測pKa位移。然而，以前使用微觀方法的一些嘗試未能準確預測pKa位移。
Simonson等人最近進行了一項分子動力學自由能（MDFE）研究，計算硫氧還蛋白中Asp26的pKa位移，證明了目前標準AMBER和CHARMM力場的問題。用兩種力場進行嚴格的MDFE模擬，不同的執行結果都高估了Asp26的pKa偏移量，誤差為4-5千卡/摩爾。Simonson等人對pKa的MDFE計算可能是迄今為止最嚴格和最廣泛的研究，其結果可以作為所研究系統的一個基準。Simonson等人的結論是，這個問題一定部分來自於所採用的力場的系統誤差。在這篇文章中，遵循Simonson等人的嚴格的MDFE程式來計算埋藏的Asp26和表面暴露的Asp20的pKa偏移，但使用了本文之前描述的新的PPC。
我們的MDFE模擬的數值細節與Simonson等人的模擬密切相關。硫氧還蛋白的核磁共振結構（Protein Data Bank，即PDB，程式碼1XOA；並見圖3）被作為初始結構。模擬是用AMBER程式進行的，使用了具有周期性邊界條件的顯式水模型（TIP3P）。首先，系統在NPT集合下平衡，在298K下執行500-ps，然後在NVE集合下繼續模擬。SHAKE[44.](44)被用來約束與氫相連的共價鍵。長程靜電相互作用是由粒子網格Ewald方法處理的。時間步長為一飛秒。我們使用的模型化合物是帶有n-乙醯基和n-甲基醯胺阻斷基的天冬氨酸。質子化和電離狀態下的電荷分佈見圖4。
表3給出了埋藏的Asp26以及表面暴露的Asp20和相應模型系統的計算結果。我們計算的Asp26去質子化的自由能轉移是5.0千卡/摩爾，這非常接近實驗值4.8千卡/摩爾。模擬一直進行到18 ns，為蛋白質的所有中間狀態提供足夠的取樣。圖5顯示了電位相對於系統引數λ=0.5的能量導數與模擬時間的關係。這相當於Asp26的質子化和d質子化狀態之間的一個虛構的中間狀態。為了驗證我們的計算，我們還計算了Asp20的pKa位移，一個表面暴露的殘基。我們的模擬結果給出了1.0千卡/摩爾的自由能差，考慮到模擬的統計不確定性，這與實驗值的零是相當合理的一致。我們的結果表明，使用正確描述蛋白質中帶電殘基周圍的極化靜電環境的PPC，可以從第一原理的嚴格的MDFE模擬中得到pKa移動的準確預測。相比之下，同樣的數值模擬但使用AMBER和CHARMM電荷產生的pKa位移是Simonson等人中實驗值的兩倍。
為了估計統計誤差，我們進行了塊平均分析，其中每個持續時間為6ns的軌跡被分為四個塊。導數的統計不確定性被估計為塊平均的標準偏差的兩倍。從表3可以看出，自由能的統計誤差為1千卡/摩爾，計算的理論值在Asp26和Asp20的實驗資料的統計誤差範圍內。由於我們的MDFE計算與Simonson等人[30.](30)的計算完全相同，只是我們用PPC取代了AMBER原子電荷，同時保持力場中的所有其他引數，所以不同的結果一定來自電荷的影響。換句話說，PPC中包含的蛋白質的極化效應是準確預測Asp26的pKa轉變的原因。這一結果真正證明了電子極化在蛋白質微妙的靜電相互作用中發揮的重要作用。
上述結果意味著電子極化的缺乏可能是早期對埋藏在蛋白質內部的可電離殘基的pKa位移的一些MD自由能計算的主要罪魁禍首。當然，需要對不同的蛋白質系統進行更多這樣的計算來證實這一結論。當然，還有其他重要的動力學效應被提出來解釋pKa移動的MD模擬的失敗。例如，Kato和Warshel提出，決定埋藏基團pKa值的一個主要因素是內部基團被電離時發生的結構變化。然而，由於有限的MD取樣所帶來的問題，大型結構重組的可能性並不容易在計算上進行研究。例如，到目前為止，葡萄球菌核酸酶中Glu66的pKa移動的計算是一個巨大的挑戰。Kato和Warshel開發了一種特殊的過度充電方法，允許更多的結構重組，並獲得了葡萄球菌核酸酶中Glu66的令人鼓舞的pKa位移。在現實中，我們預計電子極化和結構重組都應該在決定埋藏殘基的pKa移位方面發揮一些作用，而且它們很可能相互交織在一起。我們相信，特定效應的實際貢獻水平將取決於具體情況。顯然，需要進一步的計算研究來澄清這種情況。
由於PPC代表了蛋白質在給定（原生）結構中的極化靜電狀態，我們預計在使用PPC的MD模擬中，與AMBER中的平均場狀電荷相比，蛋白質的結構細節會有一些差異。我們使用PPC和AMBER電荷對硫氧還蛋白進行了進一步的MD計算，以比較結構細節。首先，我們檢查了分別使用PPC和AMBER03進行的兩組MD模擬的RMSD的特徵。如圖6所示，它們明顯表現出不同的特徵。由於兩個計算中使用的力場中的所有其他引數都是相同的，除了電荷，圖6的結果清楚地顯示了蛋白質極化對詳細的蛋白質動力學的重要影響，這在本質上是靜電的。由於蛋白質內部氫鍵的形成使供體和受體極化，這種極化效應用PPCs比用AMBER電荷更好地表示。因此，在使用PPC的MD計算過程中，蛋白質中的氫鍵通常應該更加穩定，這可能解釋了圖6所示的硫氧還蛋白在MD模擬中的RMSD的不同特徵
總結
我們提出了一種新的計算方案，通過結合蛋白質電子結構的線性縮放片段方案和溶劑的連續介質模型進行自洽處理，對溶液中的蛋白質進行量子力學計算。計算出的蛋白質原生結構中的電子密度被用來產生極化蛋白質力場（電荷）。派生的蛋白質特定的PPC正確地代表了蛋白質在其原始結構附近的極化靜電狀態，可以很容易地用於蛋白質運動的MD模擬。我們對硫氧還蛋白的MD模擬的計算研究揭示了PPC的以下特性。
1.
PPC是特定於蛋白質的，因此正確地描述了蛋白質在接近原始結構附近的極化靜電狀態。同一型別的殘基上的電荷通常是不同的，這取決於它們特定的區域性靜電環境。這與基於氨基酸的電荷相反，後者是均值場的，無論在蛋白質中的具體位置如何，都保持不變。
2.
通過採用PPCs，使用熱力學整合方法的MD自由能模擬準確地再現了硫氧還蛋白中埋藏的Asp26的實驗性pKa偏移，而採用CHARMM和AMBER的標準平均場力場的相同計算則高估了pKa偏移，大約是兩倍。
指出PPC的侷限性也很重要。由於PPC是基於一個蛋白質的給定結構（大多數情況下是原生結構），用它來描述遠離給定結構的結構可能不合適，也沒有好處。幸運的是，大多數蛋白質的MD研究都是為了模擬接近原始結構的蛋白質運動，PPC在這種蛋白質的計算研究中應該是非常有吸引力的。