CN111210867A - 基于代谢分析及酶调控的微生物生长状态估计方法 - Google Patents

Abstract

基于代谢分析及酶调控的微生物生长状态估计方法，属于工业微生物制造过程监测技术及应用领域，涉及一种将代谢分析、酶调控与先进测量方式相结合的微生物生长过程状态估计方法，由模型建立、测量方案构建、线性化环节、基于模型的滤波组成。首先进行细胞代谢途径分析和简化，并引入酶调控机制，建立描述微生物生长代谢状态的数学模型；然后通过近红外光谱、吸光度检测构建过程测量方程；在此基础上应用线性化处理，给出基于数学模型的状态估计方法；用于实现反映微生物生长状态的发酵产物浓度、底物浓度和生物量的估计。该方法适用于微生物代谢网络及化学计量学矩阵已知的酶促反应过程。

Description

基于代谢分析及酶调控的微生物生长状态估计方法

技术领域

本发明属于工业微生物制造过程监测技术及应用领域，涉及一种将代谢分析、酶调控与先进测量方式相结合的微生物生长过程状态估计方法，用于实现反映微生物生长状态的发酵过程产物浓度、底物浓度和生物量的估计。

背景技术

利用微生物细胞的繁殖和生长进行工业发酵生产，有着广泛的应用前景。在发酵过程中，很多生产状态参数不能直接测量，例如产物和底物浓度、生物量浓度等，尤其是无法观测到细胞内的代谢活动，这在很大程度上影响产品质量控制以及节能降耗，更难以实现自动化及智能化生产。一般通过建立微生物发酵过程的宏观经验动力学数学模型，反映产物和底物浓度、生物量浓度的变化规律：

其中x_s是底物浓度向量，c是生物量浓度，x_P是胞内代谢产物浓度，μ是比生长速率，π是比生成速率，Y_c、Y_p、m为动力学模型参数。

在模型(1)-(3)基础上，借助一些可检测可获得的过程变量，应用状态估计或滤波理论，可以估计出产物和底物浓度、生物量浓度等，诸如采用卡尔曼滤波、粒子滤波等方法，实现对上述发酵状态的估计。但是现有发酵过程状态估计技术所依赖的数学模型，未能考虑细胞内的代谢活动，导致发酵过程模型与实际情况有差距，难以得到理想的估计效果。事实上微生物的代谢网络和酶的调控特征对发酵工艺的优化有重要指导意义，而且很多发酵过程的代谢网络是清晰的，比如代谢通量分析、通量平衡分析等，可以在特定条件下获得胞内代谢通量的分布，其化学计量学矩阵可以获得，目前文献中已有针对聚·-羟基丁酸发酵过程，将代谢途径分析与酶的调控作用相结合，建立了细胞内外物质浓度的动力学模型，这些都为发酵过程状态估计技术的发展提供了基础。

另一方面，发酵过程状态估计还依赖于过程检测技术进步，只有当微生物生长过程有些变量能够获得，才能实施状态估计。传统的发酵过程检测多为离线进行，独立于生产环节之外，实时性差，近年来随着生化检测仪器的发展，诸如吸光度测量、质谱色谱仪，特别是近红外光谱主要反映含氢基团对近红外光的倍频和组合频吸收，从而获得含氢有机物的化学结构信息，再结合化学计量学方法建立光谱数据与有机物浓度的校正模型，可实现有机物组分的快速检测，为状态估计的应用提供了可能。

发明内容

本发明提出一种基于代谢分析及酶调控的微生物生长状态估计方法，由模型建立、测量方案构建、线性化环节、基于模型的滤波环节组成。首先进行细胞代谢途径分析和简化，并引入酶调控机制，建立描述微生物生长代谢状态的数学模型；然后通过近红外光谱、吸光度检测构建过程测量方程；在此基础上应用线性化处理，给出基于数学模型的状态估计方法；实现发酵过程产物浓度、底物浓度和生物量状态的估计。该方法适用于微生物代谢网络及化学计量学矩阵已知的酶促反应过程，包括利用微生物菌体将底物转化为产物的聚羟基丁酸、谷氨酸、赖氨酸等生产过程。

本发明的技术方案：

