布爾網(wǎng)絡(luò)動態(tài)行為研究*

王向紅， 王 欣， 劉莉莉， 許 鵬， 劉文斌

(1．溫州大學(xué)物理與電子信息學(xué)院，浙江 溫州 325035;2．溫州職業(yè)技術(shù)學(xué)院，浙江溫州 325035)

0 引言

1969年，Kauffman［1］提出了著名的布爾網(wǎng)絡(luò)(Boolean Network，BN)模型，用于研究基因調(diào)控網(wǎng)絡(luò)和細(xì)胞分化過程．他將基因的表達(dá)分為2種形式：表達(dá)(1)和不表達(dá)(0)，通過布爾函數(shù)描述基因之間的相互作用．這種網(wǎng)絡(luò)模型簡單，對數(shù)據(jù)的數(shù)量和質(zhì)量要求較低，能夠反映系統(tǒng)運(yùn)行過程中的動態(tài)特性，為網(wǎng)絡(luò)控制中出現(xiàn)的非線性動力學(xué)行為提供了很好的近似．因此，布爾網(wǎng)絡(luò)成為研究基因調(diào)控網(wǎng)絡(luò)的重要模型．目前，布爾網(wǎng)絡(luò)已經(jīng)廣泛應(yīng)用于酵母細(xì)胞周期表達(dá)、哺乳動物細(xì)胞周期表達(dá)、果蠅體節(jié)極性網(wǎng)絡(luò)、花發(fā)育形態(tài)表達(dá)等不同生物的基因調(diào)控網(wǎng)絡(luò)的研究［2-6］．本文首先介紹布爾網(wǎng)絡(luò)及動態(tài)行為的相關(guān)概念，采用計(jì)算機(jī)模擬節(jié)點(diǎn)數(shù)n=20的布爾網(wǎng)絡(luò)，用消除趨勢波動分析(DFA)方法研究其長程相關(guān)性動態(tài)行為的演化．

1 布爾網(wǎng)絡(luò)的概念

布爾網(wǎng)絡(luò)是由一個(gè)包含n個(gè)基因的節(jié)點(diǎn)集{σ1，σ2，…，σn}和一個(gè)布爾函數(shù)列表{f1，f2，…，fn}組成．每個(gè)基因的狀態(tài)σi∈{0，1}(i=1，2，…，n)是一個(gè)二值變量，它在t+1時(shí)刻的狀態(tài)完全由t時(shí)刻基因 σj1，σj2，…，σjki的值根據(jù)布爾函數(shù) fi：{0，1}ki→{0，1}確定．因此，t+1 時(shí)刻狀態(tài)可以寫為

式(1)中：ki表示第i個(gè)基因的連通度;σi代表節(jié)點(diǎn)i的狀態(tài);σi=1表示節(jié)點(diǎn)i表達(dá);σi=0表示節(jié)點(diǎn)i不表達(dá)．基因的狀態(tài)是由基因的連通度和布爾函數(shù)共同決定的．由網(wǎng)絡(luò)的n個(gè)節(jié)點(diǎn)的狀態(tài){σ1，σ2，…，σn}所構(gòu)成的向量稱為布爾網(wǎng)絡(luò)的狀態(tài)，由n個(gè)基因組成的網(wǎng)絡(luò)，其可能的狀態(tài)數(shù)為2n．布爾網(wǎng)絡(luò)的狀態(tài)分為暫態(tài)和吸引子．由于布爾網(wǎng)絡(luò)內(nèi)在規(guī)則的方向確定性以及布爾網(wǎng)絡(luò)狀態(tài)的有限性，系統(tǒng)進(jìn)入某些狀態(tài)后將會無限期地停留在這些狀態(tài)上，如果沒有擾動，系統(tǒng)就很難跳出這些0狀態(tài)，我們將這些狀態(tài)稱為吸引子．系統(tǒng)具體進(jìn)入哪個(gè)吸引子是由系統(tǒng)演化的初始狀態(tài)決定的［7］．吸引子外圍，最終能夠演化到特定吸引子的狀態(tài)，被稱為吸引子的吸引域．一般來說，一個(gè)系統(tǒng)可以有多個(gè)吸引子，如圖1所示，一個(gè)節(jié)點(diǎn)數(shù)n=5的隨機(jī)布爾網(wǎng)絡(luò)共有25=32個(gè)狀態(tài)，5個(gè)吸引子環(huán)，系統(tǒng)從任何一個(gè)狀態(tài)開始演化，最終都將進(jìn)入到吸引子00000，00100，10011，11111和吸引子環(huán)11110和11010．

