釀酒酵母代謝過程的振蕩分析

潘多濤，史洪巖，袁德成，修志龍

?

釀酒酵母代謝過程的振蕩分析

潘多濤1,2，史洪巖2，袁德成2，修志龍1

（1大連理工大學生命科學與技術學院，遼寧大連116024；2沈陽化工大學遼寧省化工過程控制技術重點實驗室，遼寧沈陽 110142)

振蕩現象是生物的固有特征，在不同生物動態過程中都起著重要作用。為探究生物代謝過程中振蕩產生的條件，以典型的釀酒酵母發酵生產乙醇過程為研究對象，采用函數連續性分析方法對其數學模型進行深入研究。首先對系統進行仿真，結合相圖確認該過程存在等幅持續振蕩現象（極限環振蕩）；以此為切入點，利用分叉分析方法考察模型參數對系統振蕩現象的影響；最后根據結果分析系統產生振蕩現象的條件。結果顯示，該代謝途徑中各步反應均有參數會對系統的振蕩現象產生影響，得出產生極限環振蕩的參數范圍，同時根據結果分析，為今后抑制振蕩或利用振蕩有利特性提供指導。

系統生物學；糖酵解；振蕩；極限環；分支分析；非線性規劃

引 言

系統生物學是通過建模分析，結合理論實踐，系統研究生物過程的學科，其首要目標是在現有實驗數據的基礎上，建立一個能夠準確預測生物過程的模型。隨著生物學研究進入“后基因組”時代，其關鍵的問題是如何解決生物過程網絡的復雜性，如生物系統經過長期進化形成的魯棒性、高度非線性以及多回路反饋等特點[1-2]。

生物系統中普遍存在的振蕩現象是生物過程復雜性的另一種表現形式。生物振蕩是生物體不斷與外界發生物質、能量交換，在特定階段會處于偏離穩態，發生非線性、非平衡的狀態周期變化的現象[3]。由外部因素引發的振蕩稱之為強迫振蕩，而源于內因的稱之為自發性振蕩。國內外學者已對生物體系各方面的振蕩現象做了大量研究，如Longo等[4]考察了誘導因子信號轉導網絡的振蕩行為。而在代謝過程方面，特別是在連續發酵過程中的振蕩行為，越來越受到重視。Williamson等[5]研究了酵母的基因控制對代謝振蕩過程的影響；Gustavsson等[6]首次觀察到游離酵母細胞的振蕩現象；Preez等[7]則對已有的糖酵解振蕩模型的參數可適范圍進行優化，并校正模型。由于振蕩過程的復雜性，僅通過數學模型很難直觀地闡明其動態特性，更為有效的方法是通過調整模型參數來分析系統的動態行為。

本文針對釀酒酵母在連續培養中的振蕩過程，在模型仿真分析的基礎上，利用非線性理論的方法對關鍵參數進行分支分析，以期通過對系統等幅持續振蕩[8]的特性分析能夠進一步為實驗研究提供指導依據。

1 振蕩過程

振蕩是生物體固有的重要性質之一，在體系內的不同水平層次都有所體現，其宏觀表現形式一般被稱為節律[9]，如心肺張弛、生殖周期以及不同表現形式的生物鐘現象。而在生物細胞內的代謝過程中出現的振蕩現象，通常是代謝物濃度出現周期性的波動。

酵母發酵產生乙醇的過程是指以葡萄糖為代表的六碳糖共同經歷的分解代謝生成丙酮酸并進一步轉化為乙醇的過程。該過程中包含的糖酵解過程是生物界最原始獲取能量（ATP）的一種方式，是生物體共同經歷的葡萄糖分解代謝的前期途徑，被認為是大多數生物細胞代謝過程中最重要的途徑。在整個過程中，ATP-ADP以及NAD-NADH形成反饋回路，在一定條件下該過程就會表現出復雜的非線性動力學特性[10]。在酵母發酵過程中，早在1957年Duysens等[11]就首次觀察到完整細胞中NADH的濃度波動現象；之后陸續有人報道在無細胞的酵母提取物中存在糖酵解振蕩現象[12]；在糖酵解途徑中，磷酸果糖激酶催化的代謝物出現振蕩現象[13]。

2 方法

2.1 模型

生物代謝過程的數學模型一般采用微分方程（ordinary differential equation, ODE）來描述。以代謝物濃度為狀態變量的一組微分方程的向量形式如下

