岳向陽(yáng) 趙忠蓋 劉 飛

（江南大學(xué)a.輕工過(guò)程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室；b.自動(dòng)化研究所）

工業(yè)發(fā)酵中的操作條件可以調(diào)節(jié)微生物的生長(zhǎng)代謝環(huán)境，如底物流加速率用于控制基質(zhì)濃度，冷/熱水流加速率則能夠調(diào)整環(huán)境溫度，當(dāng)基質(zhì)濃度過(guò)低時(shí)會(huì)導(dǎo)致菌體營(yíng)養(yǎng)不良，而濃度過(guò)高時(shí)又會(huì)使菌體耗氧量增加，溫度則會(huì)提高或抑制酶的活性［1］，因此為提高發(fā)酵生產(chǎn)的效率，對(duì)操作條件的優(yōu)化至關(guān)重要。 精確有效的發(fā)酵過(guò)程模型是實(shí)現(xiàn)優(yōu)化的前提條件。 最初對(duì)操作條件的優(yōu)化是根據(jù)人類積累的經(jīng)驗(yàn)知識(shí)［2］，隨著對(duì)過(guò)程機(jī)理的逐漸了解，通過(guò)微分方程定義的機(jī)理模型得到使用［3］，后來(lái)計(jì)算機(jī)硬件水平不斷提高，采用歷史數(shù)據(jù)驅(qū)動(dòng)的數(shù)據(jù)模型獲得廣泛應(yīng)用［4］，但這些模型只能描述發(fā)酵過(guò)程的外部特征。

微生物發(fā)酵的本質(zhì)是細(xì)胞復(fù)雜的微觀代謝反應(yīng)，其反應(yīng)速率即代謝通量能夠定量地表征菌體內(nèi)部的生長(zhǎng)代謝狀態(tài)。 代謝通量難以實(shí)際測(cè)得［5］，而 常 用 的 動(dòng) 態(tài) 通 量 分 析（Dynamic Metabolic Flux Analysis，DMFA） 方法可以利用細(xì)胞代謝網(wǎng)絡(luò)模型和細(xì)胞外時(shí)變參數(shù)來(lái)估計(jì)代謝通量［6］。 因此建立代謝通量調(diào)控模型則可以從代謝機(jī)理的角度定量描述發(fā)酵過(guò)程中操作條件與代謝通量間的聯(lián)系，為操作條件的優(yōu)化提供代謝層面的指導(dǎo)。

微生物的生命周期可分為遲滯期、對(duì)數(shù)生長(zhǎng)期、穩(wěn)定期和凋亡期，其生長(zhǎng)代謝狀態(tài)在不同階段動(dòng)態(tài)變化，具有非線性特征，因此分段建模方法得以應(yīng)用。 Gao Y等使用高斯混合模型將樣本劃分為若干個(gè)子數(shù)據(jù)集后，分別建立子模型再融合為代謝通量調(diào)控模型［7］。 然而實(shí)際上細(xì)胞外的時(shí)變參數(shù)，如生物量濃度、產(chǎn)物濃度等均需要經(jīng)過(guò)人工取樣后離線檢測(cè)， 這導(dǎo)致在用DMFA估計(jì)代謝通量時(shí)會(huì)忽略重要的過(guò)程動(dòng)態(tài)信息，為使代謝通量更真實(shí)地反映微生物生長(zhǎng)代謝狀態(tài)的變化，需要實(shí)時(shí)預(yù)測(cè)所需的細(xì)胞外時(shí)變參數(shù)。 人工神經(jīng)網(wǎng)絡(luò)（Artificial Neural Network，ANN）在發(fā)酵過(guò)程軟測(cè)量中應(yīng)用廣泛，Dach J等使用ANN對(duì)漿液發(fā)酵過(guò)程進(jìn)行建模，較好地預(yù)測(cè)了甲烷的排放水平，但使用隨機(jī)初始化深層ANN參數(shù)的策略容易陷入局部最優(yōu)或梯度消失，難以保證模型的性能［8］。 隨著計(jì)算機(jī)硬件的發(fā)展，Hinton G E等提出的深度信念網(wǎng)絡(luò)（Deep Belief Network，DBN）策略能夠有效訓(xùn)練深層ANN［9］，即先逐層進(jìn)行無(wú)監(jiān)督的預(yù)訓(xùn)練來(lái)獲得非線性的潛變量模型，然后有監(jiān)督地微調(diào)最終的回歸或分類模型。 Erhan D等利用若干基準(zhǔn)實(shí)驗(yàn)說(shuō)明無(wú)監(jiān)督預(yù)訓(xùn)練的有效性［10］。DBN的半監(jiān)督學(xué)習(xí)策略不僅能挖掘過(guò)程中所有數(shù)據(jù)的信息，還可以充分?jǐn)M合發(fā)酵過(guò)程的非線性特征，因此可將DBN用于預(yù)測(cè)所需的生物學(xué)參數(shù)和建立代謝通量調(diào)控模型。

