時(shí)滯基因調(diào)控模型的穩(wěn)定性與Hopf分支

宋永利，葉 子

（同濟(jì)大學(xué) 數(shù)學(xué)系，上海 200092）

基因表達(dá)是指遺傳物質(zhì)所攜帶的信息從DNA傳遞到蛋白質(zhì)的過程，是生命系統(tǒng)中最基本的過程，其整個(gè)過程滿足一個(gè)基本的中心法則，即遺傳信息通過DNA以自我為模板進(jìn)行復(fù)制，以DNA為模板的RNA的合成（轉(zhuǎn)錄），和以RNA為模板的蛋白質(zhì)合成［1］.并非所有細(xì)胞的全部基因總在不斷地表達(dá)，細(xì)胞在不同時(shí)間、不同空間以不同的組合表達(dá)不同的基因，所有這些差異的關(guān)鍵在于基因表達(dá)過程的調(diào)控，研究基因表達(dá)的最終目的是為了了解一個(gè)復(fù)雜基因調(diào)控網(wǎng)絡(luò)是如何工作以達(dá)到控制產(chǎn)物蛋白的濃度，以實(shí)現(xiàn)不同的生物學(xué)功能.真核細(xì)胞中常見的調(diào)控方式有：轉(zhuǎn)錄水平調(diào)控，即基因在何時(shí)以怎樣的頻率被轉(zhuǎn)錄；翻譯調(diào)控，即細(xì)胞質(zhì)中哪些mRNA被核糖體翻譯；mRNA降解調(diào)控，即細(xì)胞質(zhì)中哪些mRNA被去穩(wěn)定化；蛋白質(zhì)活性調(diào)控，即決定翻譯的蛋白質(zhì)的活性、抑制、區(qū)域化以及降解.

Novák和Tyson［2］針對(duì)類基因調(diào)控模型提出了關(guān)于一種反映mRNA和其控制合成的蛋白質(zhì)的濃度變化的數(shù)學(xué)模型.由于mRNA在轉(zhuǎn)錄剛剛完成的較短時(shí)間內(nèi)不能直接使用，mRNA和調(diào)控蛋白需要在細(xì)胞核與細(xì)胞質(zhì)之間運(yùn)輸?shù)?，這就使得在基因調(diào)控模型中時(shí)滯因素是不可避免的.本文首先研究了這個(gè)模型正平衡點(diǎn)的存在性與穩(wěn)定性；其次，探討了時(shí)滯因素對(duì)正平衡點(diǎn)的穩(wěn)定性的影響以及誘發(fā)的振蕩現(xiàn)象；最后通過數(shù)值計(jì)算驗(yàn)證了動(dòng)力學(xué)分析的理論結(jié)果.

1 正平衡點(diǎn)的穩(wěn)定性及時(shí)滯誘發(fā)的Hopf分支

Novák和Tyson［2］提出的數(shù)學(xué)模型由以下常微分方程組來描述：

式中：x表示mRNA的濃度；y表示蛋白質(zhì)濃度.蛋白質(zhì)濃度的增加會(huì)對(duì)mRNA的合成起抑制作用，這里用系數(shù)為p的Hill函數(shù)表示，而蛋白質(zhì)也會(huì)對(duì)自身的合成起抑制作用，這里用 Michaelis-Menten函數(shù)表示.另一方面，mRNA濃度的增加會(huì)促進(jìn)蛋白質(zhì)的合成，其自身還存在著降解作用.參數(shù)k1為mRNA的合成速率常數(shù)，kdx為蛋白質(zhì)的降解速率常數(shù)，ksy為蛋白質(zhì)的合成速率常數(shù)，k2為蛋白質(zhì)抑制自身合 成 速 率 常 數(shù)，S，ET，Kd，Km為 Michaelis－Menten函數(shù)與 Hill函數(shù)中相應(yīng)的常數(shù)［3－4］.

在以下研究中限定p＝2，為了方便推導(dǎo)，作一系列變量代換，轉(zhuǎn)化為更簡(jiǎn)潔的數(shù)學(xué)形式，令k1S＝m1，（Kd）2＝m2，kdx＝m3，ksy＝m4，k2ET＝m5，Km＝m6，則方程組可 轉(zhuǎn)化為以下方程組：

由于mi＞0，i＝1，…，6，方程組存在唯一正解，即方程組（2）有唯一的正平衡點(diǎn)，記為E（x＊，y＊）.下面分析該平衡點(diǎn)的局部穩(wěn)定性.考慮方程組關(guān)于此平衡點(diǎn)的的線性化系統(tǒng)，特征矩陣為

對(duì)應(yīng)的特征方程為

1.1 離散時(shí)滯

