生物大分子過渡態(tài)搜索算法及其中的機(jī)器學(xué)習(xí)*

楊建宇 席昆 竺立哲

(香港中文大學(xué)(深圳)醫(yī)學(xué)院,瓦謝爾計算生物研究院,深圳 518172)

1 引言

生物分子實(shí)現(xiàn)功能時,常伴隨著結(jié)構(gòu)的巨大轉(zhuǎn)變,即生物分子的功能性構(gòu)象變化[1-3].利用實(shí)驗(yàn)方法,往往只能獲取上述轉(zhuǎn)變過程前后重要的穩(wěn)態(tài)結(jié)構(gòu),如X射線(X-ray macromolecular crystallography)[4]、核磁共振(nuclear magnetic resonance,NMR)[5]、冷凍電子顯微鏡(cryo-electron microscopy,cryo-EM)[6]等;或者揭示分子結(jié)構(gòu)變化中的部分特征,如熒光共振能量轉(zhuǎn)移(fluorescence resonance energy transfer,FRET)可給出少數(shù)目標(biāo)殘基間的距離變化[7]等.因此,僅依賴實(shí)驗(yàn)方法難以闡明生物分子轉(zhuǎn)變過程的完整信息.

全原子(all-atom)分子動力學(xué)(molecular dynamics,MD)是從原子尺度全面描述生物分子動態(tài)行為的標(biāo)準(zhǔn)手段[8].但和化學(xué)反應(yīng)僅涉及反應(yīng)活性中心內(nèi)的數(shù)十個原子不同,構(gòu)象變化所涉及的原子數(shù)目巨大,極端情況下可包括溶質(zhì)的全部原子,甚至環(huán)境中脂類和溶劑分子的原子[9-36].眾多的原子及其三維坐標(biāo)帶來了兩個重要的瓶頸.

首先,在計算效率方面,復(fù)雜大分子百萬級的原子數(shù)量意味著需要計算萬億級數(shù)量的原子間作用力,即使在目前最優(yōu)的通用硬件上,人們所能完成的MD模擬時長也僅在微秒量級[8,37],距離生物分子的實(shí)際功能性動力學(xué)行為毫秒級的發(fā)生時間仍有巨大差距.為緩解該效率瓶頸,數(shù)十年來,人們發(fā)展了各類增強(qiáng)采樣算法,其中較有代表性的算法包括副本交換[38-45],選擇性溫度積分增強(qiáng)采樣(selective integrated tempering sampling)[46-49]、局部抬升(local elevation)[50-53]、構(gòu)象洪泛(conformational flooding)[54-56]、元動力學(xué)(metadynamics)[57-59]、高斯加速動力學(xué)[60-62]等.

更為重要的是,在數(shù)據(jù)分析層面,尤其是在提取過渡態(tài)信息這類理論化學(xué)家最關(guān)心的問題上,巨大的原子數(shù)量導(dǎo)致了維數(shù)災(zāi)難.搜尋過渡態(tài)的結(jié)構(gòu)或特征信息是準(zhǔn)確刻畫和解釋所采樣本中動力學(xué)機(jī)制的重中之重.然而,即使是在采樣數(shù)據(jù)充足的情況下,使用不恰當(dāng)?shù)姆治鍪侄?即機(jī)器學(xué)習(xí)語境下的降維算法),過渡態(tài)區(qū)域都將被扭曲以致相關(guān)信息丟失.

在已有大量模擬數(shù)據(jù)的場景中,可借助tICA(time-lagged independent component analysis)[63-65]利用已有數(shù)據(jù)中蘊(yùn)含的動力學(xué)信息進(jìn)行降維,或運(yùn)用馬爾可夫態(tài)模型(Markov state models)[66-78]等分析算法提取動力學(xué)信息來應(yīng)對維數(shù)災(zāi)難,并間接推測過渡態(tài)信息.但這類算法中并不直接含有過渡態(tài)的定義,因而超出了本文范疇.對此類算法感興趣的讀者可參看其他綜述[63,66,68,75-78].

在生物大分子模擬領(lǐng)域,因其計算效率低下,數(shù)據(jù)匱乏是常態(tài),因此人們對能高效搜尋過渡態(tài)的采樣算法需求強(qiáng)烈.但受限于維數(shù)災(zāi)難,僅有以下兩類采樣策略可供選擇.

1) 依賴CV的定向降維.在不具備先驗(yàn)數(shù)據(jù)時,依據(jù)直覺猜測少量有物理意義且可能重要的坐標(biāo),即集合變量(collective variable,CV),強(qiáng)行定向降維到該預(yù)選的低維CV空間,而后在CV空間內(nèi)搜尋過渡態(tài)[79-95].代表性方法: 溫和爬升動力學(xué)(gentlest ascent dynamics,GAD)[79-81]、有限溫度弦方法(finite temperature string,FTS)[82-87]、快速斷層掃描法(fast tomographic,FT)[88-90]、基于旅行商的路徑搜索(travelling-salesman based automated path searching,TAPS)[91-95].

