高元明

（溫州大學 物理與電子信息工程學院，浙江 溫州 325035）

概率布爾網(wǎng)絡的動態(tài)行為研究

高元明

（溫州大學 物理與電子信息工程學院，浙江 溫州 325035）

概率布爾網(wǎng)絡是布爾網(wǎng)絡的拓展。布爾網(wǎng)絡和概率布爾網(wǎng)絡已經(jīng)廣泛應用于生物系統(tǒng)模擬中。運用馬爾科夫鏈模擬概率布爾網(wǎng)絡的動態(tài)行為及采用消除趨勢波動分析（DFA）方法，分析概率布爾網(wǎng)絡的長程相關性特性發(fā)現(xiàn)，概率布爾網(wǎng)絡的動態(tài)行為具有長程相關性，大部分概率布爾網(wǎng)絡能呈現(xiàn)1/f波動，且產(chǎn)生1/f波動的噪聲值和范圍比布爾網(wǎng)絡的大。

概率布爾網(wǎng)絡；長程相關性；1/f波動；穩(wěn)態(tài)分布；DFA

0 引 言

生理過程的一個普遍特征是其產(chǎn)生的信號能呈現(xiàn)長程相關性，且偏離該相關性就被認為與疾病或老化等生理問題有關[1-2]。這些信號的相關性可以采用消除波動趨勢分析（D F A）方法分析且以一個標度指數(shù)α來衡量。1969年Kauffman提出了著名的布爾網(wǎng)絡（Boolean Network，BN）模型[3]，在此基礎上，2002年Shmulevich等提出了概率布爾網(wǎng)絡（Probabilistic Boolean Network，PBN）模型[4-5]。布爾網(wǎng)絡和概率布爾網(wǎng)絡已經(jīng)廣泛應用于基因調(diào)控網(wǎng)絡及生物系統(tǒng)的模擬。在文獻[6]和[7]中發(fā)現(xiàn)，布爾網(wǎng)絡的動態(tài)行為具有長程相關性且能呈現(xiàn)1/f波動，運用網(wǎng)絡中狀態(tài)到最大吸引子的距離刻畫網(wǎng)絡的演化軌跡；概率布爾網(wǎng)絡中沒有吸引子，因而運用馬爾科夫鏈刻畫網(wǎng)絡的動態(tài)行為。研究結(jié)果發(fā)現(xiàn)，概率布爾網(wǎng)絡的動態(tài)行為具有長程相關性，大部分概率布爾網(wǎng)絡能呈現(xiàn)1/f波動，且產(chǎn)生1/f波動的噪聲值和范圍比布爾網(wǎng)絡的大。

1 相關概念

1.1 布爾網(wǎng)絡

一個布爾網(wǎng)絡B(V,F)包括n個節(jié)點（基因）V={x1, x2,…,xn}和n個布爾函數(shù)F={f1, f2,…, fn}。每個基因可以取0或1兩個值，當xi=0時，表示該基因不表達；當xi=1時，表示該基因表達。每個基因?qū)衚i個輸入，當基因的全部ki個輸入和函數(shù)fi確定時，基因在t+1時刻的狀態(tài)可以由xi(t+1)=fi(x1(t),x2(t),…,xki(t))確定。在布爾網(wǎng)絡中，假設全部基因表達值更新是同步的，布爾網(wǎng)絡的狀態(tài)空間可以用一個二進制向量x=(x1(t),x2(t),…,xn(t))表示。給定任意一個初始狀態(tài)u0，讓其進行演化，最終會停留在一個狀態(tài)或一個狀態(tài)環(huán)，這樣的狀態(tài)（環(huán)）稱為吸引子（環(huán)），演化所經(jīng)過的狀態(tài)構(gòu)成了該吸引子（環(huán)）的吸引域。一個布爾網(wǎng)絡可以有多個吸引子，但至少必須有一個。Kauffman認為，一個吸引子對應一種細胞模型或功能[3]。

1.2 概率布爾網(wǎng)絡