由于蛋白質(zhì)的有限合成速度以及調(diào)控蛋白在細(xì)胞核與細(xì)胞質(zhì)之間運(yùn)輸需要一定的時(shí)間，因此，mRNA的濃度事實(shí)上受當(dāng)前時(shí)刻（t時(shí)刻）以前的某個(gè)時(shí)刻（比如t－τ時(shí)刻）的蛋白質(zhì)的濃度的影響.為了研究這種時(shí)間滯后對(duì)基因調(diào)控模型平衡態(tài)的穩(wěn)定性的影響，首先研究以下具有離散時(shí)滯的模型：

方程組（5）和（2）具有相同的平衡點(diǎn)E（x＊，y＊）.方程組（5）在該平衡點(diǎn)的線性化系統(tǒng)為

其中

故方程組（5）在平衡點(diǎn)E（x＊，y＊）的特征方程為

如果方程（6）有一對(duì)純虛根±iω（ω＞0），則ω滿足以下方程：

分離實(shí)部虛部得到

消去三角函數(shù)項(xiàng)，得

計(jì)算判別式得

因此，方程（9）只有一個(gè)正實(shí)根的充要條件是下式成立：

上式化簡(jiǎn)后即為

因此，當(dāng)不等式（10）滿足時(shí)，特征方程（6）ゆ在一對(duì)純虛根，假設(shè)為±iω＊.

根據(jù)方程組（8）可以得到

化簡(jiǎn)上式得到

由方程組（8）和式（12），經(jīng)過簡(jiǎn)單計(jì)算可以得到下式：

由式（10），（14）以及時(shí)滯微分方程的定性理論［6］，關(guān)于方程組（5）的穩(wěn)定性有如下定理：

定理1 （1）如果m1m2m4m6－2m3m5（y＊）3≥0，則對(duì)任意的τ≥0，方程組（5）的平衡點(diǎn)E（x＊，y＊）都是穩(wěn)定的.

（2）如果m1m2m4m6－2m3m5（y＊）3＜0，則當(dāng)τ∈0，τ］（0時(shí)，方程組（5）的平衡點(diǎn)E（x＊，y＊）是漸近穩(wěn)定的；當(dāng)τ＞τ0時(shí)，方程組（5）的平衡點(diǎn)E（x＊，y＊）是不穩(wěn)定的.

（3）τ＝τj是方程組（5）的 Hopf分支點(diǎn).

1.2 分布時(shí)滯

由于生物體內(nèi)發(fā)生的化學(xué)反應(yīng)所需時(shí)間都比較長(zhǎng)，反應(yīng)物在細(xì)胞中運(yùn)輸所需的時(shí)間也有一定的不確定性，t時(shí)刻的mRNA的濃度與t時(shí)刻以前的一個(gè)較長(zhǎng)時(shí)間段內(nèi)的蛋白質(zhì)的濃度都有關(guān)系.事實(shí)上，在生物系統(tǒng)中分布時(shí)滯的現(xiàn)象普遍存在，它更能精確刻畫時(shí)滯因素的影響作用.以下研究具有分布時(shí)滯的基因調(diào)控模型.

其中時(shí)滯核函數(shù)k：［0，∞）→［0，∞）分段連續(xù)，且滿足以下正規(guī)化條件

生物系統(tǒng)中常見的時(shí)滯弱核函數(shù)［7］為

滿足條件（16），其中β表示過去時(shí)間內(nèi)蛋白質(zhì)濃度的衰減速度.在正規(guī)化條件（16）滿足時(shí)，方程組（15）和（2）具有相同的正平衡點(diǎn)E（x＊，y＊）.以下研究β對(duì)該平衡點(diǎn)穩(wěn)定性的影響.

引入輔助變量u

由方程組（15）和式（17），可以得到以下關(guān)于x，y，u的三元常微分方程組：

方程組（18）的平衡點(diǎn)為F（x＊，y＊，y＊）.原方程組的E（x＊，y＊）的穩(wěn)定性等價(jià)于方程組（18）的平衡點(diǎn)F（x＊，y＊，y＊）的穩(wěn)定性.方程組在平衡點(diǎn)F（x＊，y＊，y＊）的特征矩陣為

相對(duì)應(yīng)的特征方程為

根據(jù)Routh-Hurwitz判別法，方程（20）所有根皆具嚴(yán)格負(fù)實(shí)部的充要條件是以下3個(gè)不等式同時(shí)成立：

由于D1＞0恒成立而且

所以，方程（20）所有根皆具有嚴(yán)格負(fù)實(shí)部，等價(jià)于D2＞0，即

當(dāng)β＝β1或β＝β2時(shí)，方程（20）具有零實(shí)部的根，由于λ＝0不是方程（20）的根，因此，假設(shè)iω（ω＞0）是方程（20）的根.將其代入方程（20），得到