基于代谢分析及酶调控的微生物生长状态估计方法，步骤包括：

第一步：建立发酵菌体细胞内外代谢物的通量平衡方程：

在一般宏观经验动力学数学模型(1)和(3)基础上，考虑细胞内外代谢物的通量平衡，引入胞外代谢物的化学计量学矩阵S_x和胞内产物的化学计量矩阵S_p：

式中r是反应速率。

第二步：在发酵代谢网络基础上选取关键代谢途径：

通常微生物发酵过程的代谢网络都较为复杂，从底物到产物的所有可能代谢途径很多。根据需要选择以底物糖分解最大化、目标代谢物生产最大化、产品产率最大化为目标，通过代谢途径分析并结合实验数据，挑选出

条关键代谢途径，并将反应速率r分解到

条代谢途径：

其中矩阵Z的每一列表示一条代谢途径，

是

条代谢途径的反应速率向量。

第三步：在反应速率中引入酶调控变量：

酶的催化作用影响反应速率，为此考虑催化酶的活性和合成作用：

式中

种催化酶的活性调控向量为

合成调控向量为

b是生物质的催化活性分数，diag(·)代表对角矩阵，

是未考虑酶调控的速率向量，一般可选用Monod形式。u和v以碳源吸收单位向量L_c为基础进行选取，u＝γdiag(L_c)和v＝δdiag(L_c)，其中γ,δ>0为标称化系数。

第四步：获得发酵过程状态方程：

根据上述式(4)-(7)和式(2)，设定x为发酵过程状态向量x＝[x_s x_p c]^T，反映发酵过程底物浓度、胞内代谢产物浓度、总生物量浓度变化的状态方程

表示如下：

第五步：构建发酵过程测量方程：

方案一是采用直接测量，建立近红外光谱数据与葡萄糖、果糖等底物以及发酵产物之间的定量模型，一般生物量也可通过氨基葡萄糖法推测，实现生物量的快速检测。利用近红外光谱构建测量方程，具体步骤1：光谱预处理和波长选择，方法包括相关系数法、方差分析法、无信息变量消除法(UVE)、连续投影算法(SPA)、区间偏最小二乘法(iPLS)、模拟退火算法(SA)、遗传算法(GA)；步骤2：光谱数据X与浓度状态x之间定量模型的建立，方法包括主成分分析(PCA)、偏最小二乘(PLS)、人工神经网络(ANN)、支持向量回归(SVR)、高斯过程回归(GPR)；步骤3：构成如下输出y＝h(x)的测量方程：

y＝h(x)＝H[x_s x_p c ]^T (9)

其中H为构造的测量矩阵，对应于x中的各分量x_s,x_p,c，一般检测一个浓度即可。

方案二是采用间接测量，考虑到培养液中菌体的生长情况和生物量浓度成正比，利用培养液的吸光度测量y来构建测量方程：

y＝h(x)＝k_px_pc+k_r(c-x_p) (10)

其中k_p、k_r为测量方程参数。

第六步：基于状态和测量方程的发酵状态递推估计：

在状态方程(8)和测量方程(9)或(10)的基础上，针对非线性系统、以及系统噪声和测量噪声ω,v是非高斯的情形，直接采用贝叶斯估计，在t_k时刻有：

x_k＝f(x_k-1)+ω_k (11)

y_k＝h(x_k)+v_k (12)

然后利用先验概率与当前的测量值来获取发酵过程状态的后验概率分布，得到状态估计值，

具体包括预测步：

p(x_k|Y_k-1)＝∫p(x_k|x_k-1)p(x_k-1|Y_k-1)dx_k-1 (13)

和更新步：

式中Y_k＝[y_k,…,y_k-N-1]是测量值序列,N为正整数。

第七步：状态方程和测量方程线性化：

对非线性系统进行线性化处理，考虑时刻t_k下的状态

k是正整数，对上述状态方程

进行线性化，其中

为线性化操作点，Δx_k是状态增量，通过一阶泰勒展开，取一次近似：

式中状态转移矩阵为：

同理对测量方程(10)线性化：

式中测量矩阵为：

第八步：基于线性化方程的发酵状态递推估计：

