蛋白質pKa預測模型研究進展*

羅方芳 蔡志濤 黃艷東

(集美大學計算機工程學院,廈門 361021)

1 引言

為保證正常的生命活動,人體細胞的細胞質、細胞核以及各個細胞器需維持在特定的pH水平.例如,線粒體和溶酶體的pH分別是8.0和4.7,偏離細胞質的7.2[1].其中,用于表征溶液的酸堿度的pH為氫離子濃度的對數取負(pH=-log[H+]),其是人體中許多重要生物過程的調控因子,例如物質跨膜轉運[2]、酶催化[3]、蛋白質折疊[4]、多肽聚集[5]、脂質分子自組裝[6]、病毒入侵細胞[7]和細胞能量代謝[8].從微觀的角度,以上生物過程均與關鍵可離子化基團的質子化(protonation)或去質子化反應(deprotonation)相關聯.可離子化基團的去質子化(正反應)和質子化反應(逆反應): AH?A-+H+,其中,AH 是一種可離子化基團的質子化態,A-是去質子化態.

以β分泌酶BACE1為例闡述蛋白質功能和可離子化基團質子化/去質子化的關系.BACE1的生物功能是裂解β淀粉樣前體蛋白APP.它與神經退行性疾病阿爾茨海默癥密切相關,是典型的結構和功能依賴于pH的蛋白質.該蛋白的催化中心含兩個天冬氨酸Asp32和Asp228 (圖1(a)).實驗指出,BACE1僅在一個狹小的pH范圍內具有活性[9].如圖1(b)所示,在最適pH條件下(約等于4.5),Asp32處于質子化態,扮演質子供體(proton donor);Asp228處于去質子化態,扮演親核試劑(nucleophile).然而,當溶液pH偏離4.5,兩個天冬氨酸同時質子化或去質子化,BACE1無法行使其生物功能[10].

圖1 BACE1催化中心質子化態和功能的關系 (a) BACE1三維結構及其催化中心酸性二分體D32和D228;(b) D32和D228質子化態和蛋白質活性隨pH的變化規律(D是Asp的縮寫)Fig.1.Relationship between protonation state of BACE1 catalytic center and the function: (a) Crystal structure of BACE1 and the acidic dyad in the catalytic center;(b) protonation states of D32 and D228 and the activity as a function of pH (D is the abbreviation of Asp).

當一個可離子化基團的質子化和去質子化達到平衡,可由以下公式計算解離常數Ka:

其中,[H+],[A-]和[AH]分別代表溶液中氫離子以及該基團去質子化和質子化態下的濃度.Ka代表一種酸(如AH)離解氫離子的能力.將方程(1)的兩邊對數取負,可得到著名的Henderson-Hasselbalch方程:

其中,pKa為解離常數Ka的對數取負,代表一種酸(如AH)去質子化的難易程度.例如,溶液中天冬氨酸的 pKa測量值是3.7[11].根據(2)式,天冬氨酸在中性(pH=7.0)水溶液中處于去質子化態(A-);在pH小于3.7的酸性溶液中,天冬氨酸質子化(AH);當pH位于 pKa附近,質子化和去質子化態共存.如上所述,pKa決定了可離子化基團在任意 pH 條件下的質子化和去質子化反應平衡.根據 pKa值,可以推斷不同 pH 條件下生物大分子質子化態的分布,進而討論結構和功能的關系.因此,pKa是研究pH相關的生物化學過程的一個核心問題.不僅如此,pKa與藥物研發中長期存在的靶向性和抗藥性問題以及蛋白質設計密切相關.然而,由于蛋白質結構的復雜性以及實驗條件的限制,人們難于通過實驗獲取蛋白質中可離子化氨基酸殘基的 pKa,需借助理論預測.

為此,將以上Henderson-Hasselbalch方程轉換為能量形式,得到游離氨基酸關于pH和 pKa的去質子化自由能 ΔGmod的表達式:

其中,kB和T分別是玻爾茲曼常數和溫度;為游離氨基酸的 pKa,是可測量值.去質子化自由能可分解為成鍵作用部分 ΔGBond和非鍵作用部分ΔGNBond.其中,成鍵作用部分描述共價鍵斷裂的能量變化,計算復雜度高,不適用于生物大分子體系[12].值得一提的是,當溶劑中的可離子化氨基酸參與蛋白質的合成,蛋白質環境對成鍵作用部分的影響可忽略不計.基于該假設,我們只需考慮非鍵作用部分.因此,可離子化氨基酸從溶劑到蛋白質的去質子化自由能改變量 ΔG-ΔGmod可表示為

根據(3)式,ΔGmod為已知量.因此,求解蛋白質中氨基酸殘基的去質子化自由能 ΔG的問題簡化為計算蛋白質環境對非鍵作用部分的自由能微擾.

基于以上框架,人們發展了基于自由能計算的蛋白質 pKa預測模型,例如恒定pH分子動力學(constant pH molecular dynamics,CpHMD)[13].許多生物大分子含有不止一個功能構象,并且構象的轉變與質子化/去質子化反應相關聯: 當活性位點質子化(pH ＜ pKa),蛋白處于構象C1;去質子化(pH ＞ pKa),構象由C1轉變到C2;當pH取pKa附近,質子化和去質子化態共存,構象C1與C2相互轉變.因此,只有考慮了構象與質子化態耦合的理論模型,才能得到和實驗相一致的宏觀 pKa(macroscopic pKa)[14].CpHMD通過分子動力學模擬實現在不同構象下對質子化態空間進行采樣.在蛋白質 pKa預測精度方面,CpHMD相對其他現有模型具有明顯的優勢[15].CpHMD的缺點是 pKa計算效率低.例如,完成一個蛋白質 pKa的計算通常需要進行幾個小時甚至幾天的分子動力學模擬,因此難以滿足工業界大批量計算的需求.目前,CpHMD多被應用于結構和功能依賴于pH的藥物靶向蛋白的分子機制研究[16].

