基于連續(xù)-離散擴(kuò)展卡爾曼建模的p53-Mdm2調(diào)控關(guān)系研究

姜 莉 李玉榕*

1(福州大學(xué)電氣工程與自動(dòng)化學(xué)院，福州 350116)2(福州大學(xué)福建省醫(yī)療器械和醫(yī)藥技術(shù)重點(diǎn)實(shí)驗(yàn)室，福州 350116)

基于連續(xù)-離散擴(kuò)展卡爾曼建模的p53-Mdm2調(diào)控關(guān)系研究

姜 莉1,2李玉榕1,2*

1(福州大學(xué)電氣工程與自動(dòng)化學(xué)院，福州 350116)2(福州大學(xué)福建省醫(yī)療器械和醫(yī)藥技術(shù)重點(diǎn)實(shí)驗(yàn)室，福州 350116)

生物醫(yī)學(xué)研究揭示，腫瘤抑制基因p53的表達(dá)水平與腫瘤的形成相關(guān)。p53相關(guān)的信號(hào)轉(zhuǎn)導(dǎo)網(wǎng)絡(luò)研究為揭示癌癥、腫瘤的發(fā)病機(jī)制和尋找其治療方法提供新思路，而p53和癌基因Mdm2的負(fù)反饋關(guān)系是該信號(hào)轉(zhuǎn)導(dǎo)網(wǎng)絡(luò)的核心部分。因此，研究p53-Mdm2的調(diào)控關(guān)系意義重大。基于電離輻射后人類白血病細(xì)胞中測(cè)得的時(shí)間序列基因表達(dá)數(shù)據(jù)，利用連續(xù)-離散的擴(kuò)展卡爾曼(EKF)算法，建立帶有延遲特性的非線性動(dòng)態(tài)連續(xù)隨機(jī)生物網(wǎng)絡(luò)模型，模擬電離輻射后的p53-Mdm2動(dòng)態(tài)調(diào)控過程。同時(shí)，對(duì)所建立網(wǎng)絡(luò)的準(zhǔn)確性進(jìn)行驗(yàn)證，可得模型的誤差率為0.85%。結(jié)果表明：所提出的算法是收斂的，可以預(yù)測(cè)基因的表達(dá)水平，準(zhǔn)確地模擬p53-Mdm2調(diào)控網(wǎng)絡(luò)受到電離輻射后的衰減震蕩響應(yīng)過程，該過程和生物實(shí)驗(yàn)結(jié)果一致。所提出的算法為生物學(xué)實(shí)驗(yàn)建模提供一種有效的方法，并可以為網(wǎng)絡(luò)系統(tǒng)動(dòng)力學(xué)特性的研究提供條件。

p53-Mdm2網(wǎng)絡(luò)；連續(xù)-離散擴(kuò)展卡爾曼(EKF)算法；時(shí)間序列數(shù)據(jù)；延遲

引言

p53是重要的腫瘤抑制因子，在細(xì)胞周期抑制、細(xì)胞凋亡、生長停滯、抑制血管生成等過程中具有重大的作用[1]，且與腫瘤的形成和治療有密切的關(guān)系。研究表明，許多癌癥中都發(fā)現(xiàn)了p53功能異常，超過50%的癌癥是由p53基因突變引起的[2]。在癌癥理療中，人工恢復(fù)p53的活性，可以使癌細(xì)胞停止增長或凋亡，來達(dá)到治療的目的[3-4]。因此，越來越多的學(xué)者開始關(guān)注p53的研究。

p53對(duì)其下游多個(gè)基因的表達(dá)進(jìn)行調(diào)控，其中p53和Mdm2的負(fù)反饋調(diào)節(jié)可以看成p53整個(gè)轉(zhuǎn)導(dǎo)通道的核心，p53和Mdm2的表達(dá)失衡會(huì)引起器官或組織中的病理變化[3]。研究p53和Mdm2蛋白功能及其之間的相互調(diào)節(jié)機(jī)制不僅可以找到新的癌癥治療方法，而且在某些疾病的預(yù)防方面也有極其重要的作用。在正常細(xì)胞中p53蛋白的含量很少，在細(xì)胞受到外界刺激時(shí)，如電離輻射(ionization radiation，IR)、紫外線輻射(ultraviolet radiation，UV)等，會(huì)引起DNA的損傷，隨后激活p53，使得p53表達(dá)增多，從而引起一系列復(fù)雜的反應(yīng)。該過程中p53-Mdm2調(diào)節(jié)過程是一個(gè)復(fù)雜且精細(xì)的過程，而生物實(shí)驗(yàn)在這方面存在著很多局限。例如：反應(yīng)過程中的微觀數(shù)據(jù)很難測(cè)量；可供實(shí)驗(yàn)的樣本數(shù)量有限；某些反應(yīng)過快或過慢，難以實(shí)現(xiàn)測(cè)量；依靠實(shí)驗(yàn)尋找靶點(diǎn)過于盲目等。而數(shù)學(xué)建模不僅可以模擬得到生物實(shí)驗(yàn)中由于技術(shù)原因不能測(cè)得的數(shù)據(jù)，而且可以從微觀上來分析生物反應(yīng)過程中的可能響應(yīng)，為生物實(shí)驗(yàn)提供更多的可能靶點(diǎn)。

