大腸桿菌蘇氨酸合成途徑動(dòng)力學(xué)模型的構(gòu)建與分析

楊雪，張彥飛，鄭陽陽，馬紅武

中國科學(xué)院天津工業(yè)生物技術(shù)研究所 中國科學(xué)院系統(tǒng)微生物工程重點(diǎn)實(shí)驗(yàn)室，天津 300308

伴隨著組學(xué)數(shù)據(jù)的大量獲得，各種組學(xué)尺度的代謝網(wǎng)絡(luò)模型不斷發(fā)展和完善。利用通量平衡分析等算法可以由代謝網(wǎng)絡(luò)模型求得某一產(chǎn)物生成的最優(yōu)途徑。獲得的結(jié)果在代謝途徑設(shè)計(jì)、敲除靶點(diǎn)預(yù)測(cè)等方面都具有重要指導(dǎo)作用[1]。與基因組規(guī)模代謝網(wǎng)絡(luò)模型相對(duì)應(yīng)，動(dòng)力學(xué)模型由于需要較多的動(dòng)力學(xué)參數(shù)數(shù)據(jù)一般只針對(duì)包含較少反應(yīng)的代謝途徑。但動(dòng)力學(xué)模型具有更強(qiáng)的預(yù)測(cè)能力，可以分析在各種擾動(dòng)下途徑速率的變化，通過模型模擬和控制分析更準(zhǔn)確地預(yù)測(cè)敲除和擴(kuò)增靶點(diǎn)，從而為高效細(xì)胞工廠的理性設(shè)計(jì)提供指導(dǎo)[2]。目前人們進(jìn)行動(dòng)力學(xué)模型構(gòu)建主要是基于已有的生化知識(shí)來選擇在模型中包括哪些反應(yīng)，而由基因組尺度代謝網(wǎng)絡(luò)模型求得的最優(yōu)途徑也可以作為動(dòng)力學(xué)模型構(gòu)建的出發(fā)點(diǎn)。與傳統(tǒng)方法相比，這種方法得到的模型自身就是滿足能量和還原力平衡的，因此不需要人為設(shè)定一些代謝物如ATP、ADP、NADH等的濃度。本文中即采用這種新的方法確定了大腸桿菌中由葡萄糖出發(fā)合成蘇氨酸的代謝途徑的動(dòng)力學(xué)模型的初始反應(yīng)集，進(jìn)而結(jié)合以前發(fā)表的相關(guān)模型中的動(dòng)力學(xué)數(shù)據(jù)確定各反應(yīng)的動(dòng)力學(xué)方程和參數(shù)，最終得到了一個(gè)完整的能量和還原力平衡的蘇氨酸合成途徑動(dòng)力學(xué)模型。

作為鳥類和哺乳類動(dòng)物的必需氨基酸之一，蘇氨酸對(duì)維持人類與動(dòng)物的營養(yǎng)和健康具有重要意義，被廣泛應(yīng)用于食品、飼料、藥品及化工領(lǐng)域，全球需求量逐年增長[3]。目前，蘇氨酸主要通過大腸桿菌等細(xì)菌的生物發(fā)酵制備[4]。為選育蘇氨酸高產(chǎn)菌株提供指導(dǎo)是人們構(gòu)建蘇氨酸合成途徑動(dòng)力學(xué)模型的主要目的。Chassagnole等曾提出一個(gè)基于實(shí)驗(yàn)數(shù)據(jù)確立的蘇氨酸合成動(dòng)力學(xué)模型，通過模型對(duì)蘇氨酸合成途徑中各步反應(yīng)對(duì)蘇氨酸合成通量的影響進(jìn)行了深入分析。但該模型中僅包括從天門冬氨酸到蘇氨酸的5步反應(yīng)，過于簡單[5-6]。2009年，該實(shí)驗(yàn)室進(jìn)一步將該模型發(fā)展為包括了從葡萄糖出發(fā)到蘇氨酸的完整合成途徑，并結(jié)合代謝途徑的通量分布的測(cè)量結(jié)果預(yù)測(cè)和驗(yàn)證了丙酮酸激酶敲除對(duì)蘇氨酸合成的促進(jìn)作用[7]。蘇氨酸合成是一個(gè)高還原力和能量需求的過程，實(shí)際蘇氨酸菌種改造過程中常常涉及到包括中心代謝內(nèi)的其他代謝途徑以滿足還原力和能量平衡的需求，但在已發(fā)表的模型中均未考慮能量和還原力的供給而是直接固定ATP、NADH等的濃度，這樣就會(huì)使得計(jì)算結(jié)果的準(zhǔn)確度受到影響。本文針對(duì)這一問題，以基因組規(guī)模代謝網(wǎng)絡(luò)分析求得的最優(yōu)途徑為出發(fā)點(diǎn)來構(gòu)建動(dòng)力學(xué)模型，不但考慮蘇氨酸合成的碳源需求，還考慮了其合成過程中消耗的ATP、還原力等對(duì)碳源的需求，可以使計(jì)算結(jié)果更為可靠。

