基于混合赛博模型的微生物代谢状态估计

潘妍如，刘飞

（江南大学轻工过程先进控制教育部重点实验室，江苏无锡214122）

引 言

在微生物发酵过程中，由于物化和生化反应的混合交叠性，组分异常复杂，大多数状态不能直接测量，例如产物和底物浓度、生物量浓度等，尤其是无法观测到细胞内的代谢活动。随着生物实验技术的不断进步，越来越多的微生物代谢网络得以精确地建立，借助代谢通量分析(MFA/DMFA)[1-3]、通量平衡分析[4-6]等，可以在特定的条件下获得胞内代谢通量的分布。一般将代谢网络中的所有可能代谢途径的集合称为基元模式，通过基元模式分析，可以得到目标产物的代谢产率、最优合成途径等，从而选取最优代谢途径[7-8]。所谓赛博模型(cybernetic model)则是基于简化的代谢网络，考虑细胞内酶的动力学方程，并引入目标函数实现细胞最优调节，确定代谢途径[9-10]。

近年来有学者将赛博模型引入到基元模式分析中，给出一种混合赛博模型(HCM)来描述微生物的代谢过程[11]，建立了细胞内外物质浓度的动力学方程，基于细胞内部拟稳态的假设，同时又考虑内部代谢产物的动态变化，将非线性系统动力学与生物体的代谢信息联系起来，更好地描述生物代谢过程[12-13]。

文献[14-15]利用HCM，给出了聚羟基丁酸(PHB)的生长代谢过程的建模分析；文献[16]则利用吸光度法来构建PHB 浓度测量方程，将HCM 模型用于状态估计，通过分析滚动时域和无迹卡尔曼两种估计，说明了HCM 模型用于实际工业的状态估计，需要寻找计算量少、权系数调整简单的算法。本文基于聚羟基丁酸生产的HCM 模型，考虑其非线性结构，采用扩展卡尔曼滤波实现PHB 代谢状态估计。众所周知，扩展卡尔曼滤波具有计算工作量小、容易实现的优点[17-18]，广泛地应用于过程工业中，在微生物发酵领域已有很多成功应用[19-21]。本文首先分析现有HCM 模型，然后对其非线性动力学进行线性化切线逼近，最后通过滤波器的递归更新，实现PHB 生产过程中产物浓度、底物浓度和生物量状态的估计。

1 模型描述

1.1 动力学模型

细胞内代谢产物PHB 主要是由外部代谢物果糖(FRU)和氯化铵(AMC)来控制产生的。根据质量守恒原则，可得细胞内外代谢物的通量平衡方程[15]，表示如下：

其中，xFRU、xAMC为胞外两个代谢物浓度向量，m为胞内代谢物浓度，c为生物量浓度，r为反应速率，μ为比生长速率，Sx(nx×nr)代表胞外代谢物的化学计量学矩阵，Sm(nm×nr)代表胞内代谢物的化学计量学矩阵。

其中，mPHB为细胞内代谢产物PHB 的浓度，Sm,s(nm,s×nr)为PHB的化学计量矩阵。

1.2 代谢网络与基元模式

通常完整的代谢网络都较为复杂，包含众多基元模式。基元模式是代谢网络中从底物到产物的所有可能代谢途径的集合，可以通过Metatool[23]或CellNetAnalyzer 软件[24-25]计算代谢网络中所有的基元模式并进行代谢途径分析。在产物的代谢途径分析中，不需要利用所有的基元模式，只需选取最关键的基元模式即能描述代谢行为[26]。因此，可以根据产物目标和产量分析对基元模式进行简化，再结合实验数据进一步挑选出活跃的基元模式。

PHB 产物的代谢网络如图1 所示，共有36 个反应，可以分解为122 个基元模式，通过代谢产率分析，最后选取5个最关键的基元模式。文献[27]详细介绍了基元模式的简化以及关键基元模式的选取方法。

1.3 HCM描述

