分子體系自由能地貌圖的變分分析及AI 算法實現*

杜泊船 田圃2)?

1) (吉林大學生命科學學院,長春 130012)

2) (吉林大學人工智能學院,長春 130012)

精確描述復雜分子體系的自由能地貌圖是理解和操控其行為,并進一步實現分子設計制造工業化的重要基礎.刻畫高維空間自由能地貌圖的主要挑戰是其往往在不同時空間尺度上具有多個層次,每個層次都可能有不止一個亞穩態被相應的自由能壘分開,且跨越路徑有可能不止一條.另外很多體系涉及非線性行為,這使得理論解析和直接使用分子模擬都有很大困難.針對這些挑戰,多年來研究者們發展了多種多樣的增強采樣方法,但往往需要很多經驗選擇和操作,從而一方面使得研究進程較為緩慢,另一方面也讓誤差控制成為困難.變分雖然在物理、統計和工程中已經被廣泛應用并取得巨大成功,但在復雜分子體系中的應用卻隨著神經網絡的發展剛剛開始.本文將對這些探索性工作的主要方向、進展和局限進行簡要總結,也對將來的可能發展給出展望,希望能夠激發更多對基于變分的分子體系自由能地貌圖人工智能算法的關注和努力,促進大分子藥物、分子生物機器等實踐應用的發展.

1 引言

大多數復雜分子,尤其是生物大分子體系,都是通過構象變化或者在一定尺度上的相變實現其功能的[1–6].和諸多分子的實驗合成與表征測試過程相比較,一方面分子模擬的代價往往更低廉;另一方面很多生物大分子復合體的大量合成非常困難甚至不可能,或者在能夠獲取的前提下動態表征很難實現.因此分子模擬被廣泛用于研究復雜分子體系[7–9].決定分子體系各種行為的基礎是對應的自由能地貌圖,因此對其準確刻畫成為必要.實現這一目標的主要挑戰是復雜分子體系一般不止一個亞穩態并且相互之間有較高的自由能壘.所以對典型的復雜分子體系(如核糖體),想要從全原子分子模擬中完成所有亞穩態的充分采樣,觀察對應的構象變化過程往往需要生成毫秒級甚至更長時間的模擬軌跡[10].這對百萬或更多原子的分子體系一方面算力需求很難滿足,另一方面在高維空間中理解所生成的軌跡也很不容易.因此人們發展了各種各樣的增強采樣方法[11–26]和軌跡降維分析方法[27].增強采樣方法大致可以分為兩大類,一類是保持分子體系的玻爾茲曼分布不變,通過改變溫度加速分子體系跨越能壘的方法[12,13].另外一類則是通過加持偏置力/勢(bias force/potential)(如元動力學方法[8]、自適應偏置力方法[28]),這類方法的主要依據是雖然一般分子體系的總自由度數目成千上萬甚至更多,但在跨越能壘的時間尺度上很多局部的原子運動都由于時間尺度的分離而成為近似白噪聲,使得體系在對應時空間尺度的運動可以用較少的反應坐標(reaction coordinates,RC)或者集合自由度(collective variable,CV)成功描述,下文中統稱集合自由度(CV).這類采樣算法的主要困難是集合自由度的構建沒有系統的方法和步驟,研究者往往依靠物理直覺選擇部分體系自由度進行組合嘗試.由于我們生活中感受到的都是三維空間中的物理存在,所以在體系維度升高后直覺判斷的準確性會大打折扣.如何準確地構建有效的CV是目前復雜分子體系模擬中尚未解決的重大挑戰之一.集合自由度空間中主要有3 類互相關聯的問題,其一是準確描述體系的集合自由度的構建;其二是繪制出該空間內主要亞穩態所在的構象空間位置和統計權重,并計算不同亞穩態之間的轉化速率;其三是構建不同亞穩態之間的過渡路徑.這幾類問題的傳統應對策略已經被多個優秀綜述覆蓋[14,29–41],本文主要簡述變分及其神經網絡實現在這些領域的應用,限于作者所熟悉研究工作的范圍,會遺漏一些優秀的研究進展,在此表示歉意.

本文的內容組織如下,首先將對CV、變分和神經網絡及自動微分進行簡要說明,其次對目前已有的針對復雜分子體系自由能地貌圖的主要變分構造方法加以討論,再次對這些基于變分的和其他CV 相關的神經網絡方法進行比較分析,最后展望將來的發展.

2 集合變量、相關神經網絡架構、自動微分和變分簡介

對一個在給定溫度T和勢能U(R) 下的分子體系,用R表示其 3N－3 維坐標,則平衡態玻爾茲曼分布為μ(R∫)=e-βU(R)/Z,其中β=(kBT)-1為逆溫度,Z=dRe-βU(R)為配分函數,kB為玻爾茲曼常數.在較長時間尺度上,這個分子體系的動力學一般可以使用比 3N維度低很多也平滑很多的d(d?N) 維自由能面描述,對應一組由原來坐標R的函數構建的新變量s(R)=(s1(R),s2(R),···,sd(R)),分子體系自由能在這個低維空間也可表示為