1 動(dòng)力學(xué)模型構(gòu)建

1.1 動(dòng)力學(xué)模型結(jié)構(gòu)的確定

我們用大腸桿菌的基因組規(guī)模代謝網(wǎng)絡(luò)模型iJO1366[8]計(jì)算了從葡萄糖出發(fā)合成蘇氨酸的最優(yōu)代謝途徑。發(fā)現(xiàn)最優(yōu)途徑中糖酵解主要通過磷酸戊糖途徑進(jìn)行，但考慮到實(shí)驗(yàn)測(cè)得的通量分布中以EMP途徑為主[7]，我們構(gòu)建的動(dòng)力學(xué)模型中將EMP途徑和磷酸戊糖途徑都包括了進(jìn)來，整體代謝途徑如圖 1所示。為了構(gòu)建該完整代謝途徑的動(dòng)力學(xué)模型，我們主要參考了兩個(gè)已發(fā)表的模型來確定動(dòng)力學(xué)方程和參數(shù)，分別為從葡萄糖到丙酮酸的中心代謝途徑動(dòng)力學(xué)模型 (圖 1中細(xì)實(shí)線部分)[9]和從天冬氨酸到蘇氨酸合成途徑的模型(圖 1中虛線框內(nèi)的反應(yīng))[5-6]。以上兩個(gè)模型已被BioModels Database收錄[10]，模型編號(hào)分別為051和066，下文簡稱為模型a和模型b。其中模型a模擬的是細(xì)胞生長過程，其中包括了大量用于生長的稀釋反應(yīng)及用于生物質(zhì)生成的合并反應(yīng)。因?yàn)槲覀兘５哪康氖悄M能量和還原力平衡的蘇氨酸合成過程，所以去掉了這些與生長相關(guān)的反應(yīng)。同樣的，模型 b中包含的兩個(gè)與生長過程相關(guān)的NADPH和ATP的內(nèi)源消耗反應(yīng)也去掉了。由于賴氨酸和蘇氨酸均由天門冬氨酸合成得到且第一步反應(yīng)相同 (同工酶)，我們?cè)谀Ｐ椭幸嗫紤]了賴氨酸合成這一分支途徑，以及賴氨酸對(duì)天冬氨酸激酶III的抑制影響，以考察酶量和調(diào)控特性的改變對(duì)兩個(gè)分支途徑間通量分配比的影響。

圖1 大腸桿菌產(chǎn)蘇氨酸模型的生物合成途徑及穩(wěn)態(tài)下的通量分布 (酶和代謝物的全名見表2、表3)Fig. 1 Metabolic pathways in the E. coli threonine biosynthesis model and the flux distribution in the metabolic pathways at steady state.

模型 a和b之間存在缺口，為了整合兩個(gè)模型，我們通過BRENDA數(shù)據(jù)庫的動(dòng)力學(xué)數(shù)據(jù)和其他文獻(xiàn)發(fā)掘工作為13個(gè)新反應(yīng)添加了相應(yīng)的動(dòng)力學(xué)方程和參數(shù)值等相關(guān)信息，模型中新增反應(yīng)及其方程詳見表1。下面對(duì)這些新增反應(yīng)的根據(jù)進(jìn)行詳細(xì)闡述。

PPC和GOT兩步反應(yīng)使模型a中的磷酸烯醇式丙酮酸可以形成草酰乙酸繼而生成天冬氨酸，與模型b實(shí)現(xiàn)連通。由于谷草轉(zhuǎn)氨酶反應(yīng) (GOT)同時(shí)將谷氨酸轉(zhuǎn)化為α-酮戊二酸，因此需要新反應(yīng) GluS由酮戊二酸生成谷氨酸以實(shí)現(xiàn)兩個(gè)代謝物的平衡。以上兩個(gè)反應(yīng)的動(dòng)力學(xué)方程均遵守Ping-Pong機(jī)制[11-12]，動(dòng)力學(xué)參數(shù)的取值參考自BRENDA數(shù)據(jù)庫，各反應(yīng)的rmax取值見表3。