考虑系统和测量噪声为高斯白噪声，在时刻t_k的协方差矩阵分别为Q_k和R_k，采用扩展卡尔曼滤波，发酵状态估计步骤1是通过状态方程

计算预测值

通过以下协方差方程：

来计算协方差的预测值

具体计算方法是利用状态方程和协方差方程，从t_k-1时刻的状态估计值

和估计协方差P_k-1，积分到当前t_k时刻获得预测值

和

步骤2是发酵状态及其协方差的估计值更新，根据时刻t_k下的测量值y_k，更新获得估计值

和协方差P_k，具体如下：

其中，K_k为卡尔曼增益。

第九步：基于线性化方程的发酵状态批处理估计：

考虑降低系统和测量噪声统计特性以及系统初始条件对估计精度的影响，采用有限脉冲响应(FIR)滤波，利用当前t_k之前的N个测量值形成批次数据，批次起点时刻t_j，j＝k-N+1，省略噪声项后线性化状态和测量方程为x_k,j＝F_k,jx_j和y_k,j＝H_k,jx_j，其中F_k,j,H_k,j是F,H的扩展后的系统和观测矩阵，而

分别为状态和测量值的N维数据矩阵；由此发酵状态的无偏FIR估计(UFIR)如下：

其中，

是F_k,j相关的系列矩阵构成的增益。

基于线性化方程的发酵状态批处理估计，还包括最优FIR估计(OFIR)、最小方差无偏FIR估计(MVUFIR)算法。

本发明的有益效果：对微生物发酵生产的产物浓度、底物浓度和生物量浓度的估计，一般在宏观经验动力学模型基础上，采用卡尔曼滤波、粒子滤波等估计方法来实现，未能考虑细胞内的代谢活动，导致难以得到理想的估计效果。本发明以细胞内外代谢物的通量平衡方程为基础，利用代谢途径分析并引入酶的调控机制，建立更精准的发酵过程状态方程，同时利用先进的检测仪器构建测量方程，最后实施的估计方法考虑到系统的非线性、噪声统计特性、状态初始值等条件。通过对聚羟基丁酸、谷氨酸、赖氨酸等生产过程中的底物浓度、产物浓度、生物量浓度的递推估计，结果表明，本发明的方法具有良好的估计效果，显示其在微生物监测、优化与控制方面的应用潜力。

附图说明

图1为基于代谢分析及酶调控的微生物生长状态估计方法流程图。

图2底物果糖(FRU)的状态估计图。

图3底物氯化铵(AMC)的状态估计图。

图4生物量的状态估计图。

具体实施方式

下面参照附图1对本发明的具体实施方式做进一步说明。

将一般宏观经验动力学的底物模型(1)和产物模型(3)用代谢物通量平衡方程代替，形成由化学计量矩阵和反应速率描述的动态模型(4)和(5)。

以发酵过程某种最大化目标为依据，分析发酵代谢网络并结合实验数据，挑选出几条关键代谢途径，按照式(6)用这几条代谢途径的反应速率来表述整体的反应速率；

对关键代谢途径的反应速率，按照式(7)引入酶调控变量，考虑催化酶的活性和合成调控作用；

根据前述式(4)-(7)和式(2)，形成反映发酵过程底物、产物和总生物量浓度变化的状态方程(8)；

采用近红外光谱直接测量方式，通过光谱预处理和波长选择、光谱数据与发酵(产物、底物、生物量)浓度之间的定量建模，建立发酵过程测量方程(9)；若无法直接测量，采用培养液的吸光度检测，间接建立测量方程(10)；

基于方程(8)-(10)对发酵状态进行估计，若已知过程噪声的统计特性，并且对方程的非线性不进行特殊处理，采用贝叶斯估计，经过预测步骤(13)和更新步(14)，获得状态估计值；

若将方程(8)-(10)线性化为(15)和(17)，并且考虑高斯白噪声，则采用扩展卡尔曼滤波(19)-(22)递推获得发酵状态估计值；进一步放宽对噪声统计特性及系统初始条件的要求，采用有限脉冲响应滤波(23)，获得状态估计值。

实施效果分析