2) 非CV依賴的高維搜索.事先不降維,堅持在高維空間內(nèi)完成采樣和過渡態(tài)搜索過程,事后再進(jìn)行降維分析[96-101].代表性方法有過渡路徑采樣(transition path sampling,TPS)[98-101].

盡管上述算法已在一定范圍內(nèi)取得成功,但在面對復(fù)雜生物分子時,仍面臨諸多限制.其中,對于依賴CV的搜索算法,最直接的問題便是如何從較高維度空間中選取合適的CV;而對于非CV依賴的路徑采樣算法,則是計算資源消耗過大和有效采樣率過低的問題.

近年來快速發(fā)展的機(jī)器學(xué)習(xí)及相關(guān)衍生算法(如強(qiáng)化學(xué)習(xí)、生成式建模等),已成功應(yīng)用于解決諸多傳統(tǒng)的復(fù)雜生物問題[102-112],如生物結(jié)構(gòu)預(yù)測及生物分子相互作用的研究[105],或基于人工智能開發(fā)蛋白質(zhì)從頭設(shè)計算法[106],或借助于機(jī)器學(xué)習(xí)實(shí)現(xiàn)蛋白質(zhì)結(jié)構(gòu)準(zhǔn)確預(yù)測的trRosetta線上服務(wù)[107],或?qū)崿F(xiàn)生物分子冷凍電鏡高分辨率結(jié)構(gòu)重建的解析算法[108]和蛋白質(zhì)間相互作用位點(diǎn)的快速預(yù)測[109],以及蛋白質(zhì)與小分子、RNA等復(fù)合物結(jié)構(gòu)性質(zhì)的預(yù)測[110,111].因此,將機(jī)器學(xué)習(xí)與現(xiàn)有過渡態(tài)搜索算法進(jìn)行有效融合,有望成為未來過渡態(tài)搜索研究實(shí)現(xiàn)進(jìn)一步突破的可行方向.

本文將首先回顧依賴CV的過渡態(tài)搜索算法的發(fā)展歷程,厘清其基本原理及潛存問題.隨后,聚焦于非CV依賴的TPS路徑采樣算法,著重介紹其融合了強(qiáng)化學(xué)習(xí)的最新版本.最后,探討一種新型的過渡態(tài)搜索策略,即結(jié)合生成模型和GAD,在保留原高維空間過渡態(tài)信息的低維空間內(nèi)實(shí)現(xiàn)過渡態(tài)搜索.完整的算法總結(jié)已展示于表1中.

表1 主要過渡態(tài)搜索算法的總結(jié)分類Table 1.Classification of the algorithms for transition state searching.

2 依賴CV的過渡態(tài)搜索算法

如前所述,為了準(zhǔn)確闡明生物分子功能性動力學(xué)的微觀機(jī)制,需要在傳統(tǒng)采樣算法的基礎(chǔ)上,發(fā)展可獲取上述轉(zhuǎn)變過程過渡態(tài)信息的過渡態(tài)搜索算法,包括依賴CV[82-95]和非CV依賴算法[96-101]兩大類.對于依賴CV的算法,需在缺乏對體系的先驗(yàn)數(shù)據(jù)和認(rèn)知的條件下,將高維相空間{x}“定向降維”至少量的依據(jù)經(jīng)驗(yàn)或直覺定義的CV上(arbitrary guess).而后續(xù)的計算采樣和過渡態(tài)搜索則發(fā)生在由這些CV構(gòu)成的低維空間(CV1,CV2,···)內(nèi)(圖1(a)).

圖1 (a)依賴集合變量的過渡態(tài)搜索示意圖,需由生物分子(以丙戊酸二肽為例)體系所在的高維相空間(phase space)選取少量集合變量CV強(qiáng)行“定向降維”,后在此低維CV空間利用非路徑類方法或路徑方法,找到過渡態(tài)(Transition State),并給出微觀機(jī)制解釋(mechanism interpretation);(b)非路徑類的GAD算法原理示意圖;(c),(d)兩類路徑類搜索算法原理示意圖Fig.1.(a) Illustration of the flow-chart of the collective variables (CVs) based transition state searching.A low dimensional space must be constructed with the CVs,which are arbitrary a priori guess about the mechanism.The transition state(s) is then determined by either the non-path or path methods.(b) The non-path method GAD.Path methods of (c) finite temperature string and(d) fast tomographic.

