JCIM｜利用模塊化細胞生物學原理進行高性能全細胞模擬

中大唯信

北京中大唯信科技有限公司

引言

全細胞建模與仿真，是21世紀的重大挑戰之一，更是系統生物學的終極目標。利用現有的實驗確定的信息，進行的詳細的全細胞模型及其模擬，可用來探索未知的、未觀測到的生物系統區域，從而進一步擴展了人類現有的生物知識的極限。

在組織建模、細胞建模、神經生物學等方面，數學建模，即生物過程的數學表示，已被證明是濕實驗室實驗的一個非常成功的替代方法。而設計和模擬一個廣泛的生物全細胞模型，是一個非常耗時的過程。當前，盡管存在一些基於隨機模擬的方法，如E-Cell、虛擬細胞、GEPASI和原始細胞的布朗動力學模擬等，但這些僅限於小的假設模型。但是，只有像Markus Schwehm在2001年預測的那樣，將問題並行化並利用現有的高性能計算(HPC)系統，大量模型和仿真數據的極端調節，才能進行全細胞仿真。因此，必須以這樣的方式構建、設計和處理全細胞模型，以便在合理的時間內，合理地利用高性能計算系統，來執行並行的全細胞模擬。

細胞功能，是由稱為功能模塊的不同相互作用的分子群來執行的。有時，多個功能模塊，共同參與完成某一細胞功能。組成每個功能模塊的相互作用分子，被分配到特定的細胞區域或隔間，它們在其中發揮功能，而它們從指定區域，穿越到細胞內其他區域的概率非常低。這一觀察幫助人們得出結論，如果能夠最小化功能模塊之間的相互依賴性，那麼每個功能模塊都可以被獨立地模擬。因此，全細胞建模，可以看作是每個單元，只包含一個功能模塊的亞細胞建模的總和。

由於全細胞模擬的計算時間，取決於許多因素，包括分子的數量、細胞的大小和模擬的持續時間，因此，正確地對整個細胞建模是很重要的，這樣才能有效地利用現有的高性能計算架構。Markus Schwehm預計，在196個CPU的開普勒簇上，模擬大腸杆菌的全細胞周期大約需要24天。在大腸杆菌的細胞周期中，在4000萬個細胞質分子中，發生了1016個生化反應。植物和動物細胞的成分，大約是大腸杆菌的1000倍。因此，與大腸杆菌細胞相比，這些細胞更複雜，它們的模擬計算也更昂貴。

全細胞模擬，需要處理大量的模型和模擬數據。因此，在小的理論模型中能夠很好執行的方法，在實際的大的全細胞模型中變得難以管理，除非利用現有的HPC系統，同時將問題並行化處理。因此，最大限度地利用高性能計算系統，是使全細胞模擬可行的絕對前提。這篇文章中，研究者提出了一種基於隨機模擬的方法，通過在合理的時間內有效地利用現有的高性能計算系統，可以模擬整個細胞的大量分子。

此文中，研究者首先描述了用於模擬的整個細胞的數學模型。隨後，研究者介紹了並行實現的計算方法和細節。最後，研究者提出了，優化全細胞模擬的方法。(由於篇幅原因，我們這裏將詳細介紹後面兩個部分，第一部分可詳見原文)

計算細節

在此，研究者選擇了單細胞細菌大腸杆菌，展示了他們的模型。與含有3748個蛋白質的大腸杆菌(K12)的蛋白質位點圖(PLG)相同，研究者設計的細胞，由3748個蛋白質分子組成。研究者主要目標，是模擬由這3748個分子單獨組成的所有功能模塊。盡管大腸杆菌細胞由大約188個蛋白質模塊組成，但研究者演示了不同數量模塊的結果，以便對模擬器進行深入的性能分析。為了可視化和易於理解，生成了一個進程監控日志(PML)文件，該文件允許用戶在使用分子可視化軟件PyMOL進行模擬之前，可視化細胞模型的整個設置。每個子單元類似於一個虛擬容器，它模擬分配給它的功能模塊。研究者提出的方法，分別執行每個虛擬子單元，邊界條件因子單元之間的不同而不同，這取決於它是駐留在單元中邊緣還是中心。在研究者的CUDA實現中，在內核模塊中設置了一些標記，可以檢測分子從一個亞細胞到另一個亞細胞的遍曆。因此，每當需要在核之間傳遞消息時，模擬數據的當前狀態，就從GPU轉移到CPU。然後CPU相應地更新數據結構，並將其傳輸回GPU進行進一步處理，主要用於下一個模擬時間瞬間。

研究者對這裏討論的所有模擬，都使用了固定的參數集。模擬的總步驟為1000步，分子之間的碰撞被認為是非彈性的，即COR(回彈系數) = 0。研究者對全細胞預聚類進行CPU模擬，串行實現了所有3748個分子。然後，進行並行仿真，研究者生成了空間定位的PLG簇。