模型a中葡萄糖的轉(zhuǎn)運(yùn)是通過PTS實(shí)現(xiàn)的，在這個(gè)過程中pep轉(zhuǎn)化為丙酮酸，并且pep還通過回補(bǔ)途徑生成草酰乙酸 (反應(yīng) PPC)用于蘇氨酸合成，因此在代謝網(wǎng)絡(luò)分析求得的最優(yōu)途徑中還包括反應(yīng)PWD以將pyr轉(zhuǎn)化回pep。同時(shí)由于該過程生成了 AMP，還需要 AdeK反應(yīng)以保證AMP的產(chǎn)生與消耗之間的平衡。以上兩個(gè)反應(yīng)均為不可逆反應(yīng)，動(dòng)力學(xué)方程采用雙底物不可逆順序反應(yīng)機(jī)制[15-16]。

表1 模型中新增反應(yīng)及其動(dòng)力學(xué)方程Table 1 New reactions and kinetic rate equations in the integrated model

續(xù)表1

Rodríguez-Prados等的實(shí)驗(yàn)分析表明蘇氨酸合成時(shí)三羧酸循環(huán)各反應(yīng)亦有較高活性, 主要是為蘇氨酸合成提供能量和還原力[7]。在我們的模型中類似于模型a采用了一個(gè)合并反應(yīng) (PDH)來表示pyr經(jīng)過三羧酸循環(huán)后分解產(chǎn)生 ATP、NADH和FADH2的過程。由于大部分的ATP都是通過氧化磷酸化途徑產(chǎn)生的，我們?cè)谀Ｐ椭刑砑恿?RCR1(Respiratory chain reaction)和RCR2兩個(gè)反應(yīng)用以體現(xiàn)NADH和FADH2經(jīng)呼吸鏈電子傳遞后生成ATP的過程。其中 ATP/NADH的轉(zhuǎn)化計(jì)量系數(shù)1.85和ATP/FADH2的轉(zhuǎn)化計(jì)量系數(shù)1都是基于文獻(xiàn)報(bào)道的呼吸鏈效率確定，均低于其最大理論值[20]。由于PDH、RCR1和RCR2這3個(gè)反應(yīng)均為多個(gè)反應(yīng)的合并反應(yīng)，無法通過酶動(dòng)力學(xué)機(jī)理確定合適的動(dòng)力學(xué)方程，因此我們補(bǔ)全PDH反應(yīng)式中的計(jì)量關(guān)系后，采用了模型a中原有的PDH反應(yīng)動(dòng)力學(xué)方程及參數(shù)。RCR1反應(yīng)的動(dòng)力學(xué)方程以呼吸鏈第一步反應(yīng)的 Ping-Pong機(jī)制及其動(dòng)力學(xué)參數(shù)值代表[18]。由于RCR2中的FADH2的消耗速率與PDH中的FADH2生成速率必須保持相同，RCR2的動(dòng)力學(xué)方程及其參數(shù)值設(shè)置與PDH中基本相同，不同之處在于其反應(yīng)速率隨FADH2的濃度變化而變化[9]。NDPK反應(yīng)實(shí)現(xiàn)了 ATP和GTP(三羧酸循環(huán)途徑中產(chǎn)生)之間的轉(zhuǎn)換，反應(yīng)動(dòng)力學(xué)機(jī)制為Ping-Pong機(jī)制[19]。三羧酸循環(huán)過程產(chǎn)生的主要是 NADH，而蘇氨酸合成需要的是NADPH。PP途徑產(chǎn)生的NADPH不足以滿足需要，需要由NADH轉(zhuǎn)化得到。因此我們添加了NT反應(yīng)以實(shí)現(xiàn)NADH和NADPH之間的轉(zhuǎn)化，反應(yīng)動(dòng)力學(xué)方程采用可逆的雙底物順序反應(yīng)機(jī)制[17]。賴氨酸合成途徑也用一個(gè)合并反應(yīng)LS表示, 其受產(chǎn)物賴氨酸的抑制，反應(yīng)動(dòng)力學(xué)方程參考 Contador等構(gòu)建的模型[14]。蘇氨酸和賴氨酸的排出反應(yīng)TD和LD均采用一級(jí)反應(yīng)動(dòng)力學(xué)方程，其反應(yīng)速率與濃度成正比 (文中所涉及的全部反應(yīng)與代謝物名稱簡寫的全名參見表2、表3)。

