羅方芳 蔡志濤 黃艷東

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

1 引言

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

以β分泌酶BACE1為例闡述蛋白質(zhì)功能和可離子化基團(tuán)質(zhì)子化/去質(zhì)子化的關(guān)系.BACE1的生物功能是裂解β淀粉樣前體蛋白APP.它與神經(jīng)退行性疾病阿爾茨海默癥密切相關(guān),是典型的結(jié)構(gòu)和功能依賴于pH的蛋白質(zhì).該蛋白的催化中心含兩個(gè)天冬氨酸Asp32和Asp228 (圖1(a)).實(shí)驗(yàn)指出,BACE1僅在一個(gè)狹小的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,兩個(gè)天冬氨酸同時(shí)質(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)一個(gè)可離子化基團(tuán)的質(zhì)子化和去質(zhì)子化達(dá)到平衡,可由以下公式計(jì)算解離常數(shù)Ka:

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

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

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

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

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

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

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

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

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

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

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

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

2.1 CpHMD

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

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

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

2.1.1 D-CpHMD

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

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

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

2.1.2 C-CpHMD

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

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

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

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

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

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

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

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

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

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

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

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

其中,Udir是實(shí)空間靜電相互作用,在庫(kù)侖勢(shì)基礎(chǔ)上增加一個(gè)補(bǔ)償函數(shù),負(fù)責(zé)截?cái)嗑嚯x以內(nèi)的短程靜電相互作用((15)式).Urec最為耗時(shí),為倒格空間(reciprocal space)下求解的長(zhǎng)程靜電能,負(fù)責(zé)截?cái)嘁酝獾拈L(zhǎng)程靜電相互作用((16)式).Ucorr是修正項(xiàng)((20)式)[39].

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

式中m是倒格矢,其表達(dá)式為,其中,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個(gè)維度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的表達(dá)式為

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

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

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

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

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

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

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

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

圖4 基于C-CpHMD的 pKa 計(jì)算 (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時(shí)所對(duì)應(yīng)的pH (圖4(b)):

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

當(dāng)兩個(gè)氨基酸的滴定動(dòng)力學(xué)存在耦合,可將二者看作一個(gè)整體,利用以下公式計(jì)算宏觀 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；植生混凝土對(duì)底部土壤的反濾攔截率達(dá)86.92%。

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

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

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

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

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

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

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

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

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

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

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

實(shí)際上,如果只考慮單個(gè)結(jié)構(gòu),可以用PB方程計(jì)算相對(duì)去質(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è),以上兩個(gè)自由能改變量的差可表示為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

上述基于AI的模型均采用PKAD中的實(shí)驗(yàn)數(shù)據(jù)來訓(xùn)練或測(cè)試模型.然而,pKa-ANI,XGBWMa和PKAI忽略了存在于PKAD的冗余數(shù)據(jù)(例如一個(gè)蛋白質(zhì)有兩組相同的 pKa值),這可能導(dǎo)致過擬合.其次,PKAD中大多數(shù) pKa處于參考值附近,因此測(cè)試結(jié)果并不能反應(yīng)模型真實(shí)的預(yù)測(cè)能力[25].值得一提的是,本課題組創(chuàng)建的測(cè)試集EXP67S不存在以上兩個(gè)問題,可較為客觀地對(duì)模型進(jìn)行評(píng)價(jià)[25].研究發(fā)現(xiàn),除了在實(shí)驗(yàn)和理論相關(guān)性方面仍舊低于CpHMD,DeepKa的預(yù)測(cè)精度明顯高于其他主流 pKa預(yù)測(cè)模型,包括PypKa,PropKa,PKAI和pKa-ANI[29].其中,PypKa代表基于PB的模型,PropKa代表基于經(jīng)驗(yàn)函數(shù)的模型,PKAI和pKa-ANI代表其他AI模型.基于樹的XGB-WMa沒有開放源代碼,所以無法利用EXP67S對(duì)其進(jìn)行測(cè)試.因此,XGB-WMa不參與下面的模型討論.同時(shí)考查精度和速度,圖6展示了5個(gè)模型的預(yù)測(cè)性能.其中,平均絕對(duì)誤差用于表征模型的精度.顯而易見,如果以PropKa的速度和CpHMD的精度作為參照,目前只有DeepKa能提供準(zhǔn)確的高通量 pKa計(jì)算[29].最近,加拿大國(guó)家研究委員會(huì)Sulea課題組[103]比較了現(xiàn)有的7種高通量 pKa預(yù)測(cè)模型,包括基于經(jīng)驗(yàn)函數(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與溫度、壓強(qiáng)一樣是基本的環(huán)境參量.傳統(tǒng)的分子動(dòng)力學(xué)假設(shè)溶劑是中性水(pH=7.0),不考慮其他pH條件;此外,傳統(tǒng)分子動(dòng)力學(xué)假設(shè)電荷是固定的,不受溶質(zhì)靜電場(chǎng)的影響.以上兩個(gè)假設(shè)限制了傳統(tǒng)分子動(dòng)力學(xué)進(jìn)一步探究細(xì)胞中許多與pH相關(guān)的生物過程,而可靠的 pKa計(jì)算將有助于解決該難題.本綜述主要介紹了4類主流的pKa預(yù)測(cè)方法.顯然,對(duì)于不同理論的 pKa預(yù)測(cè)模型,其適用范圍也存在差異.首先,不論何種特定的問題,如果不要求高通量計(jì)算,可采用預(yù)測(cè)精度較高但計(jì)算效率較低的CpHMD.當(dāng)涉及非水溶性蛋白(如膜蛋白)的 pKa計(jì)算,目前理論上可行的模型為基于雜化溶劑[37,56]或顯性溶劑[66,67]的CpHMD.另一方面,需要開發(fā)高通量的 pKa預(yù)測(cè)模型,從而滿足工業(yè)界批量的 pKa計(jì)算需求.由于隱性溶劑的理論局限性和實(shí)驗(yàn)條件的限制,上述的高通量模型僅適用于水溶性蛋白.對(duì)于水溶性蛋白質(zhì)單體的pKa計(jì)算,在所有高通量模型中DeepKa無疑是最優(yōu)的選擇[29,103].若只關(guān)心酸性氨基酸殘基(如Asp和Glu)的質(zhì)子化態(tài),也可考慮PropKa 3.0[24].而對(duì)于主要的4種可離子化氨基酸殘基(Asp,Glu,Lys和His)以外的可滴定基團(tuán)(如Cys和Tyr),可考慮基于PB的模型(如H++[18]和PypKa[21]).

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