結果討論

2.1 負載均衡

研究者將集群劃分為不同數量的核，如圖1所示，這樣每個核的工作負載都是相似的。當使用兩個GPU核時，研究者將其中一個核分配到，最大的包含1693個蛋白質的簇中，而剩下的核分配到另外三個包含1481個蛋白質的簇中。類似地，當使用四個GPU核時，研究者將兩個核分配到最大的簇中，每個核分別處理847和846個蛋白質的計算。由771個蛋白質組成的簇被分配到第三個核，而第四個核處理另外兩個共包含710個蛋白質的簇。按照相同的步驟，研究者在8、16、32、64、128和256個GPU核之間平衡工作負載。同時，研究者從最大的四個簇的3174個粒子中，移除294個隨機粒子，並將剩下的2880個粒子，分配給所有GPU的2880個CUDA核。

圖1. 繪制模擬3748個大腸杆菌分子所需的計算和通信時間圖

圖片來源於JCIM

2.2 觀察

從圖1可以看出，利用GPU的兩核模擬PLG集群系統所需的計算時間，幾乎是單核CPU所需的3.9倍。當開始使用兩個或更多的GPU內核時，內核之間的通信，會根據分子的運動而產生。因此，巨大的內存傳輸(CPU到GPU，反之亦然)開銷和核間通信，是導致兩個GPU核相對於單個CPU核計算時間，要求如此之高的兩個主要因素。對於4個GPU核，計算時間與CPU仿真時間相近。四個GPU核改進後的數據處理和計算速度隱藏了內存延遲，因此相應的時間要求低於兩個GPU核。當使用8個或8個以上的GPU CUDA核時，由於在GPU上的數據處理和計算速度上開始大幅提高，獲得了良好的性能。隨著計算單元數量的增加，分配給核的分子數量也越來越少。這導致了屬於不同核的分子之間更多的相互作用，從而導致越來越多的核間通信。雖然2核、4核和8核的通信時間較短，但通信時間從16核開始增加，但在128核左右趨於穩定。因此，對於當前的小區配置和所選的仿真參數，在128核GPU上進行仿真時得到了最優的性能。為了便於理解，研究者在每個時間步結束時，收集模擬數據或每個分子的軌跡，並使用一種化學文件格式存儲它，稱為XYZ文件格式，擴展名為a.xyz，它存儲了分子的笛卡爾坐標，可以很容易地在PyMOL的視頻中顯示出來。研究者對老鼠和人類這兩種高等生物，進行了可擴展性分析。褐家鼠(大鼠)的PLG由9554個蛋白、652738個蛋白-蛋白相互作用(protein-protein interactions, PPIs)和598個緊密連接的簇組成。同樣，對於智人(人類)，PLG由41550個蛋白，8943744個PPIs和711個簇組成。

圖2和圖3分別為模擬大鼠9554和人41550分子，所需的計算和通信時間。對於大鼠來說，計算和通信時間都穩定在256個GPU核左右，而對於人類來說，則穩定在1024個GPU核左右。正如預期的那樣，核需求隨著蛋白質數量、它們之間的相互作用以及PLG的增加，從低等生物到高等生物而增加。

圖2. 模擬褐家鼠(大鼠)9554個分子所需的計算和通訊時間圖

圖片來源於JCIM

計算優點

亞細胞內動力學所需的計算時間，總是少於亞細胞間動力學(圖1)，因為在亞細胞內動力學中，分子的運動受到模擬子空間的限制。因此，GPU核之間不需要進行信息傳遞。當PLG簇之間的連接數量非常少(可以忽略)時，就可以完全避免這樣的串擾，只限制它們各自亞細胞內分子的運動。研究者的模擬器也可以處理這種場景，這減少了通信時間要求，並提供了更高效的計算模擬。

此外，研究者的方法非常靈活，可以很容易地用於對高等生物進行並行的整體模擬。從上述例子中可以明顯看出，對於任何生物，研究者都可以計算其PLG、空間定位簇，並利用研究者提出的計算框架並行執行其細胞動力學。

對於開放系統，可能需要動態的PLG，這取決於進入/退出系統的蛋白分子濃度。研究者所提出的模擬器，可以使用基於哈希的元胞字典數據結構，來處理這種情況。首先構建整個細胞的PLG，然後將其分割成緊密相連的簇。最初，關於每個集群及其相應組件的信息，存儲在元胞字典(CD)數據結構中。在模擬過程中，每當一個蛋白質分子從系統中退出時，它對應的細節就會從CD中刪除。或者，在進入時，分子的細節在相應的簇/亞細胞的CD中更新。因此，PLG的動態特性取決於，由CD處理的進入/離開系統的蛋白分子濃度。

