生物大分子過渡態搜索算法及其中的機器學習*

楊建宇 席昆 竺立哲

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

1 引言

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

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

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

更為重要的是,在數據分析層面,尤其是在提取過渡態信息這類理論化學家最關心的問題上,巨大的原子數量導致了維數災難.搜尋過渡態的結構或特征信息是準確刻畫和解釋所采樣本中動力學機制的重中之重.然而,即使是在采樣數據充足的情況下,使用不恰當的分析手段(即機器學習語境下的降維算法),過渡態區域都將被扭曲以致相關信息丟失.

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

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

1) 依賴CV的定向降維.在不具備先驗數據時,依據直覺猜測少量有物理意義且可能重要的坐標,即集合變量(collective variable,CV),強行定向降維到該預選的低維CV空間,而后在CV空間內搜尋過渡態[79-95].代表性方法: 溫和爬升動力學(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依賴的高維搜索.事先不降維,堅持在高維空間內完成采樣和過渡態搜索過程,事后再進行降維分析[96-101].代表性方法有過渡路徑采樣(transition path sampling,TPS)[98-101].

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

近年來快速發展的機器學習及相關衍生算法(如強化學習、生成式建模等),已成功應用于解決諸多傳統的復雜生物問題[102-112],如生物結構預測及生物分子相互作用的研究[105],或基于人工智能開發蛋白質從頭設計算法[106],或借助于機器學習實現蛋白質結構準確預測的trRosetta線上服務[107],或實現生物分子冷凍電鏡高分辨率結構重建的解析算法[108]和蛋白質間相互作用位點的快速預測[109],以及蛋白質與小分子、RNA等復合物結構性質的預測[110,111].因此,將機器學習與現有過渡態搜索算法進行有效融合,有望成為未來過渡態搜索研究實現進一步突破的可行方向.

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

表1 主要過渡態搜索算法的總結分類Table 1.Classification of the algorithms for transition state searching.

2 依賴CV的過渡態搜索算法

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

圖1 (a)依賴集合變量的過渡態搜索示意圖,需由生物分子(以丙戊酸二肽為例)體系所在的高維相空間(phase space)選取少量集合變量CV強行“定向降維”,后在此低維CV空間利用非路徑類方法或路徑方法,找到過渡態(Transition State),并給出微觀機制解釋(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空間中的過渡態搜索,依照采樣開始時的已知信息可分為非路徑算法和路徑算法.非路徑算法以GAD算法為代表,而路徑算法以finite temperature string[82-87]和快速斷層掃描法[88-90]為代表.前者可在僅有一個穩定態已知時開啟過渡態搜索,而后者需事先已知至少兩個穩定態,通過尋找兩個穩定態之間的最小自由能路徑(minimum free energy path,MFEP),而后獲得沿路徑的自由能分布確定過渡態位置.此外,兩者的區別還有,前者采樣過程是主動“爬山”(即向高能區域運動,圖1(b)左紅),而后者是先通過施加外力促使分子強行翻山越嶺得到能量過高的初始路徑(圖1(b)左藍),再設法使路徑“整體下山”,落入附近的最優路徑MFEP (圖1(b)左黑).

2.1 非路徑類過渡態搜索

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

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

2.2 基于路徑優化的過渡態搜索

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

2.2.1 Finite Temperature String

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

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

2.2.2 快速斷層掃描法

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

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

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

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

圖2 (a) PCV構建[120]和TAPS Method[91-95,121]算法原理示意圖;(b)基于傘形采樣方法得到的TAPS算法確定的MEK1由Loop-Out到達Loop-In轉變過程最小自由能路徑(MFEP)的自由能圖景及相應的微觀轉變機制[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方法時,需提供目標生物系統的兩個穩態結構和連接其轉變過程的初始路徑;而后從初始路徑中確定轉變過程中變化較大的所有結構域,并以這些結構域的重原子(圖2(a)中丙戊酸二肽結構中以球形顯示的原子)為參考,通過計算構象間均方根位移偏差(root mean square distance,RMSD)來評估構象差異,并從初始路徑中在保證相鄰構象間適度的差異基礎上,均勻選取構象(即節點)來代表整個轉變過程;接著,基于此少量節點組成的轉變路徑,便可利用PCV的計算公式得到二維的路徑集合變量低維空間: 即PCV-s和PCV-z.其中,對于任意構象x,參照目標路徑計算得到的PCV-s代表其沿路徑方向的投影位置;而PCV-z表示其距離參考路徑的平均距離,見圖2(a)[120].通過在此路徑集合變量空間內,快速完成路徑搜索,將最終確定目標轉變過程的最優路徑(MFEP),如圖2(a)中基于多維度標度方法(multidimensional scaling method,MDS)[122]得到的二維路徑搜索過程展示,從黑色的初始路徑快速搜索到達綠色的最優路徑(MFEP).

此處以丙戊酸二肽由C7eq到C7ax的轉變為例,完整展示TAPS進行路徑優化的主要過程,包括以下四步(見圖2(a)中下方白色框內的TAPS迭代流程).

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

步驟2對于每個節點的采樣軌跡,通過獲取最接近軌跡PCV-z中位值的結構,并按照上輪編號連接為新的轉變路徑(藍色實線).

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

步驟4去除轉變路徑范圍外節點,并在距離較遠的相鄰節點間補入新節點.

最終,通過不斷重復迭代上述1—4步的路徑優化過程,將最終搜索到MFEP并結合傘形采樣等算法[117]得到沿MFEP的自由能景觀分布,進而給出微觀轉變機制解釋和確定相應的過渡態信息.

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

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

3 基于路徑采樣的過渡態搜索

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

圖3 路徑采樣算法的基本原理示意圖 (a)路徑采樣中生成新相空間路徑的shooting move;(b)傳統過渡路徑采樣(左側)的隨機蒙特卡羅采樣與過渡態分析原理[98-101],融合強化學習的路徑采樣(右側)在學習過程中不斷促進采樣起始點選擇向過渡態集中[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 相空間中過渡態的定義committor probability

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

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

而此條件概率在過渡態上pA=pB=0.5 時 將達到其峰值,即過渡態上的點是所有相空間點中最有可能屬于某條反應路徑的.這一點對路徑采樣算法至關重要.

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

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

路徑采樣過程就是不斷迭代選定xsel,而后進行Shooting的過程.經過迭代最終會得到從A到B轉變的路徑系綜TPE[128,129].但傳統TPS和其強化學習新版本在xsel的選擇策略上有所不同.

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

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

3.1.4 從路徑系綜中提取過渡態信息

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

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

3.2 基于強化學習的路徑采樣

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

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

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

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

值得強調的是,無論是傳統TPS還是強化學習路徑采樣,二者的理論基礎都是pA+pB=1,即不允許穩態A和B之間有第3個穩定態存在.這意味著路徑采樣只能處理單個能壘,即只能表征單個過渡態.然而,生物大分子的運動復雜,亞穩定態數量眾多,很難保證已知的兩個穩定態之間只有一個能壘.這也限制了路徑采樣在生物大分子模擬中的應用.

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

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

4.1 現有降維算法

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

圖4 物理化學家需要怎樣的降維算法 (a)現有降維算法范式不保留過渡態信息,不利于機制解析;(b)可能的替代范式,基于生成模型研發可保留過渡態信息的可逆降維算法,并與低維空間搜索過渡態的GAD聯用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 基于生成模型的可逆降維及過渡態搜索

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

當然,這一構想的實現難點是必須保證在降維過程中,低維空間保有和原高維空間一致的動力學特征以及概率密度信息,即保留過渡態信息.這里我們建議參考tICA中直接使用動力學信息進行降維的做法.此外,為保障GAD在低維空間的順利運行,該生成模型應能為低維空間自動擬合出連續可導的自由能面.

5 結論

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

在完全無需事前降維、不依賴集合變量的路徑采樣類算法中,也已出現了通過融入強化學習思想實現自適應的高效率采樣及過渡態分析優秀變體.但只能處理單個能壘和過渡態搜尋的特點限制了這類算法在生物分子模擬中的應用.

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