潘妍如，劉飛

（江南大學(xué)輕工過程先進(jìn)控制教育部重點實驗室，江蘇無錫214122）

引 言

在微生物發(fā)酵過程中，由于物化和生化反應(yīng)的混合交疊性，組分異常復(fù)雜，大多數(shù)狀態(tài)不能直接測量，例如產(chǎn)物和底物濃度、生物量濃度等，尤其是無法觀測到細(xì)胞內(nèi)的代謝活動。隨著生物實驗技術(shù)的不斷進(jìn)步，越來越多的微生物代謝網(wǎng)絡(luò)得以精確地建立，借助代謝通量分析(MFA/DMFA)[1-3]、通量平衡分析[4-6]等，可以在特定的條件下獲得胞內(nèi)代謝通量的分布。一般將代謝網(wǎng)絡(luò)中的所有可能代謝途徑的集合稱為基元模式，通過基元模式分析，可以得到目標(biāo)產(chǎn)物的代謝產(chǎn)率、最優(yōu)合成途徑等，從而選取最優(yōu)代謝途徑[7-8]。所謂賽博模型(cybernetic model)則是基于簡化的代謝網(wǎng)絡(luò)，考慮細(xì)胞內(nèi)酶的動力學(xué)方程，并引入目標(biāo)函數(shù)實現(xiàn)細(xì)胞最優(yōu)調(diào)節(jié)，確定代謝途徑[9-10]。

近年來有學(xué)者將賽博模型引入到基元模式分析中，給出一種混合賽博模型(HCM)來描述微生物的代謝過程[11]，建立了細(xì)胞內(nèi)外物質(zhì)濃度的動力學(xué)方程，基于細(xì)胞內(nèi)部擬穩(wěn)態(tài)的假設(shè)，同時又考慮內(nèi)部代謝產(chǎn)物的動態(tài)變化，將非線性系統(tǒng)動力學(xué)與生物體的代謝信息聯(lián)系起來，更好地描述生物代謝過程[12-13]。

文獻(xiàn)[14-15]利用HCM，給出了聚羥基丁酸(PHB)的生長代謝過程的建模分析；文獻(xiàn)[16]則利用吸光度法來構(gòu)建PHB 濃度測量方程，將HCM 模型用于狀態(tài)估計，通過分析滾動時域和無跡卡爾曼兩種估計，說明了HCM 模型用于實際工業(yè)的狀態(tài)估計，需要尋找計算量少、權(quán)系數(shù)調(diào)整簡單的算法。本文基于聚羥基丁酸生產(chǎn)的HCM 模型，考慮其非線性結(jié)構(gòu)，采用擴展卡爾曼濾波實現(xiàn)PHB 代謝狀態(tài)估計。眾所周知，擴展卡爾曼濾波具有計算工作量小、容易實現(xiàn)的優(yōu)點[17-18]，廣泛地應(yīng)用于過程工業(yè)中，在微生物發(fā)酵領(lǐng)域已有很多成功應(yīng)用[19-21]。本文首先分析現(xiàn)有HCM 模型，然后對其非線性動力學(xué)進(jìn)行線性化切線逼近，最后通過濾波器的遞歸更新，實現(xiàn)PHB 生產(chǎn)過程中產(chǎn)物濃度、底物濃度和生物量狀態(tài)的估計。

1 模型描述

1.1 動力學(xué)模型

細(xì)胞內(nèi)代謝產(chǎn)物PHB 主要是由外部代謝物果糖(FRU)和氯化銨(AMC)來控制產(chǎn)生的。根據(jù)質(zhì)量守恒原則，可得細(xì)胞內(nèi)外代謝物的通量平衡方程[15]，表示如下：

其中，xFRU、xAMC為胞外兩個代謝物濃度向量，m為胞內(nèi)代謝物濃度，c為生物量濃度，r為反應(yīng)速率，μ為比生長速率，Sx(nx×nr)代表胞外代謝物的化學(xué)計量學(xué)矩陣，Sm(nm×nr)代表胞內(nèi)代謝物的化學(xué)計量學(xué)矩陣。