筆者考慮發(fā)酵過(guò)程的非線性特征和過(guò)程變量間采樣率不同的特點(diǎn)，提出基于DBN-DMFA的代謝通量調(diào)控模型策略， 其中DMFA用于估計(jì)代謝通量，DBN則用來(lái)預(yù)測(cè)生物學(xué)參數(shù)和建立代謝通量調(diào)控模型。 通過(guò)青霉素仿真對(duì)比實(shí)驗(yàn)，說(shuō)明該策略構(gòu)建的調(diào)控模型能夠有效揭示工業(yè)發(fā)酵的操作條件與微觀代謝通量間的聯(lián)系，可以用于進(jìn)一步的發(fā)酵控制和優(yōu)化。

1 青霉素代謝通量調(diào)控模型問(wèn)題描述

青霉素發(fā)酵過(guò)程中發(fā)酵液的pH值會(huì)改變菌體細(xì)胞膜的通透性［11］，而菌體的代謝活動(dòng)會(huì)影響氫離子的平衡，從而導(dǎo)致pH值波動(dòng)，因此需要實(shí)時(shí)控制酸液或堿液的流加速率以使pH值保持在合適的范圍內(nèi)。 溫度會(huì)影響酶的活性，而菌體的生長(zhǎng)代謝活動(dòng)、攪拌器的工作等都會(huì)引起熱能的變化，因此要實(shí)時(shí)調(diào)節(jié)冷水或熱水的流量來(lái)保持發(fā)酵環(huán)境的溫度。 另外，溶氧濃度影響菌體的代謝途徑和產(chǎn)物產(chǎn)量，通過(guò)實(shí)時(shí)調(diào)整攪拌功率和無(wú)菌空氣流速可以滿足菌體在不同發(fā)酵階段的需氧量。

青霉素代謝通量調(diào)控模型中包含的化學(xué)信息可以用化學(xué)計(jì)量矩陣S（S∈R（M+N）×L）來(lái)表示，其中M、N分別是菌體細(xì)胞內(nèi)、 細(xì)胞外代謝物的數(shù)量，L是代謝反應(yīng)的數(shù)量，矩陣中的數(shù)值對(duì)應(yīng)于相應(yīng)反應(yīng)式的系數(shù)［12］。 使用DMFA來(lái)估計(jì)青霉素的代謝通量時(shí)，還需要生物量濃度c（bio）、底物濃度c（glu）、青霉素濃度c（peni）、溶氧濃度c（O2）和二氧化碳濃度c（CO2）的時(shí)變數(shù)據(jù)，而實(shí)際發(fā)酵過(guò)程中c（O2）和c（CO2）可以在線測(cè)得，Eun=［c（bio）c（glu） c（peni）］則需要進(jìn)行離線檢測(cè)，為使所估計(jì)的代謝通量真實(shí)反映過(guò)程的動(dòng)態(tài)，可以建立相應(yīng)的軟測(cè)量模型來(lái)實(shí)時(shí)獲得Eun。

2 基于DBN-DMFA的代謝通量調(diào)控模型

2.1 深度信念網(wǎng)絡(luò)

DBN 是由受限玻耳茲曼機(jī)（Restricted Boltzmann Machines，RBM）堆疊而成的。 RBM結(jié)構(gòu)如圖1所示，由可視層和隱含層組成，具有層內(nèi)無(wú)連接、層間全連接的特點(diǎn)。

圖1 RBM結(jié)構(gòu)

RBM的能量函數(shù)如下：

其中wij是節(jié)點(diǎn)Vj和Hi間的權(quán)重，bj和ci分別為節(jié)點(diǎn)Vj和Hi的偏置，vj和hi分別對(duì)應(yīng)節(jié)點(diǎn)Vj和Hi的狀態(tài)，根據(jù)能量函數(shù)可以定義各節(jié)點(diǎn)的概率。

RBM的訓(xùn)練目標(biāo)是擬合輸入樣本的分布［13］，目前常用對(duì)比散度方法來(lái)快速學(xué)習(xí)RBM，參數(shù)θ={w，b，c}更新規(guī)則如下：