概率布爾網(wǎng)絡中對于每一個節(jié)點xi，都有一個布爾函數(shù)集合Fi={fj(i)}j=1…l(i)。其中，fj(i)表示節(jié)點xi的第j個函數(shù)，l(i)表示xi擁有的函數(shù)個數(shù)，即任意基因xi的值的更新方式為：xi(t+1)=fj(i)(x1(t),x2(t),…,xki(t))。fj(i)通常也稱為預測器，每一個函數(shù)fj(i)都有對應的一個選擇概率 ，其中 。在任意給定的時刻，PBN中狀態(tài)轉(zhuǎn)移是由一個布爾函數(shù)向量決定的。整個概率布爾網(wǎng)絡可以看成是由 個布爾網(wǎng)絡構(gòu)成的，PBN中狀態(tài)轉(zhuǎn)移有N個可能的實現(xiàn)方式，網(wǎng)絡中就有N種諸如fk=(fk(11),fk(22),…,fkn

n)的布爾函數(shù)向量。在概率布爾網(wǎng)絡演化過程中，由一個選擇參數(shù)q決定所選擇的布爾網(wǎng)絡是否改變。當q＜1時，則稱上下文敏感PBN；當q=1時，則稱暫態(tài)隨機PBN。在每一時刻概率布爾網(wǎng)絡演化過程中都要對布爾網(wǎng)絡進行選擇。

由于基因組處于一個開放的體系中，會受到外部環(huán)境的影響，基因在外部的刺激下，如在誘變劑和熱應力等的作用下，基因能夠由激活變?yōu)橐种苹蚍催^來由抑制變?yōu)榧せ?。在概率布爾網(wǎng)絡模型中，這種擾動的作用效果是假設有一個小的概率p＞0，使基因在某一時刻的值由0變?yōu)?或由1變?yōu)?。在有擾動的情況下，網(wǎng)絡狀態(tài)各態(tài)可達。概率布爾網(wǎng)絡演化的動態(tài)行為可以運用馬爾科夫鏈描述，對于概率布爾網(wǎng)絡來說，是否存在穩(wěn)態(tài)分布是一個很重要的方面。如果存在穩(wěn)態(tài)分布，則從長期行為角度來說，網(wǎng)絡處于任意的狀態(tài)x(i)的概率與初始狀態(tài)無關。Shmulevich等已經(jīng)證明，在網(wǎng)絡存在擾動的情況下，概率布爾網(wǎng)絡存在穩(wěn)態(tài)分布π={π0,…,π2n-1}，且此時穩(wěn)態(tài)分布與平穩(wěn)分布一致，因而有π=πAr，其中Ar為r步的轉(zhuǎn)移概率[4-5，8]。

1.3 靈敏度

用總RNA提取試劑盒提取樣本總RNA。用TIANScript RT Kit對所得RNA樣本反轉(zhuǎn)錄成cDNA，反應條件：70 ℃ 5 min，42 ℃ 50 min，95 ℃ 5 min。β-actin為內(nèi)參基因，APP695和β-actin引物信息見表1。PCR反應條件：95 ℃ 10 min預變性，95 ℃ 10 s，60 ℃ 20 s，72 ℃ 30 s,進行 45個循環(huán)。熒光定量儀進行熒光定量分析。利用2-△△CT方法計算目的基因的相對表達量。

靈敏度（S）表示網(wǎng)絡中任意一個基因受到擾動后能夠影響其他基因的平均個數(shù)。當S＜1時，表示網(wǎng)絡中的擾動隨著時間的推移會消失，此時網(wǎng)絡處于有序階段；當S＞1時，表示網(wǎng)絡中的擾動隨著時間的推移會擴散，此時網(wǎng)絡處于無序階段。介于有序和無序階段的稱為臨界階段，此時對應于S=1，表示擾動的大小在網(wǎng)絡中會一直持續(xù)相同的大小。給定一個概率布爾網(wǎng)絡的穩(wěn)態(tài)分布π，網(wǎng)絡的靈敏度

，其中，xj,k=(x1,…,xj-1,k,xj+1,…,xn)，k=0,1[9-10]。