其中，=[12,…,x]T∈為系統狀態變量；=[1,2,…,k]∈為模型參數。式（1）描述了代謝物濃度隨時間的變化規律，一般以Monod或Michaelis-Menten方程描述其關系。為了盡可能提高模型的精確度，方程的形式越發復雜；代謝網絡的拓展也使需要考察的代謝物增多，對應的系統模型維度也隨之不斷提高。這些因素導致模型整體的復雜度顯著增加，需要考慮更多的因素與非線性行為的關系。

長期以來，人們試圖借助數學工具來研究發酵過程的振蕩機制，自20世紀70年代以來已有多個針對酵母振蕩過程的數學模型見諸報道[14-16]。為降低模型復雜度，本文采用Wolf等[14]給出的葡萄糖經酵母代謝生成乙醇的途徑，描述了與實驗觀測一致的酵母糖酵解途徑的最簡核心模型，途徑中包含乙酸在細胞間交換的過程。如圖1所示。

為了便于研究，圖1給出的代謝過程是一種簡化途徑，其中虛線框代表細胞膜，0和ex分別代表底物葡萄糖和代謝物乙酸在細胞間交換的流量。

Wolf等[14]依據圖1的代謝過程所建立的數學模型見式（2），其中相關參數值取自文獻[14]。這里，為了保證后續求解模型過程的順利，將NADH ? NAD、ATP ? ADP這樣的“二元環”方程用代數方程改寫并求解。

為了數學表達式描述的方便，令x=[Glc, F16p, G3P, GP, Pyr, Acd, ExAcd, ATP, NAD]。

其中

2.2 仿真分析

通過動態仿真分析，可初步了解系統的動態行為，是接下來的參數分析的切入點。應用前期工作中開發的集成優化工具GIEPT[17]，利用動態優化的方法進行求解，方法見文獻[18-20]。

2.3 參數分支分析

穩定的系統隨著一些參數連續變化，經過某臨界值時，狀態會發生改變：一種情況是出現多穩態現象；另一種情況是由穩態變為振蕩。這種由參數的變化導致的系統狀態躍遷的現象被稱為分支[21-22]。通過分支分析，能夠分析非線性系統中參數對狀態穩定性的影響。

首先只考慮一個二階自治系統[23]，如式（3）所示。系統的平衡點可以由系統方程為零求得，在平衡點處，系統的動態特性會有不同的表現。

根據系統連續性理論可知，系統平衡點處的Jacobi矩陣為式（4），可有助于分析其特性。

(4)

系統的特征值和特征向量由式（5）給出。

=v(5)

一般來說，判斷特征值的屬性，可確定平衡點的特性。除此之外，還可通過的跡式（5）和行列式值式（6）來確定相關性質。

(7)