在這裏進行的並行全細胞動力學，和圖1中提到的計算時間分析，證明了對PLG進行聚類，並利用現有HPC系統的單獨計算單元並行模擬這些聚類。事實上，研究者所提出的方案，適用於任何生物的任何類型的生物網絡，並且整個方法保持不變。並行全細胞仿真所面臨的巨大計算開銷的一個恒定因素始，終是CPU到GPU的內存傳輸開銷，當將仿真框架從CPU計算轉移到並行GPU計算時，不能忽視這個因素。在並行單元模擬過程中，研究者可以控制的另一個昂貴的計算因素，是核心間通信的數量。因此，應該盡可能減少集群間連接的數量，因為：即使盡最大努力微調集群過程，由於蜂窩網絡的高度互聯性，仍然不能完全消除集群間連接的數量。然而，在一個特定的時間步長，沒有任何集群間的移動，不需要從GPU到CPU的內存傳輸，反之亦然，這是在研究者的模擬器中使用的。

圖3. 模擬41550個智人(人類)分子所需的計算和通信時間圖

圖片來源於JCIM

範例研究

圖4中展示了用研究者的模型在16個亞細胞中模擬全細胞蛋白動力學的範例。研究者用PyMOL生成了仿真軌跡並將其可視化。在開始模擬之前，研究者將整個細胞分為16個不同大小的亞細胞。每個亞細胞由不同數量的蛋白質分子組成，它們用不同的顏色標記，如圖4所示。在HPC系統的單個計算單元上，分別對每個子單元進行了模擬。最初，所有的蛋白質分子都處於靜止狀態，它們被分配到三維亞細胞空間的隨機笛卡爾坐標中。當模擬開始時，力作用在分子上，分子開始移動。隨著模擬的進行，系統定期更新。這裏不允許亞細胞間運動，因此，蛋白質只能在它們分配的細胞間隔內運動。分子傾向於非彈性碰撞，因此，模擬中的蛋白質分子也會非彈性碰撞並粘在一起(如圖4中黃色圓圈所示)，其速度也隨之改變。隨著模擬的進行，碰撞的分子可能會再次分解成單獨的分子，這取決於作用在碰撞物體的團簇上的力。

圖4. 從細胞初始狀態開始的90次迭代後，描述平行全細胞蛋白動力學的示例模擬

圖片來源於JCIM

展望與結論

在此，研究者提出了一種有效地利用現有並行硬件架構，並利用整個細胞的動態特性的並行全細胞模擬框架，從而引導人們走向細胞動力學的新概念。研究者選擇了力場，運動方程和一個積分算法，並提出一個算法，來檢測和解決碰撞。然後，提出如何利用現有的高性能計算系統，來進行最優的並行全細胞模擬。

研究者觀察到，對於大腸杆菌，當至少使用了128個核的GPU時，得到了一個最佳的模擬時間，並且計算和通信時間都變得穩定。對於大鼠和人類，分別用256和1024個GPU核，實現了計算和通信的穩定性。此外，可以通過盡可能減少核間通信和計算時間，來達到最優的全細胞模擬時間。以空間定位的生物網絡的形式，聚集所有可能的細胞信息，通過最小化集群間連接的數量，聚集它們緊密相連的子組件，在單獨的GPU核上模擬每個集群，通過CD數據結構有效地處理核間通信，為利用高性能超級計算架構，進行並行全細胞建模和仿真，提供了新的研究視角。

在本文中，模擬是在沒有溶劑的情況下進行的。將BD(布朗動力學)應用到仿真工具中，產生了充當溶劑作用的隨機力。研究者展望，未來該方法應該進化到，支持任何種類的溶劑以及各種參數，如溶劑的粘度和溶液的溫度等。這樣，就可以利用一個合適的力場，進行全細胞模擬和特征值分析，以確定時間步長與細胞動力學之間的關系。通過所有這些循序漸進的實現，未來將能夠進行全細胞模擬，准確地模擬真實的活細胞，從而進一步拓展，現有生物學知識的極限。研究者相信，此處提出的計算框架，只要有足夠的實驗數據，對於任何生物的生物網絡都是有效的，並且可以擴展到任何CPU-GPU架構。

參考文獻

Barnali Das and Pralay Mitra. High-Performance Whole-Cell Simulation Exploiting Modular Cell Biology Principles. J. Chem. Inf. Model. 2021, ASAP. DOI: 10.1021/acs.jcim.0c01282.

相關文章

Angew Chem | 計算模擬設計大環PROTAC分子

JCIM | 從分子動力學軌跡中提取特征預測P-糖蛋白底物

Science Advances | 經典！從分子動力學到虛擬篩選發現PPI小分子抑制劑

ACS Catalysis | 基於分子動力學模擬的脫氫/還原酶的理性設計

Chemical Science | 加速動力學結合馬爾可夫模型指導蛋白設計

JMC | 反應動力學結合分子動力學研究共價抑制劑與DAGL的結合機制