本研究由加拿大滑铁卢大学化学工程系的Ittisak Promma、Marc G. Aucoin、Nasser Mohieddin Abukhdeir和Hector Budman共同完成,发表于*Computers and Chemical Engineering*期刊2024年第189卷(文章编号108806)。研究聚焦于大规模通气搅拌生物反应器中微生物代谢与流体动力学的耦合建模问题,提出了一种将计算流体力学(computational fluid dynamics, CFD)信息引导的隔室模型与动态通量平衡(dynamic flux balance, DFB)模型相结合的新方法,以降低大规模生物反应器模拟的计算成本并提高预测精度。
生物反应器在生物技术和生物制品生产中应用广泛,但从小试到工业规模的放大过程中,由于复杂多相流与微生物代谢的相互作用,反应器内会出现底物、溶氧等浓度的显著空间梯度,导致微生物群体代谢状态不均一,进而影响整体性能。例如,在大肠杆菌培养中,高葡萄糖和低氧环境会触发溢流代谢(overflow metabolism)产生乙酸。仅依靠实验难以充分理解微生物行为与胞外环境之间的动态耦合机制,因此需要发展能够准确捕捉生物反应器内多物理场耦合的建模方法。已有研究采用条件代谢模型(conditional metabolic models)结合欧拉-欧拉法或隔室模型,但此类模型参数较多,易导致过拟合。动态通量平衡分析(dynamic flux balance analysis, DFBA)基于细胞最优资源利用假设,能以较少参数描述代谢,但其线性规划(linear programming, LP)求解计算复杂度高,尤其在与流体动力学模型耦合时面临巨大计算负担。
本研究的总体目标是开发一种降低大规模生物反应器模拟计算成本的建模方法。具体目标包括:从DFB模型构建二叉搜索树(binary search tree, BST)代谢模型以降低计算复杂度;将BST代谢模型与基于流动信息的水动力隔室模型耦合;通过对比大肠杆菌补料分批发酵的实验数据验证耦合模型。
在方法上,研究首先进行多相CFD模拟以获得时间平均的流场信息。模拟采用OpenFOAM中的欧拉-欧拉双流体模型,使用多重参考系法处理叶轮运动,对一台工作体积22 m³、直径2.09 m、配备四挡板和四叶Rushton涡轮桨的30 m³生物反应器进行瞬态三相流动模拟。通过求解惰性示踪剂的纯对流传质方程估计混合时间,得到95%混合时间为226.2 s,与文献报道的165 s和250 s相近,验证了CFD流场的合理性。第二步是基于时间平均速度分量进行流动信息引导的隔室化(compartmentalization),采用Tajsoleiman等人的方法,参数阈值为Δp_i = 0.5 m/s、最小隔室体积V_min = 0.1 m³。初始得到101个隔室,但发现其混合时间预测仅为60.05 s,远低于CFD结果,原因是隔室化方法在相邻再循环区之间引入了“短路”路径。为此,研究者基于CFD流场洞察进行用户引导修正,在每对叶轮之间设置水平分界面,禁止跨平面合并隔室,最终得到108个隔室,混合时间预测提高到189.34 s,显著接近CFD值。第三步是建立隔室模型并与BST代谢模型耦合。每个隔室视为良好混合的连续搅拌釜(CSTR),其质量守恒方程包含积累项、内部对流传质项、进料项、出料项、反应项和其他项。反应速率通过DFB模型计算。本研究所用DFB模型基于Mahadevan等人提出的大肠杆菌模型,原模型包含葡萄糖、乙酸、溶解氧和生物质四种物质及四条代谢途径,经过两项关键修改以适用于补料分批操作并提高长期预测能力:一是引入生物质一级降解速率;二是引入酶与生物质比例对葡萄糖和氧摄取能力的调控函数。修改后的DFB模型以最大化生物质生长速率为目标,受氧摄取上限、葡萄糖摄取上限(Monod型动力学)以及各物质浓度非负约束的限制。研究将该LP问题通过多参数线性规划(multiparametric linear programming, mp-LP)重表述为点定位问题,识别出18个临界区域,并构建了包含67个节点、树高7的非平衡二叉搜索树。计算复杂度测试表明,BST代谢模型计算代谢通量每个时间步、每个隔室仅需约0.040 ms,而使用SciPy的单纯形算法求解原始LP需要0.466 ms;在约10^6自由度的CFD模拟中,单步代谢通量计算时间从430 s降至0.65 s。第四步是利用实验室规模(7 L)生物反应器的实验数据拟合代谢模型参数。由于小规模反应器可合理近似为完全混合,采用单隔室模型进行参数估计,使用Nelder-Mead法最小化模型预测值与实验值的误差平方和。第五步是将拟合参数应用于大规模反应器模拟并验证。模拟中采用PI控制器调节溶氧水平,通过调整k_L a控制变量维持溶氧设定值,控制参数和位置信息根据文献设定。模拟发现使用实验室规模拟合的最大葡萄糖摄取速率K_G,max为10 mmol/g/h时显著高估大规模反应器中的葡萄糖浓度,因此将K_G,max调整为20 mmol/g/h,并基于比葡萄糖摄取速率的尺度差异给出了可能的生物学解释。
模拟结果与Xu等人和Pigou与Morchain报道的实验及模拟数据进行了对比。在葡萄糖浓度方面,使用调整后的K_G,max值,模型能够更好地预测葡萄糖浓度的变化趋势,特别是在补料批次阶段。在乙酸浓度方面,模型在整个操作过程中能够准确捕捉浓度趋势,包括由于氧气控制作用导致的第二个乙酸峰值,而此前Pigou和Morchain的模型未能预测到该峰值。生物质浓度的空间梯度极小,模型预测与实验数据一致。通过校正Akaike信息准则(AICc)比较,本研究模型在实验室和大规模反应器上均优于对照模型。模拟还提供了反应速率和代谢通量的空间分布信息,表明葡萄糖主要在上部区域被消耗,上部区域同时利用好氧和厌氧途径;中部区域在指数补料阶段主要通过厌氧途径消耗葡萄糖,约10小时后由于氧气水平提高转为好氧途径;底部区域由于靠近空气供应氧气充足,但碳源已被上游消耗,生长速率较低。此外,由于探针角位置未在文献中报告,研究通过模拟评估了该不确定性对葡萄糖浓度的影响,结果表明各高度处浓度范围窄,角度变化的影响不敏感。研究还考察了流动振荡对浓度梯度的影响,通过在体积流量矩阵上引入正弦振荡(M=0.7,ω=3 h⁻¹)模拟湍流引起的局部流动波动,结果显示模拟与实验数据的定性一致性提高,尤其是顶部探针位置的葡萄糖浓度波动得到更好的再现,说明局部水动力振荡对代谢反应和浓度梯度具有重要影响。
本研究的结论表明,所开发的耦合代谢通量/隔室模型能够以显著的计算效率和精度模拟大规模通气生物反应器。BST代谢模型通过将FBA重述为点定位问题并进行多参数线性规划解析优化器映射,大幅降低了局部生长速率计算的复杂度,不仅便于与隔室模型集成,也为在优化设计和控制框架内进行动态模拟提供了可能。该方法能直接集成控制策略,更精确地预测代谢物浓度并再现实验过程中的振荡现象。更重要的是,该模型可在标准台式计算机上以几分钟的计算时间模拟超过40小时的操作过程,而耦合CFD与代谢通量模型的模拟时间估计需要10年以上,显示出该方法在实时模型预测控制和优化设计中的巨大应用潜力。研究的创新点在于:识别并解决了隔室化过程中引入跨混合区“短路”的问题,提高了混合时间预测精度;发展了基于BST的代谢模型以降低计算复杂度;实现了代谢模型与水动力隔室模型的耦合,并通过大肠杆菌补料分批发酵实验数据验证了模型性能。