1.2 模型修正及初始濃度設(shè)置

為保障模型構(gòu)建的準(zhǔn)確性，首先檢查兩個(gè)模型中的計(jì)量關(guān)系和反應(yīng)可逆性信息。兩個(gè)模型由于反應(yīng)物的代謝途徑關(guān)系不完整，往往需要固定大量的反應(yīng)物濃度才可使模型達(dá)到穩(wěn)態(tài)。整合后的模型由于添加了13個(gè)新反應(yīng)和完善了能量及還原力的計(jì)量關(guān)系，在極大程度上提高了模型的自洽能力。新模型中，ATP、NADH、NADPH和FADH2等物質(zhì)從模型 a中固定濃度的全局變量，調(diào)整為隨系統(tǒng)運(yùn)行發(fā)生濃度變化并最終可以達(dá)到動(dòng)態(tài)平衡的反應(yīng)物。因此，在固定外部葡萄糖濃度的前提下，只需固定磷酸根離子和二氧化碳的濃度為較高的數(shù)值 (視作充足供應(yīng))，對(duì)其余包括能量和還原力在內(nèi)的反應(yīng)物的初始濃度賦值不再嚴(yán)苛要求。除葡萄糖外的各反應(yīng)物的初始濃度只影響模型最初運(yùn)行階段的反應(yīng)速率，并不會(huì)影響系統(tǒng)最終的穩(wěn)定狀態(tài)和控制關(guān)系，所以模型中的初始反應(yīng)物濃度多沿用模型a和b中的數(shù)據(jù)。新增反應(yīng)物的初始濃度參考細(xì)胞組成數(shù)據(jù)范圍進(jìn)行賦值[21-22]，最終所有反應(yīng)物初始濃度值設(shè)置詳見表2。

1.3 模型調(diào)試及最終參數(shù)設(shè)置

由于兩個(gè)模型出于不同的構(gòu)建目的、創(chuàng)建于不同的反應(yīng)體系，各自的通量水平存在很大的差異。添加新反應(yīng)以及能量和還原力的計(jì)量關(guān)系后，模型初步運(yùn)行后無法達(dá)到穩(wěn)定狀態(tài)。通過觀察穩(wěn)態(tài)運(yùn)行的結(jié)果，發(fā)現(xiàn)無法得到穩(wěn)定解的原因是中心代謝部分的通量遠(yuǎn)高于蘇氨酸合成部分的通量，導(dǎo)致 asp不斷積累。這種現(xiàn)象在整合不同模型的過程中非常普遍，也是模型整合要解決的首要問題。從pep形成oaa開始，至最后生成thr的反應(yīng)過程為線性途徑，且pep生成oaa的反應(yīng)決定了糖酵解途徑向蘇氨酸合成子途徑的通量分配，因此asp的積累意味著模型a向模型b的通量分配高于模型b所能承受的通量，需將合成asp的PPC反應(yīng)速率降低，或者提高消耗asp的AK反應(yīng)的速率。調(diào)整asp的相關(guān)反應(yīng)速率后，asp的積累得到控制，但隨之會(huì)帶來其他代謝物的積累，但可采用類似思路進(jìn)行調(diào)整以最終得到穩(wěn)定解。需要說明的是由于代謝網(wǎng)絡(luò)中存在的復(fù)雜調(diào)控關(guān)系及非線性動(dòng)力學(xué)方程，這個(gè)過程常常需要不斷的試差才能得到一個(gè)可以穩(wěn)定運(yùn)行的模型。目前階段還沒有一套標(biāo)準(zhǔn)化的系統(tǒng)的調(diào)試方法，但可根據(jù)經(jīng)驗(yàn)把主要的調(diào)試方法分為3種情況：1) 對(duì)于存在于線性途徑中的反應(yīng)物，可以先將途徑外的其他反應(yīng)的速率設(shè)置為零，然后根據(jù)途徑中反應(yīng)物的積累情況，增加其消耗反應(yīng)的最大反應(yīng)速率rmax或者降低其作為生成物的最大反應(yīng)速率即可。2) 針對(duì)類似 ATP(ADP和 AMP)、NAD(NADH)以及NADP(NADPH)等這樣在很多反應(yīng)中均涉及的、但存在固定轉(zhuǎn)化關(guān)系的反應(yīng)物組，可以先將其他組反應(yīng)物的濃度固定，如固定 ATP(ADP和AMP)和NAD (NADH)兩組代謝物的濃度，單獨(dú)調(diào)整NADP和NADPH之間的平衡，就很容易發(fā)現(xiàn)導(dǎo)致崩潰的原因。同理再逐步將其他兩組反應(yīng)物改為可變的，一一調(diào)整即可。3)通過COPASI[23–24]軟件自帶的參數(shù)掃描功能對(duì)各反應(yīng)rmax進(jìn)行魯棒性分析，著重調(diào)整魯棒性較差的反應(yīng)的rmax。經(jīng)過上述過程的反復(fù)調(diào)試，我們最終使模型得到穩(wěn)態(tài)解，此時(shí)模型中的部分參數(shù)設(shè)置如表3所示。