1.4 有序布爾函數(shù)

有序布爾函數(shù)可以使網(wǎng)絡演化向有序階段進行，以下為其中兩種有序布爾函數(shù)：

（1）渠化函數(shù)。一個布爾函數(shù)f要是渠化函數(shù)需滿足的條件是：如果存在兩個變量u,v∈{0,1}，使得對所有該函數(shù)的輸入x1,x2,…,xk∈{0,1}都有當xi=u時，有f(x1,x2,…,xk)=v，則稱布爾函數(shù)f:{0,1}→{0,1}對第i個輸入變量時是渠化的，渠化值為v，渠化變量為u。一個渠化函數(shù)中至少有一個渠化變量，但可以有多個渠化變量。當渠化函數(shù)的渠化變量取渠化值時，函數(shù)的輸出與其他變量無關。

（2）post類函數(shù)。定義1：當布爾函數(shù)f所有輸出為1的向量中，有一位均為1，則函數(shù)屬于A∞類。類似地，將定義中的1換為0，則稱該函數(shù)屬于a∞類。定義2：對于μ≥2，如果任意μ個使函數(shù)f取值為1的向量，都存在一個相同的分量為1，則稱該函數(shù)屬于Aμ類。類似地，將定義中的1換為0，對應的函數(shù)為aμ[11]。

2 消除趨勢波動分析（DFA）方法

DFA方法是在分形自相似理論基礎上發(fā)展起來的，其本質(zhì)是隨機行走理論的一種修正的均方根分析方法。該方法已被證明是一種研究非線性時間序列的長程相關性的極有效的方法，能很好地去除時間序列中波動趨勢的影響，且以一個標度指數(shù)α來衡量序列的相關性。

對于一個時間序列x(i)，i=1,…,N，采用DFA方法進行分析時，其計算步驟描述如下：

（2）把累積離差Y(i)分成大小為s的不重疊的Ns個等間隔區(qū)間，取Ns=N/S的整數(shù)部分，對每個區(qū)間分別采用最小二乘法進行擬合，得到每個區(qū)間的局部趨勢函數(shù)，即多項式pv(i)，k=1,…,m。

（3）對于每個區(qū)間，消除趨勢波動后得到新時間序列，記為Ys(i)，即Ys(i)=Y(i)-pv(i)。

（5）按照一定的規(guī)律改變s值，重復上述計算過程，得到一系列s～F(s)的變化值。若時間序列是相關的，則F(s)與s為冪律關系，即F(s)∝sα。在雙對數(shù)坐標(log10(F(s)), log10(s))中采用最小二乘法擬合，其直線部分的斜率即為標度指數(shù)α。當α=0.5時，表示序列處于白噪聲階段；當α=1.5時，表示序列處于布朗噪聲階段；1/f波動對應于α=1，實際取值0.9≤α≤1.1[12]。

3 模擬模型與處理方法

3.1 網(wǎng)絡的構(gòu)造

為系統(tǒng)地研究概率布爾網(wǎng)絡的長程相關性，按照靈敏度和布爾函數(shù)的種類構(gòu)造出基因個數(shù)為n=10，連通度k=4的由一個布爾網(wǎng)絡構(gòu)成的概率布爾網(wǎng)絡。構(gòu)造網(wǎng)絡時，靈敏度范圍取值為0.7≤S≤1.3，間隔為0.1。按照函數(shù)種類分別隨機構(gòu)造出屬于隨機布爾函數(shù)、渠化函數(shù)及A2∪a2類函數(shù)三類網(wǎng)絡。對每一類函數(shù)在每個靈敏度上構(gòu)造600個網(wǎng)絡。

3.2 網(wǎng)絡的量化及處理方法