為了實現高通量的 pKa計算,人們發展了基于泊松-玻爾茲曼(Poisson-Boltzmann,PB)方程的模型,主要包括MCCE[17],H++[18],APBS[19],DelPhi-PKa[20]和PypKa[21].基于PB的模型能夠在幾分鐘內完成一個蛋白質的 pKa計算,極大地提高了計算效率.然而,基于PB的模型具有其理論局限性.例如,由于連續介質假設,PB方程不適用于非水溶性的膜蛋白.其次,蛋白質結構的復雜性增加了介電常數的不確定性,因此即便是水溶性蛋白,分子內部(例如酶的催化反應中心)的 pKa計算對介電常數敏感[22].

除了以上基于能量的模型,人們也可以用一個經驗函數描述某可離子化氨基酸殘基的蛋白質環境(如疏水環境和氫鍵)與其 pKa偏移量的映射關系.蛋白質某氨基酸殘基 pKa可表示為其游離狀態下參考值和偏移量 ΔpKa的和:

2005年,基于前期的第一性原理計算工作[12],哥本哈根大學Jensen課題組[23]提出了一個計算蛋白質 pKa的經驗函數PropKa.該模型提出一組經驗公式分別計算庫侖力、去溶劑化效應和氫鍵等關鍵因素對 pKa偏離參考值的貢獻.PropKa可在幾秒內完成一個蛋白質的 pKa計算,計算效率明顯比基于PB的模型高,近20年得到了廣泛的應用,其最新版本PropKa 3.0發表于2011年[24].

直到2021年12月,本課題組[25]發表了首個人工智能(artificial intelligence,AI)驅動的蛋白質 pKa預測模型DeepKa.隨后,美國卡內基·梅隆大學Olexandr Lsayev、美國約翰斯·霍普金斯大學Ana Damjanovic和德國拜耳公司Pedro Reis研究小組陸續提出了基于機器學習的 pKa預測模型pKa-ANI[26],XGB-WMa[27]和PKAI/PKAI+[28].其中,DeepKa和PKAI/PKAI+主要依賴于數據集,而為了在少樣本情況下建立有效模型,pKa-ANI和XGB-WMa需要一定程度的預訓練或先驗知識.值得一提的是,機器學習模型也能夠在幾秒內完成一個蛋白質的 pKa計算.

上述的CpHMD以及基于PB方程、經驗函數和機器學習的模型是目前4種主流的 pKa預測方法.最近,本課題組[29]采用CpHMD擴增了pKa數據集,進一步提高了DeepKa的預測精度.值得一提的是,DeepKa已展現出類似物理模型(如CpHMD)的高魯棒性,進一步證明了人工智能算法在蛋白質 pKa預測領域的有效性.下面將介紹這4種主流方法的理論基礎及研究進展.

2 蛋白質pKa預測方法

2.1 CpHMD

根據質子化態采樣方法的不同,恒定pH分子動力學CpHMD分為隨機采樣(discrete CpHMD,D-CpHMD)[30]和λ動力學(continuous CpHMD,C-CpHMD)[31].隨機采樣利用蒙特卡羅(Monte Carlo,MC)模擬在離散的質子化態空間(反應坐標取0或1)進行采樣[30].λ動力學則采用取值范圍0(質子化態)到1(去質子化態)的連續變量λ作為反應坐標對可離子化基團的電荷或體系哈密頓量進行標度[31].如圖2所示,先使用以上基于MC或λ動力學的采樣算法更新質子化態或者電荷.基于更新后的電荷分布,通過分子動力學模擬對構象進行采樣.更新位置坐標后,進入下一輪質子化態的采樣.模擬結束后,采用廣義Henderson-Hasselbalch方程擬合CpHMD模擬產生的不同pH條件下某可離子化基團的去質子化概率S,進而獲得其 pKa值,即S=0.5所對應的pH[31].

圖2 CpHMD模擬框架Fig.2.Framework of a CpHMD simulation.

由于滴定動力學與構象動力學相關聯,提高質子化態和構象空間的采樣是近30年CpHMD模型發展的主線.下面將分別介紹D-CpHMD和CCpHMD.

2.1.1 D-CpHMD

D-CpHMD用一個反應坐標λ表示某可離子化位點的質子化態.λ只能取0或1.其中,0和1分別表示質子化態和去質子化態.經過一定長度的分子動力學(molecular dynamics,MD)模擬,隨機選取一個可離子化基團,嘗試改變其質子化態.例如,將其λ值從0改為1.然后,計算λ值改變引起的能量變化 ΔE.將該能量變化代入Metropolis準則:

如果能量差小于或等于0,接受λ值改變的概率為1.如果能量差大于0,則接受改變的概率p小于1.在數值模擬中,通常是隨機生成一個取值范圍為[0,1]的數s.只有s小于等于p,才接受λ值改變,否則保留原值.以上為一步的MC,和開始的MD構成一個模擬周期.因此,在MC之后,便是下一個周期的MD模擬.顯性溶劑下質子化或去質子化的能量變化較大,導致較小的接受概率.起初,為了提高接受概率或質子化態的采樣效率,MC的能量計算使用隱性溶劑(implicit solvent)模型,如廣義玻恩(generalized Born,GB)[32-34]和引言提到的PB模型[31,35,36].當MC和MD均采用隱性溶劑,計算效率最高,但是犧牲了精度[32,33].為了提高構象方面的采樣精度,MD可替換成顯性溶劑,即雜化溶劑[31,34,35].其中,基于GB和PB的模型分別在分子模擬軟件Amber和GROMACS中已被實現.需要指出的是,隱性溶劑難以描述活性位點附近與功能相關的水分子或鹽離子對去質子化平衡的影響[37].