當det()>0，且tr()<0時，平衡點是穩定的；當det()>0，且tr()>0時，平衡點是不穩定的。當特征值由全負實部變為零實部的情形時，通常是系統穩定性發生變化的邊界點，在穩態譜圖中即表現為分支現象。若系統具有一對純虛數特征值且其他特征值具有負實部，或det()=0，且tr(=0，則為Hopf分支，這一般表示系統由穩定平衡點過度到極限環振蕩的臨界點[21]。

在二階系統基礎上，當系統模型的維度較高、同時包含較多的參數時，分析過程計算較為復雜[24]，分支點的數量和類型也更多，甚至會產生混沌現象[25]。

3 結果與分析

生化系統因其自身的復雜性蘊含著大量的非線性特性，振蕩是非線性動力系統和非線性理論研究的重要內容，對微生物培養過程中出現的振蕩現象進行研究有助于加深對生命活動規律的進一步認識。

3.1 動態過程仿真

使用GIEPT調用非線性規劃解題器進行求解，結果如2所示。

圖中上段曲線為1,6-二磷酸果糖F16p的曲線，下段為葡萄糖Glu的變化情況（此處略去其他代謝物的計算結果）。需要說明的是，為了使圖示效果更清晰可辨，圖2中的數據點每間隔20個繪制1次。

從圖2的動態仿真結果可以看出，該系統具有明顯的等幅振蕩特性，除圖2中的兩種代謝物外，系統式（2）中其他代謝物均產生振幅不同的振蕩現象。

對于此類周期等幅振蕩現象，一般亦稱為極限環振蕩。考察代謝物三磷酸腺苷（ATP）與輔酶I（NAD）的濃度變化關系，在二維平面中繪制兩者關系，結果如圖3所示。

結果顯示，糖酵解系統出現典型的極限環振蕩現象。由于NAD為糖酵解過程提供氧化還原當量，而ATP則為代謝過程提供能量，因此這兩種物質的變化情況與生物過程有著緊密的聯系。接下來著重分析NAD與ATP的模型參數與極限環振蕩的關系。

3.2 分支分析

首先對ATP降解為ADP的速率常數15進行參數穩態分析，使15參數值在15～45 min-1范圍內變化同時計算系統特征值，在15分別為24.89 min-1及43.78 min-1時出現兩個系統特征值變化的臨界點，判斷為Hopf分支點。其中在15=24.89 min-1處，Hopf分支點的第一李雅普諾夫系數為正值，意味著此分支點為次臨界，同時此處存在一個不穩定的極限環；而在15=43.78 min-1處分支點的第一李雅普諾夫系數為負值，這表示該分支點為超臨界且存在的極限環是穩定的。結果如圖4所示。

圖4中黑色實線為15變化的系統穩態軌跡。15在兩處Hopf分支點間取值時系統出現極限環振蕩。圖中藍色區域為系統極限環隨著參數變化，其周期、軌跡形態的變化過程，特別是在15=23.61min-1處，環軌跡出現極限分支點，這表明環軌跡在此處出現折疊，在圖4中表現為自次臨界Hopf分支點（左）開始產生極限環，在超臨界Hopf分支點（右）處結束。

由于系統涉及到乙酸在細胞間的交換過程，因此考慮乙酸降解速率常數16對NAD與ATP濃度波動的影響。結果如圖5所示。使16在大于10 min-1的范圍變化，同時計算系統特征值，16在27.03 min-1處出現一個次臨界Hopf分支點（第一李雅普諾夫系數為正值），在16=10.93 min-1處出現環軌跡極限分支點。與15不同，16在經過27.03 min-1處的Hopf分支點后，未檢測到新的分支點，而在環軌跡上表現為振幅變化不大、周期逐步放大的不穩定狀態。

在圖5中，乙酸降解速率16決定了系統能否產生極限環振蕩。只要16超過27.03 min-1，系統就將產生振蕩，并且在這之后系統振蕩的振幅、周期與參數16的增長呈弱相關性。由此可以推測，乙酸的降解量可能使得胞內代謝物產生振蕩，這為通過產物控制代謝過程的振蕩提供了依據。

限于篇幅，其他參數的分支分析結果并未列出。通過對各個參數的分支分析，可以獲知能夠影響系統振蕩產生的參數取值范圍。同時，底物流加速率也是極限環振蕩能否產生的影響因素。

3.3 驗證

經過對參數的分支分析，能夠獲得系統振蕩的產生條件，通過調節參數，可獲得穩定的代謝過程。然而，生物代謝過程中的振蕩現象分為兩類：自發振蕩和強制振蕩[26]。自發振蕩無須外部周期性變化因素來激勵，其維持動力源于發酵體系內部矛盾相互轉化；而強制振蕩一般受外部條件維持，如底物的消耗、產物的積累。強制振蕩較容易調控[27-28]。

在該模型的16個參數中，大多數參數為胞內反應速率常數，而16是胞外乙酸的降解速率常數，該參數值容易通過調節溫度等操作來改變。接下來，在分支分析的基礎上，調節16參數值，驗證系統的穩定性。由于16只有一個Hopf分支點，取其鄰近值，從仿真結果上可以看出，系統趨于穩定（如圖6所示，這里只列出ATP、NAD的變化情況）；當進一步減小16時，結果會出現發散的現象（圖略）。

有觀點[29]指出，產物乙醇的反饋抑制作用，是產生振蕩的一個內因。通過以上分析，乙酸的降解速率很有可能是產生振蕩的另一因素。

綜上，借助數學工具對模型進行參數分支分析，一方面，在連續培養過程中，可以利用這些信息為降低振蕩的不利影響提供指導；另一方面，可以利用振蕩產生后有利的影響，如利用產物濃度周期振蕩的特性，可以設計優化的產品回收及分離工藝。

4 結 論

以酵母發酵生產乙醇為代表的生物過程是復雜并高度集成化的系統，難以借助傳統的實驗手段高效準確地加以分析。對這類復雜的生物系統過程的研究，客觀上需要系統生物學和生物信息學提供方法上的支持。

針對典型的釀酒酵母發酵生成乙醇的振蕩過程，對其數學模型進行了深入分析，獲得有意義的結果如下。

（1）對過程模型進行仿真，并通過狀態關系分析確定該過程能夠產生極限環振蕩。

（2）對模型參數進行分支分析，其中15具有一個超臨界Hopf分支點和一個次臨界Hopf分支點，這是能夠引起極限環振蕩的范圍。

（3）根據分支點臨界值的分析結果，調整16參數，可使系統趨于穩定，提出乙酸降解速率可能是引起振蕩行為的一個重要因素。

從數學模型參數出發，考察其如何影響振蕩的產生，是后續研究振蕩機制的基礎。

References

[1] HEINZLE E, BIWER A P, CHARLES L C. Development of Sustainable Bioprocesses: Modeling and Assessment[M]. Sussex: Wiley, 2007.

[2] KARSENTI E. Self-organization in cell biology: a brief history[J]. Nat. Rev. Mol. Cell Biol., 2008, 9(3): 255-262.

[3] KRUSE K, JULICHER F. Oscillations in cell biology[J]. Curr. Opin. Cell Biol., 2005, 17(1): 20-26.

[4] LONGO D M, SELIMKHANOV J, KEARNS J D,. Dual delayed feedback provides sensitivity and robustness to the NF-κB signaling module [J]. PLoS Comput. Biol., 2013, 9(6): e1003112.

[5] WILLIAMSON T, ADIAMAH D, SCHWARTZ J M,. Exploring the genetic control of glycolytic oscillations in[J]. BMC Syst. Biol., 2012, 6: 108.

[6] GUSTAVSSON A K, VAN NIEKERK D D, ADIELS C B,. Heterogeneity of glycolytic oscillatory behaviour in individual yeast cells[J]. FEBS Lett., 2014, 588(1): 3-7.

[7] DU PREEZ F B, VAN NIEKERK D D, KOOI B,. From steady-state to synchronized yeast glycolytic oscillations (Ⅰ): Model construction[J]. FEBS J., 2012, 279(16): 2810-2822.

[8] GUSTAVSSON A K, VAN NIEKERK D D, ADIELS C B,. Sustained glycolytic oscillations in individual isolated yeast cells[J]. FEBS J., 2012, 279(16): 2837-2847.

[9] GOLDBETER A, GERARD C, GONZE D,. Systems biology of cellular rhythms[J]. FEBS Lett., 2012, 586(18): 2955-2965.

[10] YAMAZAKI S, MIKI K, KANO K,. Mechanistic study on the role of the NAD+-NADH ratio in the glycolytic oscillation with a pyruvate sensor[J]. Journal of Electroanalytical Chemistry, 2001, 516(1/2): 59-65.

[11] DUYSENS L N, AMESZ J. Fluorescence spectrophotometry of reduced phosphopyridine nucleotide in intact cells in the near-ultraviolet and visible region[J]. Biochim. Biophys. Acta, 1957, 24(1): 19-26.

[12] CHANCE B, ESTABROOK R W, GHOSH A. Damped sinusoidal oscillations of cytoplasmic reduced pyridine nucleotide in yeast cells[J]. Proc. Natl. Acad. Sci. USA, 1964, 51: 1244-1251.

[13] RICHARD P. The rhythm of yeast[J]. Fems. Microbiol. Rev., 2003, 27(4): 547-557.

[14] WOLF J, PASSARGE J, SOMSEN O J,. Transduction of intracellular and intercellular dynamics in yeast glycolytic oscillations[J]. Biophys. J., 2000, 78(3): 1145-1153.

[15] WOLF J, HEINRICH R. Effect of cellular interaction on glycolytic oscillations in yeast: a theoretical investigation[J]. Biochem. J., 2000, 345: 321-334.

[16] DU PREEZ F B, VAN NIEKERK D D, SNOEP J L. From steady-state to synchronized yeast glycolytic oscillations(Ⅱ): Model validation[J]. FEBS J., 2012, 279(16): 2823-2836.

[17] 潘多濤, 黃明忠, 張學軍, 等. 工程規劃建模研究及流程優化的實現[J]. 北京工業大學學報, 2012, 38(10): 1486-1490. PAN D T, HUANG M Z, ZHANG X J,. Engineering modeling programming and implementation of process optimization[J]. Journal of Beijing University of Technology, 2012, 38(10): 1486-1490.

[18] SHIVAKUMAR K, BIEGLER L T. Simultaneous dynamic optimization strategies: recent advances and challenges [J]. Computers and Chemical Engineering, 2006, 30(10/11/12): 1560-1575.

[19] 潘多濤, 史洪巖, 黃明忠, 等. 復雜生物過程的仿真分析[J]. 系統仿真學報, 2014, 26(3): 670-674. PAN D T, SHI H Y, HUANG M Z,. Simulation and analysis of complex biological processes [J]. Journal of System Simulation, 2014, 26(3): 670-674.

[20] 潘多濤, 史洪巖, 黃明忠, 等. 非線性規劃在生物代謝仿真過程中的應用[J]. 控制工程, 2014, (6): 896-899. PAN D T, SHI H Y, HUANG M Z,. Application of nonlinear programming for simulation of biological metabolic process [J]. Control Engineering of China, 2014, (6): 896-899.

[21] CRAWFORD J D. Introduction to bifurcation theory[J]. Reviews of Modern Physics, 1991, 63(4): 991-1037.

[22] KUZNETSOV Y A. Elements of Applied Bifurcation Theory[M]. New York: Springer-Verlag, 2004.

[23] GOVAERTS W. Numerical bifurcation analysis for ODEs[J]. Journal of Computational and Applied Mathematics, 2000, 125(1): 57-68.

[24] DHOOGE A, GOVAERTS W, KUZNETSOV Y A. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs[J]. ACM Transactions on Mathematical Software, 2003, 29(2): 141-164.

[25] DENG B. Food chain chaos due to junction-fold point[J]. Chaos, 2001, 11(3): 514-525.

[26] 申渝. 酵母細胞超高濃度乙醇連續發酵振蕩行為的研究[D]. 大連: 大連理工大學, 2009. SHEN Y. Exploration for oscillation in continuous VHG ethanol fermentation with[D]. Dalian: Dalian University of Technology, 2009.

[27] 申渝, 葛旭萌, 李寧, 等. 高濃度乙醇連續發酵振蕩過程中代謝通量分析及誘發機理[J]. 化工學報, 2009, 60(6): 1519-1528. SHEN Y, GE X M, LI N,. Metabolic flux analysis and mechanistic study of process oscillation in continuous VHG ethanol fermentation with[J]. CIESC Journal, 2009, 60(6): 1519-1528.

[28] 楊蕾, 陳麗杰, 白鳳武. 高濃度酒精連續發酵過程中振蕩行為的模擬及填料弱化振蕩的機理[J]. 化工學報, 2007, 58(3): 715-721. YANG L, CHEN L J, BAI F W. Dynamic models of VHG continuous ethanol fermentation and mechanisms of oscillation attenuation by packing[J]. Journal of Chemical Industry and Engineering(China), 2007, 58(3): 715-721.

[29] PORRO D, MARTEGANI E, RANZI B M,. Oscillations in continuous cultures of budding yeast: a segregated parameter analysis[J]. Biotechnol. Bioeng., 1988, 32(4): 411-417.

Analysis of metabolic oscillation processes in

PAN Duotao1,2, SHI Hongyan2, YUAN Decheng2,XIU Zhilong1

(1School of Life Science and Biotechnology, Dalian University of Technology, Dalian 116024, Liaoning, China;2Chemical Control Technology Key Laboratory of Liaoning Province, Shenyang University of Chemical Technology, Shenyang 110142, Liaoning, China)

Oscillation phenomenon is an inherent characteristic in biological systems and plays an important role in many dynamic bioprocesses. In order to explore the certain conditons that could possibly boost a oscillation, the metabolic pathway of the, glycolysis were researched, and the parameters of mathmatical model was analyzed. Firstly, the simulation results associated the phase diagrams indicated that the model exists sustained oscillations with constant amplitude (limit cycle oscillations). Next, bifurcation analysis approach was used to investigate the infulence of parameters for the system producing oscillations. The results showed that several parameters of the model could lead to oscillations and the range of parameters’ value was obtained, which could be applied to direct the manipulation of metabolic oscillations.

systems biology; glycolysis; oscillation; limit cycle; bifurcation analysis; nonlinear programming

10.11949/j.issn.0438-1157.20161628

TQ 920. 1

A

0438—1157（2017）03—0964—06

國家自然科學基金項目（21476042）；遼寧省教育廳科研一般項目（L2014168）；遼寧省博士科研啟動基金項目（201501072）。

2016-11-16收到初稿，2016-11-26收到修改稿。

聯系人：修志龍。第一作者：潘多濤（1979—），男，博士，講師。

2016-11-16.

Prof.XIU Zhilong, Zhlxiu@dlut.edu.cn

supported by the National Natural Science Foundation of China (21476042), the Educational Commission of Liaoning Province (L2014168) and the Doctoral Research Fund of Liaoning Province(201501072).