人們通常稱這組新變量s(R) 為集合變量,δ(·) 表示δ函數.

神經網絡是目前人工智能技術浪潮的核心理論方法,簡而言之是由多個神經元組成的復合函數網絡.每個神經元可以接受不同維度的輸入,經過線性組合和非線性激活函數作用后輸出.雖然原則上神經元之間的連接可以是任意的,但受視神經分層分布的啟發和隨之帶來的并行計算方便,常用的各種神經網絡架構都是層狀結構.神經網絡最有力的特點是只需要一個隱藏層,足夠多神經元組成的網絡就可以無限逼近任意函數映射,這就是著名的萬能逼近理論(universal approximation theorem)[42–44].但這個理論并沒有指出如何在有限的神經元數目的情況下有效擬合各種映射,所以其發現雖然在很大程度上增強了人們使用神經網絡擬合各種函數映射的信心,卻并沒有迅速推動其在諸多實際問題中的應用.后來多種神經網絡架構(卷積[45]、循環[45]、殘差[46]、注意力機制transformer[47]和擴散模型[48])的發展推動了神經網絡在多個學科領域應用的爆發.當然另外一個不可或缺的基礎是自動微分的發現[49]和在神經網絡中的成功應用[50],這使得理論上基于任意階導數的優化方法都能夠被有效用來訓練神經網絡參數,當然實際應用中由于算力和內存限制,人們往往限于使用基于一階和二階導數的優化方法,諸多具體實例和相關文獻可以參考 PyTorch 中的 Optim 模塊.如下所述,在眾多神經網絡架構中,分子體系自由能地貌圖刻畫中應用最為廣泛的是自編碼器(auto-encoder)架構[51](如圖1 所示),該架構把高維輸入映射到一個低維空間的降維部分被稱為編碼器(encoder),而隨后從低維逆向映射到高維(一般與輸入同維度以方便訓練)空間的部分則被稱為解碼器(decoder).這顯然與人們試圖在更低維度空間理解復雜分子體系的目標在形式上較為吻合.雖然在架構形式上非常相似,但變分自編碼器(variational auto-encoder,VAE)[52]的目標和訓練過程卻與自編碼器顯著不同,其中的隱變量(z)是個概率分布而非特定構型.如果分別用?和ψ表示編碼器和解碼器網絡中的參數,q?(z|x)和pθ(x) 表示隱變量(z)和(x)的分布,則似然函數可表述如下:

圖1 自編碼器神經網絡架構示意圖,藍色部分表示編碼器(encoder)函數 f(·),橙色部分表示解碼器(decoder)函數 g(·),維度最低的綠色表示中間隱藏層(z),對自編碼器,損失函數是輸出()與輸入 xi 的差別的函數(也可以加正則化項,如參考文獻[58] (5)式所示),每一個輸入數據點對應隱藏層空間的一個點Fig.1.Schematic representation of an auto-encoder neural network.The blue part on the left represents the encoder,the orange part on the right represents the decoder,and the middle green layer is the hidden layer (z).The loss is always a function of the difference between the input and the output vectors (xi and ),one may add some form of regularization when necessary (e.g.Eq.(5) in Ref.[58]).

其中,DKL(qφ(z|x)||pθ(z|x))≥0,所以Lθ,φ(x)=Eq?(z|x)[logpθ(x,z)-logq?(z|x)] 就是似然函數的下界,也稱為證據下界(evidence lower bound,ELBO)或變分下界,是變分優化的目標,而非自編碼器中解碼器輸出構型與數據中實際構型差別的函數.為了對隨機隱變量(z)對自動微分,Waterfall等[53]發展了二次參數化技巧(reparameterization trick).

變分的歷史非常悠久,也是諸多理工科研究生的必修課程內容.變分在物理、統計和工程領域都已經取得了非常廣泛和成功的應用[49,54],如量子力學中的 Releigh-Ritz 方法[55]也正是本文中要討論的分子體系變分計算的基礎.另外統計學中的大量應用展示了變分推斷方法同采樣計算相比高效、收斂性較好和更容易擴展的特點[56,57].在神經網絡廣泛應用之前,由于各種未知統計分布的解析和(或)參數化構造較為困難,因此基于平均場的變分成為統計變分分析中最為常用的近似[56].但在分子模擬及其增強采樣中的應用卻在最近十多年才陸續發生.原因主要有兩點,其一是和很多統計模型與工程應用不同,分子體系中的集合變量很難找到直接的方程或模型解析描述,其二是傳統數值擬合方法(如最小二乘法[58])中導數計算昂貴且精度不高,各種優化方法實現困難,而且在變量較多(大于10 個)時會收斂困難[53].不過最近十多年以來基于自動微分[49]的多個人工智能框架Pytorch[59],Tensorflow[60],PaddlePaddle[61]迅速發展成熟,與之伴隨的神經網絡架構[62]也得到了迅猛發展.這使得在擁有較為充足數據的前提下,任意函數的穩健擬合成為可能,因此增強采樣和軌跡分析的變分應用也隨之發展.傳統上人們探索復雜分子體系自由能地貌圖的主要手段是(加速)采樣,變分的突出優點是用優化取代采樣過程,從而顯著提高效率.現代神經網絡架構的強大擬合能力和基于自動微分的各種優化方法的結合為變分在復雜分子體系中的應用提供了巨大的潛力空間.這也正是本文想要討論的話題.