圖1 包含5個(gè)節(jié)點(diǎn)的隨機(jī)布爾網(wǎng)絡(luò)的狀態(tài)空間［8］

1995年，Kauffman［8］使用布爾網(wǎng)絡(luò)研究控制網(wǎng)絡(luò)時(shí)提出單個(gè)吸引子通常對應(yīng)細(xì)胞的某種顯型或細(xì)胞分化過程的不同階段．1999年，哈佛大學(xué)的Huang Sui［9］提出了布爾網(wǎng)絡(luò)的吸引子代表系統(tǒng)的一種記憶形式．隨后Huang Sui又提出了對吸引子的另外一種理解——吸引子代表細(xì)胞的狀態(tài)，如增生、凋亡、分化等．2006年，Huang Sui等［10］通過對同一白血病細(xì)胞株HL-60分別加入DMSO和atRA，7 d后2個(gè)細(xì)胞群體達(dá)到相同的狀態(tài)，從而在實(shí)驗(yàn)上證明了細(xì)胞群狀態(tài)演化中吸引子的存在．

系統(tǒng)的狀態(tài)一般分布在不同的吸引子區(qū)域．隨著時(shí)間的演化，系統(tǒng)會根據(jù)初始狀態(tài)和布爾規(guī)則，逐漸從暫態(tài)向吸引子演化．在有擾動的情況下，系統(tǒng)的狀態(tài)還會在不同吸引子區(qū)域之間動態(tài)躍遷，顯示出復(fù)雜的動態(tài)．對其進(jìn)行分析和研究可以再現(xiàn)基因調(diào)控網(wǎng)絡(luò)系統(tǒng)復(fù)雜的動力學(xué)過程．

2 布爾網(wǎng)絡(luò)動態(tài)行為特性

研究系統(tǒng)的動態(tài)主要是研究系統(tǒng)狀態(tài)在一個(gè)非常小的擾動下的演化行為，即擾動后系統(tǒng)是否會回到原來的吸引子或者一個(gè)穩(wěn)定的狀態(tài)分布．2004年，Shmulevich等［11］提出了布爾網(wǎng)絡(luò)靈敏度(Sensitivity)的概念，這一概念的提出為以后的布爾網(wǎng)絡(luò)動態(tài)特性的研究提供了極大的方便．

當(dāng)2個(gè)狀態(tài)向量其漢明距離為1，即2個(gè)向量僅有1位不同時(shí)，為漢明鄰接向量．布爾網(wǎng)絡(luò)靈敏度Sf(y)定義為：向量y的所有漢明鄰接向量的布爾函數(shù)f輸出值與向量y的布爾函數(shù)f(y)輸出值不同的鄰接向量個(gè)數(shù)．表述為

式(2)中：ei為單位向量，其第i位為1，其他位為0;χ［A］為指示函數(shù)，當(dāng)且僅當(dāng)A為真時(shí)χ［A］值為1．在向量y的某種分布下，平均靈敏度(average sensitivity)Sf表示為

平均靈敏度Sf的值介于0和K之間，通常0＜Sf＜1，Sf=1和Sf＞1分別對應(yīng)系統(tǒng)的有序(Order)、臨界(Criticality)和無序(Chaos)［7］狀態(tài)［12］．布爾網(wǎng)絡(luò)靈敏度是一種衡量網(wǎng)絡(luò)動態(tài)行為的參數(shù)，表征網(wǎng)絡(luò)對一個(gè)單位的漢明擾動的平均傳播特征，即該擾動是逐漸消失?還是保持不變?還是逐漸擴(kuò)大［13］?

目前對隨機(jī)布爾網(wǎng)絡(luò)的動態(tài)行為的研究方法主要分為2類：1)Derrida曲線法;2)長程相關(guān)性．下面將詳細(xì)介紹隨機(jī)布爾網(wǎng)絡(luò)的這2種網(wǎng)絡(luò)動態(tài)行為的研究方法．

2．1 Derrida 曲線法

設(shè) ya(t)=(σa1，σa2，…，σan)和 yb(t)=(σb1，σb2，…，σbn)是時(shí)刻 t網(wǎng)絡(luò)的2 個(gè)隨機(jī)選擇的網(wǎng)絡(luò)狀態(tài)，定義歸一化的漢明距離［14］為