低維CV空間中的過渡態(tài)搜索,依照采樣開始時的已知信息可分為非路徑算法和路徑算法.非路徑算法以GAD算法為代表,而路徑算法以finite temperature string[82-87]和快速斷層掃描法[88-90]為代表.前者可在僅有一個穩(wěn)定態(tài)已知時開啟過渡態(tài)搜索,而后者需事先已知至少兩個穩(wěn)定態(tài),通過尋找兩個穩(wěn)定態(tài)之間的最小自由能路徑(minimum free energy path,MFEP),而后獲得沿路徑的自由能分布確定過渡態(tài)位置.此外,兩者的區(qū)別還有,前者采樣過程是主動“爬山”(即向高能區(qū)域運(yùn)動,圖1(b)左紅),而后者是先通過施加外力促使分子強(qiáng)行翻山越嶺得到能量過高的初始路徑(圖1(b)左藍(lán)),再設(shè)法使路徑“整體下山”,落入附近的最優(yōu)路徑MFEP (圖1(b)左黑).

2.1 非路徑類過渡態(tài)搜索

GAD是非路徑類過渡態(tài)搜索的代表性算法,在預(yù)設(shè)的低維CV空間,從亞穩(wěn)態(tài)或任意狀態(tài)出發(fā),可在低維勢能面空間內(nèi),直接完成過渡態(tài)搜索[79-81].如圖1(b)所示,此算法的原理為由低維勢能面空間內(nèi)的任意一點(diǎn)出發(fā),根據(jù)以下規(guī)則:

來確定每輪迭代時移動至下一步的位移方向,即沿勢能函數(shù)梯度變化率的最小方向進(jìn)行小步長移動,最終收斂于鞍點(diǎn)位置(即過渡態(tài)).其中=(F(x),n)n,F(x)為分子體系在根據(jù)當(dāng)前低維CV空間內(nèi)的勢能梯度計算得到的作用力;而n被設(shè)定為趨近于勢能函數(shù)海森矩陣最小特征值對應(yīng)的特征向量,即指向曲率最小方向,其需要基于(1b)式反復(fù)迭代達(dá)到收斂,在此期間,γ則控制H對n變化的影響能力,以此消除勢能函數(shù)中的噪音.簡單而言,(1)式的規(guī)則將引導(dǎo)分子不斷沿勢能坡度最緩的方向逆勢攀登,直至收斂停滯于過渡態(tài).

2.2 基于路徑優(yōu)化的過渡態(tài)搜索

對于基于路徑優(yōu)化進(jìn)行過渡態(tài)搜索的算法,根據(jù)其輸入不同,可主要分為兩類: 1)需要高質(zhì)量預(yù)選集合變量CV的路徑優(yōu)化算法,包括finite temperature string[82-87]和快速斷層掃描法[88-90];2)基于路徑集合變量(path collective variable,PCV)的路徑優(yōu)化算法,即基于TAPS算法[91-95],此方法中避免了高質(zhì)量預(yù)選集合變量的困境,可高效且快速找到最優(yōu)轉(zhuǎn)變路徑.當(dāng)構(gòu)建完路徑優(yōu)化的低維空間后,需要從目標(biāo)系統(tǒng)的兩個穩(wěn)定態(tài)結(jié)構(gòu)出發(fā),產(chǎn)生一條較為粗糙的轉(zhuǎn)變路徑[114-116],而后對此路徑進(jìn)行迭代優(yōu)化(路徑整體下山),并最終收斂于最優(yōu)路徑(MFEP)[82-95];繼而便可通過計算MFEP的自由能圖景,準(zhǔn)確給出微觀轉(zhuǎn)變機(jī)制和過渡態(tài)信息[57-59,117].

2.2.1 Finite Temperature String

當(dāng)基于傳統(tǒng)的增強(qiáng)采樣算法(如steered MD,climber MD,targeted MD等[114-116])快速得到描述目標(biāo)生物分子過程的轉(zhuǎn)變路徑后,前人發(fā)現(xiàn)還需要通過選取合適的集合變量信息,來構(gòu)建低維空間和完成對初始轉(zhuǎn)變路徑的進(jìn)一步優(yōu)化,從而得到最優(yōu)路徑,即最小自由能路徑(minimum free energy path,MFEP).作為研究此類問題中的代表算法,finite temperature string的優(yōu)化策略[82-87]較為簡潔(以swarms-of-trajectories版本為例[87]),見圖1(c).通過對連接轉(zhuǎn)變路徑(由State A到State B)的所有節(jié)點(diǎn),依次分別完成大量(swarms)非常短時長的隨機(jī)初始速率MD采樣后,在預(yù)選的低維空間對采樣結(jié)果聚類,找到出現(xiàn)概率最高的構(gòu)象,作為代表性的采樣節(jié)點(diǎn)(圖1(c)中sampled node).這樣做是為了在路徑上各節(jié)點(diǎn)附近做非常局部的采樣,從而估計各節(jié)點(diǎn)目前所在位置的自由能梯度,等效于讓各節(jié)點(diǎn)沿著當(dāng)前所在位置的自由能梯度最大方向稍作移動(下山),類似于勢能最小化問題中的最速下降法;通過再優(yōu)化節(jié)點(diǎn)分布來保證相鄰節(jié)點(diǎn)間距離相近(equidistant nodes,圖1(c)),進(jìn)而得到新一輪的轉(zhuǎn)變路徑.

