基于最优传输理论的多细胞形状与体积约束模拟

《SCIENCE ADVANCES》：Multicellular simulations with shape and volume constraints using optimal transport

　　

许多生命与物理系统，如细胞聚集体、组织或细菌群落，表现为受体积排除（volume exclusion）与形状相互作用强烈约束的非传统粒子系统。理解这些约束如何导致宏观自组织结构是发育生物学等领域的基本问题。研究人员引入一种框架来模拟具有任意体积、动态形状及可变形特性的粒子系统。该方法基于最优传输（Optimal Transport, OT）理论及其在不可压缩流体流动、人群动力学与材料科学中的近期应用。该方法支持广泛的相互作用与变形机制，同时以最先进的数值性能自动处理体积排除约束。研究人员通过一系列实验展示该方法的通用性，证明其如何扩展与精炼先前方法的结果，特别关注生物物理学中具有挑战性的三维（3D）情形。

研究背景方面，从细胞内尺度到宏观种群尺度，形状、变形与拥堵效应在生命系统动力学中起核心作用，尤其在发育生物学中理解测量这些效应是重要且具挑战性的实验问题。过去几十年除体外与体内实验研究外，数学模型被用于研究这些问题，计算机实验可在最小廉价设置中评估推测的生物物理定律影响，其中形状与每个主体的物理属性可被单独控制。这带来建模与实现挑战，因无唯一简单数学方式表示动态形状与多细胞系统，需在物理精度、数值成本与数学分析潜力间找到权衡。现有模型包括点粒子系统（point-particle systems）、相场（Phase fields）模型、细胞自动机与格点模型（Cellular automata and lattice models）、Voronoi镶嵌（Voronoi tessellation）、顶点模型（Vertex models）等，但各有局限，如点粒子缺乏真实性，相场与顶点模型数值成本高或处理拓扑变化复杂，Voronoi镶嵌通常不考虑体积约束或控制较弱。因此研究人员开展基于最优传输理论与Laguerre镶嵌（Laguerre tessellation）的介观尺度方法研究，得出该方法能高效模拟数十至数万软体系统，严格保证体积约束并支持任意形状与变形，在3D生物物理情景中表现优越，重要意义在于连接了多种现有模型视角并为复杂组织模拟提供新工具，论文发表于《SCIENCE ADVANCES》。

主要关键技术方法包括：采用Laguerre镶嵌作为静态模型，通过半离散最优传输问题求解Kantorovich势（Kantorovich potentials）w_i以满足严格体积约束；动态模型基于梯度下降方程描述粒子位置演化，引入不可压缩力（incompressibility force）与可定制点力F^point与表面力F^surf；利用GPU加速的半符号懒惰张量框架（KeOps与GeomLoss库）求解大规模离散最优传输；采用L-BFGS-B拟牛顿法优化对偶势；在3D中通过SurfaceNets算法进行表面网格化以计算曲率与界面力；成本函数（cost function）支持参数化（如协方差矩阵Σ_i定义椭球）与水平集势函数（level-set potential functions）形式，变形参数α控制软硬度，边界条件采用周期化向量处理。

研究结果部分保留原文小标题并总结如下：

Static model: Laguerre tessellation

研究人员考虑N个粒子由位置x_i与体积v_i定义，Laguerre单元L_i由不等式c(x,x_i)?w_i≤ c(x,x_j)?w_j定义，其中c为成本函数，w_i为Kantorovich势以确保|L_i|=v_i。Laguerre镶嵌是约束最小化问题T_c的唯一解，当c为L²成本且w_i=0时退化为标准Voronoi图。通过定制c可实现任意边界形状，如幂成本产生气泡状镶嵌，扰动L²成本产生随机波动，孤立粒子边界对应于c(x,x_i)=w_i?w₀的水平集。

Dynamical model