式(4)中⊕是異或運(yùn)算．給定一個(gè)隨機(jī)布爾網(wǎng)絡(luò)，ρ(t+1)是ya(t)和yb(t)后繼狀態(tài)的漢明距離，ρ(t+1)隨ρ(t)變化的曲線稱為Derrida曲線．

運(yùn)行在有序狀態(tài)的網(wǎng)絡(luò)，其相近的狀態(tài)趨向于收斂到同一狀態(tài)，Derrida曲線位于主對角線以下(見圖2)．初始很小的漢明距離ρ(t)到ρ(t+1)數(shù)值的增大意味著Derrida曲線位于主對角線之上，其斜率在原點(diǎn)附近大于1．這意味著一個(gè)很小的基因擾動最終導(dǎo)致完全不同的狀態(tài)，即網(wǎng)絡(luò)對初始狀態(tài)非常敏感——這正是無序狀態(tài)的特征之一．臨界狀態(tài)的網(wǎng)絡(luò)，其斜率在原點(diǎn)附近等于1．對于較小的ρ(t)，Derrida曲線位于主對角線之上的面積越大，網(wǎng)絡(luò)就越無序．

對于一個(gè)基因數(shù)目為n、連通性為K、偏斜度為p的隨機(jī)布爾網(wǎng)絡(luò)，當(dāng)K2p(1-p)＞1時(shí)，系統(tǒng)為無序;當(dāng)K2p(1-p)＜時(shí)，系統(tǒng)為有序;當(dāng)K2p(1-p)=1時(shí)，系統(tǒng)為臨界，所以臨界連通度Kc為

圖2 一個(gè)與有序、臨界和無序動態(tài)特性對應(yīng)的Derrida曲線的示意圖

圖3為式(5)所確定的臨界曲線，無序狀態(tài)占據(jù)了大部分的參數(shù)空間．對大部分的偏斜度p，臨界連通度Kc都非常?。畬τ谶B通度Kc較大的網(wǎng)絡(luò)，如果希望其處于有序狀態(tài)，就必須調(diào)整偏斜度p使其遠(yuǎn)離無偏狀態(tài)．

2．2 長程相關(guān)性

1925年，Johnson在“低頻電路中的肖特基效應(yīng)”中首次提出了“1/f波動”這一術(shù)語．20世紀(jì)70年代，日本著名物理學(xué)家武者利光從混沌學(xué)的領(lǐng)域經(jīng)過多次測量，得出了如下結(jié)論：1/f波動之所以令人感覺舒適，主要因?yàn)槿嗽诎察o狀態(tài)下的心跳周期的波動規(guī)律及α腦波的波動規(guī)律與1/f波動規(guī)律相吻合［15-17］．雖然這種波動對基礎(chǔ)研究和臨床試驗(yàn)都很重要，但是產(chǎn)生這種波動的原因仍然未知［14，18］．

圖3 隨機(jī)布爾網(wǎng)絡(luò)中，由連通度K和偏斜度p確定的臨界曲線

2004 年，Amaral等［19］在 Kauffman 的隨機(jī)布爾網(wǎng)絡(luò)模型和Wolfram的細(xì)胞自動機(jī)模型的基礎(chǔ)上，研究了不同的布爾規(guī)則和連通度對一維環(huán)狀結(jié)構(gòu)網(wǎng)絡(luò)動態(tài)行為的影響并通過添加噪聲研究了網(wǎng)絡(luò)的魯棒性．Heneghan等［20］使用消除趨勢波動分析(DFA)方法量化系統(tǒng)的長程相關(guān)性．圖4為其模擬計(jì)算結(jié)果，得出了系統(tǒng)在不同的布爾規(guī)則下長程相關(guān)性不同，只有在某些特定的規(guī)則下(如232規(guī)則)才會出現(xiàn)1/f波動．他們對隨機(jī)布爾網(wǎng)絡(luò)的動態(tài)行為進(jìn)行了簡單分析，得出布爾網(wǎng)絡(luò)并不能產(chǎn)生1/f低頻功率譜衰變這一結(jié)論．布爾網(wǎng)絡(luò)的動態(tài)演化是否能產(chǎn)生1/f低頻功率譜衰變對基因調(diào)控具有重要的意義．因此，運(yùn)用計(jì)算機(jī)動態(tài)模擬布爾網(wǎng)絡(luò)的演化過程，并用DFA方法進(jìn)一步研究分析布爾網(wǎng)絡(luò)的動態(tài)特性．