出于技術(shù)、經(jīng)費(fèi)、人力的考慮，有學(xué)者提出通過建立數(shù)學(xué)模型來模擬p53-Mdm2負(fù)反饋的網(wǎng)絡(luò)模型[3-7]。對(duì)生物網(wǎng)絡(luò)的建模能夠從機(jī)制上研究生命系統(tǒng)的一些現(xiàn)象，如細(xì)胞凋亡，生長停滯等[8]。這些學(xué)者的做法都是先建立能夠符合生物特性的模型，然后通過數(shù)值分析方法，對(duì)模型的動(dòng)態(tài)特性進(jìn)行分析。這類研究在確定模型參數(shù)的時(shí)候，一般是比較盲目地去嘗試，或者根據(jù)經(jīng)驗(yàn)來取得，而且這些模型只能表示一類動(dòng)態(tài)網(wǎng)絡(luò)，并不能夠針對(duì)某一個(gè)器官或者組織給出特定的模型參數(shù)。而隨著基因芯片技術(shù)和高通量技術(shù)的發(fā)展，可以結(jié)合生物學(xué)實(shí)驗(yàn)和數(shù)學(xué)模型進(jìn)行研究。通過實(shí)驗(yàn)在有限的時(shí)間點(diǎn)上檢測(cè)細(xì)胞受到電離輻射后p53和Mdm2的表達(dá)量，而后利用這些時(shí)間序列的表達(dá)數(shù)據(jù)對(duì)模型進(jìn)行狀態(tài)和參數(shù)辨識(shí)，從而確定模型參數(shù)。

擴(kuò)展卡爾曼(extended kalman Filter，EKF)算法因其能辨識(shí)非線性模型而被廣泛應(yīng)用。Wang等建立基因調(diào)控網(wǎng)絡(luò)的狀態(tài)空間模型[9]，通過EKF算法來同時(shí)辨識(shí)網(wǎng)絡(luò)的參數(shù)和狀態(tài)，建立動(dòng)態(tài)的基因網(wǎng)絡(luò)，并對(duì)網(wǎng)絡(luò)的穩(wěn)定性能進(jìn)行分析[9]。Xiong等利用EKF線性化的特點(diǎn)，提出一種基于魯棒性狀態(tài)估計(jì)的算法，來對(duì)生物網(wǎng)絡(luò)的結(jié)構(gòu)進(jìn)行辨識(shí)[10]。Zhang等提出利用一種分?jǐn)?shù)的擴(kuò)展卡爾曼算法，對(duì)于建立的非線性基因調(diào)控網(wǎng)絡(luò)進(jìn)行建模[11]。以上建模研究中均將生物網(wǎng)絡(luò)看作離散系統(tǒng)，而生物的響應(yīng)過程本質(zhì)上是一個(gè)連續(xù)的過程，若能夠考慮到其連續(xù)的性質(zhì)，則會(huì)更加符合生物網(wǎng)絡(luò)的特性，得到更準(zhǔn)確的生化反應(yīng)動(dòng)態(tài)模型。