其中v*代表可視層v的重構(gòu)，h*是根據(jù)v*得到的隱含層，p（）為概率，ε為學(xué)習(xí)率。

DBN的訓(xùn)練過(guò)程如圖2所示， 首先對(duì)網(wǎng)絡(luò)前n-1層采用逐層貪婪學(xué)習(xí)算法進(jìn)行無(wú)監(jiān)督預(yù)訓(xùn)練，即先訓(xùn)練RBM1，然后保存該層參數(shù)，將RBM1的隱含層輸出作為RBM2的輸入， 訓(xùn)練RBM2后保存參數(shù)，以此類推，一直到RBMn-1訓(xùn)練結(jié)束。 然后將保存的參數(shù)作為網(wǎng)絡(luò)前n-1層的初始參數(shù)，最后使用反向傳播對(duì)整體網(wǎng)絡(luò)進(jìn)行有監(jiān)督地微調(diào)。

圖2 DBN的訓(xùn)練過(guò)程

2.2 動(dòng)態(tài)代謝通量分析

基于胞內(nèi)擬穩(wěn)態(tài)假設(shè)和質(zhì)量守恒原則，得到動(dòng)態(tài)通量均衡方程：

其中S∈R（M+N）×L代表菌體代謝網(wǎng)絡(luò)中的化學(xué)信 息，Sint∈RM×L和Sext∈RN×L分 別 是 細(xì) 胞 內(nèi) 和 細(xì) 胞外代謝物的化學(xué)計(jì)量學(xué)矩陣，f∈RL×1代表L個(gè)反應(yīng)的代謝通量，cext∈RN×1是細(xì)胞外N個(gè)可觀測(cè)代謝物向量。

由于細(xì)胞內(nèi)代謝物數(shù)量M通常少于代謝反應(yīng)數(shù)量L，即矩陣Sint是欠定的，其自由度為d=Lrank（Sint），因此代謝通量f可以表示為：

其中K∈RS×d是Sint的零空間，u∈Rd×1為一組自由通量。

使用DMFA時(shí)刻對(duì)整個(gè)代謝過(guò)程進(jìn)行分段線性化，DMFA起點(diǎn)和終點(diǎn)時(shí)刻與實(shí)際測(cè)量的首末時(shí)刻保持相同， 假定每?jī)蓚€(gè)DMFA時(shí)刻間的代謝通量呈線性變化，可得：

其中ND是DMFA時(shí)刻的數(shù)量，γ（t，ti）是使得代謝通量線性化的系數(shù)矩陣，ti是第i個(gè)DMFA時(shí)刻區(qū)間，ui是ti內(nèi)的自由通量。 對(duì)式（6）兩側(cè)積分后，將細(xì)胞外代謝物的估計(jì)值與實(shí)際值間的方差最小化，即可估計(jì)出代謝通量f。

2.3 算法流程

基于DBN-DMFA的代謝通量調(diào)控模型流程如下：

a. 采集發(fā)酵過(guò)程樣本數(shù)據(jù)；

b. 選擇輔助變量并使用DBN構(gòu)建軟測(cè)量模型，預(yù)測(cè)所需的細(xì)胞外時(shí)變參數(shù)；

c. 結(jié)合菌體的代謝網(wǎng)絡(luò)和細(xì)胞外時(shí)變數(shù)據(jù)，利用DMFA計(jì)算動(dòng)態(tài)代謝通量；

d. 選擇操作條件和代謝通量的數(shù)據(jù)構(gòu)成樣本，利用DBN建立代謝通量調(diào)控模型。

3 青霉素發(fā)酵過(guò)程仿真分析

以青霉素發(fā)酵過(guò)程為例，過(guò)程數(shù)據(jù)來(lái)自Pensim仿真平臺(tái)，它以Birol機(jī)理模型為內(nèi)核，能夠有效地模擬真實(shí)發(fā)酵過(guò)程［14］。根據(jù)Pensim仿真平臺(tái)中各初始條件的范圍要求，隨機(jī)生成50批青霉素發(fā)酵過(guò)程數(shù)據(jù)，發(fā)酵總時(shí)長(zhǎng)為400 h，采樣時(shí)間為0.5 h。

3.1 數(shù)據(jù)預(yù)處理

青霉素發(fā)酵過(guò)程數(shù)據(jù)具有不同量綱，會(huì)導(dǎo)致算法收斂困難，因此筆者采用Z-score方法進(jìn)行歸一化，計(jì)算式為：

