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

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

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

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

1 引言

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

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

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

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

圖1 每年結合生物分子模擬與機器學習的文獻數目隨年份的變化,數據來源于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.

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

2 基于機器學習算法的生物分子力場構建

2.1 勢能面與分子力場擬合

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

圖2 神經網絡用于生物分子構象能量面及力場的擬合Fig.2.Schematic diagram for representing the biomolecular force field by a neural network.

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

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

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

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

2.2 粗粒化力場構建

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

“自下而上”策略的基本思路是基于高精度力場模型的計算結果來確定粗粒化力場參數,主要方法有玻爾茲曼反演法 (Boltzmann inversion method)[65]、力匹配法 (force matching)[66]、漲落匹配法(fluctuating matching)[67]以及能量分解法 (energy decomposition)[68,69]等.例如,玻爾茲曼反演法主要通過全原子分子動力學模擬得到的徑向分布函數 (radial distribution function) 來提取粗粒化層次的有效相互作用參數;而力匹配法的優化目標則是使粗粒化粒子的受力與其在高精度力場中對應粒子的受力盡可能一致.需要注意的是,由于粗粒化近似,粗粒化粒子所代表的原子體系的自由度被凍結,粗粒化力場需要包含所凍結自由度構象熵對能量面的貢獻,因此是一種平均力勢 (potential of mean force).

以上基于“自下而上”方案構建粗粒化力場的策略與前述基于DFT計算結果擬合全原子力場的思路相類似,都希望基于低精度模型擬合更高精度的數據 (能量或力),從而在提升計算效率的同時,盡可能保留足夠的精確度.不同的是,從量子力學模型到全原子分子力場模型,由于原子自由度數目維持不變,因此分子力場不涉及構象熵的貢獻,原子尺度力場的擬合可以直接使用能量或力作為目標;而在構建粗粒化分子力場模型時,需要在一定程度上體現被凍結自由度的熵效應,因此對分子構象的采樣具有更高的要求,將力作為目標擬合力場參數是更常用的方法.另外,構建全原子力場模型的相關算法和構架,如神經網絡架構、體現平移與旋轉對稱性的結構特征提取方法、激活函數的選擇等,可以自然地遷移到基于力匹配的粗粒化力場擬合.近年來,基于深度學習構建粗粒化分子模型的工作越來越多地見諸于發表的論文中[34,52,70-74].例如,DeePMD團隊同時開發出與DeePMD具有相似網絡架構與結構特征提取策略的深度學習粗粒化力場方案——DeePCG[52].其中力場參數的提取使用了力匹配法和逐級擬合的辦法.同樣是基于前饋神經網絡架構和力匹配方法,Wang等[70]在2019年開發了CGNet,并展示了用于丙氨酸二肽與多肽鏈的粗粒化模擬結果,能夠很好地重現作為參考的全原子模擬得到的自由能面及其他統計性質.

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

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

圖3 基于粗粒化結構的蛋白殘基溶劑可及性表面積(SASA)計算.左圖: 蛋白分子(protein G,PDB code:1pgb)的全原子結構圖與粗粒化結構圖;右圖: 使用DeepCGSA由粗粒化結構計算得到的SASA與參考值的對比.其中參考值使用Shrake-Rupley算法由全原子結構計算得到[77].DeepCGSA能夠基于粗粒化結構給出接近參考值的SASA計算結果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 基于機器學習算法的分子模擬增強采樣與數據處理

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

3.1 基于機器學習算法提取反應坐標

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

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

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

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

圖4 用PCA (左)、t-SNE (中)和UMAP(右)對蛋白分子Protein G的基于粗粒化分子動力學的模擬軌跡[99] 降維效果對比.藍色到紅色對應表征蛋白折疊程度的Q值;Q=1 (紅色)為完全折疊結構,Q=0 (藍色)為完全解折疊結構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.

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

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

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

3.2.1 非生成模型

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

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

3.2.2 生成模型

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

圖5 不同生成模型的網絡架構.從左至右分別對應變分自編碼器、生成對抗網絡與標準化流.即便目標同為生成符合某種分布的數據,三種網絡使用了不同的架構與方法.變分自編碼器將數據降維至低維空間后,在低維空間采樣并再次變換至高維空間;生成對抗網絡則通過生成器與分類器之間的互相對抗而使生成器生成的結果符合目標分布;標準化流則是在目標分布與簡單易采樣的分布 (如高斯分布) 之間建立直接且可逆的映射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).

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

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

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

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

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

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

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

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

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

4 總結與展望

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

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

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