采用本发明提出的基于代谢分析及酶调控的微生物生长状态估计方法，应用于聚羟基丁酸发酵过程，其底物包括果糖(FRU)和氯化铵(AMC)，浓度分别为x₁和浓度x₂，即x_s＝[x₁,x₂]^T，x_P是胞内代谢物聚羟基丁酸浓度。聚羟基丁酸代谢网络可能的代谢途径超过120条，根据产量目标选取5个关键基元模式

相应的状态方程参数如下：

S_mZ＝[0.44 0.06 0.00 0.00 -1.08]，

μ＝[0.44 0.60 0.24 0.62 0.45]，

L_c＝[1.0 1.0 1.0 1.0 2.51]；

利用培养液的吸光度测量，相应的测量方程参数如下：

k_p＝4.76，k_r＝2.8；

初始状态值设为：

x₀＝[20 1.5 0.18 0.125]^T；

考虑干扰和测量噪声为高斯白噪声，其协方差矩阵分别设置为：

Q＝diag(10^-3,10^-3,10^-3,10^-3)，R＝10^-3；

对状态方程和测量方程线性化，采用扩展卡尔曼滤波，果糖(FRU)、氯化铵(AMC)和生物量的估计值如图2、图3和图4；在聚羟基丁酸合成第I阶段，培养基中加入底物果糖和氯化铵，分别提供碳源和氮源，生物量呈指数增长，约12小时后，氮源逐渐耗尽，仍有碳源；第II阶段，果糖作为合成产物的材料，生物量线性增长，约32小时后，果糖耗尽，生物体停止生长；进入第III阶段，产物的合成保持稳定。由此本发明提出的方法对底物、生物量、产物的动态变化有良好的估计效果。

Claims (8)

1.基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，步骤包括：

第一步：建立发酵菌体细胞内外代谢物的通量平衡方程：

在如下的一般宏观经验动力学数学模型基础上：

其中x_s是底物浓度向量，c是生物量浓度，x_P是胞内代谢产物浓度，μ是比生长速率，π是比生成速率，Y_c、Y_p、m为动力学模型参数；

考虑细胞内外代谢物的通量平衡，在式(1)和(3)中引入胞外代谢物的化学计量学矩阵S_x和胞内产物的化学计量矩阵S_p，得到：

式中r是反应速率；

第二步：在发酵代谢网络基础上选取关键代谢途径：

根据需要的优化目标，通过代谢途径分析并结合实验数据，选择

条关键代谢途径，并将反应速率r分解到

条代谢途径：

其中矩阵Z的每一列表示一条代谢途径，

是

条代谢途径的反应速率向量；

第三步：在反应速率中引入酶调控变量：

考虑催化酶的活性调控和合成调控作用，对

条代谢途径，分别引入活性调控向量

和合成调控向量

式中b是生物质的催化活性分数，diag(·)代表对角矩阵，

是未考虑酶调控的速率向量；u和v以碳源吸收单位向量L_c为基础进行选取，u＝γdiag(L_c)和v＝δdiag(L_c)，其中γ,δ>0为标称化系数；

第四步：获得发酵过程状态方程：

根据上述式(4)-(7)和式(2)，设定x为发酵过程状态向量x＝[x_s x_p c]^T，反映发酵过程底物浓度、胞内代谢产物浓度、总生物量浓度变化的状态方程

表示如下：

第五步：构建发酵过程测量方程：

方案一是采用直接测量，建立近红外光谱数据X与发酵过程状态向量x之间的定量模型；利用近红外光谱构建测量方程，具体过程为：

步骤1：光谱预处理和波长选择；

步骤2：光谱数据X与浓度状态x之间定量模型的建立；

步骤3：构成如下输出y＝h(x)的测量方程：

y＝h(x)＝H[x_s x_p c]^T (9)

其中H为构造的测量矩阵，对应于x中的各分量x_s,x_p,c，一般检测一个浓度即可；

方案二是采用间接测量，考虑到培养液中菌体的生长情况和生物量浓度成正比，利用培养液的吸光度测量y来构建测量方程：

y＝h(x)＝k_px_pc+k_r(c-x_p) (10)

其中k_p、k_r为测量方程参数；

第六步：基于状态和测量方程的发酵状态递推估计：

