生物分子模擬中的機器學習方法*

管星悅 黃恒焱 彭華祺 劉彥航 李文飛? 王煒?

1) (南京大學物理學院,南京 210093)

2) (國科溫州研究院,溫州生物物理重點實驗室,溫州 325000)

1 引言

以分子動力學為代表的分子模擬技術(shù)在生物大分子結(jié)構(gòu)與動力學研究中發(fā)揮著越來越重要的作用.常規(guī)分子模擬技術(shù)用于復雜生物分子體系時,不可避免地存在力場精度與構(gòu)象采樣效率瓶頸.同時,從高維分子模擬數(shù)據(jù)提取可解釋的生物大分子結(jié)構(gòu)與動力學特征也是一個挑戰(zhàn)性難題.生物分子模擬技術(shù)發(fā)展的核心任務(wù)便是解決以上難題,擴展生物分子模擬的應(yīng)用范圍.

自從20世紀70年代McCammon等[1]首次將分子動力學模擬用于生物大分子體系以來,人們在生物分子力場發(fā)展、長程靜電相互作用計算方法、增強采樣與自由能計算等方面取得了多個突破[2].分子模擬技術(shù)與高性能計算機等硬件技術(shù)的協(xié)同發(fā)展使得分子模擬能夠覆蓋的時間尺度以超過摩爾定律的速度增加,平均每10年增加約3個數(shù)量級[3].這些進展使得人們能夠直接模擬小蛋白分子毫秒時間尺度的折疊全過程[4,5],也能對固有無序蛋白(intrisically disordered protein,IDP)的構(gòu)象系綜進行合理的分子模擬表征[6,7],甚至能夠?qū)崿F(xiàn)對病毒顆粒、細胞質(zhì)等超大分子體系進行分子模擬[8,9].目前,實驗和模擬計算結(jié)合已成為生物大分子結(jié)構(gòu)與動力學研究的基本范式.另一方面,對較大的分子體系,目前的生物分子模擬能夠達到的空間和時間尺度與實驗測量仍有一定距離,從而限制了其適用范圍[10].因此,發(fā)展新的分子模擬技術(shù),擴展分子模擬技術(shù)的適用范圍,對基于生物分子模擬的基礎(chǔ)和應(yīng)用研究至關(guān)重要.

隨著計算能力的提升和海量數(shù)據(jù)的積累,機器學習算法被廣泛應(yīng)用于基礎(chǔ)與應(yīng)用科學的各個領(lǐng)域.自然地,人們也將機器學習算法應(yīng)用于計算生物學與生物信息學研究,如生物分子設(shè)計與結(jié)構(gòu)預測、分子模擬以及分子對接等.機器學習概念誕生于20世紀50年代[11],并在曲折的發(fā)展中被多次重新理解與表述.早期的機器學習算法多是對既有建模與優(yōu)化方法的重新整理與表述,如線性回歸、多項式回歸[12]以及k-近鄰算法[13]等.盡管在早期歷史中已初具雛形,目前人們廣泛使用的機器學習算法,如決策樹[14]、神經(jīng)網(wǎng)絡(luò)[15]、支持向量機[16]以及集成學習方法[17,18]等,大多成型于1980年后,并很快被應(yīng)用于蛋白質(zhì)二級結(jié)構(gòu)預測[19]、蛋白結(jié)構(gòu)與功能分類[20,21]以及藥物篩選[22]等問題.在20世紀90年代,人們也開始將神經(jīng)網(wǎng)絡(luò)用于構(gòu)建簡單分子體系 (如表面吸附氣體分子) 的勢能面并進行分子模擬[23].在這些早期的應(yīng)用中,機器學習方法往往被視為可替代的工具,且神經(jīng)網(wǎng)絡(luò)尚未表現(xiàn)出相對其他機器學習算法的顯著優(yōu)勢,因此相關(guān)算法在生物分子模擬領(lǐng)域的應(yīng)用仍非常有限.

近年來,以深度學習為代表的機器學習技術(shù)得到迅猛發(fā)展,并在多個領(lǐng)域展現(xiàn)出驚人的能力.特別是AlexNet[24]的誕生,展示了深度卷積神經(jīng)網(wǎng)絡(luò)對圖像的強大識別能力,宣布深度學習革命的到來.之后出現(xiàn)的殘差網(wǎng)絡(luò) (ResNet)[25]進一步推動了神經(jīng)網(wǎng)絡(luò)向深度發(fā)展,也出現(xiàn)了如生成對抗網(wǎng)絡(luò)(GAN)[26]與Transformer[27]等網(wǎng)絡(luò)架構(gòu)新范式.這些新的機器學習算法開始廣泛用于生物分子模擬、結(jié)構(gòu)預測與設(shè)計等領(lǐng)域.自2017年開始,機器學習與生物分子模擬相結(jié)合的研究工作大幅增加,成為勢不可擋的學科交叉趨勢.這一趨勢從近年來發(fā)表的相關(guān)研究論文數(shù)目的增長中可見一斑 (圖1).

圖1 每年結(jié)合生物分子模擬與機器學習的文獻數(shù)目隨年份的變化,數(shù)據(jù)來源于ScopusFig.1.Number of publications with the key words “molecular simulations” and “machine learning” published per year as a function of years.Data were taken from Scopus.

