微生物連續發酵建模及胞內動力學參數辨識

劉重陽， 尹 蕾， 馮恩民， 修志龍

（1.大連理工大學數學科學學院，遼寧大連 116024；2.山東工商學院數學與信息科學學院，山東煙臺 264005；3.大連理工大學生命科學與技術學院，遼寧大連 116024）

0 引 言

1，3－丙二醇（1，3－PD）是一種重要的化工原料，可用來合成許多具有優良特性的聚合物，如聚酯和聚氨酯［1］.傳統的1，3－丙二醇化學合成法生產需要高溫、高壓及貴重催化劑才能實現，此法分離提純困難、成本高，極大地限制了1，3－PD的發展.隨著現代生物技術的發展，人們開始嘗試用微生物發酵法生產1，3－PD.微生物發酵法具有條件溫和、操作簡單、副產物少、綠色環保等優點，所以這方面的研究正受到國內外越來越多的重視.特別地，通過克雷伯氏桿菌（K.pneumoniae）發酵生產1，3－PD受到眾多學者的廣泛關注［2～5］.利用克雷伯氏桿菌進行甘油發酵是一個復雜的生化過程，因為微生物的生長受到底物、各種產物和中間代謝物的抑制作用［3］，其中3－羥基丙醛（3－HPA）是一種有毒的中間代謝物，它對生產1，3－丙二醇所必需的兩種重要的酶（甘油脫水酶（GDHt）和1，3丙二醇氧化還原酶（PDOR））的活性具有抑制作用.近年來，國內外許多研究者利用數學模型來描述甘油發酵過程以實現發酵過程的優化從而提高1，3－PD的產量.Zeng等［3］提出了關于底物消耗和產物形成的過量動力學模型.修志龍等［4］改進了過量動力學模型使其適用于更大范圍的甘油注入濃度并研究了連續發酵中出現的多穩態現象.但是，上述模型都沒有考慮胞內1，3－PD、胞內甘油以及3－HPA的濃度變化.然而，3－HPA對于1，3－PD的生成起著至關重要的作用，所以上述模型不能恰當地描述甘油發酵過程.

2008年，Sun等首次提出涉及胞內1，3－PD跨膜運輸及3－HPA的數學模型來描述甘油發酵生產1，3－PD過程［6］.但是該模型中的胞內動力學參數的確定是根據大腸桿菌和其他菌種類推得到的，而沒有利用實驗數據及發酵過程的信息.顯然這種方法存在一定的偏差.為了更好地描述甘油發酵生產1，3－PD這一生化過程，本文改進文獻［6］中的數學模型，以計算值與實驗數據之間的平均相對誤差作為優化目標，建立確定胞內動力學參數的參數辨識模型.參數辨識對于確定未知動力學方程，進一步解釋和定性地描述發酵行為是非常重要的.目前，關于參數辨識問題的求解已有許多算法，其中包括非線性最小二乘法［7］、單純形法［8、9］、遺傳算法［10］和差異演化算法［11］等.本文將粒子群算法和Hooke－Jeeves算法相結合構造改進的粒子群算法求解參數辨識模型.

1 連續發酵動力學系統

在文獻［6］的基礎上，本文建立如下描述甘油連續發酵過程的動力學系統：

式中：t∈［0，T］，T∈（0，＋∞）表示連續發酵達到穩態的時刻；x1（t），x2（t），…，x8（t）分別表示生物量、胞外甘油、胞外1，3－PD、乙醇、乙酸、胞內甘油、3－HPA和胞內1，3－PD在t時刻的濃度；x i（0）（i∈I8）為發酵的初始濃度；cs0為注入甘油濃度；D表示稀釋速率；Km1、Km2分別表示GDHt和PDOR的米氏常數，且取值分別為0.53、0.14 mmol／L.基于文獻［4］和［6］，細胞的比生長速率μ、比消耗速率p1、比生成速率p2和p3，及u1、u2分別由下式給出：

對于系統（8），本文以文獻［6］中的參數值作為辨識的初始值k0（文獻［6］中無k1，根據生物意義此處取值為100；為了保證系統達到穩態，參數k3、k5、k10、k11、k15的值作了相應的修改），定義參數向量的允許范圍為另外，根據實際發酵過程，系統（8）中向量v和x分別滿足條件：

下面給出系統（8）及其解的一些性質.

性質1 對任意的k∈K和v∈V，函數f（x（t），v，k）滿足f∈C（［0，T］，R8）且f在上關于x是局部Lipschitz連續的.

性質2 對任意的k∈K和v∈V，函數f（x（t），v，k）滿足線性增長條件，即存在常數α，β＞0使得

