PhysiCell多尺度仿真集成指南:从SBML到COPASI的跨尺度建模实践
2026/9/14 21:32:46 网站建设 项目流程

PhysiCell这个系列写到现在已经是第十四篇。前面聊过它的细胞力学、细胞周期、微环境扩散这些单点功能,但隔三差五就有读者在评论区问:PhysiCell到底能不能和其他生物仿真软件配合起来用?它内置了BioFVM,那能不能接COPASI、Smoldyn、CompuCell3D这些生态,把分子、细胞、组织整个串成一条线?

这问题问到了多尺度仿真真正的核心痛点。单个软件里的功能再花哨,价值也有限;真正难的是“跨尺度怎么接”。PhysiCell的定位是细胞级代理模型,配合BioFVM能算氧气、葡萄糖这类信号的扩散和消耗,但它默认不擅长完整的细胞内信号网络,也不能直接读SBML、SBtab这类生信标准格式。真要搞一套带p53-MDM2反馈、NF-kB串扰、代谢重编程的多尺度模型,单打独斗确实吃力。

这篇来拆“集成”这件事。我把它分成四块:微环境参数怎么和实验数据对齐,分子网络怎么通过SBML挂进每个细胞,PhysiCell和Smoldyn、CompuCell3D、VCell这类工具怎么分工,以及我实际踩过的几个坑。刚入门的朋友可以照着做,已经在跑多尺度模型的也能从集成路径和排错思路里找到可复用的东西。

1. PhysiCell在多尺度仿真版图中的真实位置,以及为什么集成绕不开

1.1 PhysiCell擅长什么,不擅长什么

PhysiCell是一个开源的、C++写的、面向大规模3D细胞群体的多尺度仿真平台。它的核心设计是Agent-Based Modeling,也就是ABM。每个细胞是一个智能体,带有位置、体积、粘附、机械碰撞、细胞周期、分泌吸收能力,以及一系列表型决策规则。微环境这一层由内置的BioFVM负责,求解反应-扩散方程,模拟氧气、药物、细胞因子等底物的空间梯度。

这套设计的最大优势是,细胞群体在组织尺度上的行为涌现,比如肿瘤球生长、免疫细胞浸润、血管新生,它能跑得又快又直观。尤其适合做“细胞之间相互作用导致整体行为变化”这类问题,比如CAR-T细胞进入实体瘤后为什么容易被耗竭,T细胞在什么条件下才能穿过致密基质。

但它的短板也很明显。第一,细胞内部的分子信号网络和代谢通路不是它的主场。你可以用手写ODE的方式在细胞函数里塞进去几个分子,但一旦涉及几十个物种、几百个反应,没有标准化的格式来进行管理就是一场灾难。第二,它默认不认SBML这类生信通用格式,外面的工具导出的网络模型不能直接拖进去用。第三,实验数据导入和闭环比较吃力,比如需要根据病理切片上的细胞分布来初始化细胞位置,就需要外部脚本做图像处理。第四,再往上接药物动力学参数或者器官级模型时,PhysiCell本身不提供接口,只能靠我们自己搭桥。

所以多尺度仿真的一个现实情况是:PhysiCell适合做“细胞群体行为”这一层的主干,但分子层的通路定义、微环境层的参数标定、数据层的统计可视化,都需要跟别的工具做集成。

1.2 集成到底分成哪几个层次

我做了几年的多尺度模型,越来越倾向于把集成分层来看。每一层有它自己的代表工具和要解决的问题,别混在一起谈。

层次代表工具在PhysiCell工作流里的角色
分子通路层SBML、COPASI、BionetGen、VCell定义或转换细胞内部的生化反应网络,生成ODE后挂到每个细胞的决策逻辑上
微环境层BioFVM、实验氧/药物浓度数据定义底物的扩散系数、降解速率、边界条件,与体外实验曲线对齐
数据交换层Python、ParaView、MATLAB前处理(图像分割、初始化细胞位置)和后处理(统计、出图、动画)
外部模拟层Smoldyn、CompuCell3D、VCell处理特定尺度的建模需求,例如分子粒子级或需要精确细胞形状的场景