由于ω不為零，故式（22）可化為

消去ω2，即得

易見式（23）成立的充要條件為a2＜0且Δ＝a22－4a1a3≥0.此時(shí)，由方程組（22）得到

由方程組（22），式（24）可以化為

由式（25），得到

把β1，β2的值代入式（27）可以得到下式：

由以上的討論，關(guān)于特征方程（20）的根的分布，有以下結(jié)果：

引理1 （1）假設(shè)a2≥0或a2＜0且Δ＝a22－4a1a3＜0，則對(duì)任意β∈（0，＋∞），方程的所有根皆具有嚴(yán)格負(fù)實(shí)部.

（2）假設(shè)a2＜0且Δ＝a22－4a1a3≥0，則當(dāng)β∈（0，β1）∪（β2，＋∞）時(shí)，方程（20）的所有根皆具有嚴(yán)格負(fù)實(shí)部.

（3）當(dāng)β＝β1或β＝β2時(shí)，方程（20）有一對(duì)純虛根±iω；當(dāng)β∈（β1，β2）時(shí)，方程（20）有一對(duì)具正實(shí)部的根.

由引理1和橫截條件式（28），關(guān)于方程組（15）的正平衡點(diǎn)E（x＊，y＊）穩(wěn)定性和分支有以下定理：

定理2 （1）假設(shè)a2≥0或a2＜0且Δ＝a22－4a1a3＜0，則對(duì)任意β∈（0，＋∞），方程組（15）的正平衡點(diǎn)E（x＊，y＊）都是局部漸近穩(wěn)定的.

（2）假設(shè)a2＜0且Δ＝a22－4a1a3≥0，則當(dāng)β∈（0，β1）∪ （β2，＋ ∞）時(shí)，方 程 組 （15）的 正 平 衡 點(diǎn)E（x＊，y＊）是局部漸近穩(wěn)定的；當(dāng)β∈（β1，β2）時(shí)，方程組（15）的正平衡點(diǎn)E（x＊，y＊）是不穩(wěn)定的.

（3）β1和β2是方程組的Hopf分支值.

2 數(shù)值模擬

對(duì)取定的系統(tǒng)參數(shù)，利用Winpp軟件模擬基因調(diào)控模型的數(shù)值解.取參數(shù)：m1＝0.7，m2＝0.1，m3＝0.1，m4＝1，m5＝1，m6＝0.2.則唯一的正平衡點(diǎn)為：E（x＊，y＊），x＊＝0.813 4，y＊＝0.872 1.

取初始值為 （x，y）＝（1，3），方程（2）的數(shù)值解曲線和相圖如圖1所示.圖1表明，方程（2）的正平衡點(diǎn)是漸近穩(wěn)定的.

圖1 方程（2）的正平衡點(diǎn)漸近穩(wěn)定性Fig.1 Stability of positive equilibrium of Equation（2）

以下考察離散時(shí)滯和分布時(shí)滯對(duì)此平衡點(diǎn)的穩(wěn)定性的影響.在給定的參數(shù)值下，m1m2m4m6－2m3m5（y＊）3＝－0.118 7＜0，即式滿足.首先，對(duì)具有離散時(shí)滯的方程組（5），可以計(jì)算得到：τ0≈1.799 5，τj≈1.799 5＋16.497 9j，j＝0，1，2，….

取τ＝1.6＜τ0，具離散時(shí)滯的方程（5）的數(shù)值解曲線和相圖如圖2所示.此時(shí)正平衡點(diǎn)E（x＊，y＊）仍然是穩(wěn)定的.

圖2 當(dāng)τ＜τ0時(shí)，方程（5）的正平衡點(diǎn)漸近穩(wěn)定Fig.2 Stability of positive equilibrium of Equation（5）is stable whenτ＜τ0

取τ＝1.8＞τ0，數(shù)值解曲線和相圖如圖3所示.此時(shí)正平衡點(diǎn)Ex＊，y（）＊是不穩(wěn)定的，但在平衡點(diǎn)的附近存在小振幅的穩(wěn)定的周期解.

圖3 當(dāng)τ＞τ0時(shí)，方程（5）的正平衡點(diǎn)不穩(wěn)定Fig.3 Stability of positive equilibrium of Equation（5）is unstable whenτ＞τ0

最后，考察分布時(shí)滯對(duì)正平衡點(diǎn)Ex＊，y（）＊的穩(wěn)定性的影響.此時(shí)有：β1≈0.066 7，β2≈0.261 0，a2＝－0.089 8＜0，Δ＝a22－4a1a3＝0.002 8＞0.