機器學習與生物分子模擬的結(jié)合為推進分子生物物理學研究提供了新的機會.例如,利用機器學習技術(shù)能夠設(shè)計更準確的分子力場,開發(fā)更高效靈活的增強采樣算法,發(fā)展更具普適性的復雜生物分子體系的結(jié)構(gòu)與動力學預測算法,并輔助藥物分子的設(shè)計.這一重要的交叉領(lǐng)域正在高速發(fā)展并持續(xù)產(chǎn)生具有突破性進展的研究成果[28-35].因此對該領(lǐng)域的發(fā)展進行回顧與綜述尤為重要.關(guān)于機器學習在生物大分子結(jié)構(gòu)預測與設(shè)計方面的進展,已有非常全面的綜述可供參考[36-40],本文不再過多討論.在機器學習與生物分子模擬交叉領(lǐng)域,也有學者從不同角度進行了綜述[41-44].例如,Ramanathand等[42]在其綜述論文中介紹了使用機器學習技術(shù)表征IDP系綜以及進行多尺度模擬的方法,并提出將模擬數(shù)據(jù)集與實驗擬合的重要性及策略;Noé等[43]詳細介紹了機器學習算法在幫助解決生物分子模擬重要挑戰(zhàn)中發(fā)揮的作用,并探討了將物理學原理融入機器學習算法的必要性及相關(guān)方法;Wang等[44]詳細總結(jié)了利用機器學習算法分析分子動力學模擬軌跡的方法,以及利用機器學習與相關(guān)數(shù)據(jù)驅(qū)動方法進行增強采樣的方案.本文將在此基礎(chǔ)上,結(jié)合該領(lǐng)域的最新進展,從生物分子力場構(gòu)建、反應(yīng)坐標的選取與增強采樣、分子模擬數(shù)據(jù)處理等方面對機器學習與分子模擬交叉領(lǐng)域的代表性工作進行綜述.生物物理智識與機器學習技術(shù)迭代的融合已成為人們探索生命原理的有力手段,而結(jié)合機器學習算法的生物分子模擬是借助神經(jīng)網(wǎng)絡(luò)的強大表達性與擬合能力分析復雜生命運動密碼的重要實踐.期望本文對該領(lǐng)域的綜述有助于讀者綜合了解機器學習算法在生物分子模擬中的重要應(yīng)用,共同思考和探索基于機器學習算法解決生物分子模擬領(lǐng)域關(guān)鍵難題的可能途徑.

2 基于機器學習算法的生物分子力場構(gòu)建

2.1 勢能面與分子力場擬合

在生物分子模擬中,精度和效率通常難以兼得.不同的問題在精度和效率上有不同的偏重與要求,因此需要針對性地選擇能夠平衡精度與效率要求的折中方案.計算化學領(lǐng)域的“金標準”CCSD(T)方法能達到約1 kcal/mol的化學精度,但代價是計算效率低,通常適用于小體系的單點能計算.基于密度泛函理論 (DFT) 和Born-Oppenheimer絕熱近似的方法在精度上作出妥協(xié),從而提升了計算效率,能夠?qū)⒂嬎泱w系大小提升到數(shù)百個原子以上的規(guī)模.但是,對于絕大多數(shù)的生物大分子,計算體系通常包含上萬個原子,并涉及微秒以上的時間尺度,因此進一步提升生物分子模擬的計算效率對擴展其應(yīng)用范圍十分關(guān)鍵.分子力場模型通過參數(shù)化力場的方式在原子坐標水平近似地描述絕熱能量面,從而大幅提升計算模擬效率.這種逐級近似的框架之下,如何在提升計算模擬效率的同時盡可能減小精度的損失,成為構(gòu)建分子力場的核心問題.全原子水平的分子力場可以看作是基于原子坐標和原子類型的高維空間上的多元函數(shù).傳統(tǒng)分子力場多使用基于經(jīng)驗的結(jié)構(gòu)項和以單體、兩體勢表示的非鍵相互作用項的參數(shù)化方案[45-47].這種預先設(shè)定的具體力場函數(shù)形式不可避免地對力場精度帶來限制.盡管人們可以通過進一步引入極化和多體效應(yīng)等物理機制來提升參數(shù)化方案的表達能力[48,49],但在精度上與DFT方法仍有較大差距.深度學習算法提供了一種表達能力強大的參數(shù)化方案(圖2),可以降低對預設(shè)力場函數(shù)形式的依賴,因此原則上可以提升對分子力場的描述精度.需要注意的是,深度學習算法更強的參數(shù)化表達能力,需要由充足的計算能力和訓練數(shù)據(jù)來作為支撐.近年來,計算能力與數(shù)據(jù)規(guī)模已經(jīng)可以支持用于訓練具有足夠強表達能力的深度神經(jīng)網(wǎng)絡(luò),因此使用深度學習方法構(gòu)建生物分子力場,從而實現(xiàn)分子力場精度突破的條件已經(jīng)成熟,且在此問題上已取得重要進展[50-55].

圖2 神經(jīng)網(wǎng)絡(luò)用于生物分子構(gòu)象能量面及力場的擬合Fig.2.Schematic diagram for representing the biomolecular force field by a neural network.

機器學習算法用于生物分子力場擬合的一個典型例子是Zhang等[51,56]在2018年發(fā)表的Dee PMD工作.DeePMD使用原子尺度的構(gòu)象坐標以及量子力學精度的能量信息作為數(shù)據(jù)集,將系統(tǒng)構(gòu)象映射至其對應(yīng)的能量與力(受益于神經(jīng)網(wǎng)絡(luò)組件的求導能力).給定系統(tǒng)構(gòu)象坐標,可以通過網(wǎng)絡(luò)的前向傳播代替復雜的DFT計算,直接得到原子受力,從而在盡量保留DFT精度的前提下實現(xiàn)高效率分子動力學模擬.DeePMD的網(wǎng)絡(luò)架構(gòu)本身是深度前饋網(wǎng)絡(luò),由多個全連接網(wǎng)絡(luò)的輸出求和得到總能量.DeePMD使用分子構(gòu)型的相對坐標來保證網(wǎng)絡(luò)的輸出不依賴于生物分子體系的平移與旋轉(zhuǎn)變換.值得一提的是,DeePMD可以對接LAMMPS,Gromacs等傳統(tǒng)分子動力學模擬軟件,便于使用.

為了在神經(jīng)網(wǎng)絡(luò)訓練中保持分子構(gòu)型平移與旋轉(zhuǎn)對稱性,除使用相對坐標(或單個分子體系的內(nèi)坐標)外,另一類方法是使用Behler與Parrinello[57]在2007年提出的對稱函數(shù)方法.對稱函數(shù)方法將系統(tǒng)中每一個原子依次視為中心原子,計算其與附近原子的距離、夾角,得到對稱函數(shù)值,并作為神經(jīng)網(wǎng)絡(luò)的輸入特征量.例如,Artrith與Urban[58]發(fā)展的Aenet神經(jīng)網(wǎng)絡(luò)模型以及Smith等[59]發(fā)展的ANI-1神經(jīng)網(wǎng)絡(luò)模型均使用了該對稱函數(shù)方法,并成功用于體相TiO2等材料系統(tǒng)和有機物小分子系統(tǒng)的力場擬合.Fan等[60]在基于進化策略算法構(gòu)建用于原子模擬的機器學習勢時也采用了類似的方法.該對稱函數(shù)方法規(guī)避了笛卡爾坐標與內(nèi)坐標的相互轉(zhuǎn)換,從而提升深度網(wǎng)絡(luò)的參數(shù)表達能力和訓練效率.