在状态方程(8)和测量方程(9)或(10)的基础上，针对非线性系统、以及系统和测量噪声ω,v是非高斯的情形，直接采用贝叶斯估计，在t_k时刻有：

x_k＝f(x_k-1)+ω_k (11)

y_k＝h(x_k)+v_k (12)

然后利用先验概率与当前的测量值来获取发酵过程状态的后验概率分布，得到状态估计值，具体包括预测步：

p(x_k|Y_k-1)＝∫p(x_k|x_k-1)p(x_k-1|Y_k-1)dx_k-1 (13)

和更新步：

式中Y_k＝[y_k,…,y_k-N-1,]是测量值序列，N为正整数；

第七步：状态方程和测量方程线性化：

对非线性系统进行线性化处理，考虑时刻t_k下的状态

k是正整数，对上述状态方程

进行线性化，其中

为线性化操作点，Δx_k是状态增量，通过一阶泰勒展开，取一次近似：

式中状态转移矩阵为：

同理对测量方程(10)线性化：

式中测量矩阵为：

第八步：基于线性化方程的发酵状态递推估计：

考虑系统和测量噪声为高斯白噪声，在时刻t_k的协方差矩阵分别为Q_k和R_k，采用扩展卡尔曼滤波，发酵状态估计：

步骤1是通过状态方程

计算预测值

通过以下协方差方程：

来计算协方差的预测值

具体计算方法是利用状态方程和协方差方程，从t_k-1时刻的状态估计值

和估计协方差P_k-1，积分到当前t_k时刻获得预测值

和

步骤2是发酵状态及其协方差的估计值更新，根据时刻t_k下的测量值y_k，更新获得估计值

和协方差P_k，具体如下：

其中，K_k为卡尔曼增益；

第九步：基于线性化方程的发酵状态批处理估计：

考虑降低系统和测量噪声统计特性以及系统初始条件对估计精度的影响，采用有限脉冲响应FIR滤波，利用当前t_k之前的N个测量值形成批次数据，批次起点时刻t_j，j＝k-N+1，省略噪声项后线性化状态和测量方程为x_k,j＝F_k,jx_j和y_k,j＝H_k,jx_j，其中F_k,j,H_k,j为F,H批处理形式的系统和观测矩阵，而

分别为状态和测量值得批次数据矩阵；由此发酵状态的无偏FIR估计UFIR如下：

其中，

是F_k,j相关的系列矩阵构成的增益。

2.根据权利要求1所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，所述第二步在发酵代谢网络基础上选取关键代谢途径，具体选择方案是：以底物糖分解最大化为优化指标选择代谢途径；或者以代谢物生产最大化为优化指标选择代谢途径；或者以产品产率最大化为优化指标选择代谢途径。

3.根据权利要求1或2所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第五步方案一的步骤1中光谱预处理和波长选择，其方法包括相关系数法、方差分析法、无信息变量消除法UVE、连续投影算法SPA、区间偏最小二乘法iPLS、模拟退火算法SA、遗传算法GA。

4.根据权利要求1或2所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第五步方案一的步骤2中光谱数据X与浓度状态x之间定量模型的建立，其方法包括主成分分析PCA、偏最小二乘PLS、人工神经网络ANN、支持向量回归SVR、高斯过程回归GPR。

5.根据权利要求3所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第五步方案一的步骤2中光谱数据X与浓度状态x之间定量模型的建立，其方法包括主成分分析PCA、偏最小二乘PLS、人工神经网络ANN、支持向量回归SVR、高斯过程回归GPR。

6.根据权利要求1、2或5所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第九步中，基于线性化方程的发酵状态批处理估计，估计算法还包括最优FIR估计OFIR、最小方差无偏FIR估计MVUFIR。

7.根据权利要求3所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第九步中，基于线性化方程的发酵状态批处理估计，估计算法还包括最优FIR估计OFIR、最小方差无偏FIR估计MVUFIR。

8.根据权利要求4所述的基于代谢分析及酶调控的微生物生长状态估计方法，其特征在于，第九步中，基于线性化方程的发酵状态批处理估计，估计算法还包括最优FIR估计OFIR、最小方差无偏FIR估计MVUFIR。

Images