其中，mPHB為細(xì)胞內(nèi)代謝產(chǎn)物PHB 的濃度，Sm,s(nm,s×nr)為PHB的化學(xué)計量矩陣。

1.2 代謝網(wǎng)絡(luò)與基元模式

通常完整的代謝網(wǎng)絡(luò)都較為復(fù)雜，包含眾多基元模式?；Ｊ绞谴x網(wǎng)絡(luò)中從底物到產(chǎn)物的所有可能代謝途徑的集合，可以通過Metatool[23]或CellNetAnalyzer 軟件[24-25]計算代謝網(wǎng)絡(luò)中所有的基元模式并進(jìn)行代謝途徑分析。在產(chǎn)物的代謝途徑分析中，不需要利用所有的基元模式，只需選取最關(guān)鍵的基元模式即能描述代謝行為[26]。因此，可以根據(jù)產(chǎn)物目標(biāo)和產(chǎn)量分析對基元模式進(jìn)行簡化，再結(jié)合實驗數(shù)據(jù)進(jìn)一步挑選出活躍的基元模式。

PHB 產(chǎn)物的代謝網(wǎng)絡(luò)如圖1 所示，共有36 個反應(yīng)，可以分解為122 個基元模式，通過代謝產(chǎn)率分析，最后選取5個最關(guān)鍵的基元模式。文獻(xiàn)[27]詳細(xì)介紹了基元模式的簡化以及關(guān)鍵基元模式的選取方法。

1.3 HCM描述

賽博模型的關(guān)鍵是考慮酶的催化，假設(shè)每個基元模式都由一種酶催化，酶的含量影響反應(yīng)過程進(jìn)而影響產(chǎn)物的生成，引入控制變量u和v來調(diào)節(jié)酶的合成和活性[28]。針對前述5 個關(guān)鍵基元模式，相應(yīng)的有5種催化酶。用基元模式分解來表示反應(yīng)速率r，表示如下：

其中Z(nr×nz)為矩陣，矩陣的每一列表示一個基元模式，rM是通過5 個基元模式的反應(yīng)速率向量，通過控制變量向量v來表示：

圖1 PHB產(chǎn)物的代謝網(wǎng)絡(luò)圖[14]Fig.1 Metabolic network of PHB product [14]

其中v=[v1v2v3v4v5]是控制5 種催化酶活性的控制變量，b是生物質(zhì)的催化活性分?jǐn)?shù)，u=[u1u2u3u4u5]是控制5 種催化酶合成的控制變量向量。

酶作為催化劑，在整個代謝過程中起到關(guān)鍵作用，考慮催化相關(guān)基元模式的酶，用以下微分方程來描述[29]：

其中，e表示催化酶的濃度，α是酶的固有合成速率，rEM是酶的合成速率，diag(β)e表示酶的降解，μe表示由生長導(dǎo)致的酶的稀釋，酶合成速率rEM由控制變量向量u控制：

其中fC為碳吸收單位向量。

根據(jù)前述動力學(xué)模型并結(jié)合HCM 描述，擬對PHB 生產(chǎn)中的果糖和氯化銨濃度、產(chǎn)物PHB 濃度、總生物量濃度進(jìn)行估計，設(shè)定x=[xFRU,xAMC,mPHB,c]Τ為狀態(tài)向量，則HCM系統(tǒng)的狀態(tài)方程描述如下：

式中f(x)表示非線性函數(shù)?？紤]到培養(yǎng)液中菌體的生長情況和菌體濃度呈正比，因此可以通過測量培養(yǎng)液里菌體的生長情況，來獲得菌體的濃度，具體利用波長λ= 600nm 的培養(yǎng)液的吸光度測量[14]，構(gòu)建PHB濃度的測量方程：