以上DeepMD,Aenet,以及ANI-1均采用了深度前饋網(wǎng)絡(luò)構(gòu)架.隨著卷積神經(jīng)網(wǎng)絡(luò) (CNN) 展示出其對圖像特征提取與識別的強大能力并在機器學習領(lǐng)域帶來革命,人們也嘗試使用CNN處理圖像的范式來處理分子構(gòu)型并映射到能量面或力場.特別是殘差網(wǎng)絡(luò)構(gòu)架的引入,使得人們可以在避免過擬合的前提下,構(gòu)建足夠深度的CNN網(wǎng)絡(luò),以增強其擬合效果.一個代表性的例子是Schütt等[50]發(fā)展的SchNet.SchNet以殘差卷積網(wǎng)絡(luò)實現(xiàn)對分子構(gòu)型特征的提取.不同于處理圖像數(shù)據(jù)使用的網(wǎng)格狀離散濾波器,為了保證能量面的光滑性與精確性,SchNet采用了連續(xù)濾波器.相對于深度前饋網(wǎng)絡(luò),基于CNN架構(gòu)的SchNet能夠顯著提升在量子化學精度數(shù)據(jù)集QM9 (包含有機小分子的構(gòu)型、能量等)的預測精度,也在分子動力學數(shù)據(jù)集MD17[61]上有更好的表現(xiàn).

盡管CNN可以提取局域而抽象的特征,且相較于全連接神經(jīng)網(wǎng)絡(luò)在避免出現(xiàn)過擬合方面表現(xiàn)出色,但CNN最擅長的領(lǐng)域仍是處理規(guī)整的圖像等數(shù)據(jù).對于空間不規(guī)則且以共價鏈接為重要特征的分子構(gòu)型,圖(graph)是一種更為自然的表示.分子構(gòu)型的圖描述天然地擁有平移和旋轉(zhuǎn)不變性,并且允許將距離、化學鍵等連接信息作為“邊”數(shù)據(jù)存入圖網(wǎng)絡(luò).因為這些優(yōu)點,人們也嘗試使用圖神經(jīng)網(wǎng)絡(luò)來學習擬合分子力場.Park等[53]于2021年發(fā)表的GNNFF基于結(jié)合有向圖與消息傳遞 (message passing)的深度神經(jīng)網(wǎng)絡(luò)框架[62],構(gòu)建了神經(jīng)網(wǎng)絡(luò)分子力場模型,對有機小分子受力的預測精度超過SchNet.Wang等[63]在同年發(fā)表的sGNN,考慮了不同類型相互作用在空間尺度上的差異,對聚合物分子的主鏈共價作用和非鍵相互作用能量項分開建模,在空間尺度擴展性與對不同模擬體系的可遷移性方面表現(xiàn)良好.

2.2 粗?；鰳?gòu)建

相對于DFT等量子化學方法,基于分子力場的全原子分子動力學模型極大地擴展了計算模擬方法能夠研究的生物分子體系的空間和時間尺度.目前,人們已經(jīng)能夠?qū)崿F(xiàn)對較小蛋白體系的完整折疊過程進行全原子分子動力學模擬.另外,通過結(jié)合增強采樣算法,可以實現(xiàn)對較大生物分子體系構(gòu)象變化的全原子分子動力學模擬和自由能計算.然而,對于更大的生物分子系統(tǒng),如分子馬達、核糖體、病毒顆粒以及染色質(zhì)體系等,通常包含百萬以上原子個數(shù),并涉及毫秒以上時間尺度的動力學過程,遠超出全原子分子動力學模擬能夠達到的時間和空間尺度范圍.為了突破全原子分子動力學模擬的計算效率瓶頸,人們通常采用粗粒化的近似方法[64].在粗?；Ｐ椭?將多個原子映射為1個虛擬粒子,從而很大程度上降低了體系的自由度數(shù),實現(xiàn)分子模擬效率的提升.然而,由于采用了虛擬粒子近似,構(gòu)建具有合理精度的粗粒化分子力場是一個極具挑戰(zhàn)性的難題.已有的粗?；Ｐ偷牧鰠?shù)主要通過“自下而上”和“自上而下”兩種策略來優(yōu)化得到.

“自下而上”策略的基本思路是基于高精度力場模型的計算結(jié)果來確定粗粒化力場參數(shù),主要方法有玻爾茲曼反演法 (Boltzmann inversion method)[65]、力匹配法 (force matching)[66]、漲落匹配法(fluctuating matching)[67]以及能量分解法 (energy decomposition)[68,69]等.例如,玻爾茲曼反演法主要通過全原子分子動力學模擬得到的徑向分布函數(shù) (radial distribution function) 來提取粗?；瘜哟蔚挠行嗷プ饔脜?shù);而力匹配法的優(yōu)化目標則是使粗粒化粒子的受力與其在高精度力場中對應(yīng)粒子的受力盡可能一致.需要注意的是,由于粗?；?粗?；Ｗ铀淼脑芋w系的自由度被凍結(jié),粗?；鲂枰鶅鼋Y(jié)自由度構(gòu)象熵對能量面的貢獻,因此是一種平均力勢 (potential of mean force).