通過不斷重復(fù)上述迭代策略,路徑將最終收斂到達(dá)最小自由能路徑MFEP.最終便可通過傘形采樣等[117]方法獲取沿此MFEP的自由能景觀(free energy landscape)[82-87],進(jìn)而給出微觀機(jī)制解釋和得到相應(yīng)的過渡態(tài)信息.

2.2.2 快速斷層掃描法

快速斷層掃描法與前述的finite temperature string方法較為相似,亦需基于經(jīng)驗(yàn)或隨機(jī)預(yù)選取集合變量來構(gòu)建低維空間[88-90],而后在此低維空間進(jìn)行路徑搜索,找到MFEP,如圖1(d)所示:

首先,在選定的低維度空間內(nèi),均勻選取轉(zhuǎn)變構(gòu)象(每個構(gòu)象稱為節(jié)點(diǎn),共N個節(jié)點(diǎn))來代表初始轉(zhuǎn)變路徑(由State A到State B);隨后,對于每個節(jié)點(diǎn),都在垂直于當(dāng)前路徑的超平面空間內(nèi)進(jìn)行相同時長的MD模擬采樣,在采樣過程中還需引入SHAKE算法[118]以避免其離超平面空間過遠(yuǎn),同時,結(jié)合自適應(yīng)偏勢MD方法(adaptively biased molecular dynamics,ABMD)[119]來提高其采樣效率;接著,針對每個節(jié)點(diǎn)的采樣軌跡,直接將采樣的終態(tài)結(jié)構(gòu)進(jìn)行連接,保存為新的轉(zhuǎn)變路徑(如圖1(d)中黑色虛線代表的第i輪結(jié)果和黃色虛線代表的第i+1輪結(jié)果).按照上述流程反復(fù)迭代,將最終得到MFEP,及相應(yīng)自由能景觀分布,從而闡明其微觀轉(zhuǎn)變機(jī)制并確定目標(biāo)過渡態(tài)信息.

2.2.3 基于旅行商的自動化路徑搜索算法

在基于集合變量的搜索算法中,還存在一種基于路徑集合變量PCV的新型算法[120],即基于旅行商問題的自動路徑搜索算法(TAPS).TAPS巧妙地避開了其他路徑優(yōu)化算法中集合變量的選取問題,同時基于并行化和GPU加速,快速得到較高維度空間中的最優(yōu)路徑(MFEP),給出相應(yīng)的微觀轉(zhuǎn)變機(jī)制和過渡態(tài)信息(圖2)[91-95].

圖2 (a) PCV構(gòu)建[120]和TAPS Method[91-95,121]算法原理示意圖;(b)基于傘形采樣方法得到的TAPS算法確定的MEK1由Loop-Out到達(dá)Loop-In轉(zhuǎn)變過程最小自由能路徑(MFEP)的自由能圖景及相應(yīng)的微觀轉(zhuǎn)變機(jī)制[92]Fig.2.(a) Illustration for the construction of PCV and the flow-chart of the TAPS method;(b) TAPS revealed the free energy landscape and the transition states for the transition from the Loop-Out state of MEK1 to its Loop-In state[92].

具體來講,在使用TAPS方法時,需提供目標(biāo)生物系統(tǒng)的兩個穩(wěn)態(tài)結(jié)構(gòu)和連接其轉(zhuǎn)變過程的初始路徑;而后從初始路徑中確定轉(zhuǎn)變過程中變化較大的所有結(jié)構(gòu)域,并以這些結(jié)構(gòu)域的重原子(圖2(a)中丙戊酸二肽結(jié)構(gòu)中以球形顯示的原子)為參考,通過計算構(gòu)象間均方根位移偏差(root mean square distance,RMSD)來評估構(gòu)象差異,并從初始路徑中在保證相鄰構(gòu)象間適度的差異基礎(chǔ)上,均勻選取構(gòu)象(即節(jié)點(diǎn))來代表整個轉(zhuǎn)變過程;接著,基于此少量節(jié)點(diǎn)組成的轉(zhuǎn)變路徑,便可利用PCV的計算公式得到二維的路徑集合變量低維空間: 即PCV-s和PCV-z.其中,對于任意構(gòu)象x,參照目標(biāo)路徑計算得到的PCV-s代表其沿路徑方向的投影位置;而PCV-z表示其距離參考路徑的平均距離,見圖2(a)[120].通過在此路徑集合變量空間內(nèi),快速完成路徑搜索,將最終確定目標(biāo)轉(zhuǎn)變過程的最優(yōu)路徑(MFEP),如圖2(a)中基于多維度標(biāo)度方法(multidimensional scaling method,MDS)[122]得到的二維路徑搜索過程展示,從黑色的初始路徑快速搜索到達(dá)綠色的最優(yōu)路徑(MFEP).