研究人员提出粒子位置的一阶梯度下降方程?_i= ?τ_i?_xiT_c，源于OT理论的不相容力使粒子互斥并考虑形状与接触面。若成本仅依赖连接向量则力可表为边界?L_i上的积分，解释为内压力，总动量在无主动力与边界且τ_i≡1时守恒。通用方程引入点力F^point与表面力F^surf及噪声，允许动态体积v_i与成本c_i，通过参数化成本（如椭球Σ_i）或水平集势c_i=e^?φ_i实现形状动力学。

Soft-body simulations: First examples and benchmark

研究人员以L²成本模拟沙漏域中软球自由落体与3D聚集体指数生长至N=50000细胞（约1天计算），以及3D可变形椭球在盒域中的run-and-tumble运动。基准测试显示复杂度线性于N与M^d，GPU实现对小N次优，2000次迭代在M=512、N=10000时需约35小时。

Emergence of orientational order for rod-shape active particles

研究人员用成本函数c(x,x_i)= (|x?x_i|/r₀(θ))^α编码棒状（rod-shape）粒子如杆菌，方向由PCA（principal components analysis）主导向量定义。自推进运动因碰撞产生弯曲转向，仅体积排除交互即出现向列序参量（nematic order parameter）?的取向有序涌现，大系统t=12左右锐变且有序维持至t=25，密集小系统持久有序。

Fluid-solid phase transition

研究人员用幂成本c_α(x,y)=λ_α^?1|x?y|^α模拟软球，α为变形参数，大α硬化小α软化。在2D活跃布朗粒子（active Brownian particles）中，相位图显示归一化有效扩散D?_eff与形状指数?σ_i?随α与自推进速度c₀变化，流体相保持扩散（r²>0.9）但D?_eff<1，固体相扩散消失呈六边形密排（σ_i≈1），过渡发生于σ_i=1.015，与形状指数阈值3.81/(2√π)不严格对应。平均场近似导出含Monge-Ampère方程的PDE系统。

Deformation-driven motion

研究人员在水平集形式中取c=e^?φ，φ包含球L²项与化学吸引剂梯度δu(x,x_i)的变形项f(δ)，两种f分别产生伸长形与扇状形，仅用?_i方程因形态偏置于梯度方向诱导趋化迁移，同时保持体积与非重叠，优于传统水平集难以处理多细胞交互的局限。

Cell sorting via surface interactions

研究人员用L²成本Laguerre单元但引入界面力F_i←j=∫_Γij(γ_ij|κ|+η_ij/|x_i?x_j|)n? dσ，γ与η为表面张力（surface tension），满足Young-Dupré关系cos(θ/2)=η/(2γ)。在3D混合双类型细胞（蓝b、橙o）聚集体中，通过压实参数k_o、k_b、k_ob与相对软度γ?、排斥比η?划分六区域，复现同亲分选（homophilic sorting）的总/部分包埋、介质分离、异亲棋盘分选等构型，如η?>1时降低硬颗粒压实从包埋A经部分B到分离C，η?<1沿1/γ?轴转变，相等体积时界面曲率为零，不同体积与k显示小颗粒内吞趋势。

讨论部分总结：研究人员引入基于Laguerre镶嵌与最优传输的框架，在同一框架模拟个体可变形细胞至复杂组织聚集体，支持任意形状与严格体积约束，联系点粒子、水平集与DCM（Deformable Cell Models）等方法，并提供高效GPU实现。未来需直接生物验证如拟合3D成像数据及囊胚发育建模；理论上扩展活跃顶点理论至Laguerre镶嵌与能量最小化框架，推导软体连续介质模型；包含流体-结构交互以解决介质零成本简化的动量守恒；改进GPU性能与隐式成本及表面交互的全GPU实现。结论为该方法验证对比经典模型并具生物物理应用潜力，尤其3D情景下传统方法成本高昂时优势明显，生物动机如囊胚形态发生与分选的后续工作正在进行。