3 分子體系集合變量空間的變分方法

同物理學、工程和統計應用比較,變分在復雜分子體系自由能地貌圖應用相對較少,主要是近十多年的工作,不過目前正在迅速增長中.目前的發展大致可以分為利用轉移矩陣算子特征值和特征向量頻譜分解分析(spectral decomposition analysis)的變分構建[63–68];基于自由能壘跨越概率時間關聯函數的變分[69–71];利用偏置勢(bias potential)的變分構建[72];不受線性假設局限的可匯集性(lumpability)與可分解性(decomposobility)泛函變分構建[73];基于過去-將來信息瓶頸的變分構建[74,75];同時考慮粗粒化、集合變量和增強采樣的自適應[76];以及直接利用變分自編碼器的分析[77],這些方法的簡要總結比較見表1.具體如下所述.

表1 復雜分子體系低維隱空間的變分方法簡要總結,表中所述集合空間問題類別是指引言中提到的三類問題Table 1.A brief summary of variational methods for low-dimensional hidden spaces in complex molecular systems.The category of collective space problems mentioned in the table refers to the three types of problems defined in the introduction.

3.1 頻譜分解分析

在嚴格馬爾可夫過程和細致平衡假設下,針對給定的子態構象空間劃分,Perez-Hernandez 等[65]發展了利用演化算子P(propagator)特征函數自相關構建的變分實現了對最慢動力學過程集合變量(CV)的逼近,分子體系動力學可以被下式表述為演化算子特征函數?i(i=1,2,···,∞) 的疊加:

其中ti是和第i個特征值λi(τ)=e-τ/ti對應的時間尺度,尖括號代表標量積,〈ψi,ρt〉的結果表征概率密度ρt和ψi的重疊程度,體現了第i個特征函數對總體動力學的貢獻.因為ψi=μ-1(x)?i,也可以認為概率密度函數ρ是基于特征函數?i展開的.顯然隨著τ→∞,概率密度會趨于平衡態,(3)式中只有第1 項有貢獻,對應于λ1=1 .如果人們感興趣的時間尺度τ?td+1,則分子體系的動力學主要取決于對應于 (λ1,λ2,···,λd)的d個特征函數,也對應于前面(見方程(1)中的s(R) 定義)所說的d個集合變量.(3)式可近似為

對于分子體系坐標的任意函數f(x),其自相關函數可以表述為

顯然如果取f=ψi(x),則

顯然(8)式中轉化過程χ0和χ1一般是非線性并且未知的.神經網絡的可訓練萬能逼近能力為其構建未知轉換的提供了可能性.VAMPNets[67]正是在這種思考下構建的.對于給定的χ0和χ1變換可以構建3 個方差矩陣:

這些方差矩陣被用來構建了一個VAMP-2 打分[66]:

該分值最大化時對應在轉化后的低(d)維空間分子構象分布被準確復現.以這個分值作為損失函數的神經網絡通過訓練就有可能實現從體系原始高維坐標向低維空間較為準確的映射,實際上起到了低維空間模糊分類器的功能,消除了前述變分理論[65]中對人工聚類及以前各個步驟的需求.具體實現架構如圖2(a)所示[67].對丙氨酸二肽體系,Mardt 等[67]利用圖2(b)的特定架構,設定低維空間類別數目為6(也嘗試了從2—8 的其他類別數目),首先從250 ns 的分子動力學模擬軌跡中每皮秒提取一幀得到250000 個構型,并通過和第一幀對齊除去分子的整體平移和旋轉.使用十個重原子的三維空間(即長度為30 的向量)坐標作為神經網絡輸入,取延遲時間τ=40 ps (也嘗試了從4—32 ps的其他延遲時間),通過最大化VAMP 打分,成功實現了在二維二面角空間 (φ,ψ) 的構象聚類.他們同時使用MSM 流程聚類,當構象類別數目小于20 時得到VAMP 打分都低于VAMPnets 的結果.此外Mardt 等[67]還嘗試分析了簡單雙勢阱和NTL9蛋白折疊軌跡,均展示了和原來人工復雜流程可比擬的準確性,也說明這個思路有望在將來通過逐步發展真正實現自動分析分子模擬軌跡得到動力學特性的可能性.不過目前該方法還不夠成熟,尚不能用于多系綜組合數據[79–84],也不能有效集成模擬軌跡與相關實驗數據,另外還缺乏嚴格清楚定義的誤差估算指標[85–87].但該研究結合變分理論和非線性的神經網絡擬合,取代了原來 MSM 方法管線中一系列復雜步驟,并在簡單體系中實現了首次成功應用,是人工智能用于分析復雜分子體系軌跡的重要進展.VAMPNets 的神經網絡架構較為簡單,鑒于圖神經網絡[88–91]和注意力機制[92]在網絡型數據中的優異表現,考慮到復雜分子體系可以被視為由相互作用的單元構成的網絡圖,Brooks 等[93]構建了包含這兩種架構要素的 GraphVAMPnet,該模型實現了更高精度的構象嵌入表示,也能夠通過注意力機制給出蛋白質中對結構聚類起決定性作用的重要氨基酸.在20-氨基酸的 Trp-cage 蛋白,35-氨基酸的 Vilin 蛋白和 NTL9 蛋白軌跡上的成功應用展示了這些神經網絡構架改變的好處.