赛博模型的关键是考虑酶的催化，假设每个基元模式都由一种酶催化，酶的含量影响反应过程进而影响产物的生成，引入控制变量u和v来调节酶的合成和活性[28]。针对前述5 个关键基元模式，相应的有5种催化酶。用基元模式分解来表示反应速率r，表示如下：

其中Z(nr×nz)为矩阵，矩阵的每一列表示一个基元模式，rM是通过5 个基元模式的反应速率向量，通过控制变量向量v来表示：

图1 PHB产物的代谢网络图[14]Fig.1 Metabolic network of PHB product [14]

其中v=[v1v2v3v4v5]是控制5 种催化酶活性的控制变量，b是生物质的催化活性分数，u=[u1u2u3u4u5]是控制5 种催化酶合成的控制变量向量。

酶作为催化剂，在整个代谢过程中起到关键作用，考虑催化相关基元模式的酶，用以下微分方程来描述[29]：

其中，e表示催化酶的浓度，α是酶的固有合成速率，rEM是酶的合成速率，diag(β)e表示酶的降解，μe表示由生长导致的酶的稀释，酶合成速率rEM由控制变量向量u控制：

其中fC为碳吸收单位向量。

根据前述动力学模型并结合HCM 描述，拟对PHB 生产中的果糖和氯化铵浓度、产物PHB 浓度、总生物量浓度进行估计，设定x=[xFRU,xAMC,mPHB,c]Τ为状态向量，则HCM系统的状态方程描述如下：

式中f(x)表示非线性函数。考虑到培养液中菌体的生长情况和菌体浓度呈正比，因此可以通过测量培养液里菌体的生长情况，来获得菌体的浓度，具体利用波长λ= 600nm 的培养液的吸光度测量[14]，构建PHB浓度的测量方程：

其中，y为培养液的吸光值，kPHB、kres为参数，总生物量c由储存的PHB 含量和剩余的生物量两部分组成。通过测量方程来观察细胞大小的变化，从而反映细胞中产物PHB含量的变化。

2 底物、产物和生物量浓度估计

下面利用扩展卡尔曼滤波器对上述包含底物浓度xFRU、xAMC，产物浓度mPHB，生物量浓度c的状态向量x进行估计，基本思路是对描述PHB 生产过程的非线性函数f(x)进行泰勒展开，舍去高阶量，将非线性模型线性化，包括过程和测量模型的线性化，并利用相应的协方差矩阵，不断更新状态估计值。

考虑实际生产过程中存在干扰以及测量时的噪声，HCM 系统的状态方程式(12)和测量方程式(13)分别表达如下：

式中，h(x)是关于状态的非线性函数，w和v是均值为零的高斯白噪声。Q和R分别是w和v的协方差矩阵。

在卡尔曼滤波估计点对非线性系统进行线性化，在时刻tk有：

式中，x∗(tk)为线性化操作点状态，Δx(tk)是状态增量，对其进行一阶泰勒展开，得：

式中状态转移矩阵为：

同理，一阶线性化观测方程：

式中观测矩阵为：

基于上述线性化描述，扩展卡尔曼滤波器的预测和更新步骤如下。

(1)状态估计值及其协方差矩阵表示如下：

对式(21)、式(22)进行积分，积分过程从前一时刻的状态估计值x̂(tk-1)和协方差P(tk-1)开始，在积分结束时可得到下一时刻的状态预测值x-(tk)和协方差P-(tk)。

(2)在tk时刻，把量测值y(tk)代入状态估计和协方差估计表达式中，如下所示：

卡尔曼增益

状态更新

协方差更新

3 仿真结果

针对聚羟基丁酸的代谢模型，状态方程和测量方程参照式(12)、式(13)，式中rkinM,i和rkinEM,i采用Monod模型描述：

根据文献[15]中的描述，化学计量矩阵定义如下：

总生物量的生长率定义如下：

从果糖和氯化铵的摩尔质量及其化学计量因子计算碳吸收单位fC的数量为：

状态方程中的参数设置在表1 中给出，四个状态变量的初始值设定为：

测量方程中的参数设置为：

针对扩展卡尔曼滤波，干扰和测量噪声协方差矩阵分别设置为Q= diag(10-3,10-3,10-3,10-3)，R=10-3。