為提高顯性水溶劑下MC的接受概率,2007年Stern[38]提出了嘗試改變λ值之后,先進行一定長度的嘗試性的分子動力學模擬,再計算能量差.該嘗試性的MD使周圍水溶劑構型得到調整,可降低λ值改變前后的能量差.然而,以上嘗試性MD的長度依賴于經驗或不確定,其應用可能受到限制.盡管如此,該模型為解決顯性溶劑下質子化態空間的采樣問題提供了一條新思路.隨著高性能計算的發展,人們開始考慮將顯性溶劑應用到蛋白質D-CpHMD的MC部分.如無特別說明,以下提到的顯性溶劑均是分子動力學模擬中計算靜電相互作用的標準算法PME (particle mesh Ewald,PME)[39].2015年芝加哥大學的 Roux課題組[40]提出了顯性溶劑下的非平衡MD/MC模擬.例如,對于某可滴定位點在MC階段的去質子化(λ由0變為1)嘗試,該模型在0和1之間添加了m個中間值.對于每個λ值(m個中間值和兩個邊界值0和1),執行一定長度的非平衡MD,令可離子化基團周圍的環境根據λ值在構型上作出調整,減緩了因λ值改變而導致的能量漲落.結束λ=1的非平衡MD后,計算當λ=1和λ=0的能量差.同樣,根據Metropolis準則,如果接受該可滴定位點去質子化,繼續λ=1的MD.否則,退回到非平衡MD前的時刻,繼續λ=0的MD.通過以上的非平衡模擬,該模型提供了較合理的能量差的計算,提高了總體接受概率.Roux課題組[40,41]利用著名的Jarzynski方程將自由能變化與非平衡MD所做的功相關聯,使得以上非平衡MD的模擬時間可被量化.值得一提的是,該方法可被應用于生物大分子,目前在分子模擬軟件NAMD中已有實現.然而,可滴定氨基酸的固有 pKa(inherent pKa)是該模型的一個主要參量.為了提高預測性能,該模型要求固有 pKa盡可能接近真實值[41].因此,D-CpHMD一個潛在的研究方向是消除上述模型對固有 pKa的依賴.

2.1.2 C-CpHMD

本課題組統計了4057個蛋白質中可滴定氨基酸的個數[29].這些蛋白質來自復旦大學王任小實驗室[42]創建的蛋白質抑制劑復合物數據庫PDBbind的精細集v2016.除了半胱氨酸Cys,蛋白質中其他可滴定氨基酸類型(谷氨酸Glu、天冬氨酸Asp、賴氨酸Lys、精氨酸Arg、絡氨酸Tyr、組氨酸His)的平均個數不低于10[29].理論上,一個含有N個可滴定氨基酸殘基的蛋白質包含2N個質子化態.然而,D-CpHMD的MC每次只取一個可滴定位點來判斷是否改變其質子化態,采樣效率較低[34,43,44].

2004年,為了研究生物大分子體系(如蛋白質,DNA和RNA)的質子化和去質子化,密西根大學Brooks課題組開發了首個λ動力學框架下[45]的恒定pH分子動力學C-CpHMD[31].每個可滴定位點對應一個反應坐標λ,取值范圍同樣是0—1.和D-CpHMD不同的是,C-CpHMD的反應坐標是連續的變量.值得一提的是,C-CpHMD同時更新所有可滴定位點的質子化態.哈密頓量H代表體系的總能量,包括動能和勢能.除了真實的粒子,如模擬體系中溶劑和溶質的原子,C-CpHMD添加了虛粒子.每個可滴定基團對應一個虛粒子.這里用范圍在[0,1]的連續變量λ作為虛粒子的坐標.為了模擬虛粒子的滴定動力學,可將其質量設為10(單位是原子質量).以下是修正后的總哈密頓量:

其中,Natom是總粒子數,r是原子的位置矢量,λ是虛粒子的滴定坐標,ma和mj是原子和虛粒子的質量.第1和第4項的求和分別是原子和虛粒子的總動能.第2項Ubond是鍵相互作用能,包括鍵伸縮能、鍵角彎折能和二面角扭轉能.這里假設鍵相互作用與λ無關.第3項Unbond是非鍵相互作用能,包括靜電Uelec和范德瓦耳斯UvdW相互作用,與λ相關.最后一項U*是偏置勢,利用經驗勢描述去質子化鍵斷裂的能量變化,只和λ相關.

以下介紹如何利用λ標度非鍵相互作用能和偏置勢.對于可滴定的氫原子和周圍原子的范德瓦耳斯相互作用,直接用 1-λi標度勢能函數(這里采用6-12勒讓德瓊斯勢ULJ):

可見,當λ=1時,殘基i去質子化,殘基i的可滴定氫與j無相互作用.

對于兩個可滴定氫之間的范德瓦耳斯相互作用,采用 1-λi和1-λj進行標度:

范德瓦耳斯力是近程非鍵相互作用力,主導疏水基團間的相互作用.然而,由于原子半徑的差異,氫(半徑約1 ?)幾乎被與之成鍵(鍵長約1 ?)的重原子(半徑約2 ?)包圍,使其難以接觸到其他原子.因此,質子化和去質子化對范德瓦耳斯相互作用影響不大,相對長程靜電相互作用可以忽略不計.對于靜電相互作用,λ標度的是原子電荷[31]:

早期為了提高計算效率,Brooks課題組[31]采用隱性溶劑模型計算溶劑對溶質的平均效應.如此一來,總靜電能Uelec的溶質內靜電相互作用仍采用庫侖勢((11)式第1項),而溶質與溶劑的靜電相互作用Usolv采用GB勢能函數((12)式):

