宋永利，葉 子

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

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

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

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

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

式中：x表示mRNA的濃度；y表示蛋白質(zhì)濃度.蛋白質(zhì)濃度的增加會對mRNA的合成起抑制作用，這里用系數(shù)為p的Hill函數(shù)表示，而蛋白質(zhì)也會對自身的合成起抑制作用，這里用 Michaelis-Menten函數(shù)表示.另一方面，mRNA濃度的增加會促進蛋白質(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)化為更簡潔的數(shù)學形式，令k1S＝m1，（Kd）2＝m2，kdx＝m3，ksy＝m4，k2ET＝m5，Km＝m6，則方程組可 轉(zhuǎn)化為以下方程組：

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

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

1.1 離散時滯

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

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

其中

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

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

分離實部虛部得到

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

計算判別式得

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

上式化簡后即為

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

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

化簡上式得到

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

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

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

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

（3）τ＝τj是方程組（5）的 Hopf分支點.

1.2 分布時滯

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

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

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

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

引入輔助變量u

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

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

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

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

由于D1＞0恒成立而且

所以，方程（20）所有根皆具有嚴格負實部，等價于D2＞0，即

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

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

消去ω2，即得

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

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

由式（25），得到

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

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

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

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

（3）當β＝β1或β＝β2時，方程（20）有一對純虛根±iω；當β∈（β1，β2）時，方程（20）有一對具正實部的根.

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

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

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

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

2 數(shù)值模擬

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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