此處以丙戊酸二肽由C7eq到C7ax的轉(zhuǎn)變?yōu)槔?完整展示TAPS進(jìn)行路徑優(yōu)化的主要過程,包括以下四步(見圖2(a)中下方白色框內(nèi)的TAPS迭代流程).

步驟1基于轉(zhuǎn)變路徑節(jié)點(diǎn)間結(jié)構(gòu)差異(dx,i)和節(jié)點(diǎn)編號(i=1,2,···,N)信息,利用PCV[120]構(gòu)建路徑優(yōu)化的二維空間: 沿路徑方向,PCV-s((2a)式)和垂直于路徑方向,PCV-z((2b)式),而后從每個節(jié)點(diǎn)出發(fā)做采樣,采樣時在PCV-s方向加入限制偏勢,阻止分子在平行于當(dāng)前路徑的方向運(yùn)動,但允許其在垂直于當(dāng)前路徑的超平面內(nèi)任意運(yùn)動;同時,為了后續(xù)步驟4補(bǔ)入節(jié)點(diǎn)時能有更多候選構(gòu)象,在PCV-s進(jìn)行元動力學(xué)(well-tempered metadynamics[123])采樣.

步驟2對于每個節(jié)點(diǎn)的采樣軌跡,通過獲取最接近軌跡PCV-z中位值的結(jié)構(gòu),并按照上輪編號連接為新的轉(zhuǎn)變路徑(藍(lán)色實(shí)線).

步驟3經(jīng)步驟1非局部的垂直空間采樣后,節(jié)點(diǎn)順序很可能已發(fā)生改變需要重排.本算法將節(jié)點(diǎn)重排轉(zhuǎn)化為旅行商問題[121],并通過插入虛擬點(diǎn)(即與其他任何節(jié)點(diǎn)間的距離為零)來將旅行商問題的閉環(huán)解轉(zhuǎn)化為節(jié)點(diǎn)順序編號.

步驟4去除轉(zhuǎn)變路徑范圍外節(jié)點(diǎn),并在距離較遠(yuǎn)的相鄰節(jié)點(diǎn)間補(bǔ)入新節(jié)點(diǎn).

最終,通過不斷重復(fù)迭代上述1—4步的路徑優(yōu)化過程,將最終搜索到MFEP并結(jié)合傘形采樣等算法[117]得到沿MFEP的自由能景觀分布,進(jìn)而給出微觀轉(zhuǎn)變機(jī)制解釋和確定相應(yīng)的過渡態(tài)信息.

以TAPS對絲裂原激活蛋白激酶激酶(MEK1)由Loop-Out狀態(tài)轉(zhuǎn)變?yōu)長oop-In狀態(tài)的研究為例(圖2(b)),實(shí)驗(yàn)發(fā)現(xiàn)其在傳遞生物信號中時需經(jīng)歷Loop-Out態(tài)到Loop-In態(tài)的轉(zhuǎn)變,即兩個α螺旋(α0和α1)的局部翻轉(zhuǎn)以及連接螺旋的Loop進(jìn)入激活口袋;利用TAPS方法同時考察上述過程中涉及的所有重要?dú)埢?在較短的采樣總時間(短于32.6 ns)內(nèi)便得到了MFEP(圖2(a)最右側(cè)的MDS結(jié)果內(nèi)的綠色線)[92];沿收斂的MFEP進(jìn)一步得到了相應(yīng)的自由能圖景(圖2(b)),進(jìn)而獲得了主要轉(zhuǎn)變機(jī)制和兩個關(guān)鍵過渡態(tài)結(jié)構(gòu)(TS Ⅰ和Ⅱ).此研究所新發(fā)現(xiàn)的R227:L235及Y229:E255極性接觸作用,也被成功用于解釋實(shí)驗(yàn)關(guān)于R227或Y229的點(diǎn)突變造成MEK1無法激活的現(xiàn)象[124,125].

盡管TAPS算法巧妙地規(guī)避了預(yù)選CV空間定向降維帶來的試錯成本,但仍需選擇計算RMSD所需的原子集作為輸入信息.這意味著在復(fù)雜大分子的過渡態(tài)搜索中,即便TAPS的整體效率相比依賴CV的方法已有大幅提升,它仍在事先對所研究構(gòu)象變化的機(jī)制做出了一定假設(shè).