以上基于“自下而上”方案構(gòu)建粗?；龅牟呗耘c前述基于DFT計算結(jié)果擬合全原子力場的思路相類似,都希望基于低精度模型擬合更高精度的數(shù)據(jù) (能量或力),從而在提升計算效率的同時,盡可能保留足夠的精確度.不同的是,從量子力學模型到全原子分子力場模型,由于原子自由度數(shù)目維持不變,因此分子力場不涉及構(gòu)象熵的貢獻,原子尺度力場的擬合可以直接使用能量或力作為目標;而在構(gòu)建粗?；肿恿瞿Ｐ蜁r,需要在一定程度上體現(xiàn)被凍結(jié)自由度的熵效應(yīng),因此對分子構(gòu)象的采樣具有更高的要求,將力作為目標擬合力場參數(shù)是更常用的方法.另外,構(gòu)建全原子力場模型的相關(guān)算法和構(gòu)架,如神經(jīng)網(wǎng)絡(luò)架構(gòu)、體現(xiàn)平移與旋轉(zhuǎn)對稱性的結(jié)構(gòu)特征提取方法、激活函數(shù)的選擇等,可以自然地遷移到基于力匹配的粗?；鰯M合.近年來,基于深度學習構(gòu)建粗粒化分子模型的工作越來越多地見諸于發(fā)表的論文中[34,52,70-74].例如,DeePMD團隊同時開發(fā)出與DeePMD具有相似網(wǎng)絡(luò)架構(gòu)與結(jié)構(gòu)特征提取策略的深度學習粗?；龇桨浮狣eePCG[52].其中力場參數(shù)的提取使用了力匹配法和逐級擬合的辦法.同樣是基于前饋神經(jīng)網(wǎng)絡(luò)架構(gòu)和力匹配方法,Wang等[70]在2019年開發(fā)了CGNet,并展示了用于丙氨酸二肽與多肽鏈的粗?；M結(jié)果,能夠很好地重現(xiàn)作為參考的全原子模擬得到的自由能面及其他統(tǒng)計性質(zhì).

以上例子均采用了基于“自下而上”思路的力匹配法作為粗粒化力場擬合方案.與其相對應(yīng)的“自上而下”的思路追求粗?；瞿M結(jié)果與實驗約束或高精度模型得到的宏觀性質(zhì)的相容性.然而,因為每一步優(yōu)化都需要在當前參數(shù)下得到模擬軌跡并進行反向傳播,自上而下的方法通常會給訓練帶來較大的計算負擔,對擬合目標與參數(shù)優(yōu)化方案的選擇具有更高要求[75,76].近期Clementi和Noé等[34,71]提出了以flow-matching為例的一類新方法: 將標準化流 (normalizing flow,NF) 或去噪擴散模型 (denoising diffusion probabilistic model)等生成模型與力匹配法相結(jié)合,先利用高精度數(shù)據(jù)訓練粗粒化構(gòu)象的生成模型,再從這種生成模型中提取粗?；?這些新的方法將生成模型描述的粗粒化構(gòu)象偏好視作一種平衡采樣,從而與力場產(chǎn)生聯(lián)系.其他的生成模型,如變分自編碼器 (variational auto-encoder,VAE)[72]和使用對抗訓練思想的VADE[73]同樣可以被用于描述粗粒化坐標下的構(gòu)象分布.

另外,在基于Cα的蛋白質(zhì)粗?；Ｐ椭?由于側(cè)鏈原子位置信息的缺失,無法準確地體現(xiàn)蛋白質(zhì)分子的表面積、靜電勢分布等蛋白質(zhì)分子的基本性質(zhì).但是這些信息對理解蛋白質(zhì)分子的結(jié)構(gòu)組裝、構(gòu)象動力學以及分子識別等過程至關(guān)重要.因此,如何在粗粒化模型框架下準確地計算蛋白質(zhì)分子的表面積、靜電勢等蛋白質(zhì)分子的基本性質(zhì)是一個重要的技術(shù)挑戰(zhàn).基于深度神經(jīng)網(wǎng)絡(luò)的機器學習算法為解決這一問題提供了一個可行的方案.例如,本文作者在最近的工作中,構(gòu)建了一套深度學習網(wǎng)絡(luò)DeepCGSA,能夠基于粗粒化模型結(jié)構(gòu)高精度地估算蛋白質(zhì)、核酸等生物大分子的溶劑可及性表面積(圖3)[74].嘗試將類似的方法用于針對粗?；鞍踪|(zhì)結(jié)構(gòu)的靜電勢分布與pKa值的預測也取得了很好的效果.

圖3 基于粗?；Y(jié)構(gòu)的蛋白殘基溶劑可及性表面積(SASA)計算.左圖: 蛋白分子(protein G,PDB code:1pgb)的全原子結(jié)構(gòu)圖與粗粒化結(jié)構(gòu)圖;右圖: 使用DeepCGSA由粗?；Y(jié)構(gòu)計算得到的SASA與參考值的對比.其中參考值使用Shrake-Rupley算法由全原子結(jié)構(gòu)計算得到[77].DeepCGSA能夠基于粗?；Y(jié)構(gòu)給出接近參考值的SASA計算結(jié)果Fig.3.SASA estimation based on coarse-grained protein structure.Left: All-atom structure and coarse-grained structure of protein G (PDB code: 1 pgb).Right: Correlation plot between the SASA values from DeepCGSA based on one-bead coarse-grained structure and the reference values by Shrake-Rupley algorithm based on all-atom structure.The DeepCGSA can well reproduce the SASA values based on coarse-grained structure.

3 基于機器學習算法的分子模擬增強采樣與數(shù)據(jù)處理

由于生物大分子具有龐大的自由度數(shù)和復雜的能量面特征,全原子水平的分子模擬通常會遇到采樣困難.特別是在計算各種平衡統(tǒng)計性質(zhì)時,需要分子模擬的采樣盡可能遍歷重要的構(gòu)象空間,并在給定的系綜條件下達到平衡.盡管上述粗?；Ｐ吞峁┝艘环N解決采樣困難的有效方案,但粗粒化近似不可避免地導致計算精度的損失.特別是當特異性的氫鍵、鹽橋等原子層次的相互作用起到主導作用時,粗?；Ｐ屯ǔo法顯式地體現(xiàn)這類特異性相互作用特征,從而限制了其應(yīng)用范圍.因此,發(fā)展增強采樣算法是解決分子模擬采樣困難的另一有效方案.基于統(tǒng)計物理原理,人們已經(jīng)發(fā)展出多個有效的增強采樣算法,并廣泛應(yīng)用于生物大分子體系的蒙特卡羅模擬和分子動力學模擬[78-88].目前常見的增強采樣算法有傘形抽樣(umbrella sampling)[78]、副本交換分子動力學 (replica exchange molecular dynamics)[79]、元動力學 (metadynamics)[80]、加速分子動力學 (accelerated molecular dynamics)[81]以及溫度積分增強抽樣方法 (integrated tempering sampling,ITS)[82]等.這些增強采樣算法多已通過外部插件(如PLUMED[83]) 或直接整合到成熟的分子動力學模擬軟件.另外,人們也發(fā)展了適用于研究構(gòu)象轉(zhuǎn)變路徑的增強采樣算法,如String方法[84]與Transition path sampling方法[85]等.最近,人們將機器學習算法用于生物分子模擬的增強采樣,并取得了顯著效果,甚至還可以利用機器學習算法,基于有限的構(gòu)象采樣數(shù)據(jù)實現(xiàn)高維自由能面的構(gòu)建[89,90].