其中,星號代表排除存在鍵相互作用的原子對;rab是電荷qa和qb的距離;εp和εw是蛋白質和水的介電常數;κ是德拜長度取反((13)式);I是鹽離子強度;q是鹽離子電荷;e是基本電荷;kB是玻爾茲曼常數;T是溫度;α是有效玻恩半徑,表征某原子埋在蛋白內部的程度,為衡量GB模型精度的關鍵參數.相對PB模型,GB的計算復雜度較低,并且是解析的,適合需要對位置坐標求一階導(計算粒子所受合外力)的分子動力學模擬.GB模型的計算復雜度主要體現在有效玻恩半徑的求解.

2004和2005年Brooks課題組接連開發了CH ARMM軟件中基于隱性溶劑GBMV[31]和GBSW[46]的C-CpHMD,證明了基于GB的C-CpHMD在pKa預測方面的有效性.相對GBSW/GBMV溶劑模型,GBNeck2可提供更優的構象采樣[47].于是,馬里蘭大學Shen課題組[48]在2018年開發了Amber軟件中基于隱性溶劑GBNeck2的C-CpHMD.值得一提的是,對于實驗科學家關心的酶催化中心(如圖1活性位點Asp32和Asp228),該方法也表現較好,目前已被應用于共價抑制劑靶點的預測[49-51],蛋白質 pKa數據集的建立[25,29],以及依賴于pH的蛋白質分子機制研究[52,53].目前,基于GBSW和GBNeck2的C-CpHMD均已實現GPU加速,這進一步擴展了模型的應用范疇[54,55].

為了提高構象采樣精度以及擴展C-CpHMD的應用范圍,Shen課題組[56]提出了雜化溶劑CCpHMD: 構象動力學使用顯性溶劑;而滴定動力學保留隱性溶劑.為此,構象動力學和滴定動力學采用不同的哈密頓量.前者去掉方程(7)的最后兩項,第3項不再包含反應坐標λ,令方程(7)回歸到常規分子動力學.該方法不僅維持了質子化態空間采樣效率,而且提高了構象采樣精度.起初人們會擔心隱性溶劑GB的理論局限性(例如偏弱的疏水效應)會影響 pKa預測精度.然而,Shen課題組[56]發現,顯性溶劑PME可導致偏高的疏水效應,一定程度上抵消了隱性溶劑導致的偏弱的疏水效應.相對隱性溶劑,該雜化溶劑C-CpHMD獲得了廣泛的應用,如鈉離子質子交換蛋白[37,57],質子通道[58],類藥物分子的膜滲透[59],芬太尼激活G耦聯受體[60],糖苷水解酶[61],絡氨酸激酶藥物發現[62],以及上文提及的β分泌酶[10].

為了描述和功能相關的水分子或其他輔助因子(如金屬離子和小分子)對去質子化平衡的影響,滴定動力學部分也需采用顯性溶劑.起初,Brooks課題組和Shen課題組分別選擇了較簡單的基于截斷的顯性溶劑FSh (force shifting,FSh)[63]和GRF(generalized reaction field,GRF)[64].然而,由于截斷,這兩個模型均低估了長程靜電力對可滴定位點的影響[65].為此,Shen課題組[66]開發了基于顯性溶劑PME的C-CpHMD.最近,該模型在分子模擬軟件Amber中實現了GPU加速[67].眾所周知,PME是滿足周期性邊界條件(periodic boundary condition,PBC)的分子模擬中計算靜電相互作用的標準算法,因此基于PME的C-CpHMD是λ動力學框架下所能達到的最優版本.理論上,如果不考慮取樣問題,該模型的 pKa預測應該最接近實驗.對于一個滿足PBC的分子動力學模擬體系,PME的總靜電能是3個能量項的加和:

其中,Udir是實空間靜電相互作用,在庫侖勢基礎上增加一個補償函數,負責截斷距離以內的短程靜電相互作用((15)式).Urec最為耗時,為倒格空間(reciprocal space)下求解的長程靜電能,負責截斷以外的長程靜電相互作用((16)式).Ucorr是修正項((20)式)[39].

其中,ra和rb是中心元胞的位置矢量;n是元胞的位置矢量,其表達式為n=n1c1+n2c2+n3c3,其中c1,c2和c3代表元胞的3個正交方向矢量;星號代表被排除的原子對,包括原子自身(a=b),形成化學鍵的原子對,以及最近鄰(n的大小為1)以外的鏡像;erf是補償誤差函數;參數β決定Udir和Urec的相對收斂速度.例如,β越大,Udir計算收斂越快,而Urec計算收斂會越慢.

式中m是倒格矢,其表達式為,其中,m1,m2,m3是非零整數;是以上ci(i=1,2,3)的共軛倒格矢,二者滿足關系式=δij,這里i和j取1,2和3.另外,V=c1·c2×c3,是元胞的體積.S(m) 是結構因子:

該結構因子可近似表示為

式中通過將元胞中的電荷分布(B樣條)插值到具有相同的3個維度k1,k2,k3的網格來構造三維矩陣Q;ki/Ki是分數坐標,其中,ki(i=1,2,3)取值范圍是(1,2,3,···,Ki),正整數常數Ki代表元胞的尺寸;F(Q) 是矩陣Q的三維快速傅里葉變換.經過以上變換,Urec的表達式為

值得一提的是,Urec線性依賴于格點電荷,因此對λ求一階導和庫侖勢的一樣簡單.

Urec考慮整體的電荷分布,并未排除存在鍵相互作用的原子對,因此需采用和Udir相同的函數形式進行修正((20)式第1項).此外,Ucorr第2項的作用是排除點電荷自相互作用,第3項則是中和體系凈電荷的背景電荷(background plasma).其中,后面兩個修正只依賴于原子電荷.

