Multicellular simulations with shape and volume constraints using optimal transport

技术摘要：基于最优传输的具有形状与体积约束的多细胞模拟

问题陈述
生命系统，从细胞内组分到宏观组织及细菌菌落，均表现出由体积排斥、形状相互作用及可变形性驱动的复杂动力学。理解这些约束如何导致宏观自组织，是发育生物学和生物物理学的一项根本挑战。现有的计算模型在物理准确性、数值成本和数学可处理性之间面临权衡。常见方法包括：

• 点粒子系统： 数学上简单，但往往缺乏真实的形状表征。
• 相场和水平集模型： 能够表征任意形状，但计算成本高昂，且接触相互作用的处理具有挑战性。
• 元胞自动机和晶格模型： 将空间离散化为体素，可能限制几何灵活性。
• 沃罗诺伊（Voronoi）和顶点模型： 对类组织结构高效，但通常依赖体积的软约束，或缺乏对特定形状变形的控制。

在统一框架内对具有严格体积约束、任意动态形状和高粒子数量（数十至数万个）的系统进行建模，存在显著空白，该框架需同时支持个体可变形性和集体相互作用。

方法论
作者引入了一个基于最优传输理论的介观尺度模拟框架，具体利用拉盖尔（Laguerre）剖分。

1. 静态模型（拉盖尔剖分）：
   每个粒子 $i$ 占据的空间被建模为拉盖尔胞 $L_i$，定义为满足以下条件的点 $x$ 的集合：
   $$L_i = \{x \in \Omega, c(x, x_i) - w_i \leq c(x, x_j) - w_j \text{ for all } j\}$$
   其中，$x_i$ 是粒子位置，$v_i$ 是其目标体积，$c(x, x_i)$ 是编码粒子形状和可变形性属性的代价函数。权重 $w_i$（坎托罗维奇势）被唯一确定以满足严格体积约束 $|L_i| = v_i$。

   • 该公式将沃罗诺伊图（其中 $w_i=0$ 且 $c$ 为欧几里得距离）推广，通过选择任意代价函数（例如幂代价、各向异性代价）来允许弯曲边界和非球形形状。
   • 空空间被建模为一个具有零代价的附加粒子。

2. 动态模型：
   粒子位置 $x_i$ 的演化遵循一阶梯度下降方程：
   $$\dot{x}_i = -\tau_i \nabla_{x_i} T_c + F_i^{\text{point}} + F_i^{\text{surf}} + \dot{\xi}_i$$

   • 不可压缩力： 项 $-\nabla_{x_i} T_c$ 作为源自总传输成本 $T_c$ 的排斥力。它确保粒子不重叠，并自然地强制执行体积约束。与成对力不同，该力取决于所有接触表面的全局几何结构。
   • 外部/相互作用力（$F^{\text{point}}$）： 包括自推进、重力或吸引力，类似于点粒子模型。
   • 表面力（$F^{\text{surf}}$）： 沿界面模拟压力和表面张力效应，允许出现细胞分选和粘附等现象。
   • 形状动力学： 形状可通过参数化代价（例如由协方差矩阵定义的椭球体）或响应外部场（例如趋化性）的水平集势函数来演化。

3. 数值实现：

   • 该问题被求解为离散网格（体素）上的半离散最优传输问题。
   • 作者采用分而治之的多尺度策略和拟牛顿求解器（L-BFGS-B） 来高效计算坎托罗维奇势 $w_i$。
   • 实现利用 KEOPS 和 GEOMLOSS 库在 GPU 上进行大规模并行化，实现了在三维中对多达 50,000 个粒子的模拟，具有最先进的性能。

主要贡献

• 统一框架： 该方法在最优传输的伞盖下统一了沃罗诺伊剖分、水平集方法和顶点模型的概念。它在单一框架内支持任意形状、严格体积约束和复杂的相互作用机制。
• 严格体积执行： 与许多使用软约束或松弛力的现有模型不同，该方法通过优化坎托罗维奇势，在数学上保证了任意时刻的体积约束。
• 可扩展性： GPU 加速的求解器允许在三维中模拟大规模系统（多达 50,000 个粒子），这对于许多传统的有限元或相场方法来说在计算上是不可行的。
• 建模的多样性： 该框架成功模拟了多种现象，包括：
  • 由可变形性驱动的活动布朗粒子中的流体 - 固相变。
  • 仅由体积排斥引起的群游杆状粒子中的涌现取向有序。
  • 由形状变形驱动的趋化迁移。
  • 三维异质聚集体中的细胞分选和吞噬现象，重现了与微分界面张力假说（DITH）一致的结果。

结果
作者通过多个基准测试验证了该框架：

• 生长与分裂： 一个三维聚集体通过分裂从单个细胞生长到 50,000 个细胞，在约一天的计算时间内填满一个域。
• 相变： 在活动布朗粒子系统中，可变形性参数 $\alpha$ 驱动从流体相（扩散行为）到固体相（阻塞、六方堆积）的转变。该模型捕捉到了固体相中有效扩散的消失和均方位移的饱和。
• 杆状排列： 自推进杆状粒子的模拟表明，由于碰撞诱导的变形，长程取向有序自发涌现，这与聚合物物理理论一致。
• 三维细胞分选： 该模型重现了具有不同表面张力和软度的两种细胞类型的三维混合物中的各种分选模式（吞噬、分离、棋盘格），与理论预测和实验观察相符。
• 人群动力学： 该模型模拟了人群疏散，表明可变形粒子（流体状）能安全疏散，而硬球粒子（固体状）则形成稳定的拱门并在出口处阻塞。

意义与主张
本文声称引入了一种建模多细胞系统的新方法，弥合了详细物理描述与计算效率之间的差距。

• 验证： 主要贡献是将这一新框架与经典计算模型进行验证，证明其能够重现既定结果（例如细胞分选、相变），同时提供对体积约束和三维复杂性的优越处理能力。
• 生物学相关性： 作者强调了该方法在 in silico（计算机模拟）复杂发育过程（如囊胚发育）中的潜力，其中个体细胞变形和致密聚集体的行为都至关重要。
• 理论扩展： 这项工作为现有方法（点粒子、水平集、顶点）提供了一种新的最优传输视角，并为推导具有体积相互作用的系统的粗粒化连续模型（平均场近似）提供了一条途径，这是数学生物学中一个 largely 未解决的问题。
• 可及性： 作者提供了一个免费可用的 Python 实现，使该方法易于用于计算物理学和生物学领域的进一步研究。

作者保持谦逊，指出虽然该框架是新的，但当前的重点是对现有模型进行验证。他们确定了未来的方向，例如将活动顶点理论扩展到拉盖尔剖分，完善连续极限的推导，以及改进对复杂动态代价和表面相互作用的处理以实现实时模拟。