本研究首先根據(jù)生物調(diào)控過程的連續(xù)性、延遲性、隨機(jī)性和非線性特性，建立一個(gè)連續(xù)的帶有延遲的非線性動(dòng)態(tài)隨機(jī)生物網(wǎng)絡(luò)模型，來模擬電離輻射后p53-Mdm2的響應(yīng)過程。而后利用生物學(xué)實(shí)驗(yàn)得到的電離輻射后測(cè)得的時(shí)間序列表達(dá)數(shù)據(jù)，基于連續(xù)-離散的EKF算法[12]，同時(shí)辨識(shí)網(wǎng)絡(luò)的參數(shù)和基因表達(dá)水平。最后，對(duì)網(wǎng)絡(luò)的準(zhǔn)確性進(jìn)行估算，并對(duì)所建立模型的響應(yīng)過程進(jìn)行分析。本研究提出的算法具有通用性，可以得到特定網(wǎng)絡(luò)的更加準(zhǔn)確的模型參數(shù)，為進(jìn)一步理解p53調(diào)節(jié)通道提供了可能，從而為癌癥治療提供更加準(zhǔn)確的靶點(diǎn)。

1 材料與方法

基于時(shí)間序列的基因表達(dá)數(shù)據(jù)，綜合考慮生化反應(yīng)過程的連續(xù)性、時(shí)滯性、隨機(jī)性和非線性，首先提出符合生物特性的微分模型；然后利用連續(xù)-離散擴(kuò)展卡爾曼算法對(duì)于模型的狀態(tài)和參數(shù)進(jìn)行聯(lián)合辨識(shí)；最后利用誤差表達(dá)式對(duì)于已經(jīng)建立的模型進(jìn)行評(píng)價(jià)。

1.1 數(shù)據(jù)樣本

采用的源數(shù)據(jù)取自文獻(xiàn)[13]中：在室溫下，對(duì)人類白血病細(xì)胞(MOTL4)進(jìn)行5Gy的γ輻射后，分別獲得細(xì)胞內(nèi)0、2、4、6、8、10、12 h的基因表達(dá)水平；為了減少測(cè)量中的誤差，獨(dú)立進(jìn)行了3次數(shù)據(jù)測(cè)量過程，其均值和標(biāo)準(zhǔn)偏差如圖1、2所示。本研究選取3組數(shù)據(jù)的平均值來進(jìn)行建模。

圖1 p53的平均值和標(biāo)準(zhǔn)差Fig.1 Mean and standard deviation of p53

圖2 Mdm2的平均值和標(biāo)準(zhǔn)差Fig.2 Mean and standard deviation of Mdm2

1.2 p53-Mdm2網(wǎng)絡(luò)模型

p53-Mdm2負(fù)反饋網(wǎng)絡(luò)尚有許多細(xì)節(jié)需要完善，為了不失一般性，本研究采用簡化模型,不考慮mRNA的作用，即把轉(zhuǎn)錄和翻譯合成一項(xiàng)來建模。生物網(wǎng)絡(luò)的模型不僅要求能夠符合生物反應(yīng)機(jī)理，更要符合生物網(wǎng)絡(luò)的調(diào)節(jié)特性[3]。而基因在調(diào)控的過程中，即轉(zhuǎn)錄和翻譯過程中，都會(huì)發(fā)生延遲的現(xiàn)象[2]。p53-Mdm2動(dòng)態(tài)調(diào)控網(wǎng)絡(luò)中，p53對(duì)于Mdm2的調(diào)控也存在時(shí)滯，若不考慮到延遲的特性，必然會(huì)影響基因網(wǎng)絡(luò)的穩(wěn)定性，甚至導(dǎo)致模型錯(cuò)誤的預(yù)測(cè)響應(yīng)過程。因此，本研究在文獻(xiàn)[5]中的模型基礎(chǔ)上，考慮到生物過程具有延遲的特性[2,7]，把轉(zhuǎn)錄和翻譯過程中的延遲整合到Mdm2表達(dá)式中的一個(gè)項(xiàng)中來表示。p53、Mdm2分別表示相應(yīng)的蛋白水平，其調(diào)控網(wǎng)絡(luò)表示如下：

(1)

式中：第1項(xiàng)r1是p53蛋白本底的合成速率；第2項(xiàng)表示除了Mdm2外其他原因引起p53的降解，d1是本底的降解速率；由于p53轉(zhuǎn)導(dǎo)信號(hào)網(wǎng)絡(luò)具有非常多的正反饋回路[14]，第3項(xiàng)描述了p53系統(tǒng)中的正反饋效應(yīng)，正反饋可以被認(rèn)為是一種廣義的“自誘導(dǎo)”，模型中用一個(gè)Hill函數(shù)來描述，vp表示正反饋導(dǎo)致的本底生成速率，k1是Hill函數(shù)的常數(shù)；第4項(xiàng)表示Mdm2誘導(dǎo)的p53的降解作用；kmax表示Mdm2對(duì)p53負(fù)調(diào)控的降解速率，kp是p53和Mdm2之間的絡(luò)合物離解常數(shù)[3]。