為了避免元胞之間不真實的靜電相互作用,常規MD通過添加補償鹽離子使體系呈電中性.然而,CpHMD模擬中電荷是動態變化的.為了解決該問題,Shen課題組[64]提出了將鹽離子作為質子緩存器.然而,鹽離子如果不帶電會導致聚集,于是改使用可滴定水分子[68].酸性氨基酸(例如Asp和Glu)與水陰離子(hydroxide,TIPU)耦合(AH+OH-?A-+H2O);堿性氨基酸(例如Lys,Arg和His)與水陽離子(hydronium,TIPP)耦合(BH++H2O?H3O++B).該耦合令反應式兩端的電荷守恒.電中性的另一個好處是消除Ucorr中會導致反常 pKa偏移的背景電荷.

以上介紹了不同溶劑下靜電能的具體求解.下面介紹哈密頓量中只依賴于反應坐標λ的偏置勢[31]:

其中,第1項((22)式)和第2項((23)式)分別是游離可滴定氨基酸去質子化的非鍵相互作用能和總自由能.對于單個可滴定位點的氨基酸(如賴氨酸),Umod是一個關于λ的 一元二次函數.UpH由λ線性標度((23)式).UpH-Umod是化學能改變量的近似解.為了減少λ處于不真實的中間態(如λ=0.5)的概率,另外添加了一個二次函數勢壘Ubarr((24)式).Ubarr降低了λ的動力學,對熱力學統計沒有影響.(23)式和(24)式的參數為已知,因此,C-CpHMD的主要工作是確定Umod的參數(如(22)式中的Aj和Bj):

需要注意的是,為了將λ約束在[0,1],需定義另一個變量θ.λ和θ的關系式為λ=sin2θ.于是,數值模擬中進行迭代的是θ,而非反應坐標λ.

對于含有兩個可滴定位點的氨基酸,需要定義反應坐標x來描述處于去質子化(His)或質子化(Glu和Asp)態時質子所處的可滴定位點[46].x同樣是在0到1范圍內的連續變量.圖3展示了Asp和His側鏈3個質子化態對應的反應坐標值以及狀態間的轉化.類似變量λ,可利用插值將x加入哈密頓量的各個能量項.例如,以下分別是Asp和His電荷關于λ和x的表達式:

圖3 互變異構滴定模型的3個質子化態以及狀態間的轉化 (a)天冬氨酸Asp;(b)組氨酸HisFig.3.Three protonation states and their interconversion in the tautomeric titration model: (a) Aspartic acid;(b) histidine.

CpHMD模擬同時對構象和質子化態采樣.根據設置的輸出頻率保存每個可離子化基團的滴定坐標λ(λ∈[0,1])(圖4(a)).統計處于質子化態(0≤λ≤0.1)的次數Nprot以及去質子化態(0.9≤λ≤1)的次數Ndep,計算不同pH條件下的去質子化概率S(圖4(a))[31]:

圖4 基于C-CpHMD的 pKa 計算 (a)滴定坐標λ和去質子化概率S的軌跡;(b)采用Hill函數擬合SFig.4.The pKa calculation based on C-CpHMD: (a) Trajectories of titration coordinate λ and deprotonation fraction S;(b) fitting S to Hill function.

最后,采用如下Hill函數(廣義Henderson-Hasselbalch函數)擬合S.pKa便是S=0.5時所對應的pH (圖4(b)):

其中h是Hill系數,表征一個可離子化基團與周圍可滴定基團的滴定動力學是否存在耦合.h=1表示無耦合,如位于分子表面的殘基或游離氨基酸.h＜1表示負耦合,如形成鹽橋鍵的去質子化的Asp和質子化的Lys.h＞1 表示正耦合,如酶活性位點距離相近的兩個酸性氨基酸(質子化的Asp或Glu).h偏離1越多,耦合越強[69].

當兩個氨基酸的滴定動力學存在耦合,可將二者看作一個整體,利用以下公式計算宏觀 pK1和pK2(macroscopic sequential pKa)[64,70]:

（2）邊坡系數為1∶1的表層無覆土無植草邊坡60min內的土壤流失質量1.31kg，15～20，20～25，25～30mm植生混凝土組的底部土壤流失量分別為0.0792，0.1089，0.1584kg；植生混凝土對底部土壤的反濾攔截率達86.92%。

其中N是一定pH條件下的平均質子數.為獲得pK1和 pK2,也可以采用以下非耦合模型(31)式[71,72]:

其中S1和S2分別是兩個耦合的可滴定位點的去質子化概率.

當滴定動力學采用滿足周期性邊界條件的顯性溶劑時,需要考慮有限尺度效應[73].由于采用耦合水離子實現了電中性,有限尺度效應僅剩下和水分子模型相關的離散溶劑效應(discrete solvent effect)[66].當某個可滴定氨基酸去質子化,因離散溶劑效應引起的能量變化是

其中,κ是介電常數;ρ是水數量密度,等于水分子數N除以體積V,這里N指的是和蛋白有相互作用的水分子數,V也是這些水包絡范圍內的體積;q是可滴定氨基酸的電荷,Asp/Glu是-1e,His/Lys為+1e;γ是顯性溶劑模型范德瓦耳斯相互作用中心的電四極矩.對于溶劑模型TIP3P,γ的值為 0.764e·?2.為了估算該有限尺度效應導致的pKa偏移,需要計算相對模型分子的能量變化[66]:

其中,N和Nmod分別是蛋白質和游離氨基酸模擬體系中與溶質有相互作用的水分子數;V和Vmod是相應的周期性元胞體積.將以上表達式轉化為pKa偏移量,可得到[66]

根據N和V的定義,可以推斷有限尺寸效應對PME影響較大.PME考慮了周期性元胞內所有水分子,蛋白質體積所占比例較小,水數量密度ρ較大;另一方面,GRF和FSh僅考慮截斷以內的水,蛋白質體積所占比例較大,水數量密度可忽略不計.對于膜蛋白體系,可參考Roux課題組[74]提出的方法做相應的修正.