定理1 對任意的k∈K和v∈V，系統（8）存在唯一的解，記為x（·；v，k）.另外，x（·；v，k）在K上關于k是連續的.

證明 由性質1、2和常微分方程理論［12］知此結論成立. □

給定x0∈W，記系統（8）的解集為

由以上定理和定義，可得如下定理：

定理2 式（13）定義的可行參數向量集F是緊集.

2 參數辨識模型及優化方法

確定系統（8）中的胞內動力學參數對于明確甘油的代謝機理、甘油和1，3－PD跨膜運輸方式具有非常重要的意義.參數辨識模型通常沒有解析解，因此，構造恰當的數值優化方法來確定參數辨識模型的最優解是非常必要的.

2.1 參數辨識模型

對甘油歧化連續發酵過程，本文作如下假設.

H1：給定v∈V及任意的x0∈W，k∈F，系統（8）存在唯一的穩態.

在假設H1下，給定v j，j∈Il＝｛1，2，…，l｝，其中l表示實驗次數，已測得連續發酵達到穩態時生物量、胞外甘油、胞外1，3－PD、乙醇和乙酸的濃度分別為y j1、y j2、y j3、y j4、y j5.令y j＝（y j1y j2y j3y j4y j5）T∈，j∈Il.另外，當系統達到穩態時，系統（8）的解x（T j；v j，k）滿足

其中T j表示第j次實驗達到穩態的時刻.

在連續發酵過程中，由于堿的加入會對乙醇、乙酸的濃度產生影響，本文以生物量、胞外甘油和胞外1，3－PD三種物質的計算值與實驗值之間的平均相對誤差作為優化目標，即

這樣，以J（k）作為目標函數，本文建立如下估計連續發酵胞內動力學參數的辨識模型：

（PIM） minJ（k）

定理3 參數辨識模型（PIM）存在最優解k＊，即罷k＊∈F有下式成立：

證明 由定理1和2可證此結論成立. □

2.2 改進粒子群算法

粒子群算法（PSO）是由Kennedy和Eberhart于1995年提出的［13］.目前，PSO作為一種新的全局搜索算法已被廣泛用于智能算法、優化和其他領域［14～16］.與其他進化計算方法相比，PSO具有搜索能力強、收斂速度快、設置參數少、程序實現簡單等特點，既適合科學研究，又適合工程應用.

PSO源于對鳥群捕食行為的研究.當一群鳥在隨機搜索食物時，所有的鳥都不知道食物在哪里.但它們知道當前位置離食物還有多遠.那么找到食物最簡單有效的策略就是搜尋目前離食物最近的鳥的周圍區域.PSO就是從這種模型中得到啟發而產生，并用于解決優化問題的.PSO最初是處理連續優化問題的，目前，其應用已擴展到組合優化問題中.在PSO中，每個優化問題的解都對應于搜索空間中一只鳥的位置，這些鳥被稱之為“粒子”.所有的粒子都有一個由被優化的參數決定的適應值（目標函數），每個粒子還有一個速度，決定它們飛翔的方向和距離.在每一次迭代過程中，粒子通過跟蹤兩個最優解來更新自己的位置.一個是粒子本身所找到的最優解，即個體極值pbest；另一個是整個種群目前找到的最優解，即全局極值gbest.通過這種方式，粒子們就追隨當前的最優粒子在解空間中搜索并朝著好的區域移動.

下面給出求解（PIM）的改進粒子群算法具體步驟：

式中：分別表示第i個粒子迭代到第s步時的速度和位置的第n個分量，c1、c2是正常數，r1n、r2n、r3n是［0，1］上服從均勻分布的隨機數.如果出界，則重新執行Step 5操作直到其在界內為止.為了避免算法在迭代后期粒子出現震蕩的現象，ω（s）按下式變化：

3 數值結果

依實際發酵過程，初始濃度為x01＝0.115 g／L，x02＝495 mmol／L，x03＝x04＝x05＝x06＝x07＝x08＝0.選取連續發酵實驗測得的22組數據（其中底物過量情形的10組和底物限制情形的12組），用改進粒子群算法求解（PIM），得到最優的胞內動力學參數k＊如表1所示.特別地，本文得到最優參數對應的平均相對誤差為30.195%，而由參數k0計算所得平均相對誤差為50.830%，從而優化過程使平均相對誤差減小了20.635%.這里，在改進粒子群算法中參數m、c1、c2、Tmax、ε、ρHJ、ωmax和ωmin分別取值為50、2.0、2.0、500、0.001、0.3、0.9和0.4.另外，將改進的粒子群算法與文獻［9］和［11］中算法的收斂性進行了比較（如圖1所示）.從圖中可以看出，文獻［9］和［11］中的算法很容易收斂到局部最優解，而本文提出的算法具有較好的全局收斂性.