2 模型模擬分析

2.1 穩(wěn)態(tài)通量分布

我們利用COPASI軟件對(duì)上述構(gòu)建好的模型進(jìn)行穩(wěn)態(tài)運(yùn)算，存在賴氨酸分支時(shí)，得到 PTS反應(yīng)的通量為 0.268 mmol/(L·s)，蘇氨酸合成的通量為 0.155 mmol/(L·s)，六碳的葡萄糖轉(zhuǎn)化為四碳的蘇氨酸的碳摩爾得率為 38.56%，敲除賴氨酸分支，蘇氨酸得率提高為91.17%，接近FBA計(jì)算的途徑最大得率，分支途徑的敲除對(duì)蘇氨酸得率的提高非常重要[3]。途徑中各反應(yīng)的通量分布和主要調(diào)控關(guān)系如圖1所示。

表2 模型中設(shè)定的代謝物的初始濃度Table 2 Initial concentrations of metabolites given in the integrated model

表3 模型中最終設(shè)定的rmax值Table 3 Final rmax values in the model

2.2 代謝控制分析

通量控制系數(shù) (Flux control coefficient,FCC)表明了代謝途徑中的酶活性改變對(duì)途徑穩(wěn)態(tài)通量的影響，定義式見方程 (1)[25]。通過比較途徑中各酶反應(yīng)過程對(duì)某一途徑通量的控制系數(shù)的大小，可以確定對(duì)代謝途徑 (此模型中指對(duì)蘇氨酸合成通量)起關(guān)鍵控制作用的酶，進(jìn)而通過調(diào)節(jié)關(guān)鍵酶的濃度及活性來改變途徑中的通量分布[26]。對(duì)于FCC為正數(shù)的酶反應(yīng)，表明增加其酶量或酶活，可以提高目標(biāo)反應(yīng)的通量；對(duì)于FCC為負(fù)數(shù)的酶反應(yīng)，增加其酶量或酶活，則對(duì)提高目標(biāo)反應(yīng)的通量不利。數(shù)值越大表明該酶對(duì)通量影響越大, 是理想的改造位點(diǎn)。

2.2.1 關(guān)鍵酶預(yù)測(cè)

通過代謝控制分析，我們獲得了模型穩(wěn)態(tài)狀態(tài)下酶與反應(yīng)間的FCC矩陣。選擇FCC矩陣中各酶對(duì)蘇氨酸合成通量和賴氨酸合成通量的FCC值，分布如圖2所示。

由圖2可見FCC數(shù)值較大的反應(yīng)均為途徑中的不可逆反應(yīng)或者位于分支點(diǎn)，主要集中在糖酵解過程以及能量和還原力的生成過程，這說明前體、能量和還原力的供給是蘇氨酸合成的主要限制因素。圖2A表明對(duì)蘇氨酸合成，控制步驟主要為ASADH和HDH，其中HDH增加對(duì)蘇氨酸的合成有利而對(duì)賴氨酸合成不利，ASADH則正好相反；同時(shí)，G6PDH和GAPDH兩個(gè)脫氫酶對(duì)蘇氨酸和賴氨酸兩個(gè)競(jìng)爭(zhēng)分支的控制作用也是相反的，表明過表達(dá)這些酶可以為蘇氨酸合成提供更多還原力，使更多碳流從賴氨酸合成轉(zhuǎn)向蘇氨酸合成。PPC對(duì)蘇氨酸合成反應(yīng)的控制系數(shù)為負(fù)值，表明 PPC增加對(duì)蘇氨酸的合成并不總是有利的，這主要是由于反應(yīng)的復(fù)雜非線性特征及產(chǎn)物合成與能量和還原力供給之間的復(fù)雜協(xié)調(diào)機(jī)制造成的。這一模擬分析結(jié)果與Lee等的實(shí)驗(yàn)研究結(jié)果一致[3]。

