羅方芳 蔡志濤 黃艷東

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

1 引言

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

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

圖1 BACE1催化中心質(zhì)子化態(tài)和功能的關(guān)系 (a) BACE1三維結(jié)構(gòu)及其催化中心酸性二分體D32和D228;(b) D32和D228質(zhì)子化態(tài)和蛋白質(zhì)活性隨pH的變化規(guī)律(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).

當(dāng)一個可離子化基團的質(zhì)子化和去質(zhì)子化達到平衡,可由以下公式計算解離常數(shù)Ka:

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

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

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

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

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

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

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

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

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

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

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

2 蛋白質(zhì)pKa預(yù)測方法

2.1 CpHMD

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

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

由于滴定動力學(xué)與構(gòu)象動力學(xué)相關(guān)聯(lián),提高質(zhì)子化態(tài)和構(gòu)象空間的采樣是近30年CpHMD模型發(fā)展的主線.下面將分別介紹D-CpHMD和CCpHMD.

2.1.1 D-CpHMD

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

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

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

2.1.2 C-CpHMD

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

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

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

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

可見,當(dāng)λ=1時,殘基i去質(zhì)子化,殘基i的可滴定氫與j無相互作用.

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

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

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

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

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

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

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

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

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

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

該結(jié)構(gòu)因子可近似表示為

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

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

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

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

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

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

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

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

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

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

圖4 基于C-CpHMD的 pKa 計算 (a)滴定坐標(biāo)λ和去質(zhì)子化概率S的軌跡;(b)采用Hill函數(shù)擬合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函數(shù)(廣義Henderson-Hasselbalch函數(shù))擬合S.pKa便是S=0.5時所對應(yīng)的pH (圖4(b)):

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

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

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

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

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

當(dāng)?shù)味▌恿W(xué)采用滿足周期性邊界條件的顯性溶劑時,需要考慮有限尺度效應(yīng)[73].由于采用耦合水離子實現(xiàn)了電中性,有限尺度效應(yīng)僅剩下和水分子模型相關(guān)的離散溶劑效應(yīng)(discrete solvent effect)[66].當(dāng)某個可滴定氨基酸去質(zhì)子化,因離散溶劑效應(yīng)引起的能量變化是

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

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

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

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

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

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

除了副本交換,另一種增強采樣的方法是對生物大分子進行粗粒化(coarse graining,CG),減少模擬體系中粒子的數(shù)量,從而降低了構(gòu)象空間的自由度.該方法通常被應(yīng)用于具有較大空間和時間尺度的生物過程,如蛋白質(zhì)折疊、多肽聚集和物質(zhì)跨膜轉(zhuǎn)運等[89].近幾年,研究者們開始將CG與CpHMD結(jié)合,發(fā)展CpHMD的粗?；Ｐ蚚90-93].值得一提的是,提出Martini粗?；龅腗arrink課題組[92]已在分子模擬軟件GROMACS中實現(xiàn)了CpHMD的粗?；M.

2.2 基于PB的 pKa 預(yù)測模型

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

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

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

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

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

圖5 相對去質(zhì)子化自由能計算的熱力學(xué)循環(huán)Fig.5.Thermodynamic cycle of relative deprotonation free energy calculation.

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

2.3 基于經(jīng)驗函數(shù)的 pKa 預(yù)測模型

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

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

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

2.4 基于機器學(xué)習(xí)的 pKa 預(yù)測模型

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

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

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

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

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

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

3 結(jié)論

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

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