其中，y為培養(yǎng)液的吸光值，kPHB、kres為參數(shù)，總生物量c由儲存的PHB 含量和剩余的生物量兩部分組成。通過測量方程來觀察細(xì)胞大小的變化，從而反映細(xì)胞中產(chǎn)物PHB含量的變化。

2 底物、產(chǎn)物和生物量濃度估計

下面利用擴展卡爾曼濾波器對上述包含底物濃度xFRU、xAMC，產(chǎn)物濃度mPHB，生物量濃度c的狀態(tài)向量x進(jìn)行估計，基本思路是對描述PHB 生產(chǎn)過程的非線性函數(shù)f(x)進(jìn)行泰勒展開，舍去高階量，將非線性模型線性化，包括過程和測量模型的線性化，并利用相應(yīng)的協(xié)方差矩陣，不斷更新狀態(tài)估計值。

考慮實際生產(chǎn)過程中存在干擾以及測量時的噪聲，HCM 系統(tǒng)的狀態(tài)方程式(12)和測量方程式(13)分別表達(dá)如下：

式中，h(x)是關(guān)于狀態(tài)的非線性函數(shù)，w和v是均值為零的高斯白噪聲。Q和R分別是w和v的協(xié)方差矩陣。

在卡爾曼濾波估計點對非線性系統(tǒng)進(jìn)行線性化，在時刻tk有：

式中，x?(tk)為線性化操作點狀態(tài)，Δx(tk)是狀態(tài)增量，對其進(jìn)行一階泰勒展開，得：

式中狀態(tài)轉(zhuǎn)移矩陣為：

同理，一階線性化觀測方程：

式中觀測矩陣為：

基于上述線性化描述，擴展卡爾曼濾波器的預(yù)測和更新步驟如下。

(1)狀態(tài)估計值及其協(xié)方差矩陣表示如下：

對式(21)、式(22)進(jìn)行積分，積分過程從前一時刻的狀態(tài)估計值x?(tk-1)和協(xié)方差P(tk-1)開始，在積分結(jié)束時可得到下一時刻的狀態(tài)預(yù)測值x-(tk)和協(xié)方差P-(tk)。

(2)在tk時刻，把量測值y(tk)代入狀態(tài)估計和協(xié)方差估計表達(dá)式中，如下所示：

卡爾曼增益

狀態(tài)更新

協(xié)方差更新

3 仿真結(jié)果

針對聚羥基丁酸的代謝模型，狀態(tài)方程和測量方程參照式(12)、式(13)，式中rkinM,i和rkinEM,i采用Monod模型描述：

根據(jù)文獻(xiàn)[15]中的描述，化學(xué)計量矩陣定義如下：

總生物量的生長率定義如下：

從果糖和氯化銨的摩爾質(zhì)量及其化學(xué)計量因子計算碳吸收單位fC的數(shù)量為：

狀態(tài)方程中的參數(shù)設(shè)置在表1 中給出，四個狀態(tài)變量的初始值設(shè)定為：

測量方程中的參數(shù)設(shè)置為：

針對擴展卡爾曼濾波，干擾和測量噪聲協(xié)方差矩陣分別設(shè)置為Q= diag(10-3,10-3,10-3,10-3)，R=10-3。

根據(jù)以上所述，利用擴展卡爾曼濾波對混合賽博模型中的四個狀態(tài)變量進(jìn)行估計，狀態(tài)估計的仿真結(jié)果如圖2～圖5所示。

表1 模型參數(shù)設(shè)置［15］Table 1 Setting of model parameters［15］

圖2 底物果糖的狀態(tài)估計仿真圖Fig.2 Simulation results of state estimation of FRU

圖3 底物氯化銨的狀態(tài)估計仿真圖Fig.3 Simulation results of state estimation of AMC

圖4 生物量的狀態(tài)估計仿真圖Fig.4 Simulation results of state estimation of biomass

圖5 產(chǎn)物PHB的狀態(tài)估計仿真圖Fig.5 Simulation results of state estimation of PHB