以往的研究結(jié)果已充分證實(shí)了蘇氨酸對(duì)其合成途徑中的天冬氨酸激酶Ⅰ、高絲氨酸脫氫酶Ⅰ和蘇氨酸合成酶的反饋抑制作用[5-6]，因此我們?cè)谀Ｐ椭锌疾炝私獬K氨酸抑制對(duì)其合成通量和得率的影響，結(jié)果見圖 2B。結(jié)果表明，通過解除蘇氨酸的抑制可以有效降低ASADH和HDH對(duì)蘇氨酸合成的限制作用，更有利于碳流流向蘇氨酸的合成。在模型中將蘇氨酸抑制解除后，PTS/TS/LS的通量分布由圖 1中的 0.268/0.155/0.134變?yōu)?.248/0.296/0.027，糖轉(zhuǎn)運(yùn)系統(tǒng)PTS通量沒有明顯變化，但是蘇氨酸和賴氨酸兩個(gè)分支的合成通量比值發(fā)生明顯改變，由原來的 1.16∶1增加至 10.96∶1，蘇氨酸相對(duì)于葡萄糖的碳摩爾得率從38.56%增加至79.57%，這充分說明了解除蘇氨酸抑制對(duì)提高蘇氨酸通量和得率的必要性。

圖2 不同酶對(duì)蘇氨酸和賴氨酸合成反應(yīng)的通量控制系數(shù) (a: 存在蘇氨酸抑制；b: 解除蘇氨酸抑制)Fig. 2 Flux control coefficients of enzymes on the threonine and lysine biosynthesis. (A)FCCs of different enzymes at unrelieved threonine inhibition. (B)completely relieve the threonine inhibition.

此外，敲除賴氨酸分支途徑，重新計(jì)算后，PK對(duì) TS的 FCC為負(fù)值，這一結(jié)果與Rodríguez-Prados等不考慮賴氨酸分支途徑時(shí)提出的PK敲除策略并不矛盾[7]。添加了賴氨酸分支途徑，考慮了能量和還原力的平衡關(guān)系后，在我們的控制步驟預(yù)測(cè)結(jié)果中，PK并不是蘇氨酸合成過程的關(guān)鍵酶，敲除PK后對(duì)TS的通量影響并不明顯。對(duì)此我們還需要在后續(xù)工作中對(duì)模型進(jìn)一步驗(yàn)證和完善。

2.2.2 關(guān)鍵酶過表達(dá)分析

對(duì)考慮蘇氨酸抑制和賴氨酸分支的模型初始狀態(tài)進(jìn)行的控制分析結(jié)果表明，PTS、G6PDH和 HDH對(duì)蘇氨酸合成反應(yīng)的通量具有較大影響，即改變其酶量可以顯著提高蘇氨酸合成速率。在我們的動(dòng)力學(xué)模型中酶量體現(xiàn)在參數(shù)rmax中，因此我們通過提高PTS、G6PDH和HDH三個(gè)反應(yīng)的 rmax來考察其對(duì)提高蘇氨酸合成通量的影響，結(jié)果如圖3所示。從圖3A可以看出，隨著 PTS反應(yīng)速率的提高，代表糖轉(zhuǎn)運(yùn)系統(tǒng)的PTS反應(yīng)、蘇氨酸合成的TS反應(yīng)和賴氨酸合成的LS反應(yīng)的通量均增加，但最大反應(yīng)速率增加至1.5倍后，PTS和TS的合成通量開始下降，只有賴氨酸合成繼續(xù)增加。圖3B中，當(dāng)G6PDH的最大反應(yīng)速率增加時(shí)，PTS的通量已經(jīng)趨于平衡，蘇氨酸合成通量在上升，賴氨酸合成通量逐漸下降。圖3C中，增加HDH的最大反應(yīng)速率，PTS的通量沒有發(fā)生明顯變化，蘇氨酸和賴氨酸的通量分布發(fā)生明顯改變，通過在模型中解除蘇氨酸對(duì)HDH的抑制，也可以使蘇氨酸和賴氨酸的合成通量比值增加至10.54∶1，說明HDH對(duì)蘇氨酸合成的限制非常強(qiáng)烈。通量控制系數(shù)是系統(tǒng)中的局部概念，通過對(duì)上面3個(gè)模型初始狀態(tài)預(yù)測(cè)的關(guān)鍵酶進(jìn)行過表達(dá)，在一定范圍內(nèi)確實(shí)可以增加目標(biāo)通量，但繼續(xù)增加效果卻適得其反，因此在菌種改造過程中關(guān)鍵酶的過表達(dá)并非越多越好,而是有一個(gè)最佳范圍。我們的模型構(gòu)建和模擬分析工作都是在專用的生化反應(yīng)動(dòng)力學(xué)模型工具軟件COPASI中進(jìn)行的[23-24]。