3 基于路徑采樣的過渡態(tài)搜索

目前所有算法中,只有以TPS為代表的路徑采樣方法在事先對構(gòu)象變化機(jī)制未作任何假設(shè),因?yàn)門PS將構(gòu)象轉(zhuǎn)變路徑直接定義在了高維相空間內(nèi).傳統(tǒng)TPS通過大量隨機(jī)的不外加偏執(zhí)勢的無偏采樣,得到一個過渡路徑系綜(transition path ensemble,TPE),見圖3(a).最終通過對TPE的后處理分析,選取合適的集合變量以描述過渡態(tài)[98-101](圖3(b)左);最近,通過引入強(qiáng)化學(xué)習(xí)范式(reinforcement learning),該方法實(shí)現(xiàn)了自適應(yīng)無偏采樣(圖3(b)右),并采用符號回歸(symbolic regression)完成機(jī)制解析[113,126].

圖3 路徑采樣算法的基本原理示意圖 (a)路徑采樣中生成新相空間路徑的shooting move;(b)傳統(tǒng)過渡路徑采樣(左側(cè))的隨機(jī)蒙特卡羅采樣與過渡態(tài)分析原理[98-101],融合強(qiáng)化學(xué)習(xí)的路徑采樣(右側(cè))在學(xué)習(xí)過程中不斷促進(jìn)采樣起始點(diǎn)選擇向過渡態(tài)集中[113]Fig.3.Schematics of path sampling methods.(a) Shooting move: select a phase space point on the current path,make a small perturbation to this point (redraw random initial velocities) and perform a set of simulations.(b) Path sampling is built upon the committor probability pB.The traditional transition path sampling (left)[98-101] selects shooting points randomly and uses Monte Carlo for sampling;the transition state is characterized through post-analysis: choosing the CVs with the highest and narrowest distribution of P(TP|CV);the new reinforcement path sampling (right)[113] chooses shooting points adaptively and directly learns the committor probability pB with maximized P(TP|x).Symbolic regression of pB is used for mechanism interpretation.

3.1 過渡路徑采樣

3.1.1 相空間中過渡態(tài)的定義committor probability

由于TPS中的路徑直接定義在相空間,相應(yīng)地過渡態(tài)也無法直接套用低維空間中的鞍點(diǎn)(saddle)來具象地表征.假設(shè)我們能通過某些CV定義出兩個穩(wěn)定態(tài)A和B (并同時假設(shè)A和B中間不存在第3個穩(wěn)定態(tài)C),那么A和B之間的過渡態(tài)就能通過committor probability來定義.

對相空間中的任一點(diǎn),都可以從其出發(fā)運(yùn)行大量MD模擬并統(tǒng)計其中有多少比率分子是在抵達(dá)穩(wěn)態(tài)B之前到達(dá)了A,另有多少比率相反在到達(dá)了A之前抵達(dá)了B.這兩種比率pA和pB就是這一點(diǎn)對穩(wěn)態(tài)A和B的committor probability.顯然在不存在第3個穩(wěn)態(tài)的前提下pA+pB=1.相應(yīng)地,過渡態(tài)則可以定義為由相空間內(nèi)所有pA=pB=0.5的點(diǎn)所組成的集合.同時,依據(jù)過渡路徑理論(transition path theory)[96],我們知道對相空間中的任一點(diǎn)x而言,它是屬于連接A和B反應(yīng)路徑,即過渡路徑(transition path,TP)的其中一點(diǎn)的條件概率是

而此條件概率在過渡態(tài)上pA=pB=0.5 時 將達(dá)到其峰值,即過渡態(tài)上的點(diǎn)是所有相空間點(diǎn)中最有可能屬于某條反應(yīng)路徑的.這一點(diǎn)對路徑采樣算法至關(guān)重要.

3.1.2 Shooting move新相空間路徑的生成

假設(shè)已利用傳統(tǒng)增強(qiáng)采樣算法(如climber method/steered MD/targeted MD等[114-116])得到一條連接A到B的轉(zhuǎn)變路徑,便可以在此轉(zhuǎn)變路徑中抽選一個點(diǎn)xsel;隨 后,對xsel做出微擾Δx(典型做法為根據(jù)給定溫度的麥克斯韋-玻爾茲曼隨機(jī)重置所有分子的初始速率),而后以xnew=xsel+Δx為新的初始條件進(jìn)行多次無偏MD模擬采樣.其中,每次MD模擬采樣的終止條件為此采樣路徑到達(dá)了目標(biāo)態(tài)A或B中的一個;當(dāng)這些軌跡中既有到達(dá)過A也有到達(dá)過B態(tài)時,將到達(dá)過A態(tài)的任意路徑和到達(dá)過B態(tài)的任意路徑連接便成為由A態(tài)到達(dá)B態(tài)的轉(zhuǎn)變路徑.該過程被稱為shooting move (圖3(a))[127].