(2)

式中：第1項(xiàng)r2表示Mdm2的合成速率；第2項(xiàng)表示Mdm2的本底降解，d2表示本底降解速率；第3項(xiàng)表示p53誘導(dǎo)Mdm2的生成，vm表示p53引起的Mdm2的增長速率，τ表示調(diào)控過程存在的滯后，k2是Hill函數(shù)的常數(shù)；第4項(xiàng)表示Mdm2介導(dǎo)的自降解作用，mdeg表示Mdm2的自降解的速率，km表示Mdm2離解常數(shù)。

1.3 連續(xù)-離散的EKF算法建模

1.3.1 連續(xù)-離散的EKF算法

EKF算法在每個(gè)時(shí)間步長用一個(gè)線性化的過程來把非線性系統(tǒng)近似為時(shí)變的線性系統(tǒng)，再用傳統(tǒng)的卡爾曼算法對(duì)時(shí)變的線性系統(tǒng)進(jìn)行濾波，是卡爾曼算法在非線性系統(tǒng)中的應(yīng)用。EKF是基于最小方差的次優(yōu)算法，并不是最優(yōu)的算法，但其效果良好并成為處理非線性濾波問題中最經(jīng)典、應(yīng)用最廣泛的方法[9]。詳細(xì)介紹EKF算法可參考文獻(xiàn)[9，15]。EKF算法的收斂性討論可以參考文獻(xiàn)[16-17]。EKF算法已經(jīng)在很多時(shí)候應(yīng)用在參數(shù)識(shí)別中[9-11]，一般使用離散的EKF算法。本研究根據(jù)模型連續(xù)的特性，以及使用的時(shí)間序列基因表達(dá)數(shù)據(jù)離散的特點(diǎn)，采用連續(xù)-離散的EKF算法來對(duì)模型中的參數(shù)和p53、Mdm2蛋白的表達(dá)水平進(jìn)行辨識(shí)。此外，考慮到基因調(diào)控過程會(huì)遇到各種噪聲，如附近基因表達(dá)對(duì)于本網(wǎng)絡(luò)中基因表達(dá)信號(hào)測(cè)量的干擾、由基因芯片放置位置不同而產(chǎn)生的噪聲、測(cè)量時(shí)溫度變化引起的噪聲、基因芯片使用時(shí)產(chǎn)生的測(cè)量噪聲等，在網(wǎng)絡(luò)中加入噪聲的影響。具體過程如下：

將p53和Mdm2的表達(dá)水平分別用x1、x2表示，并在模型中加入噪聲，于是式(1)、(2)可以寫成

(3)

(4)

式中，w1(t)、w2(t)為互不相關(guān)的均值為零的高斯白噪聲。

令x(k)=[x1(t)x2(t)]T，則式(3)、(4)可以簡化寫為

(5)

為了使用EKF算法來辨識(shí)參數(shù)和基因表達(dá)水平，把需要辨識(shí)的參數(shù)擴(kuò)展到系統(tǒng)的狀態(tài)變量中。令

最終用連續(xù)-離散EKF算法辨識(shí)參數(shù)的模型寫為

式中，y(t)為基因芯片含噪聲的測(cè)量值，V(t)、W(t)為互不相關(guān)的均值為零的高斯白噪聲，其協(xié)方差矩陣分別為Q(t)、R(t)，k為采樣點(diǎn)。

1.3.2 算法的具體迭代步驟

假設(shè)采樣點(diǎn)的個(gè)數(shù)為N。

步驟1：k=0時(shí)，設(shè)置初始值。

P(0,0)=E[(X(0)-X0)(X(0)-X0)T]=PX0

步驟2：k=1，2，3，4，…，N時(shí)，計(jì)算如下：

1)時(shí)間更新(預(yù)測(cè))，其中t∈[tk-1,tk]。

狀態(tài)時(shí)間更新

(9)

(10)

誤差協(xié)方差時(shí)間更新

(11)

P(k|k-1)=P(tk)

(12)

2)測(cè)量更新(矯正) 。

計(jì)算卡爾曼增益矩陣，有

Kk=P(k|k-1)CkT[CkP(k|k-1)CkT+Rk]-1

(13)

用測(cè)量值更新估計(jì)值，有