在對概率布爾網(wǎng)絡進行刻畫過程中，假設網(wǎng)絡的穩(wěn)態(tài)分布為初始狀態(tài)，網(wǎng)絡每演化一步，用該時刻網(wǎng)絡狀態(tài)的分布與到穩(wěn)態(tài)分布的距離來描述演化過程。假設網(wǎng)絡在演化過程中每個時刻都有隨機的噪聲η，即網(wǎng)絡每一步演化都有新的噪聲作用上去，且噪聲的作用方式為對網(wǎng)絡每個時刻的狀態(tài)分布π(t)中的值產(chǎn)生微小的波動。噪聲可以表示為：η=(v0,…,v2n-1)T，且噪聲取值為0.5%0%，間隔為0.5%。

網(wǎng)絡狀態(tài)演化的軌跡可以表示為：y(t)=│π(t)-π(0)│。其中，│.│表示兩個分布之間的距離，π(t)表示網(wǎng)絡狀態(tài)在時刻t的分布，π(t)可以由π(t)=π(t-1)·P得到。而π(t-1)=(π(t-1)+η)/(1+│η│)，│η│表示η的模，即噪聲總量；P為轉(zhuǎn)移概率矩陣。網(wǎng)絡演化的軌跡可以由y(t)所構(gòu)成的時間序列表示。

為分析y(t)自相關函數(shù)的波動特性，采用DFA方法分析由y(t)生成的時間序列。分析y(t)構(gòu)成的時間序列時，時間尺度范圍取值為10≤t≤500，按照公比為100.05間隔選取時間步長，共有近25個點。對不同時間步長進行線性擬合時，選取雙對數(shù)坐標，擬合直線的斜率就是描述長程相關性的指數(shù)α值。最終得出α值為進行10次的平均值。圖1中，左圖表示運用DFA方法分析y(t)的相關性時在不同噪聲下的直線擬合情況；右圖表示在三個不同噪聲（η=1%,η=2.5%, η=5.5%）下y(t)構(gòu)成的序列波動情況，即網(wǎng)絡演化的軌跡。

圖1 在三個不同噪聲下網(wǎng)絡演化的軌跡與D FA方法分析的擬合

4 結(jié)論與討論

通過采用DF A方法研究概率布爾網(wǎng)絡的長程相關性發(fā)現(xiàn)，其長程相關性與布爾網(wǎng)絡的相似，指數(shù)α隨噪聲η的改變方式主要有兩種，即在相同的噪聲作用下，隨著噪聲的增大，長程相關性值呈現(xiàn)由高到底的過程（A類）和長程相關性值幾乎不變（B類），其中A類呈現(xiàn)1/f波動過程，如圖2所示。此外，研究表明，概率布爾網(wǎng)絡產(chǎn)生1/f波動的噪聲值和范圍比布爾網(wǎng)絡的大得多[6]，這可能是由概率布爾網(wǎng)絡自身所具有的不確定性和隨機性，導致概率布爾網(wǎng)絡在小的噪聲下呈現(xiàn)出更強的魯棒性。

考慮到1/f波動具有重要意義，運用能夠產(chǎn)生1/f波動的網(wǎng)絡作相關的相位圖，由一個BN構(gòu)成的PBN，如圖3所示。不同函數(shù)構(gòu)成的網(wǎng)絡產(chǎn)生的1/f波動所跨過的噪聲值和范圍不同，由Post類函數(shù)和渠化函數(shù)構(gòu)成的網(wǎng)絡的1/f波動所跨過的噪聲值和范圍比隨機函數(shù)的寬些。

由圖3可知，網(wǎng)絡在相同的噪聲情況下，由隨機布爾函數(shù)和渠化函數(shù)構(gòu)成的網(wǎng)絡由無序向有序過程中，α值會變大，說明網(wǎng)絡的敏感度越小，受到的影響就越小，即如果網(wǎng)絡處于有序階段時，網(wǎng)絡組織化

圖2 概率布爾網(wǎng)絡的長程相關性

程度高，受到噪聲時相對不容易改變原有結(jié)構(gòu)，因而表現(xiàn)出更強的相關性。而由A2∪a2類函數(shù)構(gòu)成的網(wǎng)絡沒有此類現(xiàn)象，但該類函數(shù)構(gòu)成的網(wǎng)絡產(chǎn)生1/f波動的噪聲要比其他兩類的大，說明由此類函數(shù)構(gòu)成的生物網(wǎng)絡具有更強的魯棒性。