圖2 (a) VAMPnets 構建VAMP 打分((10)式)的神經網絡總體架構示意圖;(b)丙氨酸二肽軌跡分析實例中的典型神經網絡架構,各層神經元數目為 32-22-16-9-6,前兩層使用10%的dropout,除最后的softmax 層外,其余各層激活函數均使用Relu[67]Fig.2.(a) Schematic illustration of VAMP score construction from VAMPnets (see Eq.(10)).(b) A typical neural network architecture for analine dipeptide analysis,with the number of neurons being 32-22-16-9-6 for five layers.The first two layers utilized a 10% dropout.Relu was selected as the activation function for all layers except the last softmax layer[67].

隨著人們使用電子顯微鏡解析生物大分子復合體的能力越來越強,如何解釋這些復合體的動力學過程變成了亟待解決的問題.為了增進處理較大分子的能力并在將來能夠有可能延伸到大復合體,Noé等[94]結合獨立馬爾可夫分解方法(independent Markov decomposition,IMD)構建了由多個獨立的VAMPNets 構成的iVAMPNets.其中不同獨立模塊的劃分由一個可訓練的MASK 實現,通過競爭訓練使每個不同的子網絡僅處理不與其他子網絡相互重疊的部分.雖然該方法在Synaptotagmin-C2A 蛋白質分子中成功應用,但顯然這種處理僅適用于不同子模塊間耦合程度較弱的情況,距離準確描述不同組成分子之間有較強關聯的復雜復合體仍然有較大距離.利用VAMPNets 輸出的子構象空間(狀態)概率,Kleiman 和Shukla[68]嘗試了結合3 種不同后續處理,包括最小計數(least count,LC),多目標強化學習(multiagent reinforcement learning-based,MA REAP)和最大熵(MaxEnt),顯著促進了構象空間搜索能力.這3 種方法的宗旨基本一致,就是利用前期生成的軌跡對VAMPNets 進行初步訓練后,在后續的采樣中按照上述不同標準聚焦前期采樣最少訪問的構象空間,從而實現更進一步的增強采樣.其中最大熵和VAMPNets 的結合在促進采樣的同時消除了聚類步驟.

3.2 自由能壘跨越概率時間關聯函數的變分

弦方法[95–99]和過渡路徑理論(transition path theory,TPT)[40]致力于尋找不同亞穩態之間過渡路徑及其過渡態的關鍵細節.不過這些方法在得到最低自由能過渡路徑的同時,卻不能直接給出人們非常感興趣的路徑上任意一點的 自由能壘跨越概率.針對此問題,文獻[71,100]基于兩個亞穩態之間的凈向前反應通量構造了自由能壘跨越概率時間關聯函數,發展了通過變分最小化該函數獲得最佳過渡路徑并同時給出自由能壘跨越概率的方法.對兩個亞穩態 A和B,集合變量空間從A 到B 在時間步長τ基于算子Pτ(s′|s) 的向前演化可表示為

其中ρ(s,t)和ρ(s′;t+τ) 分別對應于時刻t(t+τ)在路徑位置s(s′)處的概率密度.則自由能壘跨越概率q(s),即從s開始最終到達亞穩態B 并且在此前從未到達亞穩態A 的所有過渡路徑概率之和,可定義如下:

則凈向前(從 A到B)反應流為

也可以表達為自由能壘跨越概率的自相關函數:

其中二次方形式可以作為任意給定始末態時嘗試自由能壘跨越概率q(s′) 的變分優化目標.該方法使用基組展開,通過優化系數來達到變分優化的目標,在模型雙勢阱問題中展示了簡化子空間(CV 空間)中理想一維反應坐標走向沿著 自由能壘跨越概率梯度,與高維空間中的 Kramers-Langer 理論[101,102]一致.文獻[100]是針對過渡路徑變分構建的首次嘗試,并在雙勢阱問題和丙氨酸二肽中展示了應用.由于變分函數限于選定基組函數的線性組合空間,其結果顯然會受到基組選擇和線性組合的制約.Chipot 等[69]將自由能壘跨越概率時間關聯函數的變分方法延伸到了神經網絡(variational committor-based neural networks,VCN),從而可以擬合任意非線性映射.同基于特征值變分優化的 VAMPNets 相比較,在雙勢阱體系和N-acetyl-N′-methylalanylamide 異構化過程中均得到一致結果.不過顯著不同的是 VCN 需要已知始末態,針對的目標是一對始末態之間的過渡路徑,而 VAMPNets 則是從軌跡數據開始的無監督學習.另外一點是有時候人們最感興趣的慢過程可能不是分子體系中最慢的過程,這種情況下顯然VCN 更為適合.這兩類方法可以協同使用從而結合其各自優勢,當然也有可能在將來集成到更復雜的神經網絡架構中.

3.3 基于偏置勢的變分

在給定CV 的前提下,Valsson 和Parrinello[103]構建了一個基于CV 空間偏置勢V(s) 的泛函:

其中p(s) 是一個自由選擇的目標分布,這賦予人們使用該泛函的靈活性(當然也伴隨著選擇的挑戰).該泛函是一個凸函數并且不隨偏置勢任意給定的有限常數的改變而變化.用F(s) 表示體系自由能,則當V(s)=-F(s)-(1/β)logp(s) 時,泛函?[V] 取極小值,因此在選定p(s) 的前提下通過參數化的V(s),以?[V] 極小值為目標的變分優化即可求解自由能地貌圖.該方法使用線形基組組合在丙氨酸三肽分子中成功應用.另外,該泛函同Kullback-Leibler (KL)散度(DKL)的關系如下所示[46,104,105]:

其中PV和P0分別是 偏執勢為V和0 時體系的概率密度分布.由于凸函數特性,使得偏置勢與自由能面有確定關系的駐點也是其極值點.因此通過參數化偏置勢,就可以對參數實施變分優化從而求解自由能面.這在原理上比元動力學采樣方法要高效很多,不過,其表現受限于所選CV 在較長時間尺度上描述自由能面的能力.為了克服對該泛函線性展開可能出現的一些麻煩(比如自由能變化劇烈的區域需要很多項才能實現較好擬合,集合變量增大時需要變分優化的參數空間指數增長),Bonati 等[72]用神經網絡表示偏置勢泛函,在給定的集合變量定義下通過優化神經網絡參數實現,如下所示:

泛函數值微分需要統計平均((17)式中的尖括號表示系綜平均),因此需要采樣獲取.直接高精度確定最低點較為困難,因此Bonati 等[72]在實現過程選擇獲得達到一定近似程度的偏置勢,評判的標準選用了pV(s)和p(s) 在迭代次數n時的KL 散度距離:

顯然,此過程在數值實現中需要選定兩個參數,一個是選定每次迭代計算KL 散度之間的模擬更新次數,另一個是每次更新時學習率調整的幅度.為了集成CV 構建和偏置勢優化,Bonati 等[106]利用VAMPNets 的VAMP 打分作為損失函數,利用深度神經網絡和TICA (time-structure based independent component analysis)結合生成CV,隨后在更新的CV 空間采用OPES[100]增強采樣思路,實現了CV 優化和自由能地貌圖收斂的迭代.他們在丙氨酸二肽、chignolin 蛋白折疊和材料結晶過程的成功展示了該方法的應用[106].

3.4 基于可匯集性(lumpability)和可分解性(decomposability)的非線性變分描述

由于馬爾可夫假設和特征函數構建中的線性假設,基于頻譜分解分析的變分優化無法正確處理非馬爾科夫過程[40]和線性無關特征函數之間的非線性關聯,這些根本上的局限無法在后期變分優化中被消除.針對這個問題,Bittracher 等[73]通過延伸過渡流形理論(transition manifold theory)發展了不包含任何線性假設,只關注于長時間尺度分子體系行為,顯式包含誤差量且在可逆體系中互相等價的條件,lumpability和decomposability(詳見文獻[73]的definition 3.2,3.4),這兩個條件都可以作為損失函數變分.此外該變分在近似損失函數時只要求在集合變量子空間的稀疏采樣,而且損失函數的蒙特卡羅積分誤差取決于集合變量子空間而非原高維空間的方差,這會帶來巨大的算力節省.該理論和過渡路徑理論的連接仍然有待闡明.另外這些理論上的優勢在百萬級甚至更大的復雜分子體系如何得以實現也有待于進一步探索.

3.5 過去-將來信息瓶頸模型