通過對混合賽博模型的仿真，可以觀察聚羥基丁酸合成的三個階段，根據(jù)圖2～圖5 各階段進(jìn)行分析。第Ⅰ階段，培養(yǎng)基中加入底物果糖(FRU)和氯化銨(AMC)，分別提供足夠的碳源和氮源，總生物量呈指數(shù)增長，這段時間內(nèi)全部合成非PHB 生物質(zhì)，約12 h 后，氮源逐漸耗盡，但仍有碳源可用。第Ⅱ階段，氯化銨完全耗盡，果糖作為合成PHB 的內(nèi)部存儲材料，這時生物體將外部碳源儲存于細(xì)胞體內(nèi)合成PHB，由于催化活性生物量保持不變，所以生長是線性的。第Ⅲ階段，約32.5 h 后，果糖完全耗盡，生物體停止生長，產(chǎn)物的合成保持穩(wěn)定。觀察擴展卡爾曼濾波對混合賽博模型的跟蹤效果，在第Ⅰ階段，底物果糖和氯化銨的估計值與模型有微小偏差，尤其是在氯化銨逐漸耗盡的過程中；在第Ⅱ和第Ⅲ階段，擴展卡爾曼濾波對底物、產(chǎn)物的動態(tài)變化均有良好的估計效果，在10～20 h 之間對生物量濃度的估計有偏差。圖6 是吸光度測量的曲線圖，對比擴展卡爾曼濾波值和模型測量值，兩者間偏差不大，收斂效果較好。

4 結(jié) 論

圖6 模型測量值和狀態(tài)估計測量值的比較Fig.6 Comparison of model measurements and state estimation measurements

將混合賽博模型(HCM)引入微生物代謝狀態(tài)估計領(lǐng)域，利用擴展卡爾曼濾波處理HCM 中的非線性動力學(xué)，實現(xiàn)對PHB 生產(chǎn)過程中底物和產(chǎn)物濃度、生物量濃度的遞推估計，具有良好的估計效果，顯示了其在微生物代謝優(yōu)化和控制方面的應(yīng)用潛力。擴展卡爾曼濾波方法適用于具有高斯噪聲分布的非線性系統(tǒng)，即要求過程的動態(tài)噪聲和測量噪聲是高斯白噪聲，然而實際微生物發(fā)酵過程往往具有不確定性、時變性、非高斯性，另外還有初始值難以確定等問題，都將影響擴展卡爾曼估計精度，今后的研究工作可以考慮對擴展卡爾曼濾波器進(jìn)行改進(jìn)，或者采用粒子濾波、FIR濾波等狀態(tài)估計方法。

符 號 說 明

b——生物質(zhì)中催化活性部分的分?jǐn)?shù)

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

e——酶水平，g/L

fC——碳吸收單位向量

KAMC——氯化銨的飽和常數(shù)，g/L

KFRU——果糖的飽和常數(shù)，g/L

KPHB——PHB的飽和常數(shù)，g/L

ke——酶合成速率常數(shù)向量，h-1

kPHB——模型參數(shù)

kr——速率常數(shù)向量，h-1

kres——模型參數(shù)

m——胞內(nèi)代謝物濃度，g/L

mPHB——PHB產(chǎn)物濃度，g/g

P——投入回報率

r——反應(yīng)速率，h-1

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

——rEM動力學(xué)的分量

rM——通過基元模式的反應(yīng)速率，h-1

——rM動力學(xué)的分量

Sm——胞內(nèi)代謝物的化學(xué)計量矩陣

Sm,s——產(chǎn)物PHB的化學(xué)計量矩陣

Sx——胞外代謝物的化學(xué)計量矩陣

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

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

xAMC——底物氯化銨濃度，g/L

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

y——培養(yǎng)液的吸光值

Z——基元模式矩陣

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

β——酶的消耗常數(shù)向量，h-1

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