(14)

誤差協(xié)方差更新，有

P(k|k)=(I-KkCk)P(k|k-1)

(15)

其中

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

1.3.3 算法的具體實(shí)現(xiàn)步驟

不斷重復(fù)步驟1～4，直至采樣點(diǎn)全部用完。

1.4 模型評(píng)價(jià)

為了能夠定量化地衡量辨識(shí)出的模型質(zhì)量，引入誤差指標(biāo)，即觀測(cè)到的值與模型實(shí)際預(yù)測(cè)值之間的誤差百分比，模型誤差算式[9,18]具體表示為

(25)

式中：l指的是所有基因的數(shù)目；s是觀測(cè)點(diǎn)的數(shù)量，即采樣的點(diǎn)數(shù)；yck是指c基因在k時(shí)刻的實(shí)際表達(dá)水平。

2 結(jié)果

本算法是基于Matlab軟件。為探索電離輻射后網(wǎng)絡(luò)的動(dòng)態(tài)過程，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)數(shù)越多，則辨識(shí)的結(jié)果越準(zhǔn)確。而本文第1.1節(jié)中的原始數(shù)據(jù)只有7個(gè)時(shí)間點(diǎn)，使用三次樣條差值方法可以得到足夠多的點(diǎn)，以避免出現(xiàn)過度擬合的情況[19]。所以，通過三次樣條插值方法取25個(gè)采樣點(diǎn)，每30 min一次。此外，由于實(shí)驗(yàn)測(cè)得的數(shù)據(jù)過少，若辨識(shí)全部的參數(shù)難免會(huì)出現(xiàn)誤差較大。而文獻(xiàn)[5]中指出，r1(p53的合成速率)、vm(負(fù)反饋強(qiáng)度)對(duì)于系統(tǒng)動(dòng)力學(xué)影響較大，故本研究辨識(shí)式(3)、(4)中的參數(shù)r1、vm值，其余參數(shù)值根據(jù)文獻(xiàn)[5]設(shè)置(τ取為5 min)。

2.1 辨識(shí)結(jié)果

利用連續(xù)-離散擴(kuò)展卡爾曼算法和圖1、2中的實(shí)驗(yàn)數(shù)據(jù)，對(duì)p53-Mdm2模型進(jìn)行辨識(shí)，估計(jì)出模型參數(shù)和狀態(tài)的值，分別如圖3、4所示。圖5、6所示分別為連續(xù)-離散卡爾曼算法辨識(shí)的參數(shù)和狀態(tài)的估計(jì)誤差方差。

圖3 模型參數(shù)r1,vm估計(jì)值Fig.3 Estimated parameters of r1,vm

圖4 p53,Mdm2表達(dá)水平的估計(jì)值和實(shí)驗(yàn)值Fig.4 Actual and estimated expression level of p53,Mdm2

圖5 參數(shù)r1,vm的估計(jì)誤差方差Fig.5 Variances of estimated parameters r1,vm

圖6 狀態(tài)p53、Mdm2的估計(jì)誤差方差Fig.6 Variances of estimated states p53,Mdm2

從圖3可以看出，經(jīng)過迭代濾波，參數(shù)r1和vm穩(wěn)定收斂，確定出r1=0.794 7，vm=0.124 7。從圖4中看出，p53和Mdm2的預(yù)測(cè)值和實(shí)驗(yàn)值之間誤差較小。同時(shí)，圖5、6顯示參數(shù)和狀態(tài)的估計(jì)誤差方差都很小，除了開始時(shí)期，其余時(shí)刻接近于零，這意味著該方法能準(zhǔn)確地從時(shí)間序列數(shù)據(jù)中估計(jì)出模型的狀態(tài)和參數(shù)(見式(3)、(4))。

2.2 評(píng)價(jià)結(jié)果

表1列出了在不同的純滯后參數(shù)τ下所建立模型對(duì)應(yīng)的誤差率。

表1 不同τ的模型誤差率Tab.1 Error ratio of the model in different τ

由表1可看出，模型取不同純滯后時(shí)間τ,誤差率均較低，利用本算法辨識(shí)出的模型是準(zhǔn)確的。

3 討論

在研究中，基于連續(xù)-離散卡爾曼算法辨識(shí)參數(shù)和狀態(tài)，從而建立調(diào)控模型。基于該模型，可以預(yù)測(cè)12 h之后p53以及Mdm2的表達(dá)量。圖7表示了p53和Mdm2前12 h的估計(jì)值和之后的預(yù)測(cè)值，左邊的坐標(biāo)軸表示p53的表達(dá)水平，右邊坐標(biāo)系表示Mdm2的表達(dá)水平。