路徑采樣過程就是不斷迭代選定xsel,而后進(jìn)行Shooting的過程.經(jīng)過迭代最終會得到從A到B轉(zhuǎn)變的路徑系綜TPE[128,129].但傳統(tǒng)TPS和其強(qiáng)化學(xué)習(xí)新版本在xsel的選擇策略上有所不同.

3.1.3 過渡路徑采樣的shooting move策略

在原版TPS中,xsel的選擇是完全隨機(jī)的.同時,shooting move的迭代是馬爾科夫鏈蒙特卡羅的串行過程(圖3(b)左).因此,TPS天然欠缺并行化能力.

3.1.4 從路徑系綜中提取過渡態(tài)信息

經(jīng)shooting move迭代得到路徑系綜后,傳統(tǒng)TPS需要用戶自行定義CV來幫助解釋其中蘊(yùn)含的機(jī)制、提取過渡態(tài)信息.根據(jù)(3)式,如果所選的CV能夠較好地表征過渡態(tài),即無限趨近pB,那么P(TP|CV) 應(yīng)該呈現(xiàn)窄而高的分布.但由于P(TP|CV)無法直接計算,需要通過貝葉斯推測間接計算:

其中P(CV|TP) 可直接從TPE計算獲得,P(TP) 需經(jīng)額外長時間無偏采樣算出,而Peq(CV) 是CV上的平衡態(tài)分布,也需通過額外的傘形采樣獲得.在用戶選擇的CV中,以P(TP|CV) 分布最窄最高者最能表征過渡態(tài)和A到B的轉(zhuǎn)變機(jī)制[98-101].

3.2 基于強(qiáng)化學(xué)習(xí)的路徑采樣

仔細(xì)分析原版TPS的后處理分析過程,不難看出其對蒙特卡羅迭代采樣結(jié)果的要求較高,需確保所得TPE在過渡態(tài)附近有充足樣本,但由于其xsel的選擇是完全隨機(jī),這在面臨較大的生物分子體系時是難以實(shí)現(xiàn)的.

因此,Jung等[113]于近期開發(fā)了基于強(qiáng)化學(xué)習(xí)(reinforcement learning)的路徑采樣算法.與原版TPS僅在數(shù)據(jù)處理分析階段隱性地使用(4)式不同,新框架直接將P(TP|x) 用作了強(qiáng)化學(xué)習(xí)中的目標(biāo)函數(shù)(通過最大似然估計將其最大化),用以訓(xùn)練以深度神經(jīng)網(wǎng)絡(luò)表達(dá)的committor probabilitypB(圖3(b)右).因此,在此強(qiáng)化學(xué)習(xí)過程中,P(TP|x)的最大化意味著算法會自適應(yīng)地選擇xsel,自發(fā)將其聚焦至過渡態(tài)附近(即pB=0.5,圖3(b)紅線).

而后續(xù)對轉(zhuǎn)變機(jī)制的解釋,即神經(jīng)網(wǎng)絡(luò)pB物理含義的挖掘則可通過符號回歸(symbolic regression)達(dá)成,將pB(x) 的神經(jīng)網(wǎng)絡(luò)表達(dá)為容易理解的簡單解析式[125,126].

3.3 路徑采樣算法的適用場景

值得強(qiáng)調(diào)的是,無論是傳統(tǒng)TPS還是強(qiáng)化學(xué)習(xí)路徑采樣,二者的理論基礎(chǔ)都是pA+pB=1,即不允許穩(wěn)態(tài)A和B之間有第3個穩(wěn)定態(tài)存在.這意味著路徑采樣只能處理單個能壘,即只能表征單個過渡態(tài).然而,生物大分子的運(yùn)動復(fù)雜,亞穩(wěn)定態(tài)數(shù)量眾多,很難保證已知的兩個穩(wěn)定態(tài)之間只有一個能壘.這也限制了路徑采樣在生物大分子模擬中的應(yīng)用.

4 融合GAD與降維算法的可能方案

經(jīng)過對上述算法的簡單回顧,可以看出近年來依賴CV的路徑搜索算法和非CV依賴的路徑采樣算法都已呈現(xiàn)與計算機(jī)科學(xué)和機(jī)器學(xué)習(xí)算法深度融合邁向自動化的發(fā)展趨勢,但依賴CV的GAD方法尚無相似案例可循.我們推測一個可能的發(fā)展方向是將GAD在低維空間搜索過渡態(tài)的能力與降維算法結(jié)合起來.自然地,這對降維算法的性能提出了新的要求.因此,有必要先對現(xiàn)有降維算法的設(shè)計思想進(jìn)行簡要梳理.

4.1 現(xiàn)有降維算法