Wang 等[74]將分子體系中的集合變量空間視為其演化過程中的過去-將來信息瓶頸(pastfuture information bottleneck,PIB[107,108]),對 給定分子體系任意時刻坐標X和下一時刻坐標X?t,通過瓶頸變量χ(與集合變量類似的分子體系低維空間描述)分別和編碼器P(χ|X) 與解碼器P(X?t|χ) 聯系 (注意文獻[74]中結果部分第1段把坐標X誤解釋 為N個粒子 體系中的d維(1?d?N) 表示,容易引起混亂).PIB 的目標是瓶頸變量χ相對于過去應該盡量簡單但對于將來則應該有盡可能好的預測力,Wang 等[74]據此構建了如下優化目標:

其中I(χ,X?t)和I(X,χ) 分別 表示瓶頸變量與X?t和X的互信 息,常數γ∈[0,∞) 用來平衡瓶頸變量χ的復雜程度和預測力.進一步通過選擇確定性的線性編碼器,則第2 項可以忽略.他們然后利用Gibbs 不等式構建了可變分優化的PIB 下限近似:

其中H和C分別表示香農和交叉熵,Q?為隨機深度神經網絡構建的解碼器.由于選擇Pθ為確定性線性編碼器,香農熵項退出優化目標,可得更新變分下界:

其中?為隨機神經網絡中的變分優化參數.對平衡態軌跡{X1,···,XM+k}(Xn和Xn+k之間的時間間隔為 ?t),方程(20)可被離散為

其中χn從P(χn|Xn) 中采樣得到.對于有對應偏置勢{V1,V2,···,VM+k}下模擬的軌跡則可在假設偏置勢不改變解碼器的情況下近似表述為

實際計算中Wang 等[74]選擇用坐標的線性基組組合得到CV,首先對平衡態軌跡通過逐步增加 ?t觀察基組各項的權重變化,并取其趨于穩定后最小的 ?t.隨后則按照方程(24)和(25)計算偏置勢并重新估算機組系數,反復迭代:

其中w=eβVbias,Pu(χ) 是沒有偏置勢的情況下χ的平衡態分布.簡單的確定性線性編碼器在帶來方便的同時也在一定程度上限制了該方法的靈活性,但PIB 的優點之一是原則上沒有其他線形假設,不過在PIB 思路下(見(19)式)使用非線性編碼器后的變分優化方法仍有待發展.該方法在苯-溶菌酶復合體模擬中獲得了成功,在幾百納秒的加速模擬中觀察到了幾百毫秒常規模擬所觀察到的解離過程.Beyerle 等[75]后來使用該方法成功描述了雙勢阱模型和苯甲酸在雙分子層膜中擴散這兩個分別由能量和熵主導的過渡路徑,進一步展示了該方法的穩健性.

3.6 變分自適應

與前述變分方法主要關注集合變量和偏置勢不同,對有明確集合變量的情況,Zhang 等[76]結合生成式深度學習和基于能量模型[109](energy based models,EBM)發展了對抗密度估計變分,直接計算自由能地貌圖中的概率密度分布.將平衡態真實自由能對應的概率分布記為p,在集合變量空間的參數化自由能地貌圖和對應的分布分別記為Fθ(s)和pθ(s),則KL 散度DKL(p||pθ)對θ的導數可表示為

其中〈f(x)〉p(x)表示函數f(x) 在分布p(x) 下的期望值.(26)式和對抗神經網絡[76,110]高度相似,因此在原文中被稱為變分對抗密度估計 (variational adversarial density estimation,VADE).在實際操作中可以用粗粒化實驗數據PFG(s) 取代真是分布p.再通過粗粒化模擬計算〈β?θFθ(s)〉pθ(s).對于集合變量維度較高的情況,由于直接采樣計算代價過于昂貴,Zhang 等[76]通過加入可訓練生成神經網絡模塊作為神經采樣器(neural sampler)qψ,采用下式實現變分訓練:

加速分子體系自由能地貌圖統計概率分布參數的訓練.對于沒有集合變量函數的更一般情況,通過加入了強化學習模式,較好地解決了固定偏置勢在動態采樣中尷尬的同時,實現了粗細兩個不同粒度的有效采樣補充.這些方法都被集成在SPONGE[111]平臺上.

3.7 變分自編碼器的直接應用

上述構建顯式變分優化目標函數的做法能夠給出更有效的物理圖像,神經網絡主要用于擬合其中未知非線性映射過程.不過即使沒有直接顯式變分目標函數的構建,變分的思想依然可以被利用.最簡單的做法就是直接使用變分自編碼器VAE 架構[77]對自己感興趣的目標數值分布進行優化,同時在生成的隱空間(對應于分子體系集合變量空間)展開一系列增強采樣的操作,必要時再引入迭代機制.

Ribeiro 等[112]發展了重配權變分貝葉斯增強采樣(reweighted auto encoded variational Bayes for enhanced sampling,RAVE)方法,通過隱空間分布和模擬的KL 散度優化自編碼器,更新偏置勢模擬后迭代優化直至收斂,實現了獨立于傳統方法的隱空間增強采樣.針對在MSM 模型中使用過渡路徑理論方法時會得到大量子狀態之間的路徑,從而使結果難以理解的困境,Qiu 等[113,114]利用VAE的數值分布變分優化,在隱空間實現了類似過渡路徑的合并.該方法被成功應用在兩個不同的簡單體系,分別是一對疏水粒子在水溶液中的聚集和Fip35WW 結構域折疊路徑的分析中.利用VAE能夠有效預測編碼空間、隱空間和解碼空間概率密度的特性,Monroe 和Shen[115]發展了基于隱空間的蒙特卡羅移動建議方法,再通過編碼和解碼,從而實現在真實高維空間有效且高接受率的移動.該方法的突出優點是直接滿足細致平衡要求,不需要一般偏置勢加速采樣生成軌跡后的權重調整,從而避免了與之伴隨的所有潛在問題和困難.這個思路和通過粗粒化模擬促進(細粒度)全原子模擬[76],以及把低維子空間視為信息瓶頸[74,75]的具體方式雖然差別較大,但總體基本思路一致.不同粒度之間更加高效準確的構型映射和信息傳遞還有很大的方法學發展空間,這方面的新發展也大概率會顯著促進復雜分子體系高精度多尺度模型的構建.

4 其他神經網絡方法在自由能地貌圖相關研究中的應用

神經網絡網絡的萬能逼近能力使得其在自由能面探索中從多個角度被加以應用.其中很多工作都致力于獲得更好的集合變量以改善復雜體系的增強采樣.早在2005 年,Ma 和Dinner[116]就開始使用神經網絡用來尋找復雜體系的反應坐標.針對各種傳統降維方法不能直接把結果中的低維空間(集合)變量表達為原空間坐標的問題,Chen 和Ferguson[117]利用自編碼器可以實現從高維輸入空間到低維隱空間之間的可訓練映射,把通過已有軌跡數據訓練生成的隱空間自由度作為集合變量,從而實現了對集合變量偏置勢通過自編碼器對高維空間坐標的直接微分計算偏向受力,集成了集合變量的神經網絡構建和在加速采樣中的直接應用,該方法在丙氨酸二肽和TrpCage 蛋白體系中被成功使用.與此類似,Chen 等[118]也采用自編碼器進行降維訓練獲得CV,然后通過自動微分把施加于CV 上的偏置勢傳遞到分子體系中去實現模擬采樣和自由能計算.

與使用變分優化特征函數不同的另外一種思路是回歸方法.Wehmeyer 和Noé[119]嘗試了選擇對N個連續時間坐標序列(Xt,Xt+τ,t=1,2,···,N)最小化回歸誤差[120–122]:

其中D和E分別為編碼器和解碼器.在對已有軌跡數據的時間序列坐標構型按照下式進行均值歸零((28)式和(29)式)和白化((30)式和(31)式):

然后對處理后的坐標優化訓練,實現編碼器降維和解碼器對原空間的映射:

通過訓練過程中在輸出端使用相對輸入端t時刻的延后t+τ時刻坐標,也實現了演化的預測.對于在構象空間中線性可分的不同亞穩態,該方法被證明同Koopman 模型[49,68]等價.但對非線性可分的體系,與PCA 和TICA 及人工構造特征空間相比,文獻[119]的丙氨酸二肽體系顯示通過編碼器和解碼器的深度學習擬合則可以更好地處理.

Zhang 和Chen[123]針對不恰當的CV 會在其正交空間出現亞穩態簡并(degeneracy)從而導致對應方向不能加速采樣的問題,發展了利用隨機動力學嵌 入(stochastic kinetic embedding,StKE)的半監督學習方法增加對當前信息最匱乏區域(current least informative regions,CLIRs)的主動學習采樣(active enhanced sampling,AES),這與Kleiman 和Shukla[68]在VAMPNets 輸出構象類型采樣最少的部分增加后續采樣的思路類似.該方法成功在丙氨酸二肽和五肽met-enkephalin體系中從隨意給定的無效CV 開始,以較短時間實現了對自由能地貌圖的可靠采樣.Rydzewski和Valsson[124]提出的多尺度重配權重隨機嵌入(multiscale reweighted stochastic embedding,MRSE)則在此基礎上更進一步,通過高斯混合模型描述高維特征空間和重配權重,實現對平衡態和偏置勢采樣數據在訓練中的有效使用.該方法被Rydzewski和Valsson 應用到-Brown Potential以及丙氨酸二肽和四肽體系,也已被整合到開源的PLUMMD 軟件包(https://www.plumed-nest.org/eggs/21/023/).類似地,Belkacemi 等[125]發展了利用自編碼器的自由能偏置勢迭代學習 (free energy biasing and iterative learning with auto encoders,FEBILAE),該方法可以對在平衡態或者偏置勢下采樣的軌跡重配權重后作為自編碼器的輸入(既可以是原來構象空間的,也可以是某種轉換之后的構型).其中自編碼器的瓶頸層確定了CV 的維度,但顯然需要自行選擇,他們也給出了探索的建議.可能的問題是迭代收斂的CV 并不能保證自由能地貌圖的全局充分采樣.和大多數類似研究一樣,這類編碼過程不具備直接可解釋性,人們無從知道輸入構型中不同的參數對CV 的貢獻.雖然原則上可以間接從計算過程中的自動微分步驟獲取一定信息,但目前所有的方法中沒有提供這種分析.針對這個問題,Kikutsuji 等[126]利用模型無關的局部解釋(local interpretable model agnostic explanation,LIME)和沙普利加和解釋(shapley additive explanations,SHAP)框架,給出各個輸入量對RC 的貢獻,能夠在一定程度上增進我們對體系的直觀物理認知.

Sun 等[127]發展了由一個降維編碼器,構象分類器和勢能預測器組成的多任務CV 學習構架,在幾個簡單測試系統(包括5DBrown model、丙氨酸二肽和金(110)晶面重建單元反應體系)與單目標訓練優化相比較展示了一定優勢.與很多應用中系統演化過程在原有高維空間進行不同,隱空間模擬器[128](latent space simulator,LSS)在訓練產生編碼器和解碼器后,在CV 空間快速展開系統演化,然后通過解碼器生成原有高維空間的細粒度軌跡.這些在隱空間或者集合變量空間進行操作的思路是很多工作中利用自編碼器的重要方式.比大多數方法在諸如丙氨酸二肽或類似模型體系中展示更進一步的是該方法在兩個較大體系(264 殘基的PROTAC 蛋白和DNA 序列5′-GCGGTTTCCGC-3′ 對應的雙螺旋結構)獲得了較為成功的應用.

Jung 等[129]以水溶液中離子的聚集和聚合物折疊為例,集成了深度學習和過渡路徑理論實現了復雜分子體系自組織模型的構建、驗證和更新,并在此基礎上通過符號回歸總結出更容易理解的可觀測量連接,是分子復雜體系的深度學習和可解釋性方面有意義的嘗試和進步.比變分求解自由能上界更進一步,Zhao和Wang[130]用流匹配 (flow matching)同時求解上下界,從而提供更好地逼近目標體系自由能的可能途徑.

鑒于生成式模型在語言圖像繪畫等方面的巨大成功[131],Janson 等[132]基于生成對抗模型和transformer 架構訓練的構象系綜生成神經網絡成功產生了訓練數據集中沒有的內秉無序蛋白(IDP)構象,該過程與分子模擬直接采樣相比所用計算代價非常小,不過正確性依然有待進一步在更多體系中驗證.

5 結論

綜上所述,變分方法處理分子體系自由能地貌圖目前已經有了較多不同視角的嘗試,但還都限于在較為簡單的體系探索,和其他理論上不甚嚴格的刻畫分子體系自由能地貌圖的神經網絡方法相比較也還沒有展示出明顯的系統優勢.比如使用變分的VAMPNets[67]和使用回歸[119]兩種方法在丙氨酸二肽體系中就沒有明顯的表現差異.不過變分更嚴格的理論基礎有可能會讓誤差控制更加容易,也很可能會在將來較大分子體系的應用和進一步發展中體現出更多的優勢.從理論方法的角度,現有的這些不同變分目標函數都是為了更好地逼近分子體系自由能地貌圖的準確描述,如何將它們集成并能夠依據應用需求靈活選擇關注視角顯然是個有價值的任務.當前的變分和自編碼器模型中還有很多需要人工調節和嘗試的環節,最為突出的就是目前的所有方法都不能通過自主學習優化獲得自編碼器中間低維隱空間的適當維度.另外變分計算本身原則上也可以在神經網絡中數值實現,從而有可能增加靈活性和可泛化能力,不過目前尚沒有見到這類嘗試,有可能是個有價值的發展方向.

從應用的角度,目前最迫切需要解決的問題可能是將這些變分構建向更大更復雜分子體系的延伸.從自由能地貌圖構象空間的層次來看,超過兩個時空間尺度的體系顯然會帶來更多挑戰,在同一個自由能地貌圖時空間尺度層次上,多個亞穩態之間過渡路徑交匯的可能性和準確處理也有待解決.這些問題的可靠處理在較大的復合體分子機器的理解中很有必要.

目前大模型的應用如火如荼[133],不過在AI 的科學應用領域尚沒有發力.主要原因之一是作為通用大模型訓練素材的語音圖像材料非常豐富,而特定科學領域的數據一般都不夠豐富或者很多都難以理解.不過這些模型集成多模態的能力顯然對AI 在廣泛科學應用中和特定的復雜分子體系中都有參考價值.已有的這些變分構建方法,還有未來可能出現的其他新穎構建,很可能在將來被統一到一個多目標大模型中.