基于混合賽博模型的微生物代謝狀態估計

潘妍如，劉飛

（江南大學輕工過程先進控制教育部重點實驗室，江蘇無錫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