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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

2．1 Derrida 曲線法

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

4 結(jié)束語

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