以上介紹的C-CpHMD屬于對電荷插值,實現電荷對反應坐標的線性響應.實際上,由于庫侖勢對電荷線性依賴,庫侖勢和電荷兩者的線性插值是等效的.因為兩種情況下,關于插值變量(反應坐標λ)負的一階導數(作用在虛粒子上的合外力)是相等的.然而,并不是所有和靜電勢相關的能量項和電荷線性相關,如PME算法中對點電荷自相互作用和凈電荷的修正項((20)式)[66].所以,為了更好描述電荷變化對滴定動力學的影響,基于截斷的GRF和FSh較適合對靜電勢進行插值的CCpHMD,因為它們的靜電勢保留了對電荷的線性依賴.德國馬克思普朗克研究所的Grubmüller課題組[75]在分子模擬軟件GROMACS中開發的C-CpHMD便是對勢函數進行插值.最近,芬蘭的Groenhof課題組[76,77]基于該模型進行代碼優化,并實現基于CHARMM力場的CpHMD模擬.然而,該模型采用了顯性溶劑PME,而不是基于截斷的GRF或FSh.其次,該模型沒有像Shen課題組[64]一樣考慮有限尺寸效應.另一方面,同樣是對勢能進行插值,Brooks課題組[63,71]基于顯性溶劑的C-CpHMD模型合理地采用了基于截斷的FSh.除了以上正弦函數形式,Grubmüller課題組和Brooks課題組提出了其他將λ約束在區間[0,1]的方法.例如,Grubmüller課題組[75]提出了余弦形式.Brooks課題組[78]提出一個較復雜的指數形式.對于顯性溶劑C-CpHMD,體系電中性是一項重要的約束條件,Shen課題組[66]和Grubmüller課題組[79]均采用了可滴定水分子實現體系凈電荷恒等于0.然而,Brooks課題組[71]的顯性溶劑C-CpHMD還未考慮該約束.因此,為了避免溶質與其鏡像的靜電相互作用,需對FSh靜電勢設置較小截斷值.

隨著顯性溶劑CpHMD的快速發展,急需解決質子化態和構象的采樣問題.2006年Brooks課題組[82]率先將基于溫度的副本交換(replica exchange)算法應用到C-CpHMD,即將副本以一定的概率交換到較高溫度,借助熱漲落提高CpHMD模擬的采樣.受到哈密頓量副本交換算法的啟發,2011年Shen課題組[56]提出了基于pH的副本交換算法: 將副本以一定的概率p交換到較高的pH,提高去質子化態的采樣;或交換到較低的pH,提高質子化態的采樣((35)式).因為實際進行交換的pH只存在于UpH((23)式),交換前后總能量的變化Δ/β可簡化為僅含UpH的表達式((36)式).交換pH后,兩個副本將在新的pH條件下(或新的UpH)進行采樣.該算法效率極高,同時操作簡單,已被應用到其他CpHMD模型[83-86].為了增強質子化態空間采樣,美國國立衛生研究院NIH的Brooks課題組[87]提出結合包絡分布采樣(enveloping distribution sampling,EDS)和哈密頓量副本交換(Hamiltonian replica exchange,HREX).EDS通過定義一個參數s標度狀態間的能壘.較小的s對應較平滑的能壘,方便了狀態間的轉化.然而,能壘的消除促進了虛擬中間態的采樣,這將影響物理態的采樣.為了避免中間態的采樣,在EDS基礎上利用HREX提高離散的質子化態空間的采樣效率.接著,該課題組[86]加入以上基于pH的副本交換,構成二維的副本交換.從算法的角度,該方法確實提高了采樣效率,但代價是產生大量的副本以及模擬過程中副本的頻繁通訊,對計算能力要求較高.近期,為了在有限GPU顯卡數量的條件下實現基于pH的副本交換,Shen課題組[88]提出了副本同步交換.

其中,p是副本交換的概率;UpH({λj};pH) 和是兩個副本交換前的UpH.將以上兩項的 pH和pH′進行互換,得到UpH({λj};pH′) 和.

除了副本交換,另一種增強采樣的方法是對生物大分子進行粗粒化(coarse graining,CG),減少模擬體系中粒子的數量,從而降低了構象空間的自由度.該方法通常被應用于具有較大空間和時間尺度的生物過程,如蛋白質折疊、多肽聚集和物質跨膜轉運等[89].近幾年,研究者們開始將CG與CpHMD結合,發展CpHMD的粗粒化模型[90-93].值得一提的是,提出Martini粗粒化力場的Marrink課題組[92]已在分子模擬軟件GROMACS中實現了CpHMD的粗粒化模擬.

2.2 基于PB的 pKa 預測模型

實際上,如果只考慮單個結構,可以用PB方程計算相對去質子化自由能 ΔΔG=ΔG-ΔGmod.其中,ΔGmod是某可離子化氨基酸A在游離狀態下去質子化自由能改變量:

式中Gmod(A-)和Gmod(AH) 分別是去質子化(A-)和質子化(AH)狀態的自由能.同理,當該氨基酸參與蛋白質的合成,它在蛋白質中的去質子化自由能改變量 ΔG表示為

基于蛋白質環境不影響成鍵作用部分ΔGBond(見(4)式)的假設,以上兩個自由能改變量的差可表示為

其中,下標PB表示用PB方程分別計算等式右邊4個狀態下的靜電能.令ΔG(AH)=GPB(AH)-,可得到

其中,ΔGPB(AH)和ΔGPB(A-) 分別表示在水溶劑中將質子化(AH)和去質子化(A-)的氨基酸放入蛋白質的靜電能改變量.基于該等式,可以得到如圖5所示的熱力學循環(thermodynamic cycle).相對去質子化自由能 ΔΔG 可表示為

圖5 相對去質子化自由能計算的熱力學循環Fig.5.Thermodynamic cycle of relative deprotonation free energy calculation.

