Python反应-相互作用-扩散模拟器PyRID：高效模拟复杂生物分子系统的布朗动力学新工具

本研究针对生物分子系统的高复杂性、高密度及多尺度特性，开发了基于Python的布朗动力学模拟框架PyRID。该工具整合了多种计算算法，支持多分子反应、膜相关过程、多分散体系及刚性珠模型模拟，并引入分层网格数据结构提升多分散系统计算效率。验证表明PyRID能准确再现物理特性（如MSD、扩散张量），其性能优于现有工具（如ReaDDy），且支持网格分区与表面扩散，为跨膜蛋白与胞内外蛋白相互作用研究提供新途径。该工具兼具可扩展性与易集成性，显著提升了生物分子动态过程的模拟能力。

　　

在当今分子生物学研究领域，科学家们越来越关注细胞内复杂分子系统的动态行为。这些系统通常由大小不一、种类繁多的分子和蛋白质组成，以极高的密度存在于细胞中。细胞膜和细胞器结构限制了分子运动的空间，而跨膜蛋白（如离子通道）不仅需要在膜表面移动，还要与细胞内外的蛋白质发生相互作用。传统的热力学方法往往假设系统是孤立的，但生物细胞各个部分持续与周围环境交换离子、蛋白质等物质，这使得计算机模拟面临巨大挑战。

尽管已有一些优秀的模拟工具（如MCell、Smoldyn和ReaDDy）被开发出来，但它们各自存在局限性。MCell擅长在复杂三维细胞环境中模拟反应-扩散过程，但缺乏内置的力相互作用支持；Smoldyn能模拟分子扩散和化学反应，但只能近似球形分子的排除体积效应；ReaDDy虽然能模拟粒子间的力相互作用，却不支持三维网格来模拟几何复杂的区室和表面。这些工具在应对生物系统的特殊需求时显得力不从心，特别是当需要模拟多分散系统、刚性分子结构或复杂细胞几何形状时。

为了解决这些问题，Moritz Becker、Nahid Safari和Christian Tetzlaff开发了PyRID（Python反应相互作用扩散模拟器），这是一个完全用Python编写的布朗动力学模拟器，专门为高效模拟生物分子系统而设计。这项研究发表在《Cell Reports Methods》上，为研究人员提供了一个强大而灵活的工具来探索分子定位、运动和相互作用的复杂动力学。

研究人员采用了多种关键技术方法开展本研究。他们使用刚性珠模型（rigid bead models）来表示分子结构，每个分子可由一个或多个刚性连接的珠粒组成，支持各向异性的平移和旋转扩散。通过分层网格（hierarchical grids）数据结构高效处理多分散系统，利用三角网格（triangulated meshes）定义几何区室并约束分子运动。采用Numba即时（JIT）编译优化性能，结合紧密打包动态数组（tightly packed dynamic arrays）管理分子数据。表面射线行进（surface ray marching）算法处理分子在网格表面的扩散，固定浓度边界条件（fixed concentration boundary conditions）模拟与分子池的交换，并基于Doi反应方案处理双分子反应。

各向异性扩散

研究人员首先验证了PyRID在平移和旋转扩散算法上的准确性。他们使用Ilie等人的例子，比较了均方位移（MSD）和旋转时间相关性与理论预测的匹配程度。结果显示，刚性珠分子的MSD在各个维度（图1A中的彩色标记）与理论预测（黑线）高度一致。旋转弛豫（图1B）也与理论（灰线）紧密匹配，表明PyRID能准确模拟分子的各向异性扩散行为。

igG3蛋白的扩散张量

研究团队将计算扩散张量的方法直接集成到PyRID中，以便快速设置刚性珠分子系统。他们使用Hydro++文档中提供的igG3蛋白模型（图1C）进行测试，结果显示PyRID的计算结果与Hydro++高度一致（达到4位小数精度），证实了PyRID在计算复杂分子扩散特性方面的可靠性。

选择正确的反应速率和半径

在分子模拟中，准确选择反应参数（如反应速率k和反应半径R_react）对捕捉真实系统行为至关重要。研究人员发现，当反应速率较小时，系统有足够时间在反应间达到平衡，基于常微分方程（ODE）的方法与基于粒子的随机模拟算法结果匹配（图1E）。但随着反应速率增加，两种方法出现偏差，特别是在扩散限制的双分子反应中，ODE方法无法准确描述非均匀分布的体系。

刚性珠分子间的双分子反应

PyRID允许通过珠粒/斑块（beads/patches）表示复杂分子的不同反应位点。研究人员设计了一个示例：分子A和B分别由两个珠粒a₁、a₂和b₁、b₂组成，它们之间可以发生多种反应路径（图2A和2B）。当不考虑分子间排斥相互作用时，模拟结果与ODE描述高度一致（图2C）；但当引入排除体积效应后，粒子模拟与ODE结果不再匹配（图2D），突显了考虑分子间相互作用的重要性。