3.1 基于機器學習算法提取反應(yīng)坐標

常用的增強采樣算法可分為兩類: 依賴反應(yīng)坐標的增強采樣算法和不依賴反應(yīng)坐標的增強采樣算法.例如,傘形抽樣、元動力學等增強采樣算法依賴于預先定義的反應(yīng)坐標,這類算法的基本策略通常是沿預先定義的反應(yīng)坐標方向添加偏置勢,從而避免在沿反應(yīng)坐標的局部勢阱中重復采樣.因此,預先定義的反應(yīng)坐標需對應(yīng)所關(guān)注的生物分子體系最重要的運動方向,而垂直于反應(yīng)坐標方向的動力學具有更快的時間尺度.然而,定義合適的反應(yīng)坐標本身就是一項極具挑戰(zhàn)性的任務(wù).通常情況下,反應(yīng)坐標主要基于物理直覺來選取,而機器學習等數(shù)據(jù)驅(qū)動的降維方法為反應(yīng)坐標的選取給出了一個更為理性和可操作的方案.

常規(guī)的不使用神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)驅(qū)動降維方法主要基于如下思想設(shè)計: 在降維前后的空間里,盡可能維持數(shù)據(jù)的某種結(jié)構(gòu)信息不變.這種“結(jié)構(gòu)信息”可以分為全局信息和局域信息兩類.早在20世紀初就被開發(fā)的主元分析算法PCA,是一種典型的致力于維持全局結(jié)構(gòu)信息的算法[91].PCA將高維數(shù)據(jù)點相對于幾何中心的歐式距離平方和視作需要保留的“結(jié)構(gòu)信息”,在通過線性變化降維過程中最小化該結(jié)構(gòu)信息的損失,并找到承擔最大運動信息變化的反應(yīng)坐標.PCA方法的缺陷也在于此: 基于全局的歐式距離衡量信息并非總是一個合理的預設(shè);且PCA要求降維至超平面,就只允許對數(shù)據(jù)做全局的線性變換,很多時候這是一個過強的假設(shè).

更一般地,可以假設(shè)高維數(shù)據(jù)分布在一個黎曼流形(或是幾支黎曼流形)上.此時歐式距離只適用于描述數(shù)據(jù)點的局域結(jié)構(gòu),即可以構(gòu)建起離散數(shù)據(jù)點的近鄰圖,而全局結(jié)構(gòu)可視為由這些近鄰圖組合而成.基于這一思想,Isomap算法[92]和Diffussion Map算法[93]分別用測地線距離和模擬擴散距離衡量數(shù)據(jù)點的間距,并希望降維映射前后這些距離盡量保持不變,從而將流形“展平”以實現(xiàn)降維.將Isomap與Diffusion Map用于分子模擬數(shù)據(jù)分析,可以找到非線性地依賴于高維數(shù)據(jù)的反應(yīng)坐標[94-96].

在基于局域結(jié)構(gòu)信息的降維方法中,2008年提出的t-SNE算法具有突出的表現(xiàn)[97].t-SNE對數(shù)據(jù)點間的相似性做非線性變換,使得降維過程中主要維護局部團簇 (cluster) 中兩點相似性的分布不變,而對相似性低的數(shù)據(jù)點的位置關(guān)系幾乎沒有約束.因此,t-SNE的降維盡量維持了數(shù)據(jù)點基于相似性簇團的內(nèi)部結(jié)構(gòu),而對簇團間的距離朝向則幾乎沒有要求,從而帶來了降維結(jié)果的隨機性.t-SNE使用梯度下降優(yōu)化低維空間數(shù)據(jù)點的位置,通常這是一個非凸優(yōu)化,每次得到的結(jié)果會有所差別.相比于2002年提出的SNE算法[98],t-SNE構(gòu)建對稱的損失函數(shù)以代替SNE中不對稱的K-L散度,簡化了基于梯度的優(yōu)化過程;同時t-SNE以更為長尾的t-分布建立低維空間距離向概率的映射,以更好應(yīng)對高維數(shù)據(jù)點嵌入低維空間導致的擁擠問題.圖4給出了使用PCA,t-SNE以及UMAP對粗?；肿觿恿W得到的蛋白折疊軌跡[99]進行降維的效果對比: 相比于PCA,t-SNE和UMAP能更好地區(qū)分折疊態(tài)和解折疊態(tài)的結(jié)構(gòu).在分子模擬中,基于t-SNE的降維算法已被廣泛應(yīng)用于反應(yīng)坐標的定義與高維動力學軌跡的可視化[100-102].除t-SNE外,基于局域結(jié)構(gòu)信息的降維方法還有: 維持局域線性關(guān)系的LLE (locally linear embedding)[103]、維持局域鄰近圖的Laplacian Eigenmaps[104]、最小化局域曲率的Hessian LLE[105]等,然而它們在分子模擬領(lǐng)域得到的關(guān)注和應(yīng)用遠不如t-SNE.2018年McInnes等[106]提出的UMAP降維算法采用了與t-SNE類似的、基于鄰近圖提取簇團信息的策略,并同樣用梯度下降方法優(yōu)化得到低維嵌入.不同的是,相比于圍繞著“點”進行的t-SNE,UMAP采用了以“邊”為中心的優(yōu)化策略,使用交叉熵作為優(yōu)化目標,將邊存在的概率映射為低維空間的距離.在生物分子模擬中,UMAP常被用于基因組、染色質(zhì)和單細胞轉(zhuǎn)錄譜等數(shù)據(jù)[107,108].在單細胞轉(zhuǎn)錄譜數(shù)據(jù)集與蛋白質(zhì)動力學軌跡數(shù)據(jù)上的比較研究[109-111]均表明: UMAP具有不遜色于t-SNE的降維效果,但是在計算成本上遠低于t-SNE,對大規(guī)模的數(shù)據(jù)有良好的擴展性,這與UMAP原始論文中指出其計算復雜度約為N1.14一致[106].