接著,將 ΔΔG代入關系式ΔpKa=ΔΔG/(kBTln10)計算 pKa偏移量 ΔpKa.最后,利用(5)式計算 pKa.可見,熱力學循環4個狀態的靜電能計算決定了pKa的預測精度.目前,基于PB計算靜電能并預測蛋白質 pKa的方法包括MCCE[17,94],H++[18],APBS[19],DelPhiPKa[20,95,96]以及PypKa[21].其中,MCCE和PypKa利用MC對側鏈二面角進行采樣,一定程度上提高了預測精度,但總體精度仍低于CpHMD,說明了空間構象充分采樣的重要性[15].PB方程的參數主要是介電常數,原子的電荷和半徑,因此容易拓展到其他類型的體系.例如,除了蛋白質,DelPhiPKa也適用于DNA和RNA.除了蛋白質單體,H++也考慮了含有配體的復合物.

2.3 基于經驗函數的 pKa 預測模型

以上物理模型(CpHMD和基于PB的模型)需要計算體系的靜電能,計算復雜度較高.為了進一步提高 pKa計算的效率(例如將單個蛋白的pKa計算時長縮短到秒量級),2005年哥本哈根大學的Jensen課題組[23]提出了一組經驗函數PropKa分別描述點電荷相互作用(Coulomb force)、去溶劑化效應(desolvation)和氫鍵相互作用(hydrogen bonding)對 pKa偏移量的貢獻:

以上3項的函數均采用分段的一次函數,計算復雜度低,已被應用到蛋白質單體[23],蛋白質和小分子配體的復合物[97].然而,該版本的PropKa沒區分可滴定氨基酸殘基是處于蛋白質的表面還是內部.

為此,2011年Jensen課題組[24]提出了改進的PropKa 3.0.新版本考慮了相同的 ΔpKa決定因子,將(42)式的氫鍵相互作用導致的和去溶劑化效應導致的歸為自能不同的是,PropKa 3.0采取了一個折中的方案,即部分使用能量公式.例如,點電荷相互作用采用經典的庫侖勢.去溶劑化效應采用了和GB模型中求解有效波恩半徑的倒數(1/α)類似的原子體積(V)除以原子間距離的四次方(r4).此外,蛋白質表面和內部被賦予不同的介電常數.對于氫鍵相互作用,則保留了一次函數形式.該模型的參數化基于谷氨酸和天冬氨酸的 pKa實驗值,對酸性氨基酸的預測能力接近CpHMD[98].然而,該模型對堿性氨基酸(如Lys和His)的預測效果較差[25].

2.4 基于機器學習的 pKa 預測模型

上述PropKa經驗函數的提出較大程度依賴于科學家的先驗知識.理論上,如果有足夠多的pKa實驗測量值,可以結合數據和機器學習算法訓練出一個經驗函數,而不需要依靠已有的知識.2018年 波蘭華沙大學Siedlecki課題組[99]提出首個基于深度學習的蛋白質配體結合親和性(binding affinity)預測模型.這里的配體通常指具有幾何結構的小分子.我們知道,pKa表征某可滴定基團去質子的難易程度.換一種表達,pKa代表蛋白質和質子的結合親和性.可見,蛋白質配體結合親和性預測方法對 pKa預測具有參考價值[25].

由于實驗條件的限制,迄今為止蛋白質可滴定氨基酸殘基的 pKa實驗測量值不到兩千個[100,101].于是,本課題組采用基于隱性溶劑GBNeck2的CCpHMD[48]建立了一個蛋白質 pKa數據集(包含12809個 pKa)[25].2021年12月,本課題組提出了國際上首個基于機器學習的蛋白質 pKa預測模型DeepKa,證明了引入人工智能方法解決蛋白質pKa預測問題的可行性[25].本課題組對現有的pKa數據庫PKAD[100](包含1350個蛋白質 pKa實驗測量值)進行數據清洗,得到了測試集EXP67S.首先,根據氨基酸序列相似性比對排除了冗余數據.剩下的67個蛋白質的470個Asp,Glu,Lys或His的 pKa構成數據集EXP67.接著,對EXP67進行欠采樣,使得不同 ΔpKa區域分布均勻.最后剩下的167個 pKa為該模型的測試集EXP67S.該測試集的優勢將在下文的多模型比對體現出來(圖6).模型的大部分輸入特征以及三維卷積神經網絡(convolutional neural network,CNN)框架均借鑒Siedlecki課題組[99]提出的Pafnucy模型.值得一提的是,為了解決截斷導致的邊界問題,DeepKa采用格點電荷(Siedlecki課題組[99]采用原子電荷)描述對 pKa預測精度起決定性作用的靜電環境[25].雖然DeepKa第一版本的預測精度高于PropKa 3.0,但是和CpHMD還存在一定差距[25].此外,該工作只測試了DeepKa的總體性能,并未對特定的問題(如酶催化中心或無序蛋白)進行討論.

圖6 pKa預測模型性能對比Fig.6.Comparison of existing pKa predictors.

2022年1月,美國卡內基·梅隆大學Lsayev課題組[26]開發了基于神經網絡勢ANI-2X和原子環境矢量AVE的深度學習模型pKa-ANI.然而,該模型將所有的實驗數據用于模型的訓練,不利于對其性能進行客觀的評價.另外,該模型對結構敏感,需要在預處理階段對初始結構進行能量最小化,否則將得到不合理的預測結果[26].2022年3月,美國約翰斯·霍普金斯大學Damjanovic課題組[27]測試了4種基于樹的機器學習算法.其中,XGB-WMa表現最好.該小組同樣采用有限的實驗數據來訓練和測試模型.為了建立有效的模型,他們在特征描述上加入了較多的經驗知識: 首先,統計可滴定基團參與的氫鍵數量;其次,計算可滴定基團的溶劑可及表面積(solvent accessible surface area,SASA);最后,根據是否帶電或親水對可滴定基團附近氨基酸殘基進行分類.顯然,以上特征基本上覆蓋了PropKa模型中影響 pKa偏移量的3個關鍵因素:氫鍵相互作用、去溶劑化效應和點電荷相互作用.2022年7月,Reis課題組[102]利用基于PB的Pyp Ka建立了包含1200萬個 pKa值的數據集,并基于該數據集開發了深度學習模型PKAI[28].為了提高精度,在PKAI基礎上對損失函數進行正則化處理,從而得到PKAI+.然而,PKAI+在其他測試集(如EXP67S)的表現與PKAI相似,說明上述的正則化處理缺乏普適性[29].因此,如果沒有特別說明,下文只討論PKAI.