判断标准其实很直接。当单个细胞的行为规则取决于细胞内某个分子的浓度、某个信号通路的开关状态时,就必须接分子网络层。当仿真结果要跟体外实验曲线对比时,就必须确保微环境参数一致。当数据集大到用文本挨个翻看不现实时,就必须接数据交换层。

1.3 集成的本质是标准格式,不是硬凑接口

不少人一听到“集成”就以为是要在PhysiCell里封装一个什么API。实际上,绝大多数跨软件的中转站是标准格式本身。SBML负责描述生化网络,CSV和VTK负责传递空间和细胞状态,JSON/XML负责传配置。把这些标准化了,工具之间的衔接就顺了。所以下文我会反复强调SBML和VTK这两个格式,因为它们是这个生态里的通用语言。

2. 最容易上手的集成入口:先把微环境与实验数据对齐

2.1 BioFVM在PhysiCell里的实际耦合方式

PhysiCell默认内置BioFVM,很多人不知道这其实就是微环境集成的一环。BioFVM是一套独立于PhysiCell的“反应-扩散方程求解器”,被集成进来以后,每个体素里都在解一组PDE,描述各种底物的扩散、衰减和被细胞吸收/分泌的过程。氧气、葡萄糖、药物、化疗因子、细胞因子都以这种“底物”的形式存在于仿真中。

在PhysiCell的配置里,微环境变量在PhysiCell_settings.xml中定义。比如要加一个氧气变量,大概是这样的:

<microenvironment_setup> <variable name="oxygen" units="mmHg" diffusion_coefficient="100000" decay_rate="0.1" initial_condition="38"/> </microenvironment_setup>

这里每个参数背后都有物理含义。diffusion_coefficient控制氧气在组织里的扩散能力,单位通常是微米平方每分钟;decay_rate是底物的自然降解速度;initial_condition是初始浓度;如果设了dirichlet_condition,则表示边界浓度固定。调整这些参数会直接影响肿瘤球内部有没有缺氧区、坏死核心长什么样。

2.2 新增一种细胞因子,并把它耦合到细胞决策

实际做集成的时候,我们经常需要加入自定义信号分子,比如TNF-α或IL-2。做法不难,在XML里加一个变量,然后在自定义细胞函数里找到该底物在微环境中的索引,设置分泌/吸收速率。大致逻辑是这样:

void my_cell_model(Cell* pCell, Phenotype& phenotype, double dt) { int tnf_index = microenvironment.find_density_index("TNFa"); // 每个细胞都可以分泌TNF-α phenotype.secretion.secretion_rates[tnf_index] = 5.0; }

这是最朴素的一层集成:底物由细胞分泌,又反过来影响细胞行为。比如当TNF-α浓度超过阈值时,细胞进入凋亡程序,这就形成了“微环境-细胞表型”的闭环。

2.3 和实验氧分布数据对齐的具体做法

微环境集成里最容易忽略的是参数标定。很多实验数据是体外培养测得的氧浓度梯度,比如肿瘤球在特定深度会出现缺氧区。我们要做的是让仿真里的氧分布曲线和实验曲线尽量重合,这时候就需要调diffusion_coefficientdecay_rate和细胞的氧气消耗率。

我的建议是分两步。第一步先不管细胞,把纯扩散的氧分布跑出来,和实验的空白对照对齐。第二步再放入细胞,调消耗率。如果直接一上来就调所有参数,很容易过拟合,而且出了问题根本定位不到源头。这算是我做过好几个项目之后的一个经验:微环境参数对齐是一切上层集成的前提,这层不对齐,后面挂上分子网络只会更乱。

3. 分子通路集成的核心通道:用libSBML把生化网络挂进每个细胞

3.1 为什么SBML是这个场景里的标准交换格式

SBML,全称Systems Biology Markup Language,是系统生物学领域用来描述生化反应网络的XML标准格式。它定义了几类核心对象:物种(物种)、反应(反应)、速率定律(速率定律)、参数(参数)、单位(单位)等。几乎所有主流分子网络工具都支持SBML的导入导出。

对PhysiCell来说,我们需要接的是一个“从外部工具生成网络,然后进入细胞内ODE”的通道。如果不用SBML,就得把几十个反应公式一个个手抄到C++代码里,效率低不说,还特别容易抄错。用SBML等于有了一个所有工具都能认的中间语言:在COPASI里把p53网络调好参数,导出SBML;在BionetGen里写免疫受体信号规则,转成SBML;甚至在VCell里做空间反应-扩散模型,也能生成SBML。然后我再在PhysiCell这一侧统一解析。

3.2 环境准备:安装libSBML

解析SBML,建议直接用libSBML。它是对SBML标准最完整的官方解析库,支持C、C++、Java、Python等语言。PhysiCell本身是C++项目,所以我一般用C++ API来解析。

Ubuntu下的安装:

sudo apt-get install libsbml-dev

macOS下:

brew install libsbml

Windows下建议去SBML官网下载预编译库,然后把include和lib路径配到编译器里。装好以后可以先用一个简单的C++程序验证解析器是否正常。

3.3 在PhysiCell项目中解析并执行SBML模型

PhysiCell的自定义模块通常放在custom_modules/目录下。我的做法是单独写一个MyIntracellularModel类,在初始化时调用SBMLReader把模型读进来,把物种名称、初值、反应速率律都放到一张映射表里。代码结构类似这样:

#include <sbml/SBMLTypes.h> #include <map> #include <string> class MyIntracellularModel { public: std::map<std::string, double> species_values; void load_from_sbml(const std::string& filename) { SBMLReader reader; SBMLDocument* doc = reader.readSBML(filename); Model* model = doc->getModel(); if (!model) return; // 读取物种和初值 for (unsigned int i = 0; i < model->getNumSpecies(); ++i) { Species* sp = model->getSpecies(i); species_values[sp->getId()] = sp->getInitialConcentration(); } // 读取反应速率律公式 for (unsigned int i = 0; i < model->getNumReactions(); ++i) { Reaction* rx = model->getReaction(i); KineticLaw* kl = rx->getKineticLaw(); std::string formula = kl->getFormula(); // 这里需要把formula字符串解析成可执行表达式 // 可以自己写一个轻量表达式解析器,或者利用libSBML的FormulaParser } } };

拿到了速率定律字符串后,难点在“怎么把字符串变成可计算的函数”。libSBML自带的FormulaParser能解析SBML的公式语法,但如果你需要更灵活的数值积分,我建议用下面两种方式之一:要么自己实现一个轻量表达式解析器,支持加减乘除、幂、常用的函数(exp、log、sqrt),然后把函数指针存起来;要么在COPASI里直接把模型生成C代码,再手动摘出核心的速率表达式。

我在实际项目里更倾向于后者。因为SBML的速率律字符串千奇百怪,自己写表达式解析器维护成本很高。COPASI导出C代码后,公式已经是人类可读的标准C语言,移植到PhysiCell的自定义细胞函数里就非常直接。

3.4 在细胞函数里做状态更新并影响表型

当分子网络解析完成,接下来就是把它和细胞行为绑定。我在自定义细胞函数里的标准写法是:从custom_data里读取当前分子浓度,用外部解析好的速率律算一个时间步的增量,然后更新状态,再根据关键分子的浓度改变细胞表型。

void my_cell_model(Cell* pCell, Phenotype& phenotype, double dt) { double akt = pCell->custom_data["akt"]; double bad = pCell->custom_data["bad"]; // 一个非常简化的AKT-BAD网络示意 double k1 = 2.0; double k2 = 0.8; double d_akt = k1 * (1.0 - akt) - k2 * akt * bad; double d_bad = k2 * akt * bad - 0.5 * bad; // 欧拉积分,实际建议用RK4或限制单步增量 akt += d_akt * dt; bad += d_bad * dt; pCell->custom_data["akt"] = akt; pCell->custom_data["bad"] = bad; // 分子状态影响表型:AKT高表达时,细胞倾向于增殖 if (akt > 0.8) { phenotype.cycle.data.transition_rate(0, 1) = 0.9; } else { phenotype.cycle.data.transition_rate(0, 1) = 0.1; } }

关键点在于“时间步的限制”。PhysiCell的dt通常是为了解决细胞运动、扩散而设置的,分子网络的ODE尺度可能完全不一样。如果反应速率常数比较大,直接用欧拉积分容易数值发散。我的经验是先把时间步切到足够细,比如:如果反应时间尺度在分钟量级,而PhysiCell的外部时间步是0.1分钟,可以接受;但如果是毫秒级过程,就必须在细胞函数里做子循环,把dt切成更小步长,或者使用隐式求解器。

3.5 单位不归一化,仿真一定跑飞

SBML集成里最大的坑就是单位换算。SBML模型常用时间单位是秒,PhysiCell内部时间是分钟;SBML里浓度单位可能是mol/L或mmol/L,而PhysiCell微环境的浓度单位可能是mmHg或uM;空间尺度上,SBML可以是基于体积单位,不是基于微米网格。

举个我踩过的例子。从COPASI导出的一个代谢网络,速率常数是以秒为单位的,直接挂进PhysiCell后,每个真实秒的时间被当成一分钟来算,等于所有反应速率被放大了60倍。十几个时间步之后,分子浓度从1e-6直接变成了1e12,一开始我还以为是模型稳定性问题,排查很久才发现纯粹是单位没换。

所以建议在加载SBML模型时,写一个预处理脚本,把单位统一换算成PhysiCell的base units:时间用分钟,空间用微米,浓度用uM或者mmHg。libSBML里的UnitDefinition其实可以提供单位换算信息,但更省力的做法是在COPASI里就把模型单位改成min和uM再导出。做完这一步再进PhysiCell,能省掉很多麻烦。

4. 和CompuCell3D、Smoldyn、VCell这些工具怎么分工,什么时候才需要联合

4.1 三种底层模拟机制的本质差异

很多读者会拿PhysiCell和CompuCell3D、Smoldyn、VCell、COPASI做对比,其实它们并不是同一个层面的工具,底层机制完全不同。把每个工具的本质搞清楚,才知道什么时候该联合、什么时候谁替代谁。

工具底层模型核心尺度强项与PhysiCell的典型配合方式
PhysiCellAgent-Based,中心力模型细胞/组织大规模细胞群体行为、微环境扩散主干仿真
CompuCell3DCellular Potts模型细胞/组织细胞形状、细胞-细胞接触、边界张力对形态敏感的问题做交叉验证
Smoldyn分子粒子随机游走分子/细胞膜受体-配体相互作用、分子扩散轨迹提供膜表面信号或小范围梯度参数
VCellPDE/ODE空间建模分子到细胞反应-扩散的空间模型、分子通路建模用VCell生成空间反应-扩散模型,再简化为微环境参数
COPASIODE/随机/代谢网络分子通路参数拟合、稳态分析、剂量响应先在这里调节SBML网络,再导出

PhysiCell里的细胞是“球形的带半径的粒子”,通过中心力模型来处理碰撞、黏附,这种机制的好处是计算快、能跑百万细胞,坏处是细胞形状永远是圆的或者被挤成多边形,没法精确模拟上皮细胞的极化形态。如果研究的问题高度依赖细胞形状,比如集体迁移中的头尾极性、细胞重塑导致的组织折叠,那CompuCell3D的Cellular Potts模型反而更合适。

Smoldyn则是更底层的存在。它几乎是给分子做布朗运动模拟的,每个分子都是一个粒子,在空间里随机游走,遇到配体就结合。用Smoldyn去模拟跨膜受体的聚集过程,能得到受体在膜上的空间分布,这些参数可以用来修正PhysiCell里细胞对信号分子的响应阈值。

VCell是个被低估的工具。它本身就是做空间反应-扩散方程的,可以用来在连续空间里模拟一个简化的组织片段里的信号梯度,然后把梯度参数化成PhysiCell微环境里的边界条件或初始条件。这种方式比直接猜一个扩散系数要科学得多。

4.2 一个联合工作流的实例:p53-NF-kB串扰信号的跨尺度搭建

拿我之前做过的“肿瘤微环境中的p53-NF-kB串扰”的例子来说一下完整的联合路径。

第一步,在COPASI里搭建p53和NF-kB的动态网络,调好参数,做一遍稳态和动态分析,确认振荡行为正常。第二步,把网络导出为SBML。第三步,用libSBML解析SBML,在PhysiCell每个细胞里实现这套ODE,细胞状态变量存在custom_data里。第四步,使用BioFVM在微环境层模拟TNF-α的扩散。第五步,在细胞函数里读取局部的TNF-α浓度,作为NF-kB通路的输入;NF-kB激活后又促进细胞分泌更多细胞因子,反馈回微环境。

这就是一个标准的双向耦合多尺度模型:分子网络影响细胞行为,细胞行为改变微环境,微环境回过头来调节分子网络。这种闭环是单靠一个PhysiCell或者单靠一个COPASI都跑不出来的,必须做集成。

4.3 可视化与后处理层面的集成:ParaView + PhysiCell Studio

仿真做完不能只看CSV数字,可视化这一步同样属于集成。PhysiCell默认会输出带细胞位置和属性的数据文件,另外还可以输出VTK格式,可以直接拖进ParaView里渲染。ParaView能做的事情包括:显示肿瘤球的三维结构、用颜色映射表达细胞内部的p53浓度、把微环境的氧浓度用透明体渲染出来。这对于跟生物学家讨论模型结果特别有用。

PhysiCell Studio是官方出的交互式可视化工具,可以在仿真过程中暂停、拖动细胞、实时看指标。我做集成调试的时候一般先把输出间隔设得小一点,在Studio里逐步检查,确认分子浓度和表型变化符合预期,再放大规模跑正式仿真。

5. 实测踩坑记录:单位、版本、并发、可复现性

5.1 单位换算错误导致浓度指数爆炸

这个前面已经提过,但还是值得单独列出来。现象是仿真跑了几百个时间步后,某个分子浓度直接变成1e15甚至NaN,整片组织全部死掉。第一次遇到我还以为是解析器写错了,后来把每个物种的初值和速率常数逐一打印出来,才发现是SBML里k值的单位是1/s,而PhysiCell的dt单位是min,差了60倍;另外浓度单位也有mol/L和uM之间的1000倍差距。解决方案非常简单:写一个单位统一脚本,或者干脆在COPASI导出前手动把所有单位改成min和uM。以后凡是接入新模型,我第一件事就是先看单位定义,不看单位直接跑就是给自己埋坑。

5.2 libSBML版本冲突导致的链接错误

PhysiCell本身用的编译系统比较传统,如果你同时装了系统自带的和conda的libSBML,很容易出现undefined reference。我遇到过的情况是:系统里是libsbml5,但某个工具链需要libsbml6,编译时一堆符号找不到。排查思路是用ldconfig -p | grep sbml看系统库里有哪些版本,用ldd看可执行文件实际链接的是哪个so。然后统一include路径和库路径,尽量全部指向同一个小版本的libSBML。

5.3 并行仿真里的随机性和可复现性问题

PhysiCell支持用OpenMP并行,但并行之后随机数怎么分配是个大问题。如果所有线程共用同一个全局随机种子,那么每次运行结果都可能不同,而且不同机器、不同核心数下结果都不一样。这在做科研场景下非常致命,审稿人让你复现结果的时候你回一句“每次跑出来都不同”就麻烦了。

我的经验是:先在单线程模式下跑一条基线,确定要复现的标准配置;然后用显式设置随机种子的方式跑并行;如果并行逻辑改变了随机数调用顺序,最简单的方式是给每个线程单独初始化随机流,并把线程数固定下来。PhysiCell本身支持在配置里设随机种子,但如果你在细胞函数里另外用了自己的随机数生成器,一定要保证它的种子也受控。

5.4 输出文件直接把磁盘撑爆

跑大规模3D仿真,每隔0.1分钟存一次所有细胞位置和微环境场,跑几千分钟,输出几十个GB是常态。我第一次跑肿瘤免疫模型时,一个晚上醒来发现磁盘满了,所有数据都白跑。后来学乖了,输出间隔设到合理范围,比如每1分钟或每10分钟存一次,并且只保存自己关心的一组细胞属性;微环境场可以隔几帧才存,或者只输出某个切面的数据。VTK的压缩选项也开着。这类细节看着小,实际能救你一命。

问题现象对策
单位不统一浓度爆炸、NaN解析SBML前统一时间/min、浓度/uM单位
libSBML版本冲突链接时undefined reference统一include和lib路径到同一个小版本
并行不可复现每次运行结果不同单线程基线 + 固定随机种子 + 固定线程数
输出文件过大磁盘满、跑白加大输出间隔、只存关键属性、开VTK压缩

6. 再往上走:实验数据闭环、机器学习决策和更大的生态

6.1 用实验影像数据做初始化和验证

很多真实场景需要把实验切片里的细胞分布直接搬进仿真。这时集成路径通常是这样:先用ImageJ或CellProfiler对H&E染色切片做细胞分割,得到每个细胞的位置和形态;再把坐标数据转成PhysiCell的初始细胞布局;最后在仿真里复现出和切片相似的空间结构。这一步的价值在于,模型不再是从完全理想化的随机构型开始跑,而是从真实的组织状态出发,做预测的基准就完全不一样了。

6.2 把机器学习决策塞进细胞模型

再进阶一点,细胞的行为规则不一定要手写固定阈值。最近有不少工作在做“细胞级强化学习”。思路是把细胞当成一个决策智能体,观测是局部的细胞因子浓度、氧气浓度、邻居密度,动作是增殖、迁移、凋亡或者分泌某种细胞因子,奖励函数由整体肿瘤控制率来定义。

具体到PhysiCell,可以用Python离线训练一个简单的决策树或线性策略,导出成一个轻量文件,在C++细胞函数里加载并按输入特征判断动作。这种“Python训练-PhysiCell推理”的集成路径,可以把机器学习模型嵌入到大规模ABM仿真里。实测跑下来,决策树比神经网络更容易在C++里嵌入,而且可解释性强。如果你要上神经网络,可以考虑把ONNX Runtime也编译进去,但复杂度会上去不少。

6.3 和PK/PD、器官级模型的串联

还有一类集成是往上走,比如把PhysiCell放在一个更大的药物研发平台上。上游的PK模型计算药物在血液里的浓度曲线,输出给PhysiCell作为微环境里药物底物的边界条件;PhysiCell算出肿瘤细胞数量变化和耐药细胞比例,再反馈给药代动力学模型。这种串联能回答一些特别实际的问题,比如用药方案是每周一次大剂量更好,还是每天低剂量更好。

PhysiCell在这类体系里就是中间的“药效动力学”模块。它不是一个封闭的黑盒,只要输入输出格式标准化,就能嵌入到任何一种已经跑通的药研流程里。这是我认为多尺度集成最值得投入的方向。

个人经验是,无论做哪一种集成,都建议从最小闭环开始。先拿一个三到五个分子的SBML模型,在COPASI里跑通,再挂进PhysiCell,确认分子浓度曲线和COPASI趋势一致,然后再逐步扩大网络规模。不要一上来就拿一个大而全的代谢网络做集成,否则一旦跑飞,你根本分不清是解析器的问题、单位的问题还是模型本身的数值稳定性问题。最小闭环跑通了,后面加模块就是重复劳动,真正卡壳的地方大概率已经提前排掉了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询