降維是無監(jiān)督機(jī)器學(xué)習(xí)的傳統(tǒng)分支,其在生物分子模擬中的廣泛應(yīng)用已有綜述闡明[130],此處不再贅述.但在目前眾多的降維算法中,顯式利用時間序列信息,即動力學(xué)信息,進(jìn)行降維的僅有時間結(jié)構(gòu)獨(dú)立成分分析(time-lagged independent components analysis,tICA)方法[63-65].但經(jīng)tICA降維所得的低維tIC空間已被限定只能是原高維空間的線性組合,而能夠表征躍遷過程和過渡態(tài)的坐標(biāo)很可能是原高維坐標(biāo)的非線性函數(shù).其他現(xiàn)存降維算法,因在降維過程中,只關(guān)注保留高密度區(qū)域信息(即穩(wěn)定態(tài)信息),常會將高維空間過度扭曲以致過渡態(tài)信息丟失(圖4(a)).因此,現(xiàn)存降維算法都無法與GAD聯(lián)用.

圖4 物理化學(xué)家需要怎樣的降維算法 (a)現(xiàn)有降維算法范式不保留過渡態(tài)信息,不利于機(jī)制解析;(b)可能的替代范式,基于生成模型研發(fā)可保留過渡態(tài)信息的可逆降維算法,并與低維空間搜索過渡態(tài)的GAD聯(lián)用Fig.4.Requirements on dimensionality reduction algorithms by physical chemists.(a) Current paradigm for dimensionality reduction and the main difficulties for the transition state searching.(b) Proposed alternative paradigm for transition state searching:combine dimensionality reduction that preserves transition state information with GAD.

4.2 基于生成模型的可逆降維及過渡態(tài)搜索

近年來,可逆神經(jīng)網(wǎng)絡(luò)和生成模型的發(fā)展,為研發(fā)能夠保留過渡態(tài)信息的新型降維算法提供了良好契機(jī).首先,通過可逆神經(jīng)網(wǎng)絡(luò),我們可以期望利用深度學(xué)習(xí)訓(xùn)練出一個可以進(jìn)行雙向映射的生成模型,即在將高維的全原子軌跡信息映射到某一低維空間的同時,擁有把生成的低維空間樣本逆投影回原空間的能力.這樣便可利用GAD在低維空間搜得鞍點(diǎn)結(jié)構(gòu),再經(jīng)逆投影自動得到完整的高維過渡態(tài)結(jié)構(gòu).

當(dāng)然,這一構(gòu)想的實(shí)現(xiàn)難點(diǎn)是必須保證在降維過程中,低維空間保有和原高維空間一致的動力學(xué)特征以及概率密度信息,即保留過渡態(tài)信息.這里我們建議參考tICA中直接使用動力學(xué)信息進(jìn)行降維的做法.此外,為保障GAD在低維空間的順利運(yùn)行,該生成模型應(yīng)能為低維空間自動擬合出連續(xù)可導(dǎo)的自由能面.

5 結(jié)論

生物分子功能機(jī)制的有效調(diào)控有賴于對其轉(zhuǎn)變過程微觀機(jī)制的全面考察,其中以獲取其主要轉(zhuǎn)變路徑中的過渡態(tài)信息最為關(guān)鍵.當(dāng)預(yù)設(shè)靜態(tài)集合坐標(biāo)較為容易、可強(qiáng)行定向降維時,前人開發(fā)的GAD算法、finite temperature string和快速斷層掃描法,已成功闡明了諸多生物過程的微觀轉(zhuǎn)變機(jī)制,但當(dāng)面對復(fù)雜轉(zhuǎn)變過程時,仍易出現(xiàn)預(yù)設(shè)集合變量常不合理,需要消耗大量資源試錯.近年出現(xiàn)的基于旅行商的自動路徑搜索算法TAPS,則有效避免了集合變量的預(yù)設(shè)問題,還在并行化和GPU加速的基礎(chǔ)上,提升了自動化程度和過渡態(tài)搜索效率.

在完全無需事前降維、不依賴集合變量的路徑采樣類算法中,也已出現(xiàn)了通過融入強(qiáng)化學(xué)習(xí)思想實(shí)現(xiàn)自適應(yīng)的高效率采樣及過渡態(tài)分析優(yōu)秀變體.但只能處理單個能壘和過渡態(tài)搜尋的特點(diǎn)限制了這類算法在生物分子模擬中的應(yīng)用.

因此,研發(fā)可保留過渡態(tài)信息的新型降維算法或是將機(jī)器學(xué)習(xí)進(jìn)一步融入過渡態(tài)搜索的可行方向.在此,我們建議基于生成模型研發(fā)此種高質(zhì)量降維方法,并將之與GAD聯(lián)用,從而做到從任意狀態(tài)出發(fā),快速捕捉其周圍的過渡態(tài)信息.