3 結(jié)論

圖3 穩(wěn)態(tài)時(shí)糖輸入通量 (PTS)、蘇氨酸合成通量 (TS)和賴氨酸合成通量 (LS)隨關(guān)鍵酶過表達(dá)的變化Fig. 3 Effect of enzyme amounts on the PTS, threonine synthesis and lysine synthesis fluxes. (A)PTS over expression.(B)G6PDH over expression. (C)HDH over expression.

我們?cè)诨蚪M規(guī)模代謝網(wǎng)絡(luò)分析求得的蘇氨酸合成最優(yōu)途徑基礎(chǔ)上，整合了已有的中心代謝動(dòng)力學(xué)模型和從天冬氨酸出發(fā)的蘇氨酸合成動(dòng)力學(xué)模型，借助文獻(xiàn)信息對(duì)一些新反應(yīng)補(bǔ)充了動(dòng)力學(xué)方程和參數(shù)信息，并考慮了整個(gè)途徑的能量和還原力平衡關(guān)系、從而構(gòu)建了自洽的從葡萄糖到蘇氨酸的完整代謝途徑動(dòng)力學(xué)模型。對(duì)新的整合模型的代謝控制分析表明控制系數(shù)較大的反應(yīng)均為不可逆反應(yīng)或者處于途徑的分支點(diǎn)，包括前體、能量和還原力的生成和轉(zhuǎn)化過程，這表明蘇氨酸合成過程需要前體、能量和還原力的協(xié)調(diào)分配。分析結(jié)果表明解除蘇氨酸抑制確實(shí)可以顯著改變蘇氨酸和賴氨酸兩個(gè)分支間的通量分配，增加蘇氨酸的得率。對(duì)關(guān)鍵酶如PTS、G6PDH和HDH的分析表明在一定范圍內(nèi)提高其速率確實(shí)可以有效提高蘇氨酸合成途徑的通量，但關(guān)鍵酶的過表達(dá)應(yīng)在一適當(dāng)范圍內(nèi)，因?yàn)殡S著酶量的增加該酶的控制系數(shù)也會(huì)減小甚至變?yōu)樨?fù)值，繼續(xù)過表達(dá)反而對(duì)產(chǎn)物合成不利。

[1]Wang H, Ma HW, Zhao XM. Progress in genome-scale metabolic network: a review. Chin J Biotech, 2010, 26(10): 1340–1348 (in Chinese).王暉, 馬紅武, 趙學(xué)明. 基因組尺度代謝網(wǎng)絡(luò)研究進(jìn)展. 生物工程學(xué)報(bào), 2010, 26(10): 1340–1348.

[2]Kr?mer JO, Wittmann C, Schr?der H, et al.Metabolic pathway analysis for rational design of L-methionine production by Escherichia coli and Corynebacterium glutamicum. Metab Engin, 2006,8(4): 353–369.

[3]Lee KH, Park JH, Kim TY, et al. Systems metabolic engineering of Escherichia coli for L-threonine production. Mol Syst Biol, 2007, 3: 149–156.

[4]Liu S, Liang Y, Liu Q, et al. Development of a two-stage feeding strategy based on the kind and level of feeding nutrients for improving fed-batch production of L-threonine by Escherichia coli. Appl Microbiol Biotechnol, 2013, 97(2): 573–583.

[5]Chassagnole C, Ra?s B, Quentin E, et al. An integrated study of threonine-pathway enzyme kinetics in Escherichia coli. Biochem J, 2001, 356:415–423.