圖4 用PCA (左)、t-SNE (中)和UMAP(右)對蛋白分子Protein G的基于粗?；肿觿恿W的模擬軌跡[99] 降維效果對比.藍色到紅色對應(yīng)表征蛋白折疊程度的Q值;Q=1 (紅色)為完全折疊結(jié)構(gòu),Q=0 (藍色)為完全解折疊結(jié)構(gòu)Fig.4.Projection of the sampled snapshots of the coarse-grained molecular dynamics simulations for protein G [99] along the reaction coordinates constructed by PCA (left),t-SNE (middle),and UMAP (right),respectively.t-SNE and UMAP perform better than PCA in distinguishing the folded and unfolded structures.Colors from blue to red represent the structures with increasing folding extent: blue,fully unfolded;red,fully folded.

如果認為降維算法的關(guān)鍵問題在于對信息的選擇與度量,那么以上非神經(jīng)網(wǎng)絡(luò)的機器學習降維算法都是通過引入某種預設(shè) (或主觀判斷) 來解決此問題,也因此降低了對降維變換的表達能力.借助于具有強大表達能力的神經(jīng)網(wǎng)絡(luò),可以期待構(gòu)建更有效的降維算法.

在2013年被開發(fā)的VAE,通過巧妙地設(shè)計神經(jīng)網(wǎng)絡(luò)架構(gòu),將原始數(shù)據(jù)通過編碼器降維得到隱變量,再通過解碼器升維,生成與原始數(shù)據(jù)同維度的高維數(shù)據(jù)[112].如果生成數(shù)據(jù)具有和原始數(shù)據(jù)幾乎相同的分布,則說明編碼過程 (即降維過程) 幾乎沒有造成信息損失,低維的隱變量具有與原始數(shù)據(jù)相近的表達能力.就訓練過程而言,VAE通過優(yōu)化編碼器和解碼器參數(shù),以最小化生成數(shù)據(jù)與原始數(shù)據(jù)分布上的差異.其中,隱變量的“信息”通過復現(xiàn)原始數(shù)據(jù)分布的能力衡量.相比于以上非神經(jīng)網(wǎng)絡(luò)的降維算法中預設(shè)信息為數(shù)據(jù)集上的某種結(jié)構(gòu)的做法,VAE衡量信息的方式更具一般性與整體性.對于生物分子模擬系統(tǒng),這一優(yōu)勢將有利于VAE通過降維找到整體性的反應(yīng)坐標;而編碼器、解碼器所基于的深度神經(jīng)網(wǎng)絡(luò)架構(gòu)保證了VAE強大的表達能力,降低了模型對預設(shè)信息的依賴,有利于增強降維的有效性.因此VAE常被用于生物分子模擬反應(yīng)坐標的提取.另外,VAE尋找反應(yīng)坐標的思路同樣可以用于粗粒化模型的建立[72]、反應(yīng)路徑搜索[113]、甚至是藥物分子設(shè)計[114]等任務(wù).

3.2 基于機器學習算法的增強采樣

3.2.1 非生成模型

機器學習算法不僅可以用于尋找合適的反應(yīng)坐標,還可以直接用于輔助分子模擬采樣.例如,在利用Metadynamics方法進行增強采樣和自由能計算時,需要在分子體系的固有能量面添加一定形狀的高斯形偏置勢[115],而確定高斯形偏置勢的參數(shù)及其變化規(guī)律非常關(guān)鍵,直接影響采樣效率.過強的高斯形偏置勢可能會導致采樣進入非物理的區(qū)域,而過弱的高斯形偏置勢又難以遍歷感興趣的構(gòu)象空間區(qū)域.2019年,Bonati等[116]通過結(jié)合神經(jīng)網(wǎng)絡(luò)與變分增強采樣思路,靈活地以變分形式在增強采樣模擬過程中自適應(yīng)地更新偏置勢,使得反應(yīng)坐標的實際分布能夠逼近目標分布.相較于常規(guī)的Metadynamics方法,在靈活性、高效性與準確性方面得到了提升.

另一個使用神經(jīng)網(wǎng)絡(luò)給出偏置勢用以增強采樣的例子是Zhang等[117]提出的TALOS (targeted adversarial learning optimized Ssampling).類似生成對抗網(wǎng)絡(luò)GAN的思想 (見下方關(guān)于生成模型的描述),TALOS使用Wasserstein距離衡量真實分子模擬引擎生成的構(gòu)型分布與目標構(gòu)型分布的差異,將此距離的計算轉(zhuǎn)化為對一個判別器網(wǎng)絡(luò)的優(yōu)化問題.TALOS的訓練同樣類似于GAN: 對每個偏置勢,通過優(yōu)化判別器網(wǎng)絡(luò)計算兩分布Wasserstein距離的近似值,以之作為兩分布差異的數(shù)值衡量;最小化此差異以優(yōu)化偏置勢,從而使偏置勢下模擬產(chǎn)生的構(gòu)型分布盡可能接近目標分布.

3.2.2 生成模型

在以上例子中,機器學習算法僅被用作提供構(gòu)造反應(yīng)坐標或設(shè)置偏置勢的手段,即增強采樣的輔助工具,并沒有直接采樣生成分子構(gòu)型.近年來生成模型的發(fā)展,使得人們可以借助神經(jīng)網(wǎng)絡(luò)直接生成生物分子構(gòu)型,這為發(fā)展新的增強采樣算法提供了新思路.目前廣泛使用的生成模型主要包括VAE[112],GAN[26]和標準化流模型[118]等(圖5).前面已經(jīng)提到VAE通過編碼器-解碼器對原始數(shù)據(jù)進行降維再升維后,使生成數(shù)據(jù)盡可能與原數(shù)據(jù)保持相似,這一相似性的要求體現(xiàn)在損失函數(shù)中的K-L散度.但實際上VAE的損失函數(shù)是兩部分的拮抗,另一部分則是希望低維空間的隱變量盡可能接近高斯分布,于是解碼器并不是編碼器的逆變換,而是在擬合低維空間高斯分布參數(shù)后進行的采樣.生成對抗網(wǎng)絡(luò)GAN則是用一個分類器來判定生成器生成的結(jié)構(gòu)是否合理 (被生成的分子構(gòu)型與原數(shù)據(jù)是否相似這一判據(jù)由分類器訓練得到),那么由生成器生成并被分類器所選擇的分子構(gòu)型便可以作為采樣結(jié)果.標準化流模型則是用一系列變換在兩種分布間搭建可逆變換.