根据以上所述，利用扩展卡尔曼滤波对混合赛博模型中的四个状态变量进行估计，状态估计的仿真结果如图2～图5所示。

表1 模型参数设置［15］Table 1 Setting of model parameters［15］

图2 底物果糖的状态估计仿真图Fig.2 Simulation results of state estimation of FRU

图3 底物氯化铵的状态估计仿真图Fig.3 Simulation results of state estimation of AMC

图4 生物量的状态估计仿真图Fig.4 Simulation results of state estimation of biomass

图5 产物PHB的状态估计仿真图Fig.5 Simulation results of state estimation of PHB

通过对混合赛博模型的仿真，可以观察聚羟基丁酸合成的三个阶段，根据图2～图5 各阶段进行分析。第Ⅰ阶段，培养基中加入底物果糖(FRU)和氯化铵(AMC)，分别提供足够的碳源和氮源，总生物量呈指数增长，这段时间内全部合成非PHB 生物质，约12 h 后，氮源逐渐耗尽，但仍有碳源可用。第Ⅱ阶段，氯化铵完全耗尽，果糖作为合成PHB 的内部存储材料，这时生物体将外部碳源储存于细胞体内合成PHB，由于催化活性生物量保持不变，所以生长是线性的。第Ⅲ阶段，约32.5 h 后，果糖完全耗尽，生物体停止生长，产物的合成保持稳定。观察扩展卡尔曼滤波对混合赛博模型的跟踪效果，在第Ⅰ阶段，底物果糖和氯化铵的估计值与模型有微小偏差，尤其是在氯化铵逐渐耗尽的过程中；在第Ⅱ和第Ⅲ阶段，扩展卡尔曼滤波对底物、产物的动态变化均有良好的估计效果，在10～20 h 之间对生物量浓度的估计有偏差。图6 是吸光度测量的曲线图，对比扩展卡尔曼滤波值和模型测量值，两者间偏差不大，收敛效果较好。

4 结 论

图6 模型测量值和状态估计测量值的比较Fig.6 Comparison of model measurements and state estimation measurements

将混合赛博模型(HCM)引入微生物代谢状态估计领域，利用扩展卡尔曼滤波处理HCM 中的非线性动力学，实现对PHB 生产过程中底物和产物浓度、生物量浓度的递推估计，具有良好的估计效果，显示了其在微生物代谢优化和控制方面的应用潜力。扩展卡尔曼滤波方法适用于具有高斯噪声分布的非线性系统，即要求过程的动态噪声和测量噪声是高斯白噪声，然而实际微生物发酵过程往往具有不确定性、时变性、非高斯性，另外还有初始值难以确定等问题，都将影响扩展卡尔曼估计精度，今后的研究工作可以考虑对扩展卡尔曼滤波器进行改进，或者采用粒子滤波、FIR滤波等状态估计方法。

符 号 说 明

b——生物质中催化活性部分的分数

c——生物量浓度，g/L

e——酶水平，g/L

fC——碳吸收单位向量

KAMC——氯化铵的饱和常数，g/L

KFRU——果糖的饱和常数，g/L

KPHB——PHB的饱和常数，g/L

ke——酶合成速率常数向量，h-1

kPHB——模型参数

kr——速率常数向量，h-1

kres——模型参数

m——胞内代谢物浓度，g/L

mPHB——PHB产物浓度，g/g

P——投入回报率

r——反应速率，h-1

rEM——酶的合成速率，h-1

——rEM动力学的分量

rM——通过基元模式的反应速率，h-1

——rM动力学的分量

Sm——胞内代谢物的化学计量矩阵

Sm,s——产物PHB的化学计量矩阵

Sx——胞外代谢物的化学计量矩阵

u——控制酶合成的控制变量

v——控制酶活性的控制变量

xAMC——底物氯化铵浓度，g/L

xFRU——底物果糖浓度，g/L

y——培养液的吸光值

Z——基元模式矩阵

α——酶的固有合成速率向量，h-1

β——酶的消耗常数向量，h-1

μ——比生长速率，h-1