圖7 p53和Mdm2的響應(yīng)過程Fig.7 The response process of p53 and Mdm2

從p53和Mdm2的響應(yīng)過程中可以看出，白血病細(xì)胞在受到電離輻射后，p53比Mdm2先達(dá)到峰值。p53第一個(gè)脈沖寬度大約為250 min，達(dá)到第一個(gè)脈沖峰值為200 min左右；而Mdm2的第一個(gè)脈沖寬度大約為600 min，達(dá)到第一個(gè)脈沖峰值為400 min。生物實(shí)驗(yàn)的結(jié)果表明[20]，細(xì)胞受到電離輻射時(shí)，會(huì)進(jìn)行振蕩，振蕩脈沖的寬度為350±160 min，第一個(gè)脈沖峰值對(duì)應(yīng)的時(shí)間相差很大，大約在損傷后的360±240 min?？梢姳灸Ｐ头治鏊玫降慕Y(jié)果很好地符合了生物學(xué)實(shí)驗(yàn)結(jié)果。

此外，從圖7中可以看出，電離輻射后，p53的濃度增加，促進(jìn)Mdm2的生成，當(dāng)Mdm2的濃度增加到一定程度后，反過來抑制Mdm2的生成。并且可以看出，p53和Mdm2的響應(yīng)過程是一個(gè)衰減震蕩的過程，最后p53和Mdm2的濃度基本不變，達(dá)到了一個(gè)平衡的狀態(tài)。已有研究表明[21-23]，電離輻射會(huì)引起DNA損傷，在輕度的DNA損傷中，p53表現(xiàn)出脈沖動(dòng)力學(xué)特性，這一模式具有很好的靈活性，既有利于細(xì)胞周期停滯期間損傷的有效恢復(fù)，又可避免輕度損傷細(xì)胞的過度凋亡。在達(dá)到平衡狀態(tài)后，說明此時(shí)DNA的損傷已經(jīng)被修復(fù)。此外，震蕩的過程即為DNA的修復(fù)過程，p53發(fā)出第一個(gè)脈沖后，等待看損傷是否被正確修復(fù)，若沒有就會(huì)產(chǎn)生第二個(gè)脈沖，直到損傷有效修復(fù)，p53的濃度不再變化[5,22]。從圖7中可以看出，大約800 min細(xì)胞完成DNA的修復(fù)過程。

本算法去除了離散擴(kuò)展卡爾曼算法不符合生物過程連續(xù)性的缺點(diǎn)，不僅考慮到了實(shí)驗(yàn)得到的微陣列數(shù)據(jù)的離散性，而且考慮到了生物調(diào)控的連續(xù)性。建立模型后，通過圖5、6和表1對(duì)辨識(shí)結(jié)果進(jìn)行評(píng)價(jià)，說明建立模型的準(zhǔn)確性較高。而圖7進(jìn)一步表明，該辨識(shí)的模型可以對(duì)之后的基因表達(dá)進(jìn)行預(yù)測(cè)，且辨識(shí)得到的結(jié)果和預(yù)測(cè)結(jié)果均具有實(shí)際意義，最后得到的調(diào)控過程基本符合生物學(xué)實(shí)驗(yàn)結(jié)果，為廣大系統(tǒng)動(dòng)力學(xué)研究工作者提供了一種新的建模方式。

4 結(jié)論

本研究討論了p53-MDM2網(wǎng)絡(luò)在受到電離輻射后的響應(yīng)過程。首先，根據(jù)生物過程延遲和非線性的特性，并依據(jù)生物學(xué)實(shí)驗(yàn)，建立生物網(wǎng)絡(luò)的非線性隨機(jī)動(dòng)態(tài)模型。一般模型中的參數(shù)依據(jù)經(jīng)驗(yàn)試湊確定，而且無法針對(duì)某個(gè)個(gè)體或某個(gè)組織、器官的特異性以及不同疾病模型來確定。本研究根據(jù)生物響應(yīng)連續(xù)的特性，提出連續(xù)-離散的EKF算法，利用實(shí)驗(yàn)獲得的電離輻射后測(cè)的時(shí)間序列基因表達(dá)數(shù)據(jù)，對(duì)模型中的參數(shù)進(jìn)行辨識(shí)，進(jìn)而得到該個(gè)體或者該組織、器官的模型參數(shù)。然后，根據(jù)得到的模型，不僅可以得到p53-Mdm2網(wǎng)絡(luò)在測(cè)量時(shí)間內(nèi)的表達(dá)水平及動(dòng)態(tài)響應(yīng)，而且可以定量地預(yù)測(cè)之后任意時(shí)刻的基因表達(dá)量及動(dòng)態(tài)響應(yīng)。