3.2 仿真分析

采用均方根誤差σRMSE和最大誤差絕對(duì)值σMAXE來(lái)量化模型的預(yù)測(cè)性能，表達(dá)式如下：

其中Ns是樣本數(shù)，Yi和Y^i分別是真實(shí)值和預(yù)測(cè)值，σRMSE和σMAXE值越小，說(shuō)明模型性能越好。

3.2.1 預(yù)測(cè)細(xì)胞外時(shí)變參數(shù)Eun

將50批原始數(shù)據(jù)按7∶2∶1劃分為訓(xùn)練集、驗(yàn)證集和測(cè)試集，然后隨機(jī)去掉訓(xùn)練集中40%數(shù)據(jù)點(diǎn)中的Eun向量，從而把訓(xùn)練集分為預(yù)訓(xùn)練集和微調(diào)集。

選取發(fā)酵罐溫度、pH值、底物流加速率、底物流加溫度、發(fā)酵液體積和發(fā)酵產(chǎn)生的熱量作為預(yù)測(cè)Eun的輔助變量， 通過(guò)實(shí)驗(yàn)確定網(wǎng)絡(luò)結(jié)構(gòu)為6-5-4-4-3， 將DBN與傳統(tǒng)ANN進(jìn)行比較，DBN和ANN預(yù)測(cè)性能見(jiàn)圖3和表1。

圖3 兩種模型對(duì)c（glu）的預(yù)測(cè)誤差

表1 測(cè)試集性能指標(biāo)對(duì)比

由仿真結(jié)果可知，DBN相比ANN能更好地預(yù)測(cè)Eun，有利于進(jìn)行動(dòng)態(tài)代謝通量分析。

3.2.2 計(jì)算動(dòng)態(tài)代謝通量青霉素代謝網(wǎng)絡(luò)［15］包含66個(gè)代謝反應(yīng)、49個(gè)胞內(nèi)代謝物和5個(gè)細(xì)胞外代謝物（L=66，M=49，N=5）。 基于DBN可對(duì)細(xì)胞外時(shí)變參數(shù)進(jìn)行有效預(yù)測(cè)，根據(jù)式（6）可估計(jì)出66個(gè)代謝通量，其中3個(gè)胞內(nèi)外交 換 反 應(yīng) 的 代 謝 通 量F=［f （glu） f （peni）f（O2）］與對(duì)應(yīng)的濃度E=［c（glu） c（peni） c（O2）］對(duì)比如圖4所示。

圖4 交換反應(yīng)代謝通量

以圖4中f（glu）為例，f（glu）的變化會(huì)引起c（glu）發(fā)生相應(yīng)的改變。由此可說(shuō)明DMFA所計(jì)算的代謝通量能夠有效地揭示菌體內(nèi)部的生長(zhǎng)代謝狀態(tài)。

3.2.3 代謝通量調(diào)控模型

基于對(duì)細(xì)胞外時(shí)變參數(shù)Eun的預(yù)測(cè)和對(duì)代謝通量的估計(jì)， 獲得建立調(diào)控模型所需的樣本數(shù)據(jù)。 操作條件變量（輸入）為通風(fēng)率、攪拌速率、底物流加速率、底物流加溫度、酸/堿流加速率、冷/熱水流加速率。 代謝通量變量（輸出）為f（peni）和f（glu）。 將樣本按7∶2∶1分為訓(xùn)練集、驗(yàn)證集和測(cè)試集， 由實(shí)驗(yàn)確定網(wǎng)絡(luò)結(jié)構(gòu)為8-7-6-6-4-2，ANN與DBN的預(yù)測(cè)性能見(jiàn)圖5與表2， 觀察可知DBN比ANN更能有效建立代謝通量調(diào)控模型。

圖5 兩種模型對(duì)f（peni）的預(yù)測(cè)誤差

表2 測(cè)試集性能指標(biāo)對(duì)比

4 結(jié)束語(yǔ)

操作條件對(duì)發(fā)酵生產(chǎn)效率影響顯著，代謝通量模型則可以結(jié)合微觀代謝機(jī)理來(lái)為操作條件的優(yōu)化提供指導(dǎo)。 考慮實(shí)際發(fā)酵中存在過(guò)程非線性和數(shù)據(jù)的多采樣率，筆者提出基于DBN-DMFA的建模策略，并利用Pensim仿真數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，說(shuō)明了該策略的有效性，所得模型可用于代謝層面的在線檢測(cè)、控制和優(yōu)化。