圖5 不同生成模型的網(wǎng)絡(luò)架構(gòu).從左至右分別對應(yīng)變分自編碼器、生成對抗網(wǎng)絡(luò)與標準化流.即便目標同為生成符合某種分布的數(shù)據(jù),三種網(wǎng)絡(luò)使用了不同的架構(gòu)與方法.變分自編碼器將數(shù)據(jù)降維至低維空間后,在低維空間采樣并再次變換至高維空間;生成對抗網(wǎng)絡(luò)則通過生成器與分類器之間的互相對抗而使生成器生成的結(jié)果符合目標分布;標準化流則是在目標分布與簡單易采樣的分布 (如高斯分布) 之間建立直接且可逆的映射Fig.5.Network architecture of different generative models:Variational autoencoder (VAE,left),generative adversarial network (GAN,middle),and normalizing flow (NF,right).Three networks have different architectures.VAE first reduces data to a low-dimensional space,samples in the lowdimensional space,and then transforms back to a high-dimensional space.GAN generates target distribution by combining a generator and the discriminator.Normalizing flow model establishes a direct and reversible mapping between the target distribution and a simple and easy-tosample distribution (such as Gaussian distribution).

用神經(jīng)網(wǎng)絡(luò)本身作為采樣核心的一個代表性例子是Noé等[28]2019年發(fā)表于Science的Boltzmann Generator.該工作完全展示了標準化流模型作為數(shù)學優(yōu)美的可逆生成模型在兩種分布之間建立聯(lián)系的能力.真實的分子構(gòu)型的分布符合玻爾茲曼分布,通常這一分布難以直接采樣.Boltzmann Generator的標準化流模型通過構(gòu)建可逆的坐標變換,將簡單且容易采樣的高斯分布映射到玻爾茲曼分布,從而實現(xiàn)滿足玻爾茲曼分布的采樣.

假設(shè)變量x取自高斯分布,概率密度為q(x),z代表生物大分子結(jié)構(gòu)的原子坐標,概率密度為r(z),應(yīng)當符合玻爾茲曼分布.需要注意的是,即使通過訓練建立起變量之間的映射函數(shù)z=f(x),變量z的概率密度并不是簡單的r(z)=r(f(x)),而是需要考慮體積元變化:

在保持變換可逆的同時計算出雅可比因子是一個難點.而標準化流模型使用了如下精巧設(shè)計: 用一系列小變換一步步將x變換到z,并且每一個小變換都是可逆且雅可比因子易計算的.每一步變換將x劃分為兩部分x=(x1,x2),:

該變換的可逆性也顯而易見.如果變換使得雅可比矩陣呈三角矩陣,則可非常容易地計算得到雅可比因子:

在Boltzmann Generator的實際訓練中,可以選擇多種訓練方式: 1)將實際軌跡作為玻爾茲曼分布端的輸入z,并通過訓練將經(jīng)變換輸出的x優(yōu)化為高斯分布;2)將根據(jù)高斯分布采樣得到的x變換得到構(gòu)型z,并根據(jù)分子力場計算其能量,通過訓練使z的分布優(yōu)化到玻爾茲曼分布;3)結(jié)合前兩種方法訓練;4)在前三種選項的基礎(chǔ)上加入額外的依賴于反應(yīng)坐標的損失項,使得變換后的分子構(gòu)型盡可能在反應(yīng)坐標空間均勻分布 (類似Metadynamics的思想),此后再進行reweighting操作,得到正確的分布.在將Boltzmann Generator用于BPTI蛋白的構(gòu)象采樣時,成功得到了其“X”態(tài)到開放的“O”態(tài)之間的構(gòu)象轉(zhuǎn)變,即使這種轉(zhuǎn)變的過渡態(tài)并不存在于訓練集中,展現(xiàn)了Boltzmann Generator在用于生物分子構(gòu)象采樣的強大能力.

另外一類常見的增強采樣算法采用了強化學習方法.強化學習使用獎懲機制,在不同的環(huán)境條件下強化學習器采取不同的動作時將給出一定的獎勵或懲罰,而訓練的目標是使得強化學習的動作能夠?qū)ⅹ剟钭畲蠡?Shamsi等[119]提出了基于強化學習的REAP算法,將獎懲機制與分子構(gòu)象空間的探索綁定在一起,尋找最利于在構(gòu)象空間擴散的反應(yīng)坐標.該方法用于丙氨酸二肽和Src激酶體系時展示了出色的增強采樣的效果.基于類似的思想,人們也可以基于強化學習,在沿所設(shè)定反應(yīng)坐標采樣的不確定度(uncertainty)上施加獎懲來鼓勵體系在未遍歷的構(gòu)象區(qū)域采樣(在已經(jīng)遍歷的方向上反應(yīng)坐標的不確定度較低),因此能對增強采樣模擬施加一個自適應(yīng)的靈活偏置勢[120],達到增強采樣的效果.

綜上,機器學習在增強采樣領(lǐng)域表現(xiàn)出強大的功能和應(yīng)用前景,既可以在傳統(tǒng)增強采樣算法框架下通過構(gòu)建反應(yīng)坐標發(fā)揮作用,也可以通過自適應(yīng)的方式提供更高效靈活的偏置勢,還可以直接利用生成模型作為采樣核心.隨著新的機器學習算法的開發(fā),將機器學習用于輔助生物分子模擬增強采樣是未來生物分子模擬領(lǐng)域的重要課題.

3.3 基于機器學習算法的生物分子模擬數(shù)據(jù)處理