本算法可以較準(zhǔn)確地估計(jì)出模型的參數(shù)，這樣在研究某一疾病或者腫瘤治療的時(shí)候，就可以準(zhǔn)確得到該疾病患者癌細(xì)胞或腫瘤細(xì)胞中p53的含量以及響應(yīng)過程，為進(jìn)一步?jīng)Q定電離輻射治療的量提供參考，為醫(yī)療工作者提供指導(dǎo)。存在的問題是，由于技術(shù)的限制，研究中使用的原始數(shù)據(jù)的點(diǎn)數(shù)太少。這樣，為了考慮參數(shù)的準(zhǔn)確性，不得不只選擇比較重要的參數(shù)進(jìn)行辨識(shí)，其他參數(shù)用經(jīng)驗(yàn)取值。隨著基因芯片技術(shù)的不斷發(fā)展，肯定可以得到更多時(shí)間點(diǎn)的數(shù)據(jù)，這樣，把模型中的所有參數(shù)辨識(shí)出來成為可能，并且提高其準(zhǔn)確性，為研究網(wǎng)絡(luò)的動(dòng)力學(xué)特性提供基礎(chǔ)，進(jìn)而可更精確地選擇特定癌癥化療中輻射使用量以及更準(zhǔn)確的治療靶點(diǎn)。

[1] Taberlay PC, Jones PA. DNA methylation and cancer [J]. Oncogene, 2012, 21(35): 89-95.

[2] Yang Y, Lee KS, Xiang C, et al. Biological mechanisms revealed by a mathematical model for p53-Mdm2 core regulation [J]. Iet Systems Biology, 2009, 3(4): 229-238.

[3] Voropaeva OF, Shokin YI, Nepomnyaschchikh LM, et al. Mathematical modeling of the tumor markers network [C] // International Conference on Biomedical Engineering and Computational Technologies. Novosibirsk: IEEE, 2015: 225-229.

[4] Wang Shaomeng, Sun Wei, Zhao Yujun, et al. SAR405838: an optimized inhibitor of MDM2-p53 interaction that induces complete and durable tumor regression [J]. Cancer Research, 2014, 74(20): 5855-5865.

[5] 孫延哲. P53信號(hào)轉(zhuǎn)導(dǎo)網(wǎng)絡(luò)動(dòng)力學(xué)模型研究 [D]. 南京: 南京大學(xué), 2014.

[6] 畢遠(yuǎn)宏, 楊卓琴, 何小燕. Mdm2生成速率調(diào)控的p53-Mdm2振子的全局動(dòng)力學(xué)和穩(wěn)定性 [J]. 物理學(xué)報(bào), 2016, 65(2): 360-367.

[7] Mihalas GI, Simon Z, Balea G, et al. Possible oscillatory behavior in p53-Mdm2 interaction computer simulation [J]. Journal of Biological Systems, 2011, 08(1): 21-29.

[8] Lavoie H, Hogues H, Mallick J, et al. Evolutionary tinkering with conserved components of a transcriptional regulatory network [J]. PLoS Biology, 2010, 8(3): e1000329.

[9] Wang Zidong, Liu Xiaohui, Liu Yurong, et al. An extended Kalman filtering approach to modeling nonlinear dynamic gene regulatory networks via short gene expression time series [J]. IEEE/ACM Transactions on Computational Biology & Bioinformatics, 2009, 6(3): 410-419.

[10] Xiong Jie, Zhou Tong. Structure identification for gene regulatory networks via linearization and robust state estimation [J]. Automatica, 2013, 50(11): 2765-2776.

[11] Zhang Yongqing, Pu Yifei, Zhang Haisen, et al. An extended fractional Kalman filter for inferring gene regulatory networks using time-series data [J]. Chemometrics & Intelligent Laboratory Systems, 2014, 138: 57-63.

