引言
全细胞建模与仿真，是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.