圖4 系統(tǒng)在不同布爾規(guī)劃下演化的長程相關(guān)性

3 結(jié)果與討論

在本次布爾網(wǎng)絡(luò)計(jì)算機(jī)模擬中，筆者選取的節(jié)點(diǎn)數(shù)n=20，時(shí)間尺度范圍是10≤t≤103．圖5中縱坐標(biāo)表示靈敏度，其范圍是0．7≤S≤1．3，最小分度為0．02，共有31個(gè)數(shù)據(jù)點(diǎn);橫坐標(biāo)表示添加噪聲范圍0．000 5≤p≤0．02，最小分度為0．000 5，共有40個(gè)數(shù)據(jù)點(diǎn)．即圖5的面積子圖由40×31=1 240個(gè)數(shù)據(jù)點(diǎn)構(gòu)成．對每一個(gè)靈敏度，筆者建立20個(gè)網(wǎng)絡(luò)，每個(gè)網(wǎng)絡(luò)在同一噪聲下演化20次．同一靈敏度、同一噪聲下表征網(wǎng)絡(luò)長程相關(guān)性值α是20×20=400個(gè)數(shù)值的平均．2個(gè)子圖分別表示連通度k=3和k=4的布爾網(wǎng)絡(luò)長程相關(guān)性的變化情況．

從圖5可以看出，布爾網(wǎng)絡(luò)模型的長程相關(guān)性的動態(tài)變化是比較有規(guī)律的．首先，隨著噪聲的增加，絕大多數(shù)網(wǎng)絡(luò)的長程相關(guān)性值逐漸減小;其次，絕大多數(shù)布爾網(wǎng)絡(luò)演化過程能夠出現(xiàn)我們所期望的1/f波動范圍0．9≤α≤1．1，雖然與Amaral等模擬的網(wǎng)絡(luò)一樣均出現(xiàn)了1/f波動，但是噪聲范圍不同;最后，從圖5還可以看出，不同連通度呈現(xiàn)出來的1/f波動區(qū)域?qū)挾炔煌琸=3時(shí)的寬度要比k=4時(shí)的寬，說明連通度低的網(wǎng)絡(luò)更具有穩(wěn)定性，這與圖3中反映的網(wǎng)絡(luò)狀態(tài)與平均連通度和偏斜度的關(guān)系是一致的．

圖5 布爾網(wǎng)絡(luò)長度相關(guān)性值α與噪聲及靈敏度的關(guān)系

圖6 表示靈敏度對布爾網(wǎng)絡(luò)動態(tài)行為的影響，可以看到，隨著布爾函數(shù)靈敏度數(shù)值的增加，網(wǎng)絡(luò)的長程相關(guān)性數(shù)值整體往下偏移，網(wǎng)絡(luò)呈現(xiàn)的1/f波動區(qū)域向右移動，即由無序狀態(tài)漸漸進(jìn)入有序狀態(tài)．這是由于靈敏度是一種衡量網(wǎng)絡(luò)動態(tài)行為的參數(shù)，表征布爾網(wǎng)絡(luò)對微小擾動的敏感程度，而長程相關(guān)性值表征的是系統(tǒng)對其動態(tài)演化狀態(tài)的記憶行為．所以，靈敏度數(shù)值較低的網(wǎng)絡(luò)，其穩(wěn)定性相對較高，外界的擾動較難影響網(wǎng)絡(luò)的演化，網(wǎng)絡(luò)較易記憶原有的演化軌跡，故長程相關(guān)性較強(qiáng)．靈敏度數(shù)值較高的網(wǎng)絡(luò)，網(wǎng)絡(luò)穩(wěn)定性相對較弱，外界微小的擾動就能引起系統(tǒng)狀態(tài)較大的波動，網(wǎng)絡(luò)較易脫離原有的演化軌跡，所以長程相關(guān)性較弱．也就是說有序網(wǎng)絡(luò)產(chǎn)生的1/f波動范圍比無序網(wǎng)絡(luò)要寬，抗干擾能力強(qiáng)．

圖6 布爾網(wǎng)絡(luò)的長程相關(guān)性與靈敏度的關(guān)系

4 結(jié)束語

作為一種數(shù)學(xué)模型，布爾網(wǎng)絡(luò)已經(jīng)廣泛應(yīng)用于基因調(diào)控網(wǎng)絡(luò)的建模、動態(tài)行為分析及網(wǎng)絡(luò)干預(yù)的研究．本文研究了在噪聲擾動下，布爾網(wǎng)絡(luò)的長程相關(guān)性．結(jié)果表明：隨著噪聲的增大，大多數(shù)的網(wǎng)絡(luò)都能經(jīng)過從布朗噪聲到1/f，最后到白噪聲的相變過程．由于已有很多研究生命系統(tǒng)具有1/f特征，這一研究表明在微觀層次上基因調(diào)控網(wǎng)絡(luò)的運(yùn)行也具有1/f特征．當(dāng)然，不同網(wǎng)絡(luò)的1/f的噪聲范圍各不相同，筆者下一步將主要研究影響1/f噪聲范圍的原因及1/f這一特征與網(wǎng)絡(luò)干預(yù)方面的關(guān)系．

［1］Kauffman S A．Metabolic stability and epigenesis in randomly constructed genetic nets［J］．J Theoret Biol，1969，22(3)：437-467．

［2］Li Fangting，Long Tao，Lu Ying，et al．The yeast cell-cycle network is robustly designed［J］．Proceedings of the National Academy of Sciences of the United States of America，2004，101(14)：4781-4786．

［3］Faure A，Naldi A，Chaouiya C，et al．Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle［J］．Bioinformatics，2006，22：124-131．

［4］Albert R，Othmer H G．The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster［J］．Journal of Theoretical Biology，2003，223(1)：1-18．

［5］Espinosa-Soto C，Padilla-Longoria P，Alvarez-Buylla E R．A gene regulatory network model for cell-fate determination during Arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles［J］．The Plant Cell，2004，16(11)：2923-2939．

［6］Chen Lan，Cai Lun，Geir S，et al．Analysis of transcriptional regulation collaboration networks in saccharomyces cerevisiae［J］．Progress in Biochemistry and Biophysics，2008，35(1)：35-42．

［7］Shmulevich I，Dougherty E R．基因組信號處理［M］．劉文斌，高琳，譯．北京：科學(xué)出版社，2010．

［8］Kauffman S A．Search strategies for applied molecular evolution［J］．Journal of Theoretical Biology，1995，173(4)：427-440．

［9］Huang Sui．Gene expression profiling，genetic networks，and cellular states：an integrating concept for tumorigenesis and drug discovery［J］．Journal of Molecular Medicine，1999，77(6)：469-480．

［10］Huang Sui，Ingber D E．A non-genetic basis for cancer progression and metastasis：self-organizing attractors in cell regulatory networks［J］．Breast Disease，2006，26：27-54．

［11］Shmulevich I，Kauffman S A．Activities and sensitivities in Boolean network models［J］．Physical Review Letters，2004，93(4)：048701．

［12］Goldberger A L，Amaral L A N，Haussdorf 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．

［13］Shmulevich I，Kauffman S A．Activities and sensitivities in Boolean network models［J］．Physical Review Letters，2004，93(4)：08701．

［14］Derrida B，Pomeau Y．Random networks of automata：a simple annealed approximation［J］．Europhys Letters，1986，1(2)：45-49．

［15］Peng C K，Havlin S，Stanley H E，et al．Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time-series［J］．Chaos，1995，5(1)：82-87．

［16］Ivanov P C，Nunes-Amaral L A，Goldberger A L，et al．Multifractality in human heartbeat dynamics［J］．Nature，1999，399(6735)：461-465．

［17］Amaral L A N，Ivanov P C，Aoyagi N，et al．Behavioral-independence features of complex heartbeat dynamics［J］．Physical Review Letters，2001，86(26)：6026-6029．

［18］Lipsitz L A．Dynamics of stability：The physiologic basis of functional health and frailty［J］．J Gerontol A Biol Sci Med Sci，2002，57(3)：B115-B125．

［19］Amaral L A N，Moreira A A，Goldberger A L，et al．Emergence of complex dynamics in a simple model of signaling networks［J］．Proceedings of the National Academy of Sciences of the United States of America，2004，101(44)：15551-15555．

［20］Heneghan C，McDarby G．Establishing the relation between detrended fluctuation analysis and power spectral density analysis for stochastic processes［J］．Physical Review Letters，2000，62(5)：6103-6110．