[6]Chassagnole C, Fell DA, Ra?s B, et al. Control of the threonine-synthesis pathway in Escherichia coli:a theoretical and experimental approach. Biochem J,2001, 356: 433–444.

[7]Rodríguez-Prados JC, De Atauri P, Maury J, et al. In silico strategy to rationally engineer metabolite production: a case study for threonine in Escherichia coli. Biotechnol Bioeng, 2009, 103(3): 609–620.

[8]Orth JD, Conrad TM, Na J, et al. A comprehensive genome-scale reconstruction of Escherichia coli metabolism——2011. Mol Syst Biol, 2011, 7: 535.

[9]Chassagnole C, Noisommit-Rizzi N, Schmid JW, et al. Dynamic modeling of the central carbon metabolism of Escherichia coli. Biotechnol Bioeng,2002, 79(1): 53–73.

[10]Le Novère N, Bornstein B, Broicher A, et al.BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Oxford J, 2005,34(suppl 1): 689–691.

[11]Barden RE, Fung CH, Utter MF, et al. Pyruvate carboxylase from chicken liver: steady state kinetic studies indicate a “Two Site” ping-pong mechanism. J Biol Chem, 1972, 247: 1323–1333.

[12]Sung HAC, Tachiki T, Kumagai H, et al. Properties of glutamate synthase from Brevibachterium flavum. J Ferment Technol, 1984, 62: 569–575.

[13]Vanoni MA, Nuzzi L, Rescigno M, et al. The kinetic mechanism of the reactions catalyzed by the glutamate synthase from Azospirillum brasilense.Eur J Biochem, 1991, 202(1): 181–189.

[14]Contador CA, Rizk ML, Asenjo JA, et al. Ensemble modeling for strain development of L-lysine producing Escherichia coli. Metab Eng, 2009, 11:221–233.

[15]Tjaden B, Plagens A, D?rr C, et al.Phosphoenolpyruvate synthetase and pyruvate,phosphate dikinase of Thermoproteus tenax: key pieces in the puzzle of archaeal carbohydrate metabolism. Mol Microbiol, 2006, 60(2): 287–298.

[16]Narindrasorasak S, Bridger WA.Phosphoenolypyruvate synthetase of Escherichia coli: molecular weight, subunit composition, and identification of phosphohistidine in phosphoenzyme intermediate. J Biol Chem, 1977,252(10): 3121–3127.

[17]Freeman H, Shimomura K, Horner E, et al.Nicotinamide nucleotide transhydrogenase: a key role in insulin secretion. Cell Metab, 2006, 3(1):35–45.

[18]Velázquez I, Pardo JP. Kinetic characterization of the rotenone-insensitive internal NADH: ubiquinone oxidoreductase of mitochondria from Saccharomyces cerevisiae. Arch Biochem Biophys,2001, 389(1): 7–14.

[19]Morin-Leisk J, Lee TH. Nucleotide-dependent self-assembly of Nucleoside Diphosphate Kinase(NDPK)in vitro. BBA Proteins Proteom, 2008,1784(12): 2045–2051.

[20]Oberhardt MA, Pucha?ka J, Martins dos Santos VAP, et al. Reconciliation of genome-scale metabolic reconstructions for comparative systems analysis. PLoS Comput Biol, 2011, 7(3): e1001116.

[21]Pramanik J, Keasling JD. Effect of Escherichia coli biomass composition on central metabolic fluxes predicted by a stoichiometric model.Biotechnol Bioeng, 1998, 60(2): 230–238.

[22]Pramanik J, Keasling JD. Stoichiometric model of Escherichia coli metabolism: incorporation of growth-rate dependent biomass composition and mechanistic energy requirements. Biotechnol Bioeng, 1997, 56(4): 398–421.

[23]Hoops S, Sahle S, Gauges R, et al. COPASI–a complex pathway simulator. Bioinformatics, 2006,22(24): 3067–3074.

[24]Mendes P, Hoops S, Sahle S, et al. Computational modeling of biochemical networks using COPASI.Syst Biology, 2009, 500: 17–59.

[25]Pryor HJ, Smyth JE, Quinlan PT, et al. Evidence that the flux control coefficient of the respiratory chain is high during gluconeogenesis from lactate in hepatocytes from starved rats. Biochem J, 1987,247: 449–457.

[26]Nikolaev EV. The elucidation of metabolic pathways and their improvements using stable optimization of large-scale kinetic models of cellular systems. Metab Eng, 2010, 12: 26–38.