[12] Kulikova MV, Kulikov GY. Adaptive ODE solvers in extended Kalman filtering algorithms [J]. Journal of Computational & Applied Mathematics, 2014, 262(10): 205-216.

[13] Barenco M, Tomescu D, Brewer D, et al. Ranked prediction of p53 targets using hidden variable dynamic modeling [J]. Genome Biology, 2006, 7(3): R25.

[14] Vilborg A, Glahder JA, Wilhelm MT, et al. The p53 target Wig-1 regulates p53 mRNA stability through an AU-rich element [J]. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(37): 15756-15761.

[15] Hatami E, Salarieh H, Vosoughi N. Design of a fault tolerated intelligent control system for a nuclear reactor power control: using extended Kalman filter [J]. Journal of Process Control, 2014, 24(7): 1076-1084.

[16] Huang S, Dissanayake G. Convergence and consistency analysis for extended Kalman filter based SLAM [J]. IEEE Transactions on Robotics, 2007, 23(5): 1036-1049.

[17] Krener AJ. The convergence of the extended Kalman filter [J]. Mathematics, 2003, 286: 173-182.

[18] Zeng Nianyin, Wang Zidong, Li Yurong, et al. Inference of nonlinear state-space models for Sandwich-Type lateral flow immunoassay using extended Kalman filtering [J]. IEEE Transactions on Biomedical Engineering, 2011, 58(7): 1959-1966.

[19] Tu CT, Chen BS. On the increase in network robustness and decrease in network response ability during the aging process: a systems biology approach via microarray data [J]. IEEE/ACM Transactions on Computational Biology & Bioinformatics, 2013, 10(2): 468-480.

[20] Lahav G, Rosenfeld N, Sigal A, et al. Dynamics of the p53-Mdm2 feedback loop in individual cells [J]. Nature Genetics, 2004, 36(2):147-150.

[21] Bar-Or RL, Maya R, Segel LA, et al. Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study [J]. Proceedings of the National Academy of Sciences of the United States of America, 2000, 97(21): 11250-11255.

[22] 張小鵬, 劉鋒, 王煒. p53信號(hào)網(wǎng)絡(luò)的非線性動(dòng)力學(xué)研究 [J]. 中國科學(xué):物理學(xué)力學(xué)天文學(xué), 2014(12): 1311-1318.

[23] Lahav G. The strength of indecisiveness: oscillatory behavior for better cell fate determination [J]. Journal of Health Services Research & Policy, 2014, 19(1): 1-2.

Research of p53-Mdm2 Regulatory Networks Based on Modeling by Continuous-Discrete Extended Kalman Filter

Jiang Li1,2Li Yurong1,2*

1(CollegeofElectricalEngineeringandAutomation,FuzhouUniversity,Fuzhou350116,China)2(FujianKeyLaboratoryofMedicalInstrumentation&PharmaceuticalTechnology,FuzhouUniversity,Fuzhou350116,China)

The studies of biomedicine show that the expression level of p53 (tumor suppressor gene) is related with the tumor forming. Researches of the p53 related signal transduction networks would provide new ideas for revealing the pathogenesis of cancer or tumor and looking for the treatment method. As the p53-Mdm2 negative feedback network plays core roles, the study of the p53-Mdm2 regulation net work is of great significance. In this paper, based on the gene expression time series data of human leukemia cells after ionizing radiation, the nonlinear dynamic continuous random biological networks with time delay were established, then continuous-discrete extended Kalman filter (EKF) algorithm was used to simulate the dynamic regulation of p53 and Mdm2. Meanwhile, the accuracy of the model established in this paper was validated and we found that the error rate of the model was only 0.85%. Results showed that the algorithm thatis convergent could predict gene expression level at any time and simulate accurately the damped oscillation process of p53-Mdm2 regulatory networks after ionizing radiation. The response process was consistent with biological experiment results. The algorithm proposed in this paper provided an effective method for biological experiment modeling, and supplied the foundation for the research of system dynamics.

p53-Mdm2 network; continuous-discrete extended Kalman filter (EKF) algorithm; time-series data; delay

10.3969/j.issn.0258-8021. 2017. 02.008

2016-08-05， 錄用日期:2016-10-13

國家自然科學(xué)基金(61403319);福建省科技廳國際合作項(xiàng)目(2015I0003);福建省教育廳科技項(xiàng)目(JK2014001)

R318

A

0258-8021(2017) 02-0180-07

*通信作者(Corresponding author)， E-mail: liyurong@fzu.edu.cn