5 結(jié)束語

本文研究了概率布爾網(wǎng)絡的長程相關性，結(jié)果表

圖3 概率布爾網(wǎng)絡的長程相關性相位圖

[1]Goldberger A L，Amaral L A N，Hausdorff J M，et al.Fractal dynamics in physiology:alterations with disease and aging[J]. Proceedings of the National Academy of Sciences of the United States of America，2002，99(1):2466-2472.

[2]Lipsitz L A.Dynamics of stability the physiologic basis of functional health and frailty[J].The Journals of Gerontology，Series A:Biological Sciences，2002，57(3):115-125.

[3]Kauffman S A.The origin of order:self-organization and selection in evolution[M].Oxford：Oxford University Press，1993. [4]Shmulevich I，Dougherty E R，Kim S C，et al.Probabilistic Boolean networks:a rule-based uncertainty model for gene regulatory networks[J].Bioinformatics，2002，18(2):261-274.

[5]Shmulevich I，Dougherty E R，Zhang W.Gene perturbation and intervention in probabilistic Boolean networks[J].Bioinformatics，2002，18(10):1319-1331.

[6]Xu P，Wang X H，Liu W B.The influence of the basin structure of Boolean networks on their long range correlated dynamics [C].IEEE 6th International Conference on Systems Biology(ISB)，2012:17-21.

[7]Liu L L，Xu P，Wang X，et al.The Study of the long-range correlated behaviors of Boolean networks[J].Information，2013.

[8]Shmulevich I，Dougherty E R.Genomic signal processing，princeton series in applied mathematics[M].Princeton:Princeton University Press，2007.

[9]Vahedi G，Shmulevich I，Liu W B，et al.Critical Boolean networks require minimal intervention rate[C].IEEE International Workshop on GENSIPS，2009:1-2.

[10]Shmulevich I，Kauffman S A.Activities and sensitivities in Boolean network models[J].Physical Review Letter，2004(93): 48071-48074.

[11]Shmulevich I，Lahdesmaki H，Dougherty E R，et al.The role of certain post classes in Boolean network models of genetic networks[J].Proceedings of the National Academy of Sciences of the United States of America，2003，100(19):10734-10739.

[12]Hu K，Ivanov P C，Chen Z，et al.Effect of trends on detrended fluctuation analysis[J].Physical Review E，2001(64): 11114-11133.

[責任編輯：王向紅]

Study on Dynamic Behavior of Probabilistic Boolean Networks

GAO Yuanming
(College of Physics and Electronic Information Engineering, Wenzhou University, Wenzhou, 325035, China)

The probabilistic Boolean networks are the extension of Boolean network. Both of them are widely used in biological system simulation. By simulating the dynamic behavior of probabilistic Boolean network by applying Markov chain and analyzing the long range correlation of probabilistic Boolean networks by applying detrended fluctuation (DFA), it shows that the dynamic behavior of probabilistic Boolean networks has long range correlation with most of the probabilistic Boolean network presenting 1/f fluctuation and the noise and the range produced bigger than those of the Boolean network.

Probabilistic Boolean networks; Long range correlation; 1/f fluctuation; Steady-state distribution; DFA

O631

A

1671-4326(2013)02-0058-04

2013-04-02

高元明（1988—），男，江西撫州人，溫州大學物理與電子信息工程學院碩士研究生.

明，網(wǎng)絡演化過程具有長程相關性，大部分概率布爾網(wǎng)絡能呈現(xiàn)1/f波動，且產(chǎn)生1/f波動的噪聲值和范圍比布爾網(wǎng)絡的大。由于生命活動能產(chǎn)生1/f波動，這一研究表明，在微觀層次上基因調(diào)控網(wǎng)絡的運行也具有1/f波動特征，且驗證了概率布爾網(wǎng)絡模型用于基因調(diào)控網(wǎng)絡的可行性,為基因調(diào)控網(wǎng)絡的建模和模擬等提供參考平臺。