硬球流体

为了验证PyRID的准确性，研究人员模拟了一个硬球流体系统。他们使用谐波排斥势而非阶跃函数来近似硬球相互作用，结果显示径向分布函数（图3B蓝线）与Trokhymchuk等人的理论预测（红线）高度一致。压力计算结果（图3C）也与理论预期相符，且在使用Berendsen恒压器时，系统能够达到目标压力（图3D）。

性能基准测试

在多分散性测试中，PyRID使用分层网格方法处理多分散系统，性能表现优异。与经典链接细胞列表算法相比，分层网格将粒子更新速度从约80,000 pu/s提升到超过500,000 pu/s（图4A），使多分散系统的性能与单分散系统相当。

在分子相互作用和反应基准测试中，PyRID的表现优于ReaDDy（图4B）。对于粒子数在1,000到10,000之间的系统，每个粒子更新的计算时间保持在大约1.25微秒，相当于约800,000 pu/s。即使粒子数量超过10,000，性能也只有轻微下降。

与其他模拟器（ReaDDy、MCell和Smoldyn）的性能评估显示，所有模拟器产生一致的反应动力学（图4C），验证了它们的准确性。在内存使用方面，PyRID因支持更广泛的功能而具有较高的基线内存占用，但随着系统规模增大，内存需求增长较慢（图4D）。在运行时间方面，当包括Numba JIT编译时间时，MCell和Smoldyn表现最快；排除编译时间后，PyRID在所有模拟器中达到最佳运行时性能（图4E）。

网格区室和表面分子

PyRID使用三角网格表示区室（图5A），采用基于重心坐标的高效算法来检查点是否位于三角形内或位移向量是否与三角形边相交。通过表面射线行进方法（图5C），分子在网格表面上遵循测地线路径（图5D和5E）。扩散表面分子的MSD与理论一致（在二维中MSD=4Dt，其中D=43nm²/μs）（图5F）。

模拟边界

当区室与模拟盒相交时，相交的三角形被分配到透明类别（黄色），相应的边界相交边（紫色线）也是如此（图5G）。如果边界条件设置为"固定浓度"，透明三角形和边充当吸收边界，但同时将新分子释放到模拟体积中。用户还可以将网格三角形组合成面组，例如在网格表面的特定区域分布分子（图5G左侧蓝色区域）。对于周期性边界条件，PyRID遵循最小镜像约定（图5H）。

多分散球的泊松盘采样

PyRID支持多分散粒子的泊松盘采样，包括三种不同大小的粒子类型限制在网格区室体积内的分布（图5I）、表面分子的泊松盘采样（图5J）以及限制在由三角形面组定义的表面区域的表面分子采样（图5K）。

研究结论表明，PyRID成功整合了现有工具的多种特性，包括用于简化粒子模拟的刚性珠模型、通过实验和理论扩散张量实现准确扩散运动、珠模型表面斑块支持多价蛋白质-蛋白质相互作用等。该模拟器通过分层网格数据结构考虑粒子大小的多分散性，支持三角网格几何形状以考虑分区化的3D环境，固定浓度边界条件允许模拟较大系统中的子区域，并允许实施大量单分子和双分子反应以及表面扩散。

PyRID的研究意义在于它提供了一个强大而灵活的工具，能够模拟复杂生物分子系统的动态行为。与现有工具相比，PyRID在多个方面表现出优势：它支持多分散系统和刚性分子结构，能够高效处理几何复杂的区室和表面，性能与现有工具相当甚至更优，且完全用Python编写使得它更容易被广大科研社区接受和修改。

然而，研究也指出了PyRID目前的一些局限性。它尚未纳入分子间的流体动力相互作用，因为每时间步重新计算6N×6N扩散张量的计算成本会使大规模模拟不切实际。与ReaDDy不同，PyRID目前仅支持键约束的 pairwise 相互作用势，尚不包括粒子三元组的角约束或粒子四元组的扭转势。此外，PyRID目前无法模拟几秒或几分钟长时间尺度的过程，这对于研究细胞信号传导、网格蛋白自组装和蛋白质运输等过程是一个重要限制。

尽管存在这些局限性，PyRID仍然是粒子基反应扩散模拟的强大工具，特别适合需要模拟分子结构、相互作用和复杂几何形状的研究场景。它的开发为生物分子动态过程的研究提供了新的可能性，有望在计算生物学、生物物理学和药物设计等领域发挥重要作用。

研究人员已经将完整的PyRID代码连同广泛的在线文档和示例发布在GitHub存储库上，这将会促进更广泛的研究社区采用和发展这一工具。随着未来功能的进一步扩展和优化，PyRID有望成为研究生物分子系统动态行为的重要计算平台。