表1 系統（8）中最優胞內動力學參數Tab.1 The optimal intracellular kinetic parameters in the System（8）

圖1 3種算法的收斂性曲線Fig.1 Convergence curves of three algorithms

4 結 論

本文給出了描述連續發酵甘油生產1，3－PD過程的動力學系統，并討論了該系統的一些性質.為了精確估計模型中胞內動力學參數，建立了參數辨識模型，并證明了最優參數的存在性.最后，構造了改進粒子群算法求解參數辨識模型.數值結果表明：計算值與實驗數據的平均相對誤差減小了20.635%，且該算法比其他算法具有更好的全局收斂性.

［1］修志龍.1，3－丙二醇的微生物法生產分析［J］.現代化工，1999（3）：33－35

［2］BIEBL H，MENZEL K，ZENG An－ping，etal.Microbial production of 1，3－propanediol［J］.Applied Microbiology and Biotechnology，1999，52（3）：289－297

［3］ZENG An－ping，ROSE A，BIEBL H，etal.Multiple product inhibition and growth modeling ofClostridiumbutyricumandKiebsiellapneumoniaein glycerol fermentation［J］.Biotechnology and Bioengineering，1994，44（8）：902－911

［4］修志龍，曾安平，安利佳.甘油生物歧化為1，3－丙二醇過程的動力學數學模擬和多穩態研究［J］.大連理工大學學報，2000，40（4）：428－433

（XIU Zhi－long，ZENG An－ping，AN Li－jia.Mathematical modeling of kinetics and research on multiplicity of glycerol bioconversion to 1，3－propanediol［J］.Journal of Dalian University of Technology，2000，40（4）：428－433）

［5］高彩霞，馮恩民，王宗濤，等.微生物間歇發酵生產1，3－丙二醇過程辨識與優化［J］.大連理工大學學報，2006，46（5）：771－774

（GAO Cai－xia，FENG En－min，WANG Zong－tao，etal.Parameter identification and optimization of process for bio－dissimilation of glycerol to 1，3－propanediol in batch culture［J］.Journal of Dalian University of Technology，2006，46（5）：771－774）

［6］SUN Y Q，QI W T，TENG H，etal.Mathematical modeling of glycerol fermentation byKlebsiella pneumoniae：Concerning enzyme－catalytic reductive pathway and transport of glycerol and 1，3－propanediol across cell membrane［J］.Biochemical Engineering Journal，2008，38（1）：22－32

［7］CHEN Q H，HE G Q，SCHWARZ P.Studies on cultivation kinetics for elastase production byBacillus sp.EL31410［J］.Journal of Agricultural and Food Chemistry，2004，52（11）：3356－3359

［8］FENG Y Y，HE Z M，SONG L F，etal.Kinetics of β－mannanase fermentation by Bacillus licheniformis［J］.Biotechnology Letters，2003，25（14）：1143－1146

［9］TAXARWS A P M，COELHO M A Z，COUTINHO J A P，etal.Laccase improvement in submerged cultivation：induced production and kinetic modeling［J］.Journal of Chemical Technology and Biotechnology，2005，80：669－676

［10］PINCHUK R J，BROWN W A，HUGHES S M，etal.Modeling of biological processes using selfcycling fermentation and genetic algorithms［J］.Biotechnology and Bioengineering，2000，67（1）：19－24

［11］WANG F S，SU T L，JANG H J.Hybrid differential evolution for problems of kinetic parameter estimation and dynamic optimization of an ethanol fermentation process［J］.Industrial ＆Engineering Chemistry Research，2001，40（13）：2876－2885

［12］龐特里亞金ЛC.常微分方程［M］.6版.林武中，倪明康，譯.北京：高等教育出版社，2006

［13］KENNEDY J，EBERHART R.Particle swarm optimization［C］／／Proceedings of the IEEE International Conference on Neural Networks.Piscataway：IEEE，1995：1942－1948

［14］CLERC M，KENNEDY J.The particle swarm：explosion，stability，and convergence in a multidimensional complex space［J］.IEEE Transactions on Evolutionary Computation，2002，6（1）：58－73

［15］YU J B，XI L F，WANG S J.An improved particle swarm optimization for evolving feedforward artificial neural networks［J］.Neural Process Letter，2007，26（3）：217－231

［16］LI H Q，LI L，KIM T H，etal.An improved PSO－based of harmony search for complicated optimization problems［J］.International Journal of Hybrid Information Technology，2008，1（1）：57－64