取β＝0.03＜β1，具分布時(shí)滯的方程組（15）的數(shù)值解曲線和相圖如圖4所示.此時(shí)正平衡點(diǎn)E（x＊，y＊）仍然是穩(wěn)定的.取β1＜β＝0.25＜β2，具分布時(shí)滯的方程組（15）的數(shù)值解曲線和相圖如圖5所示.此時(shí)正平衡點(diǎn)Ex＊，y（）＊是不穩(wěn)定的，類似于離散時(shí)滯誘發(fā)的Hopf分支，此時(shí)在平衡點(diǎn)Ex＊，y（）＊的附近存在小振幅的穩(wěn)定的周期解.取β＝0.45＞β2，具分布時(shí)滯的方程組（15）的數(shù)值解曲線和相圖如圖6所示.此時(shí)正平衡點(diǎn)Ex＊，y（）＊又是穩(wěn)定的.

圖4 當(dāng)β＜β1時(shí)，方程（15）的正平衡點(diǎn)漸近穩(wěn)定Fig.4 Stability of positive equilibrium of Equation（15）is asymptotically stable whenβ＜β1

3 結(jié)論

從本文的分析可以看出，蛋白質(zhì)的有限合成速度以及調(diào)控蛋白在細(xì)胞核與細(xì)胞質(zhì)之間運(yùn)輸所引起的時(shí)滯對(duì)基因調(diào)控模型的正平衡態(tài)的穩(wěn)定性有很大的影響.在離散時(shí)滯的情況下，出現(xiàn)了無限個(gè)Hopf分支點(diǎn)，當(dāng)時(shí)滯大于第一個(gè)分支點(diǎn)，正平衡點(diǎn)失去穩(wěn)定性，在平衡點(diǎn)附近存在穩(wěn)定的小振幅周期振蕩.在分布時(shí)滯的情況下，出現(xiàn)了兩個(gè)Hopf分支點(diǎn)，在這兩個(gè)分支點(diǎn)所確定的區(qū)間外邊，正平衡點(diǎn)是穩(wěn)定的，而在這兩個(gè)分支點(diǎn)所確定的區(qū)間內(nèi)正平衡點(diǎn)是不穩(wěn)定的，但有穩(wěn)定的小振幅周期解存在.

對(duì)于有關(guān)生物振蕩現(xiàn)象出現(xiàn)的條件：① 有負(fù)調(diào)控作用將反應(yīng)中的物質(zhì)不斷拉向起始水平；② 負(fù)調(diào)控作用必須有足夠的延遲，使化學(xué)反應(yīng)不會(huì)趨于一個(gè)絕對(duì)穩(wěn)定的狀態(tài)；③ 動(dòng)力學(xué)作用規(guī)則必須足夠的非線性，使物質(zhì)的穩(wěn)定狀態(tài)瓦解；④ 化學(xué)反應(yīng)的反應(yīng)和生成速率必須在一個(gè)合適的時(shí)間尺度內(nèi)進(jìn)行，使整個(gè)調(diào)控網(wǎng)絡(luò)共同出現(xiàn)振蕩現(xiàn)象.結(jié)合本文模型的具體環(huán)境，可以發(fā)現(xiàn)條件①由于mRNA和蛋白質(zhì)之間存在的負(fù)調(diào)控作用顯然是成立的.離散和分布時(shí)滯的引入恰好使條件②得到滿足.米氏函數(shù)和希爾函數(shù)呈現(xiàn)明顯的非線性，因此條件③滿足.模型中推導(dǎo)的，對(duì)于周期解的出現(xiàn)，系數(shù)所必須滿足的條件則對(duì)應(yīng)條件④.

［1］ Goodwin B C. Oscillatory behavior in enzymatic control processes［J］.Advances in Enzyme Regulation，1965，3：425.

［2］ Novák B，Tyson J J.Design principles of biochemical oscillators［J］.Nature Reviews Molecular Cell Biology，2008，9：981.

［3］ 雷錦志.系統(tǒng)生物學(xué)——建模，分析，模擬 ［M］.上海：上海科學(xué)技術(shù)出版社，2010.

LEI Jinzhi.Systems biology—modeling，analysis，simulation［M］.Shanghai：Shanghai Science and Technology Press，2010.

［4］ Klipp E，Herwig R，Kowald A，et al.Systems biology in practice：concepts，implementation and application ［M］.Berlin：Wiley-VCH，2005.

［5］ Chow S N，Hale J K.Methods of bifurcation theory［M］.New York：Springer，1982.

［6］ Hale J K，Verduyn Lunel S M.Introduction to functional differential equations［M］.New York：Springer，1993.

［7］ Murray J D.Mathematical biology［M］.Berlin：Springer，1993.