2023年5月,本課題組發布了DeepKa的最新版本[29].該版本的輸入特征和模型框架與舊版本相同,僅僅是增加了訓練和驗證集的 pKa樣本量.這些樣本出自549個蛋白質的26552個Asp,Glu,Lys和His.相對舊版本,該版本預測性能更接近CpHMD.此外,在這個工作中特定的蛋白質體系被用于進一步評估DeepKa的可靠性.例如,酶催化中心具有復雜的靜電環境,是 pKa預測的一個重要挑戰.新版本通過 pKa計算準確預測了5個酶催化中心的質子供體.除了具有穩定三維結構的蛋白,該模型也可被應用于無序蛋白.理論預測pKa偏移量較小的滴定位點往往容易做到預測精確,但難以做到預測相關,而即使在 pKa偏移量小于1.0的情況下,理論和實驗仍然表現出較高的相關性,證明了該模型的高魯棒性[29].如無特別說明,下文的DeepKa代表該新版本.

上述基于AI的模型均采用PKAD中的實驗數據來訓練或測試模型.然而,pKa-ANI,XGBWMa和PKAI忽略了存在于PKAD的冗余數據(例如一個蛋白質有兩組相同的 pKa值),這可能導致過擬合.其次,PKAD中大多數 pKa處于參考值附近,因此測試結果并不能反應模型真實的預測能力[25].值得一提的是,本課題組創建的測試集EXP67S不存在以上兩個問題,可較為客觀地對模型進行評價[25].研究發現,除了在實驗和理論相關性方面仍舊低于CpHMD,DeepKa的預測精度明顯高于其他主流 pKa預測模型,包括PypKa,PropKa,PKAI和pKa-ANI[29].其中,PypKa代表基于PB的模型,PropKa代表基于經驗函數的模型,PKAI和pKa-ANI代表其他AI模型.基于樹的XGB-WMa沒有開放源代碼,所以無法利用EXP67S對其進行測試.因此,XGB-WMa不參與下面的模型討論.同時考查精度和速度,圖6展示了5個模型的預測性能.其中,平均絕對誤差用于表征模型的精度.顯而易見,如果以PropKa的速度和CpHMD的精度作為參照,目前只有DeepKa能提供準確的高通量 pKa計算[29].最近,加拿大國家研究委員會Sulea課題組[103]比較了現有的7種高通量 pKa預測模型,包括基于經驗函數的PropKa 3.0[24],基于深度學習的DeepKa[29]、PKAI和PKAI+[28]以及基于PB方程的DelPhiPKa[95]、MCCE2[94]和H++[18].該研究指出在以上高通量模型中DeepKa的精度最高,與圖6的結論一致.

3 結論

pH與溫度、壓強一樣是基本的環境參量.傳統的分子動力學假設溶劑是中性水(pH=7.0),不考慮其他pH條件;此外,傳統分子動力學假設電荷是固定的,不受溶質靜電場的影響.以上兩個假設限制了傳統分子動力學進一步探究細胞中許多與pH相關的生物過程,而可靠的 pKa計算將有助于解決該難題.本綜述主要介紹了4類主流的pKa預測方法.顯然,對于不同理論的 pKa預測模型,其適用范圍也存在差異.首先,不論何種特定的問題,如果不要求高通量計算,可采用預測精度較高但計算效率較低的CpHMD.當涉及非水溶性蛋白(如膜蛋白)的 pKa計算,目前理論上可行的模型為基于雜化溶劑[37,56]或顯性溶劑[66,67]的CpHMD.另一方面,需要開發高通量的 pKa預測模型,從而滿足工業界批量的 pKa計算需求.由于隱性溶劑的理論局限性和實驗條件的限制,上述的高通量模型僅適用于水溶性蛋白.對于水溶性蛋白質單體的pKa計算,在所有高通量模型中DeepKa無疑是最優的選擇[29,103].若只關心酸性氨基酸殘基(如Asp和Glu)的質子化態,也可考慮PropKa 3.0[24].而對于主要的4種可離子化氨基酸殘基(Asp,Glu,Lys和His)以外的可滴定基團(如Cys和Tyr),可考慮基于PB的模型(如H++[18]和PypKa[21]).

隨著計算機軟件和硬件的快速發展,國際著名的美國藥物設計公司薛定諤(Schr?dinger)開始嘗試利用自由能微擾(free energy perturbation,FEP)方法計算 pKa,說明蛋白質 pKa理論計算開始引起工業界的關注[104].值得一提的是,基于機器學習的pKa預測模型雖處于起步的階段(2021年至今),卻已表現出和物理模型同水平的預測精度,例如本課題組開發的DeepKa.我們相信: AI模型有可能突破先驗知識,在不久的將來提供更為高效的預測;利用物理模型CpHMD建立的 pKa數據集PHMD549和基于 pKa數據庫PKAD建立的測試集EXP67S將為基于機器學習的 pKa預測工具的研發奠定基礎[29].最近,基于DeepKa本課題組開發了國內首個蛋白質 pKa在線計算平臺(http://www.comput biophys.com/DeepKa/main),這對未來參與到人工智能驅動的新藥研發產業具有重要意義[105,106].