生物分子模擬通常在高維空間中進行,所得到的分子模擬軌跡包含了豐富的結(jié)構(gòu)與動力學信息,如何從這些高維的分子模擬軌跡提取出可解釋的熱力學與動力學數(shù)據(jù),并實現(xiàn)與實驗結(jié)果的定量比較是分子模擬領(lǐng)域的另一個挑戰(zhàn)性難題.分子動力學模擬數(shù)據(jù)處理主要包括以下幾個方面: 高維分子模擬數(shù)據(jù)特征提取、分子模擬軌跡降維與反應(yīng)坐標構(gòu)建、分子模擬微觀狀態(tài)粗粒化與馬爾可夫狀態(tài)模型構(gòu)建,以及低維自由能面構(gòu)建等.顯然,適合于處理復雜數(shù)據(jù)的各類機器學習算法在生物分子模擬數(shù)據(jù)處理中扮演著越來越重要的角色.事實上,前述關(guān)于增強采樣算法部分介紹的基于機器學習的反應(yīng)坐標構(gòu)建是機器學習用于分子模擬數(shù)據(jù)處理的重要方面.除此之外,人們也發(fā)展了深度網(wǎng)絡(luò)模型,用于提取生物分子體系的動力學與自由能信息.例如,Mardt等[121]設(shè)計了VAMPNet,能夠端到端地直接實現(xiàn)從分子模擬數(shù)據(jù)軌跡得到馬爾科夫狀態(tài)模型(Markov state model)的映射.以馬爾可夫過程的變分法(VAMP)為基礎(chǔ),深度網(wǎng)絡(luò)用于表達特征變換的形式,通過變換后的特征空間內(nèi)近似得到弛豫時間τ范圍內(nèi)的狀態(tài)轉(zhuǎn)移矩陣,從而用于提取生物分子的動力學信息.另外,Schneider等[90]通過訓練神經(jīng)網(wǎng)絡(luò),實現(xiàn)了高維自由能的計算以及典型系綜平均性質(zhì)的計算.

4 總結(jié)與展望

本文對機器學習方法在生物分子模擬領(lǐng)域的應(yīng)用進行了綜述.借助其突出的特征提取和參數(shù)擬合能力,機器學習方法(特別是神經(jīng)網(wǎng)絡(luò)算法)在全原子/粗粒化分子力場構(gòu)建、分子模擬數(shù)據(jù)降維與反應(yīng)坐標提取、以及生物分子構(gòu)象采樣等方面已經(jīng)開始發(fā)揮重要作用.隨著以深度神經(jīng)網(wǎng)絡(luò)為代表的機器學習算法的迭代更新,結(jié)合機器學習算法的生物分子模擬技術(shù)將成為人們在分子層次探索生命原理的重要研究范式.需要指出的是,目前機器學習算法大多作為輔助工具在生物分子模擬中發(fā)揮作用.即使整合了機器學習算法,對較大的生物分子體系能夠達到的分子模擬時間尺度仍與真實生物學相關(guān)時間尺度有較大差距.完全解決生物分子模擬精度與效率瓶頸,實現(xiàn)生物分子模擬與實驗測量的定量比較,需要在分子模擬的理論框架與算法方面同時進行探索.近年來,整合全原子模型和粗?；Ｐ偷亩喑叨壬锓肿幽M技術(shù)越來越受到人們的重視[67,122-124],是解決生物分子模擬精度與效率瓶頸的一個值得重點嘗試的思路.神經(jīng)網(wǎng)絡(luò)等機器學習算法的發(fā)展將成為進一步推動多尺度分子模擬技術(shù)發(fā)展的新突破口.

盡管本文將機器學習用于生物分子模擬的工作分為力場構(gòu)建、增強采樣以及數(shù)據(jù)處理等不同的主題來進行綜述,近年來突破性的工作通常打破了主題分類的邊界,并依賴于多個步驟的集成耦合.因此,實現(xiàn)機器學習在生物分子模擬多方面的融合應(yīng)用,需要開發(fā)能夠集成機器學習算法與生物分子模擬的軟件平臺.例如機器學習與生物物理交叉領(lǐng)域代表工作——AlphaFold2與ESM大語言模型,均得益于對多模態(tài)數(shù)據(jù)與算法的集成整合能力[30,125].國內(nèi)在相關(guān)領(lǐng)域的集成軟件平臺開發(fā)方面也取得了很大進展.由深勢科技開發(fā)的RiDYMO平臺集成了神經(jīng)網(wǎng)絡(luò)、分子動力學引擎、增強采樣方法,不僅能進行分子動力學模擬,分析蛋白質(zhì)構(gòu)象空間、還能探索藥物結(jié)合位點并計算藥效相關(guān)動力學參數(shù),適合藥物的設(shè)計與開發(fā)工作[126].北京大學與華為等團隊開發(fā)的MindSPONGE[127]在華為昇思MindSpore框架下整合了多種分子模擬、結(jié)構(gòu)預測設(shè)計以及全面的神經(jīng)網(wǎng)絡(luò)支持.這些集成平臺將降低新算法的開發(fā)和使用門檻,促進生物分子模擬技術(shù)的應(yīng)用范圍.

關(guān)于機器學習與生物分子模擬融合應(yīng)用的研究進展給我們帶來一個重要的啟示: 生物物理智識與機器學習發(fā)展是相輔相成的.例如,AlphaFold2的架構(gòu)借鑒了由序列比對得到的共進化信息,而AlphaFold2的成功又是機器學習推進生物分子結(jié)構(gòu)預測領(lǐng)域的代表例子.生物分子模擬與神經(jīng)網(wǎng)絡(luò)結(jié)合的需求也同樣在推進機器學習領(lǐng)域的發(fā)展.SchNet為了擬合光滑連續(xù)力場而在卷積神經(jīng)網(wǎng)絡(luò)架構(gòu)下提出的連續(xù)濾波器,可以被推廣到其他機器學習任務(wù)情景;而主要受分子結(jié)構(gòu)拓撲相關(guān)研究驅(qū)動而發(fā)展的圖神經(jīng)網(wǎng)絡(luò),也被推廣到諸如社交網(wǎng)絡(luò)等應(yīng)用情景中.機器學習架構(gòu)的每一次突破性進展都會為生物分子研究領(lǐng)域帶來難以估量的靈感與啟發(fā).如何借助神經(jīng)網(wǎng)絡(luò)的成功進一步反哺生物物理智識與經(jīng)驗將是未來生物物理與人工智能交叉領(lǐng)域的重點研究課題.