復(fù)雜網(wǎng)絡(luò)上流行病傳播動力學(xué)的爆發(fā)閾值解析綜述

李睿琪, 王　偉,舒盼盼, 楊　慧,潘黎明,崔愛香,唐　明

(1.電子科技大學(xué)互聯(lián)網(wǎng)科學(xué)中心,成都611731; 2.北京師范大學(xué)系統(tǒng)科學(xué)學(xué)院,北京 100875)

?

復(fù)雜網(wǎng)絡(luò)上流行病傳播動力學(xué)的爆發(fā)閾值解析綜述

李睿琪1, 2, 王偉1,舒盼盼1, 楊慧1,潘黎明1,崔愛香1,唐明1

(1.電子科技大學(xué)互聯(lián)網(wǎng)科學(xué)中心,成都611731; 2.北京師范大學(xué)系統(tǒng)科學(xué)學(xué)院,北京 100875)

摘要：對流行病傳播爆發(fā)閾值的理論解析方法進行總結(jié)，主要介紹平均場、點對近似、主方程、邊滲流、空穴理論、邊劃分以及譜分析這7種常用的動力學(xué)解析方法的前提假設(shè)、具體思路、步驟及其應(yīng)用局限，并且梳理總結(jié)了SIS與SIR模型爆發(fā)閾值的異同。

關(guān)鍵詞：復(fù)雜網(wǎng)絡(luò)；流行病傳播；爆發(fā)閾值；理論解析

0引言

流行病的預(yù)警與干預(yù)作為公共安全的重中之重，已受到社會各界的廣泛關(guān)注[1-2]。然而想要進行流行病傳播規(guī)律的實體研究是不可接受、不可行的，數(shù)理模型自然成為流行病傳播研究方向的希望與重點，通過對數(shù)理模型的分析及與實證的對比，讓我們對于流行病傳播的機制與防控的效率能夠有更為深入而客觀的認識。早在1760年，Daniel Bernoulli就提出了第一個流行病傳播模型[3]；1927年，Kermack和McKendrick確立了流行病傳播的現(xiàn)代數(shù)理模型框架[4]，激勵了無數(shù)科學(xué)家的后續(xù)工作；1992年，Anderson和May的重要著作[1]則是對過去流行病傳播研究工作的精到總結(jié)。自從1998年復(fù)雜網(wǎng)絡(luò)理論誕生以來[5-7]，網(wǎng)絡(luò)思維逐步深入人心。網(wǎng)絡(luò)化的思考方式，可以將流行病傳播過程中的重要因素(例如個體交互的結(jié)構(gòu)、移動性、交互模式)進行很好的建模[7-11]，復(fù)雜網(wǎng)絡(luò)上的傳播動力學(xué)研究成為了流行病傳播研究新的方向和希望。十幾年間，復(fù)雜網(wǎng)絡(luò)上的流行病傳播研究受到了極大關(guān)注，取得了飛躍性的發(fā)展[12-13]。早在年，Pastor-Satorras等率先研究了因特網(wǎng)上的病毒傳播問題[14-15]，通過理論解析發(fā)現(xiàn)在熱力學(xué)極限情況下病毒爆發(fā)的閾值接近于零，這暗示無論多小的傳染概率都將導(dǎo)致病毒長期存在于網(wǎng)絡(luò)中。2004年德國馬普所的Hufnagel等在全球性流行病傳播方面做出了開創(chuàng)性工作[16]，他們利用全球航空網(wǎng)絡(luò)及相關(guān)航班數(shù)據(jù)重現(xiàn)了2003年SARS的時空演化斑圖。美國印第安納大學(xué)的Vespignani等主持了一項名為全球性流行病傳播建模(Global Epidemic and Mobility Modeler)的項目。這個項目基于復(fù)雜網(wǎng)絡(luò)理論，根據(jù)病毒傳播的特點和世界交通的數(shù)據(jù)，計算病毒可能的傳播情況，從而使我們能夠?qū)砜赡馨l(fā)生的疫情進行預(yù)警[17]。令人振奮的是，他們較為準(zhǔn)確地預(yù)測出了2009年甲型H1N1流感的時空演化斑圖[18]。2013年，Brockman等利用航空網(wǎng)絡(luò)的流量數(shù)據(jù)，定義了地點之間的有效距離，將原先觀察到的較為復(fù)雜的流行病傳播演化過程重新映射回反應(yīng)擴散過程，并且可以通過這一方法有效地反推出流行病的源頭[19]。

在流行病研究中，上述這類模型通常被稱作復(fù)合種群模型，它們的基本假設(shè)是將地點看作節(jié)點，將個體看作節(jié)點中的粒子。目前很多研究已經(jīng)足夠精細、足夠定量化，已然能夠為政府決策提供一定建議與支持。然而更為基礎(chǔ)的流行病研究則是將個體看作網(wǎng)絡(luò)中的節(jié)點，分析其上傳播動力學(xué)特征。真實的網(wǎng)絡(luò)其結(jié)構(gòu)通常具有時效性、動態(tài)演化性，然而靜態(tài)網(wǎng)絡(luò)是這些更為復(fù)雜情況的基礎(chǔ)，所以如何能夠?qū)㈧o態(tài)網(wǎng)絡(luò)傳播動力學(xué)有效地分析出來，是理論研究的第一步，也是本文關(guān)注的重點。

為了能夠解釋并預(yù)測流行病傳播中的一些重要問題，許多研究者利用計算機模擬進行了許多相應(yīng)的嘗試；但科學(xué)家們對這樣的結(jié)果往往并不滿意，想要找尋出其內(nèi)在的規(guī)律，用更為準(zhǔn)確的方式來更加精確地進行預(yù)測。另一方面，即便對于靜態(tài)網(wǎng)絡(luò)上傳播動力學(xué)的解析，想要得出一般的解析解也是極其困難甚至不可能的，網(wǎng)絡(luò)層級性、聚集性等復(fù)雜性因素也迫使研究者們發(fā)展相適應(yīng)的數(shù)理解析方法、計算模型。而且值得一提的是，雖然不同的動力學(xué)具體機制不盡相同，流行病傳播動力學(xué)模型及其解析方法還被應(yīng)用、借鑒到了更為廣泛、更為一般的動力學(xué)行為研究上，例如信息傳播、觀點傳播、文化規(guī)范和社會行為、電網(wǎng)級聯(lián)失效[11]等等。

在流行病傳播動力學(xué)的研究中，流行病的爆發(fā)閾值無疑是非常重要的，它對于流行病的評估、預(yù)警、干預(yù)策略的選擇都有極為重大的影響。研究者發(fā)現(xiàn)網(wǎng)絡(luò)上的流行病爆發(fā)這類宏觀涌現(xiàn)行為的研究方法路線與統(tǒng)計物理中的非平衡相變非常接近[20]。受此啟發(fā)，在過去十年，誕生了眾多流行病傳播解析方法：從經(jīng)典的平均場方法到更為嚴(yán)格的定量化數(shù)值解析方法[21-22]。然而目前復(fù)雜網(wǎng)絡(luò)研究中，閾值解析的方法繁多、魚龍混雜，部分文章往往只是對某些細節(jié)進行了考慮、可以適用于某些略為特殊的情況，但卻未必會有較好的拓展性與普適性。我們進行這樣的解析方法綜述，一個重要的目的就是對目前各類繁多的理論方法進行一個梳理，探尋不同方法的解析思路、適用條件、準(zhǔn)確性等重要問題，讓我們對于這些相應(yīng)的方法能有一個更為直觀而全面的認識和理解。

本文將以流行病傳播的解析方法為主要研究對象，力求對傳播動力學(xué)爆發(fā)閾值的解析方法進行較為全面的分析。按照各類方法使用的數(shù)學(xué)工具，將其分為3大類：微分方程法，生成函數(shù)法以及譜分析法。具體介紹平均場(也會在這部分簡介SIS與SIR模型以及流行病傳播對應(yīng)的更為廣義的相變)、點對近似、主方程、邊滲流、空穴理論、邊劃分以及譜分析這7種常用的動力學(xué)解析方法，力求較為全面地覆蓋傳播動力學(xué)研究中重要的解析方法，并探尋各類解析方法內(nèi)在的機制原理、前提假設(shè)、解析思路、物理含義及其局限性,以期為初學(xué)者提供一個較為清晰的學(xué)習(xí)脈絡(luò)，為廣大復(fù)雜網(wǎng)絡(luò)研究人員提供一定的參考。

1平均場方法

1.1平均場理論概述

在復(fù)雜網(wǎng)絡(luò)研究中，進行解析時最為簡潔的方法當(dāng)數(shù)平均場近似。平均場方法計算過程簡單且原理易于理解，基本上只要熟悉常微分方程的人都具備運用平均場方法的能力，這也大大降低了理論解析的門檻。平均場方法被廣泛地應(yīng)用于流行病傳播動力學(xué)解析及預(yù)測[15,23-29]、投票模型[30]以及同步現(xiàn)象[31]的解析當(dāng)中，其解析結(jié)果可較好地反應(yīng)實際情況。概言之，平均場是把環(huán)境對對象的作用進行集體處理，以平均作用效果替代單個作用效果的加和的方法。雖然其精確適用的限制條件非常嚴(yán)格，但是由于其解析思想簡潔、過程簡明，平均場近似被廣泛地應(yīng)用。然而不能完全滿足平均場應(yīng)用的假設(shè)條件時，其理論解析結(jié)果的準(zhǔn)確性和如何進行相應(yīng)的改進等問題就受到了眾人的關(guān)注，也得出了許多重要的結(jié)論[24-29]。本節(jié)兼顧最新結(jié)論對于原始平均場方法的指導(dǎo)與更正，但并不主要討論這些最新的發(fā)現(xiàn)，只在后面部分大略涉及一些關(guān)于平均場最新進展的討論。我們力爭對平均場方法進行一個完整而清晰的描述并給出具體的解析案例，使初學(xué)者對于平均場有一個整體的把握。我們將主要關(guān)注平均場方法在全接觸SIR模型閾值解析中的應(yīng)用，并以小世界網(wǎng)絡(luò)與無標(biāo)度網(wǎng)絡(luò)為例給出具體的解析過程及相應(yīng)分析結(jié)論。

1.2平均場方法提出的背景

平均場理論最初是由朗道為解釋二級相變而提出[32-35]。對于所有的連續(xù)相變，都可以定義一個序參量，使它在臨界點的一邊為零，表示相對的無序狀態(tài)，而在臨界點的另一邊不為零，表示相對的有序狀態(tài)[32-35]。與一級相變(即不連續(xù)相變，如物質(zhì)氣、液、固三態(tài)間的轉(zhuǎn)變)不同，二級相變(連續(xù)相變，如鐵磁相變、超導(dǎo)相變、超流相變、以及合金的有序無序相變等)不具有兩相共存的臨界變化階段和清晰的界面(例如在冰融化這一一級相變過程中，存在冰與水兩種相位共存的情形)，相變在物質(zhì)各處同時地突然發(fā)生。在二級相變過程中，各處同時出現(xiàn)新相與舊相互相嵌套(而非清晰的相位變化界面，此時很難判定它究竟處于哪一個相位)，形成所謂“花斑”[32]。這說明二級相變中，物質(zhì)分子與其近鄰的其它分子如何互相作用、結(jié)合成什么相，主要取決于跨越一切尺度的所有分子之間的相互作用的總體效果，亦即“在連續(xù)相變的臨界點，關(guān)聯(lián)長度趨于無窮”；相比之下，對于一級相變至關(guān)重要的分子的周圍環(huán)境的“局部信息”(如晶格常數(shù)、晶格結(jié)構(gòu)等)對二級相變的影響并不重要。因此，在二級相變中，可以把跨越一切尺度的所有分子之間的相互作用的總體效果等價于一個“平均場”，不去計算局部的、處處不同的相互作用情況——而這正是平均場方法的核心思想。例如對于二維Ising模型的相變求解，在低溫時，自旋關(guān)聯(lián)很大，熱強度漲落較大(亦即偏離平均值較多)，平均場方法很難正確描述相應(yīng)現(xiàn)象[36]。

雖然連續(xù)相變是一個理想的平衡過程(即不存在能量與物質(zhì)的流動，或者說流動極其緩慢、在局部可認為無變化的物理過程)，然而，以廣義的序參量隨驅(qū)動量變化的冪律為標(biāo)志的臨界現(xiàn)象涉及更為廣泛的領(lǐng)域。對于流行病的傳播過程，可以將其等價為一種廣義的相變。在流行病停止傳播的狀態(tài)，節(jié)點的感染密度ρ為零時，系統(tǒng)相對無序；而在流行病傳播的狀態(tài)，ρ大于零，系統(tǒng)相對有序。因此可以把ρ定義為廣義的序參量。對于流行病的傳播，大量的實驗和理論結(jié)果都證明對于有效感染率λ存在一個閾值λc，只有在λ-λc>0時，流行病才能全局傳播，因此可以把λ-λc定義為熱力學(xué)驅(qū)動量。本節(jié)將介紹如何運用平均場近似方法的思想，寫出流行病傳播的動力學(xué)方程，并在一定限制條件下進行求解分析。

1.3SIR模型的假設(shè)及其微分方程表達形式

經(jīng)典的SIR(Susceptible-Infected-Recovered/Removed)與SIS(Susceptible-Infected-Susceptible)模型是Reed和Frost在1920年的一篇未發(fā)表論文中首先提出的，這也是首次開創(chuàng)性地利用微分方程來進行動力學(xué)的描述[1,37]。他為個體定義了3種狀態(tài)：易感態(tài)(Susceptible)、感染態(tài)(Infected)、恢復(fù)態(tài)(Recovered/Removed)，任意個體必處于這3種狀態(tài)中的某一種。SIR模型假設(shè)流行病傳播的周期遠遠小于個體的壽命，所以不考慮節(jié)點的死亡，整個流行病傳播的過程中，網(wǎng)絡(luò)節(jié)點總數(shù)恒定，亦即網(wǎng)絡(luò)為靜態(tài)網(wǎng)絡(luò)。一個易感態(tài)節(jié)點(S)以某一概率λ被任一個感染態(tài)節(jié)點(I)感染, 一個感染態(tài)節(jié)點(I)以概率μ康復(fù)且變?yōu)榻K生免疫的恢復(fù)態(tài)(R)。SIR模型中，節(jié)點狀態(tài)演化過程為S→I→R，或一直保持S狀態(tài)；當(dāng)傳播過程結(jié)束后，只會有S與R兩種狀態(tài)存在。

對于隨機規(guī)則網(wǎng)絡(luò)(節(jié)點度相同)和均勻網(wǎng)絡(luò)(度分布服從正態(tài)或者泊松分布)，設(shè)網(wǎng)絡(luò)中節(jié)點的平均度為〈k〉(它實則是對任意節(jié)點在某一時刻接觸能力的描述，對于單點接觸(Contact Process, 又譯作接觸過程, 但筆者認為單點接觸這一說法更加形象直觀)或者有限接觸能力的傳播情形，只須將此處修改為相應(yīng)數(shù)值即可，稍后在1.6節(jié)中還會詳細介紹)，流行病的感染率為λ，節(jié)點的恢復(fù)概率為μ，用S(t)表示t時刻處于S態(tài)節(jié)點的密度；ρ(t)、R(t)依次對應(yīng)I態(tài)與R態(tài)的節(jié)點密度。則SIR模型的微分方程組如式(1)～(4)[15]。

(1)

(2)

(3)

1=S(t)+ρ(t)+R(t)

(4)

第一個方程其含義是t時刻S態(tài)節(jié)點被感染為I態(tài)的概率，亦即S態(tài)節(jié)點密度減少的變化率。為了得到這樣一個變化率，就必須知道一個S態(tài)節(jié)點會連接到多少個處于感染態(tài)的節(jié)點，以確定其被感染的概率。平均地說來，一個節(jié)點會有〈k〉個鄰居，其鄰居在t時刻處于感染態(tài)的概率為ρ(t)，進而可知其處于感染態(tài)的鄰居平均說來有〈k〉ρ(t)個，所以第一個方程的精確的形式應(yīng)當(dāng)是dS(t)/dt=-S(t)(1-(1-λ)〈k〉ρ(t))；但是這樣一個方程較難處理，將S(t)(1-(1-λ)〈k〉ρ(t))使用泰勒展開后，舍掉λ的二階及以上的高次項而得到的線性近似即可得到S(t)λ〈k〉ρ(t)這一線性表達，而這樣的線性近似就自然而然地會引入誤差，只有在其近似適用范圍內(nèi)(λ〈k〉ρ(t)?1)，才能認為較為精確。

第二個方程是表征t時刻處于I態(tài)節(jié)點密度變化的概率，這當(dāng)中有原先處于S態(tài)被感染后新增的個體，也有由I態(tài)以概率μ恢復(fù)為R態(tài)減少的個體，故由兩部分組成。而R態(tài)t時刻的變化率則只可能是由I態(tài)節(jié)點轉(zhuǎn)化而來，故第三個方程僅有μρ(t)一項。

分析SIR模型對應(yīng)的微分方程，可以看出它隱含了許多近似假設(shè)以及平均場的思想?！磌〉的引入是平均場思想的重要體現(xiàn)，也是平均場方法最大的假設(shè)所在：亦即認為波動可以忽略，每一個節(jié)點的度都可以用平均度〈k〉來代替。上面提到的線性近似在處理閾值時適用的原因，是由于在低于傳播閾值λc的情況下，當(dāng)網(wǎng)絡(luò)中節(jié)點感染密度遠遠小于1時，對S(t)(1-(1-λ)〈k〉ρ(t))進行泰勒展開，得到S(t)λ〈k〉ρ(t)是較為精確的；然而當(dāng)感染密度較大(大略為ρ(t)>1/〈k〉)時，上述的近似結(jié)果從理論上講將不再十分精確(偏差將增大,而且通常解析值會高于實際值)，但已有的工作大多都表明它仍能較好地描述其趨勢、定性解釋許多現(xiàn)象，甚至定量地來看也可以得到不錯的結(jié)果[12,13,28,38-41]。

1.4平均場理論應(yīng)用的約束條件及準(zhǔn)確性分析

平均場方法從平均的、概率的角度去考量問題，免去了對于細節(jié)的計算，故而易于使用；但原則上來講，其適用條件是非常嚴(yán)格的，歷史上著名的科學(xué)案例(如二維Ising模型的相變問題)明確警示了濫用平均場方法的不良后果。數(shù)學(xué)上，平均場方法的適用范圍只能是完全圖，或者說系統(tǒng)結(jié)構(gòu)是充分混合的：在這種情況下，系統(tǒng)中的任何一個個體能夠以等概率接觸其他個體。只有滿足這一條件，平均場方法得到的解析結(jié)果才是完全精確的。

在復(fù)雜網(wǎng)絡(luò)研究中，與其說平均場理論是一種技術(shù)手段，把它看作一種思想方法或許更為恰當(dāng)，即從平均的、概率的角度去考量問題，進行宏觀的把握。其技術(shù)手段無外乎常微分方程、泰勒級數(shù)等等。從平均場方法自身發(fā)展的歷史來看，它也是不斷地縮小近似的范圍、進行更為合理的近似，著名的異質(zhì)平均場方法(Heterogeneous Mean-Field Method)[15,23]就是將度相同的節(jié)點進行平均近似，而不再把所有節(jié)點看作相似的。平均場的這種思想方法也滲透到了點對近似(Pairwise Appropriation)、主方程(Master equation)[28]等更為精確的方法中，這些方法近似范圍更小；當(dāng)然這些方法超越平均場近似之處在于，考慮了點對情況以捕捉動力學(xué)相關(guān)性[28]、關(guān)注度的演化[42]等等，使得相應(yīng)微分方程更為精確。

在復(fù)雜網(wǎng)絡(luò)研究中，標(biāo)準(zhǔn)的平均場方法一般隱含了3個前提假設(shè)[28](我們會在1.5中結(jié)合具體的平均場方程來依次具體闡述這些假設(shè)是如何被體現(xiàn)的)：

1)無局部聚類：對于任意度為k的節(jié)點A，其鄰居的狀態(tài)互相獨立：亦即A某一鄰居B1的狀態(tài)只受A的影響，而不受A其它鄰居(B2-Bk)的影響。故而當(dāng)鄰居之間存在連邊時，這一假設(shè)并不完全滿足；然而應(yīng)當(dāng)注意的是，實證發(fā)現(xiàn)在感染率小于或等于閾值時，流行病的傳播過程基本上是樹形分支過程，產(chǎn)生的回路極少[27]，而這也能部分地說明為何平均場方法在預(yù)測閾值上較為準(zhǔn)確。

2)無模體：所有度相同的節(jié)點，其動力學(xué)行為可用其平均情況代替。網(wǎng)絡(luò)中如果存在大量模體，即便兩個節(jié)點度相同，但由于可能處于不同模體的不同位置中，它們對于傳播所產(chǎn)生的影響將非常不同。例如三元模體中，度相同的節(jié)點在不同的模體中可能處于不同的位置、發(fā)揮不同的作用：有些可能是接收者，而有些則是作用者[43]；有些是純?nèi)坏膫鞑フ撸行﹦t是完全的被感染者。這類區(qū)分在有向網(wǎng)絡(luò)中更為明顯，例如圖1a中被雙環(huán)選中的度為2的兩個節(jié)點，在此種情形下，雖然二者度相同但其動力學(xué)行為必將大不相同。姑且不論有向網(wǎng)絡(luò)，即便是在無向網(wǎng)絡(luò)中，圖1b中4個被選中的節(jié)點雖然度均為3，然后它們在傳播過程中，第二個時間步所導(dǎo)致的影響將極不相同。在許多文獻中通常會看到平均場方法的一個基本假設(shè)是波動可忽略[23,41,44]，它實則是指所有的情形與平均情況相近，其波動不大：具體說來即網(wǎng)絡(luò)中所有節(jié)點度及動力學(xué)行為相近；或者起碼在度相同的情形下，其動力學(xué)行為相近。

3)無動力學(xué)相關(guān)性：更新某一節(jié)點A的狀態(tài)時，A的狀態(tài)與其鄰居(B1-Bk)的狀態(tài)相互獨立，即A狀態(tài)的更新變化不會作用于Bi的狀態(tài)。這一假設(shè)實則是忽略了鄰居節(jié)點間的動力學(xué)相關(guān)性，有理論研究發(fā)現(xiàn)，假設(shè) 3)正是平均場理論誤差產(chǎn)生的主要來源[28](點對近似恰恰是捕捉到了這種動力學(xué)相關(guān)性；對于這一假設(shè)的包含與否，恰是平均場方法與點對近似最大的區(qū)別之所在)。

但是當(dāng)聚類及模體作用較強時，點對近似同樣無法精確描述相應(yīng)動力學(xué)行為[28]，退火網(wǎng)絡(luò)鄰接矩陣中存儲的為其連接概率)由于每個時間步都會以相應(yīng)概率進行邊的重連，故而它取消了動力學(xué)相關(guān)性[12,45](但是應(yīng)當(dāng)注意的是，退火網(wǎng)絡(luò)仍然有拓撲相關(guān)性，對于任意節(jié)點其連接情況還是明確的，只是從0-1二元值變?yōu)榱烁怕手?，所以同配網(wǎng)絡(luò)和異配網(wǎng)絡(luò)的退火鄰接矩陣的整體數(shù)值分布情況還是很不同的)；另外，對于動態(tài)網(wǎng)絡(luò)，當(dāng)節(jié)點移動的速度較大時，也能夠較好地滿足平均場的假設(shè)[46]。故而一般認為網(wǎng)絡(luò)中節(jié)點平均度較大時解析結(jié)果相對準(zhǔn)確：因為當(dāng)平均度較大時，網(wǎng)絡(luò)更靠近于全連通圖，所以這種情況下動力學(xué)相關(guān)性反而較弱，不再起到?jīng)Q定性作用；當(dāng)平均度較小時，聚類、模體、動力學(xué)相關(guān)性的作用往往更為明顯而無法被忽略。

1.5SIR模型閾值解析

1.5.1無網(wǎng)絡(luò)結(jié)構(gòu)時的模型求解

首先對于方程(1)～(4)進行初步分析。由于流行病的有效傳播率等于λ/μ，這一分式可化為x/1的形式，μ值的變化只會影響感染演化的時間尺度[1]，故而假設(shè)μ=1可不失一般性。λ/μ也是單點接觸中的基本再生數(shù)R0，基本再生數(shù)R0表示一個感染個體在充分混合的S態(tài)群體中一次接觸平均能感染的個體數(shù)；當(dāng)有網(wǎng)絡(luò)結(jié)構(gòu)時，SIS過程的平均基本再生數(shù)R0=〈k〉λ/μ。通過這一概念可以初步定義流行病爆發(fā)閾值：只有當(dāng)R0>1時流行病才可能傳播開來，R0<1時流行病會在有限時間步內(nèi)消亡。

(5)

要考查流行病最終是否在全局爆發(fā)(或者說流行病能否長時間的存在于網(wǎng)絡(luò)中)，可直接分析傳播過程結(jié)束達到穩(wěn)態(tài)時的情況

R∞=1-S∞=1-e-λ〈k〉R∞

(6)

λc=1/〈k〉

(7)

但是稍后會發(fā)現(xiàn)它并不是SIR模型閾值的精確解，而是其上界。

除閾值這一問題外，我們往往還很關(guān)心最終的感染密度。由函數(shù)導(dǎo)數(shù)性質(zhì)可知，f′(0)的值越大，其函數(shù)曲線在此處斜率越大，最終與x軸交點的值也可能更大，亦即R∞值較大，且此時隨λ〈k〉的增大f(1)=-e-λ〈k〉→0。故而從平均場角度來看，λ〈k〉數(shù)值越大，則最終感染的比例R∞的數(shù)值也越大。

在λ=λc時，R∞的值較小趨近于0，|-λ〈k〉R∞|<1，故進行泰勒展開，可以得到

R∞=1-e-λ〈k〉R∞

(8)

(9)

又因為λc=1/〈k〉，可得

(10)

然而當(dāng)|-λ〈k〉R∞|>1時，式中的泰勒展開將不再適用；雖然方程R∞=1-e-λ〈k〉R∞形式非常簡單，然而不幸的是對于R∞并沒有一個簡單的閉合形式解即可由有限個初等函數(shù)與算術(shù)運算、指數(shù)、對數(shù)、根式表示，一定不包含無窮級數(shù)、連分?jǐn)?shù)、積分、微分或極限)。但是仍然可以得到相對精確的數(shù)值解：將方程看作y=R∞與y=1-e-λ〈k〉R∞聯(lián)立，在[0,1]區(qū)間內(nèi)求解即可，這時通過程序模擬便可以較為容易地獲得較為精確的結(jié)果，而且這一方法具有較好的普適性，應(yīng)用范圍也更廣。

1.5.2復(fù)雜網(wǎng)絡(luò)上的求解

方程組(1)～(4)從平均的角度出發(fā)，認為每個節(jié)點度均為〈k〉；然而當(dāng)網(wǎng)絡(luò)度分布取值較廣時，這樣的假設(shè)將無法刻畫真實情況。在無標(biāo)度網(wǎng)絡(luò)中，節(jié)點度分布呈現(xiàn)冪律，且往往存在胖尾效應(yīng)，〈k〉的近似與實際情況大相徑庭，無法用平均度較好地描述節(jié)點的度值，而在ER(Erd?s-Rényi)、WS(Watts-Strogatz)網(wǎng)絡(luò)中，度分布往往呈現(xiàn)指數(shù)有界特性，概率密度在均值〈k〉兩側(cè)迅速減小。

為了更準(zhǔn)確地刻畫異質(zhì)網(wǎng)絡(luò)上的流行病傳播過程，Pastor-Satorras和Vespignani提出了異質(zhì)平均場方法，假定度為k的所有節(jié)點所處的作用環(huán)境相近、遵從相同的動力學(xué)規(guī)律(這正是在1.4節(jié)中提到的假設(shè)1)[14-15]，這樣就不能在宏觀的層面而須要在介觀層面單獨考慮度不同的節(jié)點相應(yīng)的行為，此時引入變量Xk(t)表示t時刻度為k的節(jié)點中處于X(可為S、I、R)狀態(tài)節(jié)點的密度。相應(yīng)的SIR模型微分方程如下[23]：

(11)

(12)

(13)

1=Sk(t)+ρk(t)+Rk(t)

(14)

對于方程組(11)～(14)含義的解釋與方程組(11)～(14)相同，不同之處在于對網(wǎng)絡(luò)中節(jié)點按度分類后，考慮不同的度對于傳播的影響時，ρ(t)不再能有效描述任意度為k的給定節(jié)點的一條邊連接到一個感染態(tài)節(jié)點的概率；所以引入一個新的變量Θ(t)來表征這一概率。根據(jù)相應(yīng)物理意義可得其表達式：

(15)

當(dāng)網(wǎng)絡(luò)無度關(guān)聯(lián)性時，即度為k的節(jié)點并不以特定的概率而是會隨機地連向度為k′的節(jié)點。

(16)

對于同配與異配網(wǎng)絡(luò)條件概率P(k′|k)表達式將有很大不同。文獻[24]詳細討論了網(wǎng)絡(luò)存在度關(guān)聯(lián)時的情況。

同時，式(16)還包含了平均場近似的假設(shè)3：亦即認為節(jié)點A的k個鄰居(B1-Bk)之間狀態(tài)互相獨立，A與其任何兩個鄰居Bi、Bj不會構(gòu)成三角形，故而B1-Bk感染與否只與A有關(guān)，是相互獨立事件，所以式(16)中才能夠?qū)λ笑裬′進行加和。

由式(11)～(16)可以解得未感染節(jié)點的密度函數(shù)：

(17)

其中輔助函數(shù)

(18)

Φ(t)的物理意義是到t時刻、網(wǎng)絡(luò)中任意一條邊的另一端節(jié)點處于R態(tài)的概率。當(dāng)網(wǎng)絡(luò)上的動力學(xué)在到達穩(wěn)態(tài)時，流行病究竟有沒有在全局爆發(fā)，就可以通過Φ(t)得到清晰的認識。為了得到感染節(jié)點密度的表達式，我們發(fā)現(xiàn)觀察Φ(t)的均值隨時間的變化率會更為方便一些，結(jié)合式(13)可以得到：

(19)

(20)

(21)

1.5.3同質(zhì)網(wǎng)絡(luò)的求解

當(dāng)考慮具體的網(wǎng)絡(luò)結(jié)構(gòu)之后，可以得到更為具體的解析情況：

先來看隨機網(wǎng)，當(dāng)其規(guī)模N足夠大(N?k)時，將其度分布近似為泊松分布是較為嚴(yán)格、精確的：由于ER隨機網(wǎng)絡(luò)中任意兩邊都會以p的概率進行重連，所以某一節(jié)點其度為k的概率

(22)

通過這樣一個簡單的推導(dǎo)，在熱力學(xué)極限情況下，ER隨機網(wǎng)絡(luò)的度分布就是服從均值為〈k〉的泊松分布。

將P(k)帶入到等式(21)中，進而可知

(23)

網(wǎng)絡(luò)最終感染密度

R∞=1-S∞=1-e-λ〈k〉Φ∞≈λ〈k〉Φ∞

(24)

方程中省略了二次及以上的高階項。為了得到較高的精確度，保留最相關(guān)的項，在dΦ(t)/dt中將泰勒級數(shù)展開到二階：

(25)

(26)

其中A為形如ec的常數(shù)。將(10)代入(8)即得

(27)

通過計算機模擬也驗證了這一冪律關(guān)系(見圖2)。在臨界點附近λc存在冪律行為，表明這一過程是連續(xù)相變過程。

圖1　模體示意圖[43]

圖2　閾值點附近的感染密度變化情況[23]

1.5.4異質(zhì)網(wǎng)絡(luò)的求解

(28)

其中kmax～N1/(γ-1)[51]，在熱力學(xué)極限下(亦即當(dāng)N→∞)，kmax→∞，分母發(fā)散，進而有λc→0。在這里也可以明顯看出，在熱力學(xué)極限情況下，其閾值會趨近于0。雖然現(xiàn)在有文章指出在SIS模型下SF網(wǎng)絡(luò)的閾值為0，對于γ>2.5的情形閾值消失的原因主要是由于度最大的中心節(jié)點在任何傳染概率之下均會處于活躍狀態(tài)，而當(dāng)γ<2.5時則主要是由于進行k-core分解后，k-core值最高的團簇處于活躍態(tài)所致[25,52]；但原有的解釋在SIR模型中的解釋仍是準(zhǔn)確的[25]。

相應(yīng)地

(29)

Φ∞在穩(wěn)態(tài)時被看作一個常數(shù)，所以可將其提到積分號外部，而在下面Φ(t)的求解中，則無法采用如此假設(shè)，因為它是與k直接相關(guān)的函數(shù)：

(30)

(31)

其中,γE是相應(yīng)的Γ函數(shù)。對其進行積分后可得：

(32)

其中A為一積分常數(shù)。當(dāng)r→∞時，

(33)

因而結(jié)合式(29)和(33)，可以得到

(34)

對于非零的λ值，最終的R態(tài)節(jié)點密度都不為0；而這一分析也與λc的解析相一致(見圖3)。

圖3　感染密度與1/λ的函數(shù)關(guān)系[23]

1.6平均場方法的應(yīng)用

可以發(fā)現(xiàn)，應(yīng)用平均場方法的求解過程充滿了這樣或那樣的近似與假設(shè)，但這也大大簡化了解析的難度，讓我們得以抓住最主要的影響因素，得到相對準(zhǔn)確的結(jié)論。在上述的解析過程中，大多都是進行線性假設(shè)，舍掉了相應(yīng)的高階項；在對于閾值的解析中，由于感染率一般很小，故而舍棄高階項產(chǎn)生的誤差會相對較小；然而在最終感染密度的解析中，這樣的線性近似會使其解析結(jié)果偏大：改進過的非微擾平均場(NonPerturbative Heterogeneous Mean-field method)保留了這些高階項，使得解析的準(zhǔn)確度大大提高[38]，但由于它還是假設(shè)度相同的節(jié)點動力學(xué)行為相同，進行了相應(yīng)的平均代換，所以其數(shù)值結(jié)果仍與真實模擬值略有差異，而且通常比實際值略高。

普遍認為平均場方法主要有本文介紹的3個假設(shè)，但在文中也可以看到標(biāo)準(zhǔn)的異質(zhì)平均場方法也沒有明確地去考慮度關(guān)聯(lián)性和節(jié)點的接觸能力(即感染節(jié)點會接觸自己鄰居的多大比例，是全部抑或只是一小部分)；然而這兩個因素并不是其應(yīng)用的前提，因為這兩種因素都可以被平均場方法捕捉到，改進后較為完整的微分方程組如下[53]：

(35)

(36)

(37)

1=Sk(t)+ρk(t)+Rk(t)

(38)

方程組(35)～(38)中表征節(jié)點的接觸能力，亦即某一節(jié)點只能接觸到目標(biāo)節(jié)點比例的邊；當(dāng)A=1時即為單點接觸(Contact Process：每次只選擇某一鄰居進行感染，若鄰居已然處于感染態(tài)，則跳過感染過程)。這樣一來，就將度關(guān)聯(lián)性與接觸能力這兩個因素納入了考慮范疇，屆時只須根據(jù)具體情況給出P(k′|k)以及相應(yīng)的A值即可。對于有權(quán)網(wǎng)絡(luò)，其上單點接觸往往存在依賴權(quán)重進行偏好選擇，對于這一情況，只須將接觸概率與感染概率分開考慮并將其相乘即可[38]。

另外對于SIS模型，由于I態(tài)節(jié)點不會變?yōu)镽態(tài)，而會以一定概率μ恢復(fù)為S態(tài)，其微分方程如下[15]：

(39)

(40)

1=Sk(t)+ρk(t)

(41)

此外，當(dāng)SIR傳播模型中感染節(jié)點不以某一概率μ而是經(jīng)過固定時間步τ后恢復(fù)時，其傳播過程是可以嚴(yán)格對應(yīng)于滲流過程的：流行病的傳播過程相當(dāng)于以概率p(p=1-(1-λ)τ)從源節(jié)點不斷進行加邊[27,54]；對于時間步為變量(或者以概率恢復(fù)，二者效果基本相同)的情況，由于引入了更強的隨機性和時序關(guān)聯(lián)性，而無法被滲流方法所準(zhǔn)確描述[55]。滲流方法會在后面的章節(jié)詳細闡述。

1.7平均場方法小節(jié)

縱觀平均場方法，它更多的是一種思想方法而非技術(shù)手段：亦即從平均的、概率的角度去考慮問題，進行宏觀的把握，把環(huán)境對物體的作用進行集體處理，以平均作用效果替代局域作用效果的加和，其具體的數(shù)學(xué)技術(shù)手段無外乎率方程、常微分方程、泰勒級數(shù)等等。

其特點正是只考慮最為主要的影響因素，考查其速率變化，略去其它相比之下影響較小的因素：如不同度之間的轉(zhuǎn)換[42]、點對選取后的情況[29]、子節(jié)點對父節(jié)點的重感染概率[27]等相應(yīng)因素。平均場方法正是在一系列假設(shè)與近似基礎(chǔ)上，犧牲了一定的準(zhǔn)確性來換取其極大的簡潔性。然而是否存在整合上述多種影響因素的精確方法仍然是一個極大的挑戰(zhàn)；是否一定是考慮的因素越多其結(jié)果就越精確，同樣是一個重要的問題；探究不同因素之間的互相影響以及不同因素對于傳播動力學(xué)作用的影響程度，將是非常有意義與挑戰(zhàn)性的工作。

平均場理論被廣泛地應(yīng)用到流行病傳播閾值及最終感染密度的解析中[15,23]，它出現(xiàn)之后也成為其它理論方法爭相參照的對象,將平均場方法看作復(fù)雜網(wǎng)絡(luò)動力學(xué)解析中最基礎(chǔ)的方法應(yīng)該也會被大多數(shù)人所接受。通過對于平均場方法更深入的研究和應(yīng)用，也得出了許多重要結(jié)論，諸如：對于SIS模型，無標(biāo)度網(wǎng)絡(luò)在二階矩發(fā)散的情況下，其閾值均趨近于0[24]。而對于SIR模型，則無論標(biāo)度因子的大小，都是k-core值最高的團簇在影響著流行病的爆發(fā)閾值和傳播過程[52]。另外高集群系數(shù)和負度關(guān)聯(lián)會保護相應(yīng)的無標(biāo)度網(wǎng)絡(luò)(抑制流行病的傳播)[26]。

平均場方法在穩(wěn)態(tài)時感染密度的解析中也占有主導(dǎo)地位[15]；同時在謠言傳播解析中也被廣泛應(yīng)用[44,56]，文獻[30]中在信息傳播中引入了信息之間的互相競爭，其數(shù)學(xué)解析也是通過平均場方法完成的。在真實流行病的預(yù)測應(yīng)用中，許田[57]應(yīng)用平均場方法下的SIR模型有效地模擬了北京2003年4月18日到6月16日每天新增SARS病例統(tǒng)計數(shù)據(jù)的研究結(jié)果，作者得到的結(jié)果清晰地顯示出數(shù)據(jù)模擬結(jié)果相當(dāng)于實際數(shù)據(jù)的一種“平滑化”，意即其結(jié)果相當(dāng)于全局的、平均的描述。而且在新興的多層網(wǎng)絡(luò)研究中，平均場方法也為進一步理解其上傳播閾值減小這一動力學(xué)行為提供了有力解釋[58]。

2點對近似

2.1點對近似方法概述

在上一節(jié)中，我們已經(jīng)指出平均場方法誤差產(chǎn)生的主要原因即是忽略了節(jié)點間的動力學(xué)相關(guān)性；為了更準(zhǔn)確地描述流行病傳播動力學(xué)過程，許多研究者不再假設(shè)節(jié)點間動力學(xué)無關(guān)(即節(jié)點狀態(tài)與鄰居節(jié)點狀態(tài)無關(guān)，只考慮處于不同態(tài)節(jié)點數(shù)目的變化)，開始考慮節(jié)點狀態(tài)之間的關(guān)聯(lián)性[59-61]，也就是將處于不同狀態(tài)節(jié)點的連邊數(shù)量或比例作為一個動態(tài)參量，考慮不同類型連邊數(shù)目的變化。想要得出不同類型節(jié)點對數(shù)目變化的微分方程，需要依賴于網(wǎng)絡(luò)中不同類型的三元組數(shù)目(三元模體只是三元組的一種特例：三元組只要保證AB、BC分別相連即可，無論AC相連與否均可；三元模體則必須保證三者均兩兩相連)。為使相應(yīng)的微分方程閉合，需要用節(jié)點對的數(shù)目來近似網(wǎng)絡(luò)中三元組的數(shù)目，這種方法被稱做點對近似(Pairwise Approximation)[7,61]。之所以用節(jié)點對數(shù)目來進行相應(yīng)的近似，原因在于在非線性模型中，描述第n階矩的微分方程會依賴于第n+1階矩，這樣得出的一系列方程是不可解的；為了使方程可解，必須使方程在某一階矩上閉合。假設(shè)要使方程在第n階閉合，必須用n階及n階以下的矩來估計第n+1階矩，這個過程被稱為矩閉合近似。

點對近似方法同樣是用一組常微分方程來描述簡單空間模型的動力學(xué)行為，例如簡單的宿主—寄生蟲模型[60,72-73]、疾病傳播模型[7,42,61,64]、投票模型[65]、空間博弈[74-76]等等。它不僅描述網(wǎng)絡(luò)中節(jié)點的狀態(tài)變化，而且還同時考慮“存在連接的節(jié)點對”(下文中為簡便起見,將簡稱“節(jié)點對”)中兩個節(jié)點狀態(tài)之間的相關(guān)性。這樣就解決了困擾平均場方法的動力學(xué)相關(guān)性，其假設(shè)條件相對平均場方法更弱，進而可以更為準(zhǔn)確地描述網(wǎng)絡(luò)中的動力學(xué)過程[42,64-65]。

點對近似方法不但提高了解析的精度，而且它還具有更廣的適用范圍和拓展性：它既適用于靜態(tài)網(wǎng)絡(luò)，包括同質(zhì)網(wǎng)絡(luò)[42,62]和異質(zhì)網(wǎng)絡(luò)[42,64-65]；又適用于動態(tài)網(wǎng)絡(luò)[66-68]。當(dāng)網(wǎng)絡(luò)結(jié)構(gòu)并非固定不變時，例如節(jié)點或邊存在增加或刪除的情形時，節(jié)點與其鄰居的連邊會不斷發(fā)生變化，此時網(wǎng)絡(luò)結(jié)構(gòu)的變化對傳播過程的影響很大，所以必須考慮網(wǎng)絡(luò)中不同狀態(tài)的連邊數(shù)量隨時間的變化，因此點對近似方法很適于描述動態(tài)網(wǎng)絡(luò)上的動力學(xué)過程。

下面將分別討論在不同類型的網(wǎng)絡(luò)上，如何使用點對近似方法來解析不同傳播模型的動力學(xué)過程。

2.2靜態(tài)網(wǎng)絡(luò)上的求解

2.2.1點對近似方法預(yù)備知識

在真正給出描述動力學(xué)過程相應(yīng)的微分方程之前，需要一些有別于平均場方法的預(yù)備知識。

與平均場方法相似，考慮處于不同態(tài)節(jié)點數(shù)量的變化率時，仍然使用微分方程進行描述。定義抽象函數(shù)f用于描述網(wǎng)絡(luò)中節(jié)點狀態(tài)(例如用f表示處于某種狀態(tài)的節(jié)點數(shù)目(如[A])或某種類型的連邊數(shù)目(如[AB]))，當(dāng)網(wǎng)絡(luò)趨近于無窮大時，可將f近似看作是連續(xù)的，那么其變化率可用如下微分方程表示為

(42)

接下來，問題就集中在如何求解[ABC]和[ABA]上。

當(dāng)A≠C時，

(43)

其中,

(44)

當(dāng)A=C時，

(45)

其中,

(46)

接下來分別求解Γσ(A|B|A)和Γσ(A|B|C)[62]：

=EBj=1{[Qj(A|B)-EBj=1(Qj(A|B))]2}=varBj=1[Qj(A|B)]

(47)

其中EBj=1[]表示狀態(tài)為B的節(jié)點上相應(yīng)狀態(tài)節(jié)點數(shù)目的平均值，varBj=1[]表示方差。也就是說，Γσ(A|B|A)是網(wǎng)絡(luò)中所有B態(tài)節(jié)點所連接的A態(tài)鄰居節(jié)點的個數(shù)的方差。

下面分別來分析Qj(A|B)服從不同分布時適用的情況：

(48)

2)Qj(A|B)服從二項分布，這種情況主要適用于規(guī)則網(wǎng)絡(luò)(每個節(jié)點有m個固定的鄰居)：

(49)

與求解Γσ(A|B|A)相類似，

(50)

由此可知Γσ(A|B|C)是Qj(A|B)和Qj(C|B)的協(xié)方差。下面來考慮在B態(tài)節(jié)點其鄰居中C態(tài)A和態(tài)節(jié)點的分布情況：

[ABC]=[AB][BC]/[B]

(51)

(52)

(53)

到這里，就求解出了網(wǎng)絡(luò)中的三元組相應(yīng)的數(shù)目，然而上述這些計算仍然沒有涉及到網(wǎng)絡(luò)上的動力學(xué)。下面針對不同的傳播模型(SIS/SIR)來給出具體的解析。假設(shè)流行病的感染概率為β，恢復(fù)概率為γ。

2.2.1.1SIS模型上的求解

首先來看節(jié)點數(shù)目變化率的微分方程組,依然應(yīng)用類似于平均場方法的線性近似：假設(shè)一個S態(tài)節(jié)點有n個I態(tài)鄰居，則其被感染概率的準(zhǔn)確值是1-(1-β)n，進而在滿足約束條件(nβ?1)時近似為nβ,進而有：

(54)

(55)

[S]+[I]=N

(56)

接下來再來看點對近似方法所特有的關(guān)于不同連邊類型的變化率方程：

(57)

(58)

(59)

[SS]+2[SI]+[II]=M

(60)

其實從上述的方程組中，仍然可以發(fā)現(xiàn)這里面仍然有許多必要的假設(shè)以及近似：首先仍然是把每個感染或者恢復(fù)看作獨立事件，這也正是每個式子中都有求和項的原因；而且由于這一假設(shè)，真實情況下可能存在的II→SS這一類的轉(zhuǎn)化就沒有被考慮其中。

[S]+[I]=N

(56)

[SS]+2[SI]+[II]=M

(60)

(61)

(62)

(63)

接下來，如若求解流行病爆發(fā)閾值，只需設(shè)定好相應(yīng)的初始值，然后不斷地迭代方程，當(dāng)方程的最終解由0變?yōu)榉?時，就可以數(shù)值地確定出相應(yīng)的閾值。

對于SI模型，只需要將SIS模型中的恢復(fù)概率置為0即可，此處不作贅述。

2.2.1.2SIR模型上的求解

SIR模型解析思路與SIS模型情形相似，只是多了一個R態(tài)，下面只給出相應(yīng)的微分方程組，不再展開后續(xù)的求解：

(64)

(65)

(66)

(67)

(68)

(69)

(70)

(71)

(72)

2.2.2異質(zhì)網(wǎng)絡(luò)上的求解

這里靜態(tài)的異質(zhì)網(wǎng)絡(luò)[42,64-65]是指網(wǎng)絡(luò)度分布較廣、大多數(shù)節(jié)點的度值較明顯偏離于網(wǎng)絡(luò)平均度的網(wǎng)絡(luò)，例如無標(biāo)度網(wǎng)絡(luò)；在此種情況下，用平均度來代表所有節(jié)點的度往往無法得出正確的估計，仍須在介觀的層次上來對不同的度分別進行考慮，寫出各自相應(yīng)的微分方程，在這一解析思路上它與異質(zhì)平均場方法是相似的。

2.2.2.1SIS模型上的求解

對于感染個體的動力學(xué)，包括感染個體以概率γ恢復(fù)為易感個體和易感個體被感染態(tài)鄰居感染；這就須要把鄰居節(jié)點的狀態(tài)信息包含到節(jié)點的狀態(tài)方程中：

(73)

對于易感個體的動力學(xué)變化，相類似地可構(gòu)建出方程(76)。

(74)

為了能夠求解方程(73)～(74)，接下來需要構(gòu)建相應(yīng)的方程來模擬節(jié)點對的動力學(xué)：

(75)

式(75)中的各項表示[SnIm]的增加量(由[SnSm]中Sm的被感染、[InIm]中In的恢復(fù)所導(dǎo)致)和減少量(由[SnIm]中的Sn被Im或其他與此Sn相連的I態(tài)節(jié)點感染、[SnIm]中Im的恢復(fù)所導(dǎo)致)。

類似地可得其他類型點對的動力學(xué)方程：

(76)

(77)

理論上來說，方程中的三元組可以用更完整的四元組來表示，更高階的處理可依此類推。但是，這樣會導(dǎo)致方程變得非常復(fù)雜甚至不可解，所以為了簡化，類似于均勻網(wǎng)絡(luò)的分析，進行點對近似，也就是用不同類型邊(二元組)的數(shù)量來近似估計三元組的數(shù)目[61,69]：

(78)

進而依照同質(zhì)網(wǎng)絡(luò)中相應(yīng)的求法，即可解出網(wǎng)絡(luò)上流行病爆發(fā)的閾值。

2.2.2.2SIR模型求解

其相應(yīng)的微分方程為

(79)

(80)

(81)

(82)

(83)

(84)

(85)

(86)

(87)

其求解過程與SIS模型相似，此處不再贅述。

2.3動態(tài)同質(zhì)網(wǎng)絡(luò)上的求解

相比于靜態(tài)網(wǎng)絡(luò)，動態(tài)網(wǎng)絡(luò)最大的特點在于網(wǎng)絡(luò)的拓撲結(jié)構(gòu)不再一成不變，而是隨著時間的演化、周圍環(huán)境的變化，可能會發(fā)生相應(yīng)的改變。Gross等[67]使用點對近似方法來模擬動態(tài)網(wǎng)絡(luò)上的傳播，其模型為：隨機圖的規(guī)模為N，平均度為〈k〉，每個節(jié)點處于S態(tài)或I態(tài)之一，在每個時間步對于每條SI類型的邊，S節(jié)點以概率β被感染為I態(tài)節(jié)點，感染節(jié)點以概率γ恢復(fù)為易感態(tài)，此外對每條SI邊，S節(jié)點以概率w(w∈[0,1]自適應(yīng)地斷開與I態(tài)節(jié)點的連邊，并將斷邊隨機重連到一個S態(tài)節(jié)點上。在這個模型中，零階矩是感染態(tài)和易感態(tài)節(jié)點的數(shù)目，即[I]和[S]；一階矩是平均到每個節(jié)點上的SS、SI和II邊數(shù)目，即[SS]、[SI]和[II]；二階矩是三元組的數(shù)目，即[ABC](A,B,C∈{I,S})，這里的邊和三元組都是無向的(仍可將其看作雙向邊，這里只會影響到[M]的數(shù)值以及相應(yīng)方程前面的系數(shù)“2”)。[I]及[SS]和[II]的動力學(xué)微分方程加上兩個守恒條件[I]+[S]=N和[SS]+[SI]+[II]=N〈k〉/2就可以完整地刻畫零階矩和一階矩的動力學(xué)。

由此得出的方程組(88)～(90)。

(88)

每個S態(tài)節(jié)點的鄰居中I態(tài)節(jié)點的平均數(shù)量為[SI]/[S]，所以每個S態(tài)節(jié)點被感染的概率為β[SI]/[S]，S態(tài)節(jié)點密度為[S]，因此由感染引起的I態(tài)節(jié)點增加的速率為β[SI]；因為每個I態(tài)節(jié)點都會以概率γ恢復(fù)為健康態(tài)，所以由于感染態(tài)節(jié)點恢復(fù)引起的I態(tài)節(jié)點減少的速率為γ[I]。

(89)

在一次感染事件中，流行病通過一條SI邊來傳播，會將SI邊變?yōu)镮I邊，因此一次成功感染事件會產(chǎn)生至少一條II邊。但是如果除了感染它的那個I態(tài)節(jié)點之外，這個新被感染節(jié)點的鄰居中還有其他I態(tài)節(jié)點，那么這次感染事件還會產(chǎn)生其他的II邊；因此，一次感染事件總共產(chǎn)生的II邊的數(shù)目是1+[ISI]/[SI]，其中“1”表示感染發(fā)生的那條SI邊，第二項表示由這條SI邊形成的ISI三元組的數(shù)目。由此可得，II邊產(chǎn)生的速率為β[SI](1+[ISI/[SI]=β([SI]+[ISI])。如果II邊上一點以概率β恢復(fù)，那么II邊減少；一個I態(tài)節(jié)點形成的II邊的平均數(shù)量為2[II]/[I]，其中的“2”是由于一條II邊連接了兩個I態(tài)節(jié)點、在統(tǒng)計邊數(shù)時須要將一條II邊重復(fù)計數(shù)一次,因此II邊減少的速率為2γ[II]。

(90)

SI邊中的I態(tài)節(jié)點恢復(fù)會使SS邊增加，相應(yīng)的速率為γ[SI]，同時，每個重連事件都會增加一條SS邊，所以由重連引起的SS邊增加的速率為w[SI]；因為每條SS邊中的S節(jié)點可能同其他I態(tài)節(jié)點相連，所以若這條SS邊中的S節(jié)點被感染就會造成SS邊減少，每條SS邊形成的SSI類型的三元組數(shù)目為[SSI]/[SS]，由此可知，SS邊減少的速率為β[SSI]。

此時，這個方程組仍不是一個閉合的模型，它們還依賴于兩個未知二階矩[SSI]和[ISI]。而方程中的一階矩已經(jīng)能夠很好地反映重連的效果，因此使用點對近似方法來使系統(tǒng)閉合：用零階和一階矩來近似地估計二階矩。

ISI三元組的一半是一條SI邊，其密度為[SI]。為了近似一條給定SI邊形成的ISI三元組數(shù)量，必須計算出連向此S節(jié)點的其余I態(tài)節(jié)點的平均數(shù)量。假設(shè)除去給定SI邊之外，S節(jié)點平均還有〈q〉條邊，其中〈q〉表示平均剩余度，即通過一條隨機選擇的邊到達的節(jié)點它除去此邊之外擁有的平均邊數(shù)。又因為這〈q〉條邊中每條邊是一條SI邊的概率為[SI]/〈k〉[S]。已知SI邊的密度、SI邊中S節(jié)點有其他SI邊的概率，所以[ISI]=κ[SI]2/[S]，其中κ=〈q〉/〈k〉。初始網(wǎng)絡(luò)是隨機網(wǎng)絡(luò)，剩余度分布q(k)=(k+1)P(k+1)/〈k〉，平均剩余度

(91)

點狀圖對應(yīng)數(shù)值模擬結(jié)果，實線為方程

又因為P(k)服從均值為λ的泊松分布(P(k)=e-λλk/k!)，因此相應(yīng)的均值E(k)和方差D(k)=E(k2)=[E(k)]2=〈k2〉-〈k〉2的關(guān)系為E(k)=D(k)=λ=〈k〉，由以上兩式可知，〈k2〉=〈k〉2=〈k〉，因此〈q〉=〈k〉；在此情況下κ=1，兼由式(48)可得[ISI]=[SI]2/[S]。類似的，由式(51)可得，[SSI]=2[SS][SI]/[S]，其中的“2”表示一條SS邊中的兩個S節(jié)點(在有向圖中將沒有這個系數(shù)“2”)。由此可得到一個閉合系統(tǒng)的微分方程組：

(88)

(92)

(93)

其中圓圈表示實際模擬的結(jié)果，方框表示用點對近似解出的結(jié)果，此處初始感染比例i(0)=0.3，感染概率p=0.008，重連概率w=0.3

2.4點對近似方法適用條件及準(zhǔn)確性概析

點對近似方法適于對沒有較短回路的網(wǎng)絡(luò)結(jié)構(gòu)進行解析，它對拓撲多為局域樹形結(jié)構(gòu)的網(wǎng)絡(luò)動力學(xué)過程近似非常準(zhǔn)確。當(dāng)網(wǎng)絡(luò)中沒有強連通模體時，若i和j直接相連，同時節(jié)點i與j存在另一條非直接相連的路徑，即i和j之間存在的回路中較長的那條路徑，這條路徑越長，由它導(dǎo)致的節(jié)點i和j的狀態(tài)之間的關(guān)聯(lián)性就越小，越可以忽略，因此點對近似方法的近似較準(zhǔn)確。而且點對近似方法擴展性極強，可以處理幾乎各種類型的網(wǎng)絡(luò)。但是對于動態(tài)網(wǎng)絡(luò)中，當(dāng)在演化過程中會出現(xiàn)較強的關(guān)聯(lián)性時，點對近似就不能比較準(zhǔn)確地刻畫其演化情況，如圖5對比了自適應(yīng)模型[67]在雙穩(wěn)態(tài)參數(shù)下感染i(t)比例隨時間演化的模擬結(jié)果和點對近似分析的結(jié)果，由圖可知點對近似此時解析不準(zhǔn)確，因為網(wǎng)絡(luò)的重連導(dǎo)致網(wǎng)絡(luò)中出現(xiàn)了很強的社區(qū)性，關(guān)聯(lián)性也隨之變強。

然而許多實際社會網(wǎng)絡(luò)中都存在很多回路，存在大量模體時[70-71]，點對近似方法同樣無法精確刻畫，其誤差仍然存在；對于模體這一因素，平均場方法與點對近似都無能為力。點對近似對于平均場方法最大的改進就在于其加入了對動力學(xué)相關(guān)性的考慮，使得分析結(jié)果更為準(zhǔn)確；同時，點對近似方法具有很好的拓展性，適用于許多模型(例如簡單的宿主-寄生蟲模型[60,72-73]、流行病傳播模型[7,42,61,68]、投票模型[65]、空間博弈[74,76]等等)。

3主方程方法

3.1主方程方法概述

到目前為止，應(yīng)用微分方程的解析方法中，最為精確的當(dāng)數(shù)主方程方法，這是馬爾可夫過程概率分布的演化方程，研究概率分布隨時間的演化。主方程是統(tǒng)計物理里最重要的方法之一，它幾乎是普遍適用的，廣泛地并應(yīng)用于化學(xué)、生物學(xué)、金融、人口動力學(xué)、復(fù)雜網(wǎng)絡(luò)動力學(xué)等諸多學(xué)科領(lǐng)域，以及布朗運動、流體、半導(dǎo)體等具體問題的解析當(dāng)中。

在前面介紹的平均場和點對近似方法都假設(shè)度相同的節(jié)點動力學(xué)行為相似，而點對近似卻更近一步考慮了節(jié)點對間的動力學(xué)相關(guān)性。為了更精確地描述傳播動力學(xué)過程，Gleeson更深入一步，不僅考慮度為k的節(jié)點的狀態(tài)變化，而且還考慮其鄰居的狀態(tài)變化，提出了相應(yīng)的主方程方法來對復(fù)雜網(wǎng)絡(luò)上SIS模型傳播動力學(xué)進行解析[42]；該方法基于Marceau等[77]提出的思路，并進行了相應(yīng)改進。實驗結(jié)果表明主方程方法比平均場和點對近似方法得到的解析結(jié)果更準(zhǔn)確；而且可以通過適當(dāng)?shù)慕疲嘶癁辄c對近似方法和平均場方法的演化方程；另外此方法還可用于對自旋模型等其他只包含兩個狀態(tài)轉(zhuǎn)換的動力學(xué)過程進行解析。

3.2對SIS模型的求解

給定網(wǎng)絡(luò)的度分布，利用配置模型[78]可以生成相應(yīng)的靜態(tài)無向網(wǎng)絡(luò)。在SIS模型中，每個節(jié)點要么處于易感狀態(tài)S，要么處于感染狀態(tài)I。

Gleeson提出的主方程方法假定節(jié)點的感染概率和康復(fù)概率依賴于節(jié)點的度k和感染鄰居數(shù)m，令Fk，mdt表示一個度為k且在t時刻具有m個感染鄰居的易感節(jié)點，在t+dt時刻變?yōu)楦腥緫B(tài)的概率，令Rk，mdt表示一個度為k且具有m個感染鄰居的感染節(jié)點，在dt時間間隔內(nèi)恢復(fù)為易感態(tài)的概率。對于SIS模型，有Fk，m=λm和Rk，m=μ[1]，其中λ和μ分別表示感染率和康復(fù)率。

對于每個m值(0≤m≤k)，sk,m(t)和ik,m(t)演化的主方程為

(94)

(95)

其中，每個方程右邊的前兩項分別表示度為k的易感節(jié)點被感染和感染節(jié)點康復(fù)為易感態(tài)；其余4項考慮其鄰居的狀態(tài)變化(S態(tài)鄰居節(jié)點被感染或者I態(tài)鄰居節(jié)點康復(fù))：其中βs(k-m)sk，m表示度為k的易感節(jié)點其k-m個易感鄰居節(jié)點變?yōu)楦腥緫B(tài)的個數(shù)，βs表示在dt時間間隔內(nèi)易感節(jié)點的易感鄰居變?yōu)楦腥緫B(tài)的比例；相應(yīng)地，γs表示dt時間間隔內(nèi)易感節(jié)點的感染態(tài)鄰居康復(fù)為易感態(tài)的比例；βi表示dt時間間隔內(nèi)感染節(jié)點的易感鄰居變?yōu)楦腥緫B(tài)的比例；γi表示dt時間間隔內(nèi)感染節(jié)點的感染態(tài)鄰居康復(fù)為易感態(tài)的比例。其余的4項旨在考慮度相同但鄰居狀況不同的節(jié)點，它們產(chǎn)生哪些變化后會對sk，m的數(shù)目產(chǎn)生影響。在解析相應(yīng)鄰居節(jié)點變化率的過程中，研究者實際上應(yīng)用了點對近似和平均場的思路：式中βs表示在dt時間間隔內(nèi)整個網(wǎng)絡(luò)中所有sk，m中度為k的易感節(jié)點被感染后，所導(dǎo)致的原有S-S邊變?yōu)镾-I邊占原有S-S邊的比例，亦即對于網(wǎng)絡(luò)中任意一條S-S邊其在dt時間間隔后變?yōu)镾-I邊的概率；相應(yīng)地，γs表示dt時間間隔內(nèi)整個網(wǎng)絡(luò)中所有ik，m中度為k的感染態(tài)節(jié)點恢復(fù)后，所導(dǎo)致的原有S-I邊變?yōu)镾-S邊占原有S-I邊的比例；βi表示dt時間間隔內(nèi)整個網(wǎng)絡(luò)中所有sk,m中度為k的易感節(jié)點被感染后，所導(dǎo)致的原有I-S邊變?yōu)镮-I邊占原有I-S邊的比例；γi表示dt時間間隔內(nèi)整個網(wǎng)絡(luò)中所有ik,m中度為k的感染態(tài)節(jié)點恢復(fù)后，所導(dǎo)致的原有I-I邊變?yōu)镮-S邊占原有I-I邊的比例。

我們發(fā)現(xiàn)，βs雖然是關(guān)注sk,m中那個度為k的s態(tài)節(jié)點被感染，但是由于是對于全網(wǎng)中所有S態(tài)節(jié)點進行了考慮，所以可以用這樣一個概率來描述任意一條S-S邊中某一端S節(jié)點被感染的概率。

而且通過這樣一個式子，也可以發(fā)現(xiàn)為什么系統(tǒng)(1)～(2)中在考慮sk,m的變化率時只考慮了sk,m-1,sk,m+1，而沒有考慮sk,m-2,sk,m+2及之后的項，因為后面的項它們本身的轉(zhuǎn)換概率階數(shù)更高，將其忽略所產(chǎn)生的影響將非常小。該方法對節(jié)點的分類較點對方法更為細致，而且其所進行的近似也非常合理，這也是它高精確度的原因所在。

上述這4個參數(shù)的值可以通過跟蹤網(wǎng)絡(luò)中每條邊的兩個節(jié)點的狀態(tài)變化計算得到：例如對于βs，首先計算t時刻兩個節(jié)點均是易感節(jié)點的邊數(shù)(即S-S邊數(shù))，其次計算在dt時間內(nèi)由S-S邊變?yōu)镾-I邊的數(shù)目，并且對于所有不同的度都進行計算后求取平均，后者與前者的比值即是βs；其它3個參數(shù)γs、βi、γi通過相同的計算方法得到。βs、γs、βi、γi這4個參數(shù)的計算公式為

(96)

(97)

(98)

(99)

給出一定的初始條件，數(shù)值迭代上述方程即可得出其感染密度，當(dāng)感染密度從0變?yōu)榉?時即可求得其爆發(fā)閾值。

3.3主方程方法近似條件

主方程方法與點對近似、平均場理論都使用微分方程描述相應(yīng)的動力學(xué)變化，而且它也應(yīng)用了前兩者的一些解析思想，通過適當(dāng)?shù)慕?，就可以得到?biāo)準(zhǔn)點對近似以及平均場的方程。

我們?nèi)匀豢梢酝ㄟ^考慮點對的情況，亦即pk(t)(對于任意度為的易感節(jié)點,隨機選擇的一個鄰居在t時刻處于感染態(tài)的概率，亦即某一點度為k的S-I邊的比例)和qk(t)(對于度為k的感染節(jié)點,隨機選擇的一個鄰居在t時刻處于感染態(tài)的概率，亦即某一點度為k的I-I邊的比例)，來對系統(tǒng)(94)，(95)進行近似推導(dǎo)出點對近似方法。

(100)

(101)

(102)

更進一步，可以從平均場的角度，利用先前1.5小節(jié)中所介紹過的概率Φ(在SIS模型中的含義為：在任意時刻,隨機選擇一條邊、其另一端節(jié)點處于感染態(tài)的概率)對系統(tǒng)(94)，(95)進行近似，推導(dǎo)出平均場近似方法；當(dāng)網(wǎng)絡(luò)不存在度關(guān)聯(lián)性時，Φ=P(k)kρk/k=〈ρkk/〈k〉〉。利用概率Φ代替系統(tǒng)(100)～(102)中的pk(t)和qk(t)，亦即將sk，m≈(1-ρk)Bk，m(Φ)和ik，m≈ρkBk，m(Φ)代入系統(tǒng)(94)，(95)即可推導(dǎo)出平均場近似方法[79]：

(103)

系統(tǒng)(103)的初始條件為ρk(0)=ρ(0)。當(dāng)k的變化范圍是[0,kmax]時，系統(tǒng)(103)至多包括kmax+1個方程。進行進一步展開，即可得到

(104)

這與1.5節(jié)中式(12)是一致的(SIS與SIR模型中ρk變化率的方程是相同的)。

3.4主方程方法小節(jié)

Gleeson提出的對SIS傳播模型進行解析的主方程方法，是目前基于微分方程組解法中最為精確的方法，它不僅考慮了節(jié)點狀態(tài)的變化，而且還考慮節(jié)點了的鄰居節(jié)點狀態(tài)的變化，將動力學(xué)演化過程考慮得更為細致。主方程方法比點對近似方法和平均場方法的解析結(jié)果更為準(zhǔn)確(如圖6所示)，它解析出的感染密度隨時間的演化情況和模擬結(jié)果非常接近。而且通過適當(dāng)?shù)慕疲捎芍鞣匠掏茖?dǎo)出點對近似方法和平均場方法。另外該方法不僅可以解析SIS傳播動力學(xué)，它還適合于其它只包括兩個狀態(tài)的動力學(xué)，例如自旋模型、投票模型[79]等。

主方程方法的不足是系統(tǒng)中演化方程的個數(shù)隨最大度的平方規(guī)模增長，當(dāng)最大度的值較大時，此方法的計算復(fù)雜度較高。

4邊滲流方法

4.1邊滲流方法概述

從圖論的角度來看，滲流理論描述的是在一個隨機圖中的連通簇的行為。在物理、化學(xué)和材料科學(xué)中，滲流理論被用來描述有孔材料中液體的流動，例如巖石中的石油等。同時滲流理論還可以用來研究路由器網(wǎng)絡(luò)故障、森林火災(zāi)、網(wǎng)絡(luò)上流行病的免疫、流行病傳播。在復(fù)雜網(wǎng)絡(luò)動力學(xué)研究中，滲流方法通常分作邊滲流與點滲流，點滲流即為刪除網(wǎng)絡(luò)中的節(jié)點(同時刪除連接到這個節(jié)點的邊)的過程，例如網(wǎng)絡(luò)中路由器的故障、在流行病傳播過程中個體的免疫，都可以用點滲流來描述。而邊滲流即為刪除節(jié)點之間連邊的過程。本節(jié)將主要介紹邊滲流的方法在求解SIR及SIS模型閾值上的應(yīng)用，邊滲流方法從網(wǎng)絡(luò)中邊的角度去考慮網(wǎng)絡(luò)中流行病的傳播。

對于SIR模型，系統(tǒng)的長期行為(例如網(wǎng)絡(luò)中最終的感染密度等不包含時間演化過程的性質(zhì))可以近似等價地映射為網(wǎng)絡(luò)中的邊滲流過程。用滲流來研究流行病傳播的方法最早由Mollison[80]和Grassberger[81]提出，并由Newman等[82-84]發(fā)展完整。本節(jié)會簡要回顧Newman在任意度序列的隨機圖上得到的一些基本結(jié)論[7]。另外，先前的理論認為只有SIR模型傳播動力學(xué)才可以映射為滲流過程，然而有研究者發(fā)現(xiàn)只需要引入邊的重感染概率[27]就可以解析SIS模型的傳播動力學(xué)。

4.2對SIR模型的求解

在SIR模型中，感染節(jié)點以λ的傳播率傳播流行病，并且假設(shè)節(jié)點的染病周期為τ，即節(jié)點經(jīng)歷確定的τ時間步后從I狀態(tài)恢復(fù)為R態(tài)(這與節(jié)點以1/τ的概率恢復(fù)有較大的差異，以概率恢復(fù)會導(dǎo)致更大的隨機性，一定程度上破壞了平均假設(shè)，會使得邊的滲流概率呈現(xiàn)出一定的異質(zhì)性)。假設(shè)動力學(xué)演化過程中時間步是連續(xù)的(對于離散情況引入一個無窮小量Δt→0即可)，那么對于任意一條邊，流行病在τ時間內(nèi)沒有通過其成功傳播的概率為

(105)

那么，流行病通過任意一條邊成功傳播的概率為

φ=1-e-λτ

(106)

如果令網(wǎng)絡(luò)中的邊分別以概率φ和1-φ被占據(jù)和不占據(jù)，那么流行病傳播過程就相當(dāng)于一個邊滲流概率為φ的過程：被占據(jù)即意味著這條邊足以傳播流行病，不被占據(jù)則反之。當(dāng)一條邊被占據(jù)，意味著當(dāng)流行病達到這條邊任意某一端的時候，這條邊足以傳播流行病，而非一定會傳播流行病。因此，被占據(jù)的邊意味著潛在的傳播可能性。當(dāng)流行病從初始感染節(jié)點開始傳播，可以看到流行病最終可以傳播到的節(jié)點恰好就是那些與初始節(jié)點之間有被占據(jù)的邊組成通路的那些節(jié)點。因此，傳播的最終狀態(tài)就是初始節(jié)點所在的邊滲流簇中的節(jié)點處于恢復(fù)態(tài)(R)。由此，流行病爆發(fā)的閾值與相對應(yīng)的邊滲流的閾值是相同的。

以任意度序列的隨機網(wǎng)絡(luò)為例。設(shè)網(wǎng)絡(luò)的度分布為pk，u為一個節(jié)點不通過某條邊連接到最大簇的平均概率。這里假定每條邊的占有概率是相等的，那么，計算u時有兩種情形須要考慮：一是這條邊沒有被占據(jù)；二是這條邊被占據(jù)，但是這條邊另一端的節(jié)點同樣不屬于最大連通簇。而第二種情況當(dāng)且僅當(dāng)另一端的節(jié)點不通過任意一條其他的邊連接到最大連通簇才成立。如果另一端的節(jié)點有k條剩余邊(亦即除了連接到節(jié)點的這條邊之外其它邊的數(shù)目)，那么發(fā)生的概率為uk。因此不通過某條邊連接到最大簇的總概率為1-φ+φuk。因此，針對k取平均，就可以得到一個自洽方程。

(107)

qk=(k+1)pk+1/〈k〉

(108)

根據(jù)式(106),(108)及式(107)導(dǎo)函數(shù)可以圖像地解出傳播相變的閾值。式(107)的解為y=u與y=1-φ+φg1(u)這兩個函數(shù)的交點。顯然u=1恒成立為一個平凡解。只有當(dāng)存在u=1之外的解時，網(wǎng)絡(luò)中才存在一個最大連通簇使得流行病爆發(fā)。因為qk為概率必有qk大于0，所以對于u≥0，g1(u)及其各階導(dǎo)數(shù)必然非負。因此，1-φ+φg1(u)與u在u=1處恰相切(此即為相變點)，對式(107)等號兩邊進行求導(dǎo)可得

(109)

再將式(106)帶入到式(109)中，即可得到

(110)

因此當(dāng)λτ高于這個閾值ln(〈k2〉-〈k〉)/(〈k2〉-2〈k〉)的時候，即很有可能發(fā)生流行病的爆發(fā)，因為初始感染者有可能沒有落在最大連通簇里。而當(dāng)λτ小于這個值的時候，流行病則一定不會爆發(fā)。

當(dāng)然，在上面介紹的經(jīng)典滲流方法中，將SIR模型映射為邊滲流只能反映系統(tǒng)長期的行為，而無法包含隨時間實時演化的信息，這也是上述方法最大的局限所在[83,85]。

4.3對SIS模型的求解

先前鮮有理論將SIS傳播動力學(xué)映射到滲流過程，原因在于節(jié)點被感染后會以一定概率恢復(fù)，這樣流行病的傳播就不再能精確地等效于一個簡單的加邊過程。但是Parshani和Havlin等[27]創(chuàng)造性地運用子節(jié)點對父節(jié)點的重感染概率，有效地克服了這一問題，對于SIS模型傳播閾值的預(yù)測也比一般的平均場方法更加準(zhǔn)確。具體做法為

首先計算每一條邊被占據(jù)的概率

p=1-(1-λ)τ

(111)

然后從一個相對微觀的角度來考慮流行病傳播的過程：當(dāng)網(wǎng)絡(luò)不存在度關(guān)聯(lián)性時，對于任意一個節(jié)點j，它被度為k的感染節(jié)點i(節(jié)點i會有1條入邊，k-1條出邊；此處入邊與出邊是針對動力學(xué)過程而言的，網(wǎng)絡(luò)仍然是無向網(wǎng)絡(luò))觸及到的概率是kP(k)/〈k〉。只要網(wǎng)絡(luò)是樹形結(jié)構(gòu)(或者環(huán)極少進而可以被忽略時)，在SIR模型中，平均說來被節(jié)點i所感染的節(jié)點數(shù)目[85]

(112)

為了簡便起見，在后面的公式中沿用這一代換κ-1=〈k2〉/〈k〉-1。

然而對于SIS模型，上述的方程不再準(zhǔn)確。在SIS傳播模型中，流行病傳播過程不再是一個有向無環(huán)的樹形分支過程：網(wǎng)絡(luò)中的連邊可能會多次地傳播流行病，更為重要的是先前滲流的邊也會恢復(fù)成非占有態(tài)。因此，直觀地說來這應(yīng)當(dāng)會導(dǎo)致最終穩(wěn)態(tài)下某條邊滲流的難度增大，使SIS模型的閾值在同等參數(shù)條件下相較SIR模型要小些，也就是說SIS模型中的流行病更容易大范圍地傳播開。對應(yīng)到具體的數(shù)學(xué)解析上，研究者引入子節(jié)點對于父節(jié)點重感染概率π來克服上述的問題，這樣相應(yīng)的滲流方程為

〈ni〉SIR=p(κ-1)+π

(113)

這里對于重感染概率的求解是關(guān)鍵點。大道至簡，復(fù)雜的問題往往能夠溯源到一個最簡單的模型：

以圖7中的i節(jié)點為例，它在感染j1節(jié)點之后，會在之后一定時間步之后恢復(fù)，具體值取決于i感染j1發(fā)生在i被感染后的第幾個時間步，設(shè)這一變量為s，那么在接下來的(τ-s)步時間內(nèi)，i仍處于感染態(tài)，而i恢復(fù)之后j1尚處于感染態(tài)，因此i、j1在感染態(tài)上存在一個時間差τ-(τ-s)=s。在這s時間步內(nèi)，j1節(jié)點在每一時間步都會以概率λ去感染i節(jié)點。據(jù)此寫出子節(jié)點j1感染其父節(jié)點i的概率的方程

圖6　感染密度的時序演化[42]

圖7　感染過程示意圖[27]

(114)

從式(114)可以看出其過程，由于i節(jié)點已將j1節(jié)點感染所以整個式子首先是一個條件概率P(j→i|i→j)，節(jié)點i感染j1節(jié)點的概率為1-(1-λ)τ，故而在條件概率中這一項會出現(xiàn)在分母上；再來看分子，既然i在第s步感染了節(jié)點j1，那么在前s-1步，都未感染成功：(1-λ)s-1λ；再接下來，就是j1在余下的s步時間內(nèi)感染節(jié)點i：1-(1-λ)s。式子展開后是一個典型的等比數(shù)列，化簡后即得右式。

然而這一過程仍然略為理想，因為除去j1節(jié)點，i節(jié)點的其它子節(jié)點j2—jk以及其父節(jié)點a也會以πt的概率去感染i節(jié)點，所以在i的鄰居節(jié)點之間存在一個競爭

(115)

其中πt表示了j1感染其父節(jié)點i的概率，求和符號中每一項表示除了節(jié)點j1外，i感染了它的k′個鄰居并被它們重感染的概率。對于這k′+1個重感染節(jié)點，僅有第一個感染i的節(jié)點會產(chǎn)生作用，其他k′個鄰居便不會再有任何影響。既然這些節(jié)點都處于等同地位，那么將相應(yīng)值除以k’+1就可以得到先前被i感染的子節(jié)點重復(fù)感染i的概率。

當(dāng)〈ni〉SIS=1時，由式(111),(113)～(115)即可解得閾值λc。并且在圖8中依次比較了不考慮重感染概率(p對應(yīng)的曲線)、考慮重感染但不考慮競爭(πt對應(yīng)的曲線)以及考慮競爭效應(yīng)的重感染情況(π對應(yīng)的曲線)。另外Parshani和Havlin等[27]還研究了不同網(wǎng)絡(luò)結(jié)構(gòu)、以及SIS與SIR模型不同的動力學(xué)對于傳播閾值的影響(見圖9)。

圖8　3種近似方法對應(yīng)的理論觀測值(相應(yīng)曲線)與模擬結(jié)果(點狀圖)[27]

圖9　不同情況下流行病爆發(fā)閾值與的關(guān)系[27]

4.4邊滲流方法的適用性分析及最新進展

在原有的理論框架下，當(dāng)節(jié)點感染后不以一定概率而是經(jīng)過特定時間步后恢復(fù)為健康態(tài)，則SIR模型下的流行病傳播可以映射為邊滲流[55]；當(dāng)節(jié)點以一定概率恢復(fù)時，其情況較為復(fù)雜；在4.3中介紹的方法認為重感染率π是此種概率恢復(fù)情形中的上界，故而仍可較準(zhǔn)確地描述此類的傳播問題[27]。

總體來說邊滲流方法相較于傳統(tǒng)的平均場方法在解析SIR模型下傳播閾值時更為精確。但是Newman在文獻[82]中的經(jīng)典方法也有一定的局限性：例如對于“網(wǎng)絡(luò)大小有限、網(wǎng)絡(luò)中存在簇會導(dǎo)致較強的動力學(xué)相關(guān)性)、流行病傳播過程隨時間演化、感染和恢復(fù)概率異質(zhì)不恒定”等問題都不能求解。然而Newman在文獻[82]中的方法中的這些局限性在研究者們的不斷努力下，一定程度上都得到了克服：Marder[85]將邊滲流的方法做了適當(dāng)?shù)母淖?，使得能觀察到流行病隨時間的演化過程；Pierre-André No?l等在文獻[86]的基礎(chǔ)上研究了有限大小網(wǎng)絡(luò)上流行病的演化過程；Miller[88]研究了對于異質(zhì)的感染和恢復(fù)概率的問題，求出了感染比例的上界和下界。Newman[89-90]給出了網(wǎng)絡(luò)中存在聚類時的解析方法，文獻[83],[91]研究了兩種流行病在單個網(wǎng)絡(luò)上傳播的情形。Allard等[92]則研究了多重屬性網(wǎng)絡(luò)上的流行病傳播，即使網(wǎng)絡(luò)中每條邊的占有概率都不同，也能精確求解其閾值。Funk等[93]將此方法推廣到重疊網(wǎng)絡(luò)上，研究了兩種流行病在網(wǎng)絡(luò)上的傳播過程。Dickison, Wang Y等[94-95]則將滲流的方法應(yīng)用到理解雙層互連接網(wǎng)絡(luò)上的流行病動力學(xué)?？梢哉f滲流方法自身的擴展性較強，能夠克服自身先前的各種局限，在傳播動力學(xué)的解析方面也發(fā)揮著越來越重要的作用。

5空穴理論

空穴理論最早由Perl[96]提出，爾后由Mézard等[97-98]將其應(yīng)用到統(tǒng)計力學(xué)中。物理學(xué)中的空穴是指電子掙脫價鍵的束縛成為自由電子后，其價鍵中所留下的空位?？昭ɡ碚撛谖锢怼⒒瘜W(xué)、生物、計算機科學(xué)中有著廣泛的應(yīng)用，僅在網(wǎng)絡(luò)科學(xué)領(lǐng)域中，就有用于解析流行病傳播動力學(xué)[99]、計算有向網(wǎng)絡(luò)中環(huán)的個數(shù)[100]、求解圖的二分問題[101]、計算最大連通子圖[102]、分析網(wǎng)絡(luò)中節(jié)點的作用[103-104]等等。

利用空穴理論的思想求解SIR流行病傳播動力學(xué)最早是由Karrer等[99]提出的。他們在求解過程中，會將原先的無向網(wǎng)絡(luò)轉(zhuǎn)換成有向網(wǎng)絡(luò)，其中邊的方向表示流行病的傳播方向。由于處于S態(tài)節(jié)點不能感染它的鄰居節(jié)點、只能被它的I態(tài)鄰居節(jié)點感染，也就是說S態(tài)節(jié)點只有入邊而沒有出邊。借用空穴理論的思想該S態(tài)節(jié)點處于空穴態(tài)，對應(yīng)到物理過程中處于空穴態(tài)的節(jié)點只能作為受體接觸電子。因此，求解一個節(jié)點t在時刻處于S態(tài)的概率，也就是計算它在t時刻處于空穴態(tài)的概率，再根據(jù)SIR模型中3種狀態(tài)的轉(zhuǎn)換關(guān)系就可以求解出節(jié)點處于每種狀態(tài)的概率以及流行病傳播的爆發(fā)閾值。

此前研究SIR模型的傳播動力學(xué)時，為了簡化解析過程，通常會假設(shè)恢復(fù)時間服從指數(shù)分布，恢復(fù)過程為泊松過程[23,105]，然而這一假設(shè)與現(xiàn)實生活很不相符。舉例來說，季節(jié)性流行病恢復(fù)時間服從泊松分布而并非指數(shù)分布，其均值通常在一周左右[99]。這使得人們對流行病傳播過程的描述無論從定性理解還是定量分析上都有一定的偏差，Karrer等[99]利用空穴理論的思想解決了這一難題。

圖10　空穴理論網(wǎng)絡(luò)示意圖

5.1SIR流行病模型求解過程

本小節(jié)介紹Karrer等[99]如何利用空穴理論求解無向網(wǎng)絡(luò)上SIR模型流行病動力學(xué)。如圖10a所示，流行病在傳播過程中，節(jié)點既可能被鄰居節(jié)點感染，也可能感染鄰居節(jié)點。為了描述流行病的傳播方向，他們將網(wǎng)絡(luò)中的無向邊變?yōu)閮蓷l有向邊(如圖10b)，邊的方向表示流行病的傳播方向。由于處于S態(tài)節(jié)點不能感染鄰居節(jié)點，只能被鄰居節(jié)點感染。因此，刪除從它出發(fā)的所有邊，則該節(jié)點處于空穴態(tài)(如圖10c中節(jié)點i為空穴態(tài))；此時，原來的無向網(wǎng)絡(luò)已變?yōu)橛邢蚓W(wǎng)絡(luò)。求解任意節(jié)點在t時刻處于S態(tài)的概率時，只需要考慮它在初始時刻處于S態(tài)(空穴態(tài))的概率且在t時間步內(nèi)沒有被鄰居節(jié)點感染的概率即可。接下來將考察兩種網(wǎng)絡(luò)上的流行病傳播動力學(xué)：樹形網(wǎng)絡(luò)和非樹形網(wǎng)絡(luò)。

5.1.1樹形網(wǎng)絡(luò)上的求解

所謂樹形網(wǎng)絡(luò)，是指網(wǎng)絡(luò)中不存在有限環(huán)即環(huán)的邊數(shù)為有限大小)。例如配置網(wǎng)絡(luò)中節(jié)點之間隨機連邊，且網(wǎng)絡(luò)規(guī)模趨于無窮時，就可以近似地看作樹形結(jié)構(gòu)或者局部樹形結(jié)構(gòu)[7,82]。相對于非樹形網(wǎng)絡(luò)而言，樹形網(wǎng)絡(luò)上的流行病傳播研究過程要簡單一些。如果節(jié)點i在t時刻處于S態(tài)(空穴態(tài))，那么只可能是因為它在初始時刻處于S態(tài)并且在之后t時間步內(nèi)沒有被鄰居感染，因而相應(yīng)概率為

(116)

其中，Si(t)表示節(jié)點i在t時刻處于S態(tài)，P(Si(t))表示相應(yīng)概率，簡記為P(Si)；z表示節(jié)點在初始時刻處于S態(tài)的概率；N(i)表示節(jié)點i的鄰居構(gòu)成的集合；Hi→j(t)表示節(jié)點在j時間t內(nèi)沒有成功感染節(jié)點i的概率。等式(116)給出了節(jié)點i在t時刻處于S態(tài)的概率，但是式子含有未知量Hi→j(t)。

根據(jù)動力學(xué)機制可知，只有I態(tài)節(jié)與它的S態(tài)鄰居接觸才可能傳播流行病。因此，一個節(jié)點在被感染后在t到dt內(nèi)，感染它的一個鄰居的概率為

(117)

對式(117)左右兩邊對t積分，可得一個I態(tài)節(jié)點在恢復(fù)之前將流行病傳給它的S態(tài)鄰居的概率，即Newman等[82]所說的有效傳播率。其中，s(t)dt表示被節(jié)點感染后，t到t+dt內(nèi)感染一個S態(tài)鄰居的概率；r(t)dt表示節(jié)點被感染后，t到t+dt內(nèi)恢復(fù)的概率。

(118)

(119)

對等式(119)進行數(shù)值求解，就可以得到節(jié)點i在t時刻處于I態(tài)的概率。由于節(jié)點在任意時刻只能處于3種狀態(tài)中的某一種，即P(Si)+P(Ii)+P(Ri)=1成立，再聯(lián)立等式(116)和等式(119)就可求出節(jié)點i在t時刻處于各種狀態(tài)的概率。

5.1.2非樹形網(wǎng)絡(luò)

在實際網(wǎng)絡(luò)中，由于邊的連接并非隨機，樹形網(wǎng)絡(luò)幾乎不存在[7,99]。即使在配置網(wǎng)絡(luò)中，當(dāng)網(wǎng)絡(luò)規(guī)模有限時也會導(dǎo)致網(wǎng)絡(luò)中出現(xiàn)有限環(huán)[108]?？昭ɡ碚撝荒芫_求解樹形網(wǎng)絡(luò)上SIR模型，對于非樹形網(wǎng)絡(luò)只能給出流行病傳播范圍的上界。下面將介紹利用空穴理論求解非樹形網(wǎng)絡(luò)上SIR模型的傳播動力學(xué)。

類似于樹形網(wǎng)絡(luò)上流行病傳播的求解過程，首先將原來的無向網(wǎng)絡(luò)變?yōu)橛邢蚓W(wǎng)絡(luò)，邊的方向表示流行病傳播方向，S態(tài)節(jié)點在傳播過程中任意時刻都為空穴態(tài)。此時，還需要對節(jié)點和邊都賦予一定的權(quán)值。其中，節(jié)點的權(quán)值τi表示它被感染后所需要的恢復(fù)時間為τi；有向邊的權(quán)值wij表示j節(jié)點被感染后，感染它所指向的節(jié)點i所需要耗費的時間為wij。根據(jù)恢復(fù)時間分布函數(shù)r(τ)對每個節(jié)點賦予權(quán)值，同理根據(jù)分布函數(shù)s(τ)對每條邊賦予權(quán)值。在對邊賦予權(quán)值時需要注意，若wij<τj表明節(jié)點i被節(jié)點j感染時，節(jié)點j還沒有恢復(fù)為R態(tài)，這表明節(jié)點j可以成功感染節(jié)點i，此時為有向邊j→i賦值wij；否則，表示節(jié)點j在其處于感染態(tài)的時間步內(nèi)無法將節(jié)點i感染，此時令wij=∞。

假設(shè)網(wǎng)絡(luò)中任意兩點l，k存在一條從l到k的有向路徑，其路徑長度Lkl表示節(jié)點l被感染后節(jié)點被k感染需要等待的時間為Lkl。因此，求解節(jié)點i在t時刻處于S態(tài)的概率就變得相當(dāng)容易，即只需要滿足：它在初始時刻沒被感染，相應(yīng)概率為z；不存在初始時刻被感染且它到節(jié)點i的有向路徑長度小于t的節(jié)點j。

不妨以節(jié)點i為中心，將所有到節(jié)點i有向路徑小于t的節(jié)點構(gòu)成集合記為Ω(ni表示集合Ω中節(jié)點數(shù)目)。因此，集合Ω中任意一個節(jié)點k到節(jié)點i的有向路徑都小于t。若集合Ω中所有節(jié)點在初始時刻都沒有被感染，那么節(jié)點i在t時刻也不會被Ω中的節(jié)點感染。此時，節(jié)點i在t時刻處于S態(tài)的概率為zni+1。對所有的wij和τj進行平均，則節(jié)點i在t時刻處于S態(tài)的概率為

P(Si)=z〈zni〉

(120)

其中，〈…〉表示對所有wij和τj進行平均。

(121)

圖11　網(wǎng)絡(luò)上的流行病傳播示意圖

由切比雪夫不等式[99]，式(121)可化為

(122)

其中Hi←j(t)=〈znij〉，表示時間t內(nèi)節(jié)點j沒有將流行病傳給節(jié)點i的概率。

在計算P[Sj(t′)|i處于空穴態(tài)]時，可以利用不等式(122)。但是，此時節(jié)點i和節(jié)點j都處于空穴態(tài)，需要做適當(dāng)?shù)奶幚?，即添加從?jié)點i出發(fā)到除節(jié)點j之外的所有鄰居的向邊(圖10d)，相當(dāng)于原網(wǎng)絡(luò)中只刪除了從節(jié)點i到節(jié)點j的有向邊。在做此處理后，勢必增加了到達節(jié)點j的路徑數(shù)，從而增加了節(jié)點j被感染的概率。于是有

P[Sj(t′)|i處于空穴態(tài)]≥P[Sj(t′)|刪除從i到j(luò)的有向邊]

(123)

再聯(lián)立不等式(122)，有

(124)

成立。不等式(124)和等式(118)形式很相像，只是這里是不等式。因此，無法精確計算出網(wǎng)絡(luò)中S態(tài)節(jié)點的比例。

(125)

經(jīng)過一次迭代后，令不等式(125)右邊為

(126)

它對于任意的τ都有f(τ)≥0成立，根據(jù)不等式(126)有表達式

(127)

成立。再根據(jù)不等式(126)有表達式

(128)

成立。經(jīng)過m次迭代，得到關(guān)系式

(129)

(130)

利用等式(122)，有

(131)

表達式(131)給出了在非樹形網(wǎng)絡(luò)中P(Si)滿足的關(guān)系，進而可以求解出P(Si)最小值和P((Ii)+P(Ri)的最大值，即網(wǎng)絡(luò)中被感染節(jié)點的最大比例。

5.1.3配置網(wǎng)絡(luò)上的流行病傳播

配置網(wǎng)絡(luò)的度分布P(k)的生成函數(shù)為

(132)

剩余度分布q(k)的生成函數(shù)為

(133)

因此，在任意節(jié)點在時刻t，節(jié)點沒有被它的某個鄰居感染的平均概率

(134)

進一步可知，t時刻網(wǎng)絡(luò)中處于S態(tài)節(jié)點的比例為

(135)

(136)

即為沒有被鄰居感染的概率。根據(jù)表達式(131)t時刻節(jié)點i處于S態(tài)的概率為

(137)

當(dāng)t→∞時，網(wǎng)絡(luò)中沒有節(jié)點處于I態(tài)，只有S態(tài)和R態(tài)節(jié)點，即P(Ri)=1-P(Si)成立。對于配置網(wǎng)絡(luò)，當(dāng)網(wǎng)絡(luò)規(guī)模很大時可以近似看成樹形網(wǎng)絡(luò)，此時Fi←j和等式(116)的Hi←j含義相同，而非樹形網(wǎng)絡(luò)根據(jù)表達式(130)都有Fi←j≤Hi←j成立，即Fi←j是Hi←j的下界。因此對于樹形網(wǎng)絡(luò)，可精確地求出節(jié)點在穩(wěn)態(tài)時處于S態(tài)的概率；然而對于非樹形網(wǎng)絡(luò)而言，只能近似求出節(jié)點i在穩(wěn)態(tài)時處于S態(tài)的最小概率，即處于R態(tài)的最大概率。

當(dāng)t→∞時，等式(134)變?yōu)?/p>

H1=1-p+pzG1(H1)

(138)

等式(135)變?yōu)?/p>

P(S)=zG0(H1)

(139)

其中，H1=H1(∞)表示一條邊在整個傳播過程中沒有被傳播流行病的概率。根據(jù)等式(138)和等式(139)就可得穩(wěn)態(tài)時處于S態(tài)和R態(tài)的比例。根據(jù)等式(138)，可用圖解法求出pc=〈k〉/(〈k2〉-〈k〉)，這和Newman等[82]利用滲流理論得到的閾值相同。當(dāng)p>pc時，網(wǎng)絡(luò)中就會有一定比例的節(jié)點被感染；反之，則幾乎沒有更多的節(jié)點被感染。

最后，實驗?zāi)M結(jié)果也非常精準(zhǔn)地驗證了之前理論分析結(jié)果的正確性與準(zhǔn)確性。圖12給出了在ER網(wǎng)絡(luò)上感染密度隨時間變化的理論值和實驗結(jié)果。其中γ和β分別表示恢復(fù)率和感染率，因此r(t)=γe-γ t，s(t)=e-β t[7]。根據(jù)式(117)有f(t)=βe-(β+γ)t，將t′=t-τ帶入等式(134)可得

(140)

在等式(140)兩邊同時對t微分有

(141)

網(wǎng)絡(luò)節(jié)點數(shù)量N=105，平均度為3，模擬重復(fù)100次以上。紅色實線對應(yīng)根據(jù)等式(135)在每個t時刻求得的理論值，點狀圖對應(yīng)模擬值

初始條件H1(0)=1，因此等式(141)的解

所以數(shù)值地求解出H1的值，將其帶入等式(138)即可求解出網(wǎng)絡(luò)穩(wěn)態(tài)時S態(tài)節(jié)點的比例，進而也能求出I態(tài)和R態(tài)節(jié)點的比例。

5.2空穴理論總結(jié)及展望

對于流行病動力學(xué)的研究，最為重要的問題便是考察在什么條件下疾病會成為流行病，以及流行病的影響范圍。此節(jié)介紹的Karrer等[99]利用空穴理論求解SIR模型的方法，可以求解感染時間和恢復(fù)時間服從任意分布的情況，且能精確求解樹形網(wǎng)絡(luò)上的流行病傳播范圍以及動力學(xué)演化過程?？昭ɡ碚摕o疑成為人們認識流行病動力學(xué)的一個新工具。當(dāng)然也應(yīng)該清楚地認識到，對于非樹形網(wǎng)絡(luò)空穴理論只能求解出最大感染比例，而無法給出準(zhǔn)確的結(jié)果。

張昊等[107]在Karrer等[99]和Shiraki等[102]工作基礎(chǔ)之上，對節(jié)點度之間存在關(guān)聯(lián)性時的流行病控制問題做了研究。他們利用空穴理論得到的解析結(jié)果和仿真結(jié)果非常吻合，在一定程度上加深了我們對于流行病傳播的認識。Altarelli等[108]最近利用空穴理論解決了傳播源優(yōu)化問題。當(dāng)然，空穴理論對于SI、SIR、SEIR流行病傳播模型也都適用[99]。

6邊劃分方法

在研究復(fù)雜網(wǎng)絡(luò)上的流行病傳播時，人們往往會關(guān)注動力學(xué)過程、穩(wěn)態(tài)時感染密度和傳播閾值。異質(zhì)平均場、馬爾可夫過程、點對近似等建模方法已成為研究流行病傳播的主流工具。然而，這些方法在描述動力學(xué)過程時，微分方程組維數(shù)很高，難以求解和找到對應(yīng)的解析表達式。Volz[109]利用邊劃分的方法很好地解決了這一難題。他用低維非線性方程組描述了復(fù)雜網(wǎng)絡(luò)上SIR流行病傳播動力學(xué)模型，可求解在任意時刻流行病感染傳播范圍。

6.1配置網(wǎng)絡(luò)上SIR流行病傳播

在研究流行病傳播模型時，非常重要的問題之一就是隨機選擇一節(jié)點u，它處于S、I和R態(tài)的概率分別是多少。由于節(jié)點u是隨機選取的，當(dāng)網(wǎng)絡(luò)規(guī)模很大時它處于S態(tài)的概率就為網(wǎng)絡(luò)中S節(jié)點所占的比例S(t)，處于I態(tài)和R態(tài)的概率也是類似的。如果知道S(t)，就能很容易求出I(t)和R(t)。

節(jié)點u處于S態(tài)只有它的所有鄰居節(jié)點到此時刻都沒有成功感染它，其概率大小取決于它的鄰居數(shù)量、I態(tài)鄰居的變化率和鄰居處于I態(tài)的概率。由于節(jié)點u的二層鄰居數(shù)通常大于它的平均度，因此即使知道網(wǎng)絡(luò)中I態(tài)節(jié)點比例也不知道節(jié)點u的鄰居處于I態(tài)的概率。因此，求解節(jié)點u的鄰居節(jié)點處于I態(tài)的概率也就成為重點之一。一旦求出此概率，就很容易得到節(jié)點u處于S態(tài)的概率。

假設(shè)節(jié)點u的每個鄰居相互獨立。然而，節(jié)點u的一個鄰居被感染，它可能感染節(jié)點u，而節(jié)點u又可能會感染鄰居節(jié)點w，使得假設(shè)不成立。這里只需要假設(shè)節(jié)點u不傳播流行病即可保證它的鄰居之間相互獨立。雖然做了如此假設(shè)這并不影響節(jié)點處于各種狀態(tài)的概率。因為，節(jié)點u處于S態(tài)時，它的鄰居之間沒有關(guān)聯(lián)性，因此不影響它處于S態(tài)的概率。當(dāng)它處于I態(tài)時，鄰居節(jié)點不會相互獨立，但是這也不影響它處于I態(tài)的時間和恢復(fù)概率，因此不影響它處于I態(tài)和R態(tài)的概率。因此，不允許節(jié)點u傳播流行病，并假設(shè)它的鄰居相互獨立并不會影響它處于各種狀態(tài)的概率。

在介紹邊劃分求解方法時，還有必要說明只有當(dāng)初始時刻有足夠多的I態(tài)節(jié)點使得流行病在傳播的過程中不會受到隨機擾動而受影響時才會很準(zhǔn)確。當(dāng)然，這里足夠多的I態(tài)節(jié)點相對于整個網(wǎng)絡(luò)而言也是非常少的。

(142)

(143)

再聯(lián)立

(144)

為在配置網(wǎng)絡(luò)上SIR流行病的演化方程。根據(jù)等式(143)和等式(144)就可以求出在任意時刻節(jié)點處于各種狀態(tài)的概率。在這里只需要求解一個方程，其求解時間相對于異質(zhì)平均場、點對近似和馬爾可夫等方法更少。

(145)

(146)

當(dāng)R0>1時，等式(146)有兩個解，其中有一個解為1，另一個解小于1，小于才是此時尋求的解。再根據(jù)等式(144)可求出穩(wěn)態(tài)時流行病傳播范圍R(∞)=1-G(θ(∞))。

初始時刻選取1%的節(jié)點作為傳播源。理論值為虛線，模擬值為實線

6.2最新研究進展

在Volz 提出邊劃分方法后，Miller[110-111]對此方法做了進一步簡化，得到方程和。最近，Miller、Slim和Volz在綜述文獻[112],[114]中詳細介紹了利用邊劃分所解決的流行病傳播問題。他們在文獻[112]中主要介紹了利用邊劃分如何求解不同度分布、不同接觸模式下的SIR流行病傳播。文獻[113]中主要介紹對于不同的情況如何選取不同邊劃分模型，而文獻[114]重點在于介紹在其他情況下的應(yīng)用。比如，有向網(wǎng)絡(luò)上、異質(zhì)的感染率和恢復(fù)率[88]、節(jié)點類型多樣性網(wǎng)絡(luò)[115]、動態(tài)網(wǎng)絡(luò)[116]上的流行病傳播。Volz等[117]還利用邊劃分方法研究了異質(zhì)網(wǎng)絡(luò)和時空聚集模式，以及初始感染比例對動力學(xué)的影響[118]。不僅如此，Valdez等還利用此方法對自適應(yīng)流行病傳播模型[119]、動態(tài)滲流[120-121]做了相應(yīng)的研究。然而，這里介紹的邊劃分方法并不適用于SIS模型，楊紫陌與周濤[48]利用另外的邊劃分方法求解權(quán)重網(wǎng)絡(luò)上SIS單點接觸流行病傳播模型。

7馬爾可夫鏈方法

7.1馬爾可夫鏈概述

馬爾可夫鏈方法因安德烈·馬爾可夫(Markov A. A.)得名，在數(shù)學(xué)中是指具有馬爾可夫性質(zhì)的離散時間隨機過程：亦即在這一隨機過程中，預(yù)測未來的狀態(tài)只依賴于先前的部分狀態(tài)信息，而與更早之前的歷史狀態(tài)無關(guān)。通常可分為一步條件轉(zhuǎn)移概率與N步(或稱多步)條件轉(zhuǎn)移概率，前者下一時刻狀態(tài)只與當(dāng)前狀態(tài)和前一時間步有關(guān)；后者則與之前N個時間步狀態(tài)都有關(guān),例如物理學(xué)中矩磁的相變過程。對于服從一步條件轉(zhuǎn)移概率的指定個體，已知當(dāng)前狀態(tài)信息，根據(jù)馬爾可夫過程便能預(yù)測下一時刻的狀態(tài)信息。另外特別值得一提的是，馬爾可夫方法可以處理離散時間，不需要像其它方法假設(shè)時間連續(xù)。

馬爾可夫鏈?zhǔn)敲枋鰧ο蟮囊环N狀態(tài)序列，其每個時刻的狀態(tài)值取決于前面的狀態(tài)，是具有馬爾可夫性質(zhì)的隨機變量{X1，X2，…，Xt}的一個數(shù)列。這些變量的取值范圍(即它們所有可能取值的集合)，被稱為“狀態(tài)空間”(或“相空間”)；Xn的值表示相應(yīng)變量在t時刻所處的狀態(tài)，如果Xt+1對于過去狀態(tài)的條件概率分布僅是Xt的函數(shù)：亦即P(Xt+1=x|X1=x1，X2=x2，…，Xt=xt)=P(Xt+1=x|Xt=xt)，其中x的不同取值對應(yīng)過程中的不同狀態(tài)；那么上面這個恒等式可以被稱作具有馬爾可夫性質(zhì)，更為具體的是指“一步條件轉(zhuǎn)移概率”。

馬爾可夫鏈通常用于排隊理論和統(tǒng)計學(xué)中的建模，還可作為信號模型用于熵編碼技術(shù)(例如算術(shù)編碼)。馬爾可夫鏈還有眾多生物學(xué)方面的應(yīng)用，模擬生物人口過程是最為常見的。隱蔽馬爾可夫模型還被用于生物信息學(xué)，用以編碼區(qū)域或進行基因預(yù)測。此外，馬爾可夫鏈因其特有的馬爾可夫性質(zhì)，被廣泛應(yīng)用在病毒傳播的預(yù)測研究上，例如模擬植物病毒傳播[135]和解析流行病的爆發(fā)[136-137]。

P(Xu=xu|{Xv},v≠u}=P(Xu=xu|Xv,v∈Vu)

(147)

其中Vu表示節(jié)點u的鄰居節(jié)點所構(gòu)成的集合。在馬爾可夫隨機網(wǎng)絡(luò)中，只有節(jié)點對中兩節(jié)點間存在關(guān)聯(lián)性，因而可以用度分布P(k)和條件概率分布P(k′|k)這兩個參量來刻畫其上的動力學(xué)。這類網(wǎng)絡(luò)的馬爾可夫特征意味著，所有高階相關(guān)性都可以表示為P(k)和P(k′|k)的函數(shù)。

7.2基于馬爾可夫鏈的SIS傳播動力學(xué)建模

考慮無向連通網(wǎng)絡(luò)G=(V,E)，感染概率為β，感染節(jié)點的恢復(fù)概率為u。定義pi，t是節(jié)點i在t時刻處于感染態(tài)的概率,ξi，t是節(jié)點i在t時刻沒有被其鄰居感染的概率(在SIS模型中，將所有節(jié)點的狀態(tài)看成是一個狀態(tài)組合，由于每個節(jié)點狀態(tài)可能為易感態(tài)或感染態(tài)，所以N個節(jié)點的狀態(tài)組合共有2N種。換言之，可以將傳播過程看成是包含2N種不同狀態(tài)的馬爾可夫鏈：節(jié)點在t時刻的狀態(tài)只取決于t-1時刻的狀態(tài))，

(148)

其中Vi表示節(jié)點i的鄰居節(jié)點集合，節(jié)點i在t時刻沒有收到來自鄰居節(jié)點的感染包括兩種可能：1)鄰居節(jié)點j處于感染態(tài)，但是在t時刻沒能感染節(jié)點i，即式中pj，t-1(1-β)；2)鄰居j在t-1時刻處于易感態(tài)，即式中(1-pj，t-1)。

進一步，若滿足以下3種情況之一，可推知節(jié)點i在t時刻處于健康態(tài)：

1)i在t時刻之前是健康態(tài)，并且在t時刻沒有被其鄰居感染；

2)i在t時刻之前已被感染，并且在t時刻恢復(fù)同時沒有再次被其鄰居感染；

3)i在t時刻之前已被感染，在t時刻接觸到了感染態(tài)鄰居，在這些無效的接觸(因發(fā)生接觸時處于感染態(tài),因而不會被其鄰居感染,稱這些接觸是無效的)之后恢復(fù)為健康態(tài)。

因此節(jié)點i在t時刻處于健康狀態(tài)的概率為

1-pi，t=(1-pi，t-1)ξi，t+μpi，t-1ξi，t+μpi，t-1(1-ξi，t)/2，(i=1，2，…，N)

(149)

本節(jié)將主要關(guān)注基于上述馬爾可夫鏈模型的傳播動力學(xué)(我們前面各個章節(jié)所涉及的傳播動力學(xué)模型均可用馬爾可夫鏈來進行描述)解析方法：特征值分析法和微觀尺度上的離散時間馬爾可夫方法。兩種方法都基于節(jié)中的數(shù)學(xué)模型，其差別在于后續(xù)的解析過程中：前者主要通過非線性系統(tǒng)解的穩(wěn)定性進行相應(yīng)分析；而后者則是一種更為數(shù)值化的非線性系統(tǒng)迭代解析過程。

7.2.1特征值分析

7.2.1.1特征值分析方法

定義為網(wǎng)絡(luò)的鄰接矩陣，AT是矩陣A的轉(zhuǎn)置，At是矩陣A的t次冪，λi，A表示矩陣A的第i大特征值。下面不妨以一種更為數(shù)學(xué)化、更為嚴(yán)謹(jǐn)?shù)男问絹斫o出對于相應(yīng)結(jié)論的證明。

定理: 若一種流行病消亡，則有β/μ<τ=1/λi，A，其中τ表示流行病傳播閾值。

證明：由公式(149)可知，

(150)

基于近似公式，當(dāng)時a?1，b?1時，省略二階項(1-a)(1-b)≈1-a-b,上述證明中進行了兩處近似：

由公式(150)可得到

(151)

將式(151)轉(zhuǎn)換為矩陣表示法

Pt=((1-μ)I+βA)Pt-1=SPt-1=StP0

(152)

其中t時刻的各節(jié)點狀態(tài)空間Pt=(p1，t，…，pN，t)T，S=(1-μ)I+βA稱為系統(tǒng)矩陣。

引理：矩陣S的第i大特征值λi，S=1-μ+βλi，A，而且S的特征項量與A的完全相同。

引理證明：vi,A為矩陣A的特征值λi，A對應(yīng)的特征向量，由定義可知Avi,A=λi，Avi,A進而有

Svi,A=(1-μ)vi,A+βAvi,A=(1-μ)vi,A+β λi，Avi,A=(1-μ+β λi，A)vi,A

7.2.1.2穩(wěn)定性分析

在均勻網(wǎng)絡(luò)中，節(jié)點的度都比較集中在平均度附近(即k≈〈k〉，?k)，ρ(t)表示感染個體的平均密度。其平均場速率方程應(yīng)用平均場理論可寫作

dρ(t)/dt=-ρ(t)+λ〈k〉ρ(t)(1-ρ(t))

(153)

此近似方程在ρ(t)→0時(即處于傳播的初始階段)成立，原因已在平均場部分進行了解釋。當(dāng)t→∞時，由?ρ(t)=0即可解得傳播閾值λc=〈k〉-1[14]。當(dāng)時λ<λc時，ρ=0；若λ≥λc，則有ρ～(λ-λc)。

對于異質(zhì)網(wǎng)絡(luò)，節(jié)點度值的波動性很大，無法再使用全局的平均近似；此時定義度為k的感染個體密度為ρk(t)，相應(yīng)的反應(yīng)速率方程為[14-15]

dρk(t)/dt=-ρk(t)+λk(1-ρk(t))Θk(t)dt

(154)

然而對于節(jié)點之間有關(guān)聯(lián)性的網(wǎng)絡(luò)，上面的推導(dǎo)是不正確的，因為在Θk(t)中沒有考慮邊的連出節(jié)點的度k。對于馬爾可夫網(wǎng)絡(luò)，網(wǎng)絡(luò)的關(guān)聯(lián)性是完全由條件概率P(k′|k)來定義的，其中

(155)

kP(k′|k)P(k)=k′P(k|k′)P(k′)≡〈k〉/P(k,k′)

(156)

很明顯ρk=0是一個平凡解；在閾值附近，只有少許節(jié)點被感染(即ρk→0)，省略掉ρkρk′這類高階項，我們可以將式(149)線性化為

(157)

如果雅可比矩陣L至少有一個正特征值，那么ρk=0這個解就不穩(wěn)定[7]。依據(jù)相應(yīng)的定義，連接矩陣C中元素Ckk′=kP(k′|k)，由式(156)可得，若vk是C某一特征值對應(yīng)的特征向量，那么P(k)vk是C的轉(zhuǎn)置矩陣CT里對應(yīng)相同特征值λi，C的特征向量，由此可知C的所有特征值都為實數(shù)。假設(shè)λ1，A是C的最大特征值，那么只要雅可比矩陣L的最大特征值-1+λλ1，C>0，即λ>1/λ1,C時，ρk=0這個解就不穩(wěn)定，就會有另一個非零解成為實際的穩(wěn)態(tài)，即網(wǎng)絡(luò)處于染病態(tài)，由此可知傳播閾值λc=1/λ1，C。

可以發(fā)現(xiàn)，無論使用網(wǎng)絡(luò)的何種矩陣表示方式，都得到了定性相同且定量相同的結(jié)論。然而真實網(wǎng)絡(luò)往往都存在一定程度的度關(guān)聯(lián)性(例如同配性或異配性導(dǎo)致的度關(guān)聯(lián))，而且通常很難用一個閉合形式的條件概率P(k′|k)來刻畫。因此Wang等[123]提出了一個普適的病毒傳播分析模型，能夠刻畫底層網(wǎng)絡(luò)拓撲對傳播的影響，并且適合于任何網(wǎng)絡(luò)結(jié)構(gòu)。

對于均勻網(wǎng)絡(luò)或隨機網(wǎng)絡(luò)，其鄰接矩陣的最大特征值是平均度，因此傳播閾值為λc=1/〈k〉，這與平均場得出的結(jié)果一致；對于服從冪律分布且大小有限的網(wǎng)絡(luò)，文獻[7]通過對實際網(wǎng)絡(luò)和模型網(wǎng)絡(luò)中的閾值模擬來對比鄰接矩陣特征值方法和平均場方法確定閾值的準(zhǔn)確性，模擬結(jié)果明顯表明λc=1/λ1，A更準(zhǔn)確。Chung等[124]更從理論上計算出了服從冪律度分布且大小有限的網(wǎng)絡(luò)的鄰接矩陣最大特征值

(158)

(159)

具體計算中需要注意以下3點：1)如果初始向量x(0)剛好與最大特征值對應(yīng)的特征向量正交，那么此方法無效。由佩龍—弗洛比尼斯定理[7]可知，一個非負實矩陣的最大特征值非負，且最大特征值對應(yīng)的特征向量的元素都為非負，因此若初始向量x(0)的元素都為正數(shù)，則x(0)不會和最大特征值對應(yīng)的特征向量正交，所以選擇初始向量時應(yīng)該使其元素都是正數(shù)；2)向量的元素每次迭代之后都會增長，因此最終這些元素在計算機中會因為過大而溢出，因此需要周期性地對向量進行歸一化；3)判斷x(t)是否收斂為最大特征值對應(yīng)的特征向量時可以針對不同的初始向量計算兩次，由此在允許的誤差范圍內(nèi)判斷其何時收斂。

圖14　流行病爆發(fā)閾值隨傳播率的演化情況[14-15]

模擬結(jié)果均是在俄勒岡州自治系統(tǒng)網(wǎng)絡(luò)上進行的，其度分布服從冪律N=10 900,

在得到最大特征值對應(yīng)的特征向量后，再次乘以鄰接矩陣，由矩陣的特征值和特征向量的關(guān)系可知，兩個向量的任意元素的比值即為最大特征值，通常取幾個較大的元素的比值的平均值，這樣可以減小誤差、提高準(zhǔn)確性。

7.2.1.3特征值方法的準(zhǔn)確性及適用性分析

通過觀察對于閾值和感染密度的具體理論預(yù)測與數(shù)值模擬情況(詳見圖14與圖15)，通過求網(wǎng)絡(luò)鄰接矩陣的最大特征值來確定傳播閾值的方法，其結(jié)果比平均場的解析更為準(zhǔn)確，并且此方法既適用于實際網(wǎng)絡(luò)也適用于模型網(wǎng)絡(luò)。

通過以上分析，網(wǎng)絡(luò)鄰接矩陣的最大特征值的大小決定了閾值，此外最大特征值對于諸如信息傳播[127]、同步[128-129]、滲流[130-132]等其他動力學(xué)解析的影響也非常大。通過求解網(wǎng)絡(luò)鄰接矩陣的最大特征值來解析傳播閾值雖然普遍適用于各種網(wǎng)絡(luò)結(jié)構(gòu)，但當(dāng)網(wǎng)絡(luò)較大時很難直接快速求出其最大特征值，因此只能用近似的方法，例如上面提到的冪乘法，這會使結(jié)果的準(zhǔn)確性有所降低。

7.2.2離散時間的非線性系統(tǒng)數(shù)值迭代解析法

這一方法同樣基于7.2中的數(shù)學(xué)模型，但與7.3.1中分析解的穩(wěn)定性完全不同。

對于二狀態(tài)SIS傳播模型，目前的理論研究多集中于全接觸又譯作反應(yīng)過程,每個時間步,感染節(jié)點會接觸其全部鄰居節(jié)點)和單點接觸(每個時間步,感染節(jié)點只隨機接觸它的一個鄰居節(jié)點進行試圖的感染,若鄰居已處于感染態(tài)則跳過感染過程)兩種特殊的接觸過程；但實際上，節(jié)點之間的接觸過程則與具體問題本身密切相關(guān)，具有不確定性。文獻[134]采用基于接觸過程的概率離散馬爾可夫鏈模型，并將其分別應(yīng)用到含權(quán)和無權(quán)網(wǎng)絡(luò)，研究不同的接觸過程對流行病動態(tài)傳播的影響。

pi(t+1)=(1-qi(t))(1-pi(t))+(1-μ)pi(t)+μ(1-qi(t))pi(t)

(160)

對于由N個式(151)組成的非線性系統(tǒng)，應(yīng)用不動點迭代(或牛頓迭代)就可以得到對應(yīng)的數(shù)值解(若為便捷,可使用Matlab中的fsolve函數(shù)進行相應(yīng)的求解)，其與蒙特卡洛模擬的對比結(jié)果如圖16和圖17所示。圖16和圖17表明，不論是在微觀尺度(插圖中的具體節(jié)點概率)上，還是在不同接觸情形下，基于離散馬爾可夫鏈方法得到的數(shù)值解都與模擬結(jié)果吻合得很好，充分驗證了該方法的準(zhǔn)確性。

其中N=10 000，μ=1，λ→∞，

圖17　接觸個數(shù)取不同值時,航空網(wǎng)絡(luò)上

7.2.2.1上述方法相應(yīng)的拓展變形

為了能夠更準(zhǔn)確地描述流行病傳播動力學(xué)以及自適應(yīng)網(wǎng)絡(luò)的相應(yīng)情況，文獻[29]將退火鄰接矩陣(Annealed Adjacency Matrix, AAM,該矩陣中的每個元素Aij都是0-1之間的實數(shù),表示節(jié)點i-j相連的概率)引入到離散馬爾可夫鏈方法[126]中,得到節(jié)點i在t+1刻為健康態(tài)的退火方程：

si(t+1)=si(t)+μ(1-si(t))-

si(t)(1-qi(t))

(161)

將鄰接矩陣換成是與傳播過程一起隨時間演化的退火矩陣，文獻[29]進一步考慮了自適應(yīng)網(wǎng)絡(luò)上的馬爾可夫鏈SIS傳播模型，如圖18所示，結(jié)果再一次驗證了退火近似的準(zhǔn)確性，而且對于度較小的節(jié)點預(yù)測性也很高。

如果應(yīng)用平均場假設(shè)，認為度相同的節(jié)點具有相同的感染概率，那么?i，當(dāng)ki=k時，有si=ρk。從而可對方程(161)進行一個簡單代換，得出異質(zhì)平均場方程的一個變形(姑且用NewHMF與先前的方法進行一個區(qū)分)：

(162)

其中Nk=NP(k)。

通過對比發(fā)現(xiàn)，相較于一般的異質(zhì)平均場方程(HMF)，退火方程(Annealed,由式(161)迭代解析可得)和上述異質(zhì)平均場方法(New HMF,由式(162)迭代解析可得)的數(shù)值解與蒙特卡羅模擬結(jié)果(MC)有著很好的吻合(見圖19)。

一般來講，P(k，k′)=kP(k′|k)/(NP(k′))。從而，上述異質(zhì)平均場方程可以寫成如下形式：

(163)

當(dāng)對于任意的k，當(dāng)λ(1-ρk)?1時，將式(154)線性化后可得出與經(jīng)典HMF方法相同的解析式：

(164)

換句話說，HMF方程只在流行病爆發(fā)閾值附近才是精確的，在λ較大時它所使用的近似都會使得誤差產(chǎn)生。

相應(yīng)地，如果式(162)中使用的鄰接矩陣，對應(yīng)得到的方法被稱作非微擾異質(zhì)平均場方法(Non-Perturbative HMF)[38]。

7.2.3小結(jié)

回顧前面介紹的7種方法，可以發(fā)現(xiàn)大多數(shù)方法是在從網(wǎng)絡(luò)層面、在較宏觀的尺度分析問題，或者也只是精確到度層階(以度對網(wǎng)絡(luò)中的節(jié)點進行分類)的介觀尺度，只有點對近似和主方程精確到了點對的狀態(tài)這一較為微觀的視角；那么有沒有方法可以從微觀的角度去研究宏觀尺度上的一些性質(zhì)呢？馬爾可夫鏈方法帶給我們一個新的視角。它與常見的微分方程(以平均場方法為代表)有較大差異。從群體角度來看，經(jīng)典的平均場方法的確能夠較好地描述、預(yù)測傳播動力學(xué)的臨界現(xiàn)象，但往往不能給出指定個體的感染信息，然而在許多實際應(yīng)用和理論分析中, 特定個體的信息往往是人們所需要和關(guān)注的。馬爾可夫方法恰好能解決這一問題。

N=104，μ=1

N=5 000，〈k〉=4

7.3.2小節(jié)中介紹的數(shù)值迭代方法準(zhǔn)確性非常之高，其優(yōu)點在于幾乎沒有進行近似，保留了網(wǎng)絡(luò)結(jié)構(gòu)的幾乎所有信息，這也正是它明顯優(yōu)于其它方法的特點。然而，也正是由于這樣的特性，使得計算的時間復(fù)雜度較高，比如馬爾可夫鏈方法的狀態(tài)空間隨著節(jié)點個數(shù)呈指數(shù)增長,很難用于處理實際的大型網(wǎng)絡(luò)。這或許也是它最大的問題所在。另外，它的另一顯著優(yōu)點是可以處理離散時間，無須連續(xù)時間假設(shè)；這也克服了之前許多方法的一些局限，使得結(jié)果更為準(zhǔn)確。而基于退火矩陣的馬爾可夫鏈方法具有更好的拓展性，可以對基于節(jié)點狀態(tài)的動態(tài)演化網(wǎng)絡(luò)實行動力學(xué)分析，且對于大多數(shù)節(jié)點預(yù)測十分準(zhǔn)確；雖然所做的改動并不大，只是將鄰接矩陣變?yōu)榱送嘶鹁仃?，但其收效卻是頗有些“意想不到”的：基于退火鄰接矩陣(矩陣的每個元素值代表節(jié)點對之間的連接概率)的馬爾可夫鏈方法還能夠很好地用于動態(tài)網(wǎng)絡(luò)中的傳播動力學(xué)分析，如預(yù)測自適應(yīng)網(wǎng)絡(luò)中的感染密度分析。

8總結(jié)與討論

本文前述各節(jié)分析了如何具體求解具體網(wǎng)絡(luò)在特定動力學(xué)規(guī)則下的爆發(fā)閾值，細心的讀者或許已經(jīng)發(fā)現(xiàn)SIS和SIR的臨界感染分布是否存在本質(zhì)區(qū)別是一個聚訟已久的話題。接下來，我們將盡量清晰地梳理一下這一問題的發(fā)展脈絡(luò)并給出確定閾值的方法。

8.1SIS與SIR模型閾值求解的演進

(165)

(166)

很接近，而QMF閾值只對γ<5/2有效。局域化分析的相關(guān)工作[143-144]則證明度分布指數(shù)γ>5/2時SIS真實閾值有限，與異質(zhì)平均場結(jié)論近似。文獻[145]指出任何度分布衰減小于指數(shù)的小世界網(wǎng)絡(luò)上都不存在流行病。

較早的SIR閾值結(jié)論等同于SIS[1]，即無關(guān)聯(lián)的均勻網(wǎng)絡(luò)上的SIR傳播閾值為λc=1/〈k〉。Moreno等[147]考慮節(jié)點度的異質(zhì)性，推得復(fù)雜網(wǎng)絡(luò)上SIR傳播閾值與先前SIS相同λc=〈k〉/〈k2〉。如此得出，二階矩〈k2〉對傳播閾值起著決定性作用。在無標(biāo)度網(wǎng)絡(luò)結(jié)構(gòu)中，當(dāng)N趨于無窮大時二階矩發(fā)散，傳播閾值為0。Boguna等[47]對此進行了進一步補充，發(fā)現(xiàn)考慮節(jié)點度的異質(zhì)性時，在無關(guān)聯(lián)網(wǎng)絡(luò)上的SIR閾值為

(167)

關(guān)聯(lián)網(wǎng)絡(luò)上的SIR閾值為最大特征值的倒數(shù)

(168)

上述研究都是針對全接觸模型，有研究者單獨研究了單點接觸情況下的閾值行為。Zanette等[151]考慮基于單點接觸的SIR，發(fā)現(xiàn)當(dāng)小世界網(wǎng)絡(luò)的隨機性p比較大時，SIR感染范圍分布形式為冪律分布尾部截斷。冪律指數(shù)約為-1.5，分布的尾部呈現(xiàn)二次凸起，且直觀認為最大凸起對應(yīng)的感染范圍為Rmax～0.25N，隨著網(wǎng)絡(luò)規(guī)模N線性增長。Ben-Naim等[152-153]基于均勻混合假設(shè)考慮SIR閾值點的爆發(fā)規(guī)模。得出閾值點的感染范圍服從冪律分布P(n)～n-3/2。對于規(guī)模為N的網(wǎng)絡(luò)，SIR的最大感染范圍n～N2/3，平均感染范圍〈n〉～N1/3，流行病最長生存時間t*～N1/3，平均生存時間〈t〉～lnN?；陔S機游走的思想，文獻[154]再次驗證了上述結(jié)論。Khaleque[155]考慮空間網(wǎng)絡(luò)(節(jié)點距離異質(zhì)分布)上的SIR傳播，模擬得出在閾值之下感染范圍呈指數(shù)衰減；閾值點為冪律分布，閾值之上呈雙峰分布。Petter Holme[156]發(fā)現(xiàn)，SIR的流行病最長生存時間對應(yīng)的感染概率大于且獨立于傳播閾值。

8.2對于SIS和SIR模型傳播閾值的模擬判定方法

鑒于理論傳播閾值的重要性，對理論模擬閾值進行準(zhǔn)確檢測便顯得尤為重要。然而如何確定一個傳播模型的模擬閾值，這一問題研究者們一直爭論不休，提出了許多不同的指標(biāo)。我們同樣對其進行一個梳理。

鑒于易感性對估計SIS傳播閾值的有效性，我們嘗試將其進一步應(yīng)用于SIR傳播模型。基于隨機規(guī)則網(wǎng)絡(luò)，我們對SIS和SIR傳播過程進行了大量模擬，并根據(jù)易感性

χ=N(〈ρ2〉-〈ρ〉2)/〈ρ〉[141]

(169)

的峰值估計傳播閾值?；谝赘行阅軌驕?zhǔn)確估計SIS的傳播閾值[12]，但對估計SIR傳播閾值失效。

在磁化系統(tǒng)中，通常采用序參量的比值來確定平衡相變的臨界點，其一般定義[159]為

(170)

顯然，當(dāng)q=s=1時即對應(yīng)可變性Δ：

(171)

另外我們的研究小組發(fā)現(xiàn)，實際網(wǎng)絡(luò)的理論閾值與模擬閾值存在明顯區(qū)別。絕大多數(shù)網(wǎng)絡(luò)中，由可變性和易感性兩種方法確定的SIS模擬閾值差不多(具有明顯負關(guān)聯(lián)性的網(wǎng)絡(luò)除外)，但二者確定的SIR模擬閾值存在很大差別。特別是，由易感性確定的SIR閾值總是偏離理論閾值很多，由可變性確定的模擬閾值則更偏向于異質(zhì)平均場(整體看來，實際網(wǎng)絡(luò)的特征值閾值要小于異質(zhì)平均場閾值)。當(dāng)特征值閾值與模擬閾值差很多時，對應(yīng)的異質(zhì)平均場和特征值閾值也差很多。通過計算IPR參數(shù)[143]可以發(fā)現(xiàn)，模擬閾值和特征值閾值較吻合時IPR值較小，二者不一致時IPR值較大[161]。也就是說，特征值閾值和模擬閾值不一致應(yīng)該是特征向量局域化所致。另外，當(dāng)網(wǎng)絡(luò)呈負關(guān)聯(lián)(如ego-facebook,AS2000)時，平均場和特征值閾值都與模擬閾值不一致，這與無標(biāo)度網(wǎng)絡(luò)上的結(jié)果分析基本一致。

8.3討論

網(wǎng)絡(luò)動力學(xué)作為復(fù)雜網(wǎng)絡(luò)研究中的重要課題，具有極其重要的意義。我們在文章中介紹的方法大多都可以應(yīng)用于其它動力學(xué)過程，如投票者模型、網(wǎng)絡(luò)同步等等，而不只局限于流行病的傳播。

在這篇綜述中，縱觀我們介紹的這七種方法，其所應(yīng)用的數(shù)理方法大致可分為3類：微分方程、生成函數(shù)、基于馬爾可夫鏈方法。這些看似不同的方法，卻都能比較準(zhǔn)確地描述現(xiàn)實網(wǎng)絡(luò)中動力學(xué)的行為特征；然而也應(yīng)當(dāng)看到不同的方法所對應(yīng)的解析難度以及分析角度著實有較大不同。

基于微分方程的3種方法(平均場近似、點對近似和主方程)里，主方程方法考慮的動力學(xué)過程、網(wǎng)絡(luò)拓撲等情況最為全面，也最為精確；然而它的計算機復(fù)雜性也是最高的，其維數(shù)很高，無法像平均場方法那樣在人力可解的情況下給出一個還算不錯的預(yù)測。這3種方法都可以描述動力學(xué)的時間演化過程，相比之下，經(jīng)典的滲流理論和矩陣分析在這一方面略遜一籌。然而這3種方法它們都假設(shè)時間是連續(xù)的，而馬爾可夫方法則很好地克服了這一問題，而且通過它的框架也可以衍生出時間離散的平均場方法。

另外，可以發(fā)現(xiàn)影響傳播過程的網(wǎng)絡(luò)拓撲因素有很多：例如度相關(guān)性、動力學(xué)相關(guān)性、局域效應(yīng)(如聚類和模體等)、有限尺度效應(yīng)、網(wǎng)絡(luò)異質(zhì)性、傳播率的異質(zhì)性、感染的非對稱性、拓撲隨時間的演化、鄰居節(jié)點度轉(zhuǎn)化等等；另外網(wǎng)絡(luò)上個體的自適應(yīng)行為、人類動力學(xué)與網(wǎng)絡(luò)拓撲的共演化[162]、網(wǎng)絡(luò)多層耦合性與連邊的多關(guān)系特性也是極其重要、近年研究較熱的問題。那么是否能夠有一種考慮所有重要因素的綜合方法是一個有趣的問題。想要達成這樣一個遠大的目標(biāo)，有很多問題依然擺在我們面前：在傳播動力學(xué)中如此繁復(fù)的影響因素里，它們各自對動力學(xué)產(chǎn)生的影響究竟有多大？哪些對于傳播動力學(xué)來說是至關(guān)重要不可或缺、哪些又是相對無足輕重的？甚至更深一步，基于網(wǎng)絡(luò)的流行病傳播動力學(xué)解析方法其根本局限在哪里？這些都是我們未來最終須要回答的問題。

我們運用和拓展本文介紹的7種常用理論方法來描述網(wǎng)絡(luò)傳播動力學(xué)，取得了一系列研究成果。對于異質(zhì)平均場，我們將其拓展后用于描述信息-疾病非對稱耦合傳播動力學(xué)，發(fā)現(xiàn)信息擴散有利于抑制流行病傳播[165-166]。楊慧等[167]拓展點對近似方法來分析個體異質(zhì)性對自適應(yīng)傳播所帶來的影響。Xu等[164]利用主方程方法準(zhǔn)確地刻畫了多關(guān)系社會網(wǎng)絡(luò)中的信息傳播。王偉等[168-170]將邊劃分理論拓展用于描述非馬爾可夫社會傳播動力學(xué)，發(fā)現(xiàn)網(wǎng)絡(luò)結(jié)構(gòu)和動力學(xué)參量都對相變類型有很大的影響。舒盼盼等[160]利用可變性方法對常用理論方法得到的閾值有效性進行了分析。

我們?nèi)ピu判一種方法的優(yōu)劣與否，不只要看它在某些特定情景下的解析是否準(zhǔn)確，也須要考慮這一方法是否具有較好的可用性與擴展性，能夠加入對于更多重要因素的考慮。一個理論方法所基于的假設(shè)越少，其適用性往往越廣；同時，還要看到，理論解析不同于數(shù)值模擬，解析一定是建立在一定的假設(shè)上的，它對于數(shù)據(jù)的需求性不應(yīng)當(dāng)過強，最好只需要基于較少的初始數(shù)據(jù)(如度分布、度關(guān)聯(lián)性)，其更多的應(yīng)當(dāng)強調(diào)理論解析性，應(yīng)較少的加入過多的數(shù)值模擬；有時從某種意義上講，數(shù)值求解未必會比蒙特卡羅模擬更具優(yōu)勢?？山忉屝砸约耙锥裕抢碚摻馕龇椒ㄋ鶓?yīng)具有的；如果一種解析方法自身就很復(fù)雜、很難為我們提供一些直觀的定性理解，那么即使它十分準(zhǔn)確也未必能稱作好的方法。理論解析事關(guān)對于機制的理解，哪些假設(shè)和近似是合理的關(guān)系到是否能夠描述出動力學(xué)過程中的關(guān)鍵因素，這是非常本質(zhì)而且極其重要的方面。

參考文獻:

[1]AndersonRM,MayRM.InfectiousDiseasesofHumans:dynamicsandcontrol[M].Oxford:OxfordUniversityPress, 1992.

[2]DaileyDJ,GaniJ.EpidemicModeling:AnIntroduction[M].Cambridge:CambridgeUniversityPress, 2001.

[3]BernoulliD.Essaid'unenouvelleanalysedelamortalitecauseeparlapetiteveroleetdesavantagesdel'inoculationpourlaprevenir[J].MemMathPhysAcadRoySci, 1760, 1-45.

[4]KermackWO,McKendrickAG.Acontributiontothemathematicaltheoryofepidemics[J].ProcRSocLondA, 1927, 115, 772: 700-721.

[5]AlbertR,BarabsiAL.Statisticalmechanicsofcomplexnetworks[J].RevModPhys, 2002, 74:47.

[6]BoccalettiS,LatoraV,MorenoY,etal.Complexnetworks:structureanddynamics[J].PhysRep, 2006,424: 175.

[7]NewmanMEJ.Networks:AnIntroduction[M].Oxford:OxfordUniversityPress, 2010.

[8]ButtsCT.Revisitingthefoundationsofnetworkanalysis[J].Science, 2009, 325: 414.

[9]JacksonM.SocialandEconomicNetworks[M].Princeton:PrincetonUniversityPress, 2010.

[10] Vespignani A. Predicting the behavior of techno-Social systems[J]. Science, 2009, 325 (5939): 425.

[11] Vespignani A. Modeling dynamical processes in complex socio-technical systems[J]. Nature Physics, 2012, 8: 32-39.

[12] Dorogovtsev S N,Goltsev A V. Critical phenomena in complex networks[J]. Rev Mod Phys, 2008, 80: 1275.

[13] Barrat A, Barthelmy M, Vespignani A.Dynamical Processes on Complex Networks[M]. New York: Cambridge University Press, 2008.

[14] Pastor-Satorras R, VespignaniA. Epidemic spreading in scale-free networks[J]. Phys Rev Lett, 2001, 86: 32001.

[15] Pastor-Satorras R, Vespignani A. Epidemic dynamics and endemic states in complex networks[J]. Phys Rev E, 2001, 63(6): 066117.

[16] Hufnagel L, Brockmann D, Geisel T. Forecast and control of epidemics in aglobalized world[J]. Proc Natl Acad Sci, 2004, 101: 15124.

[17] Colizza V, Barrat A, Barthélemy M, et al. The role of the airline transportation network in the prediction and predictability of global epidemics[J]. Proc Natl Acad Sci, 2006, 103: 2015.

[18] Balcan D, Hu H, Goncalves B,et al. Seasonal transmission potential and activity peaks of the new influenza A(H1N1): a Monte Carlo likelihood analysis based on human mobility[J]. BMC Medicine, 2009, 7: 45.

[19] Brockmann D, Helbing D.The hidden geometry of complex, network-driven contagion phenomena[J]. Science, 2013, 342 (6164): 1337-1342.

[20] Henkel M, Hinrichsen H, Lubeck S. Non-equilibrium phase transition: Absorbing Phase Transitions[M]. Netherlands: Springer Verlag, 2008.

[21] Danon L,Ford A P,House T,et al. Networks and the epidemiology of infectious disease[J]. Interdisciplinary Perspectives on Infectious Diseases, 2011:284909.

[22] Keeling M, Rohani P. Modeling Infectious Diseases in Humans and Animals[M]. Princeton: Princeton University Press, 2010.

[23] Moreno Y, Pastor-Satorras R, Vespignani A. Epidemic outbreaks in complex heterogeneous networks[J]. Eur Phys J B, 2002, 26: 521-529.

[24] Boguna M, Pastor-Satorras R, Vespignani A.Absence of epidemic threshold in scale-free networks with degree correlation[J]. Phys Rev Lett, 2003, 90: 2.

[25] Castellano C, Pastor-Satorras R. Thresholds for epidemic spreading in networks[J]. Phys Rev Lett, 2010, 105(21): 218701.

[26] Eguiluz V M, Klemm K. Epidemic threshold in structured scale-free networks[J]. Phys Rev Lett, 2002, 89: 10.

[27] Parshani R, Carmi S, Havlin S. Epidemic threshold for the susceptible-infectious-susceptible model on random networks[J]. Phys Rev Lett, 2010, 104(25): 258701.

[28] Gleeson J P, Melnik S, Ward J A,et al. Accuracy of mean-field theory for dynamics on real-world networks[J]. Phys Rev E, 2012, 85(2): 026106.

[29] Guerra B, Gómez-Gardees J. Annealed and mean-field formulations of disease dynamics on static and adaptive networks[J]. Phys Rev E, 2010, 82(3): 035101.

[30] Sood V,Redner S. Voter model on heterogeneous graphs[J]. Phys Rev Lett, 2005, 94(17): 178701.

[31] Pikovsky A S, Rosenblum M G, Kurths J. Synchronization in a population of globally coupled chaotic oscillators[J]. EPL, 1996, 34: 165.

[32] 于淥, 郝柏林. 相變和臨界現(xiàn)象[M]. 北京: 科學(xué)出版社, 1984.

[33] 林宗涵. 熱力學(xué)與統(tǒng)計物理學(xué)[M]. 北京: 北京大學(xué)出版社, 2007.

[34] Landau L D, Lifshitz. Statistical Physics[M]. 3rd Edition London: Pergamon Press, 1986.

[35] Stanley H E. Introduction to Phase transition and Critical Phenomena[M]. New York: Oxford Univ Press, 1971.

[36] 蘇汝鏗. 統(tǒng)計物理學(xué)[M].第2版,上海:高等教育出版社, 2004.

[37] Bailey N T J. The Mathematical Theory of Infectious Diseases and Its Applications[M]. New York: Hafner Press, 1975.

[38] Gomez S, Gómez-Gardenes J, Moreno Y,et al. Nonperturbative heterogeneous mean-field approach to epidemic spreading in complex networks[J]. Phys Rev E, 2011, 84(3): 036105.

[39] Durrett R. Some features of the spread of epidemics and information on a random graph[J]. Proc Natl Acad Sci, 2010, 107: 4491.

[40] Baronchelli A, Loreto V. Ring structures and mean first passage time in networks[J]. Phys Rev E, 2006, 73(2): 026103.

[41] Baronchelli A, Pastor-Satorras R. Mean-field diffusive dynamics on weighted networks[J]. Phys Rev E, 2010, 82(1): 011111.

[42] Gleeson J P. High-accuracy approximation of binary-state dynamics on networks[J]. Phys Rev Lett, 2011, 107(6): 068701.

[43] Milo R, Shen-Orr S,Itzkovitz S,et al. Network motifs: simple building blocks of complex networks[J]. Science, 2002, 298: 824.

[44] Zanette D H, Negro R, Argentina. Dynamics of rumor propagation on small-world networks[J]. Phys Rev E, 2001, 65(4): 041908.

[45] Baronchelli A, Castellano C, Pastor-Satorras R. Voter models on weighted networks[J]. Phys Rev E, 2011, 83(6): 066117.

[46] Liu Z, Wang X, Wang M. Inhomogeneity of epidemic spreading[J]. Chaos, 2010, 20(2): 023128.

[48] Yang Z, Zhou T. Epidemic spreading in weighted networks: an edge-based mean-field solution[J]. Phys Rev E, 2012, 85(5), 056106.

[49] Watts D J, Strogatz S H. Collective dynamics of ‘small-world’ networks[J]. Nature, 1998, 393: 4.

[50] Barabási A L, Albert R. Emergence of scaling in random networks[J]. Science, 1999, 286: 509-512.

[51] Newman M E J. Power laws, Pareto distributions and Zipf's law[J]. Contemporary Physics, 2005, 46: 5.

[52] Castellano C, Pastor-Satorras R. Competing activation mechanisms in epidemics on networks[J]. Sci Rep, 2012, 2: 371.

[53] Yang R, Wang B H, Ren J,et al. Epidemic spreading on heterogeneous networks with identical infectivity[J]. Phys Lett A, 2007, 364: 189-193.

[54] Alexandrowicz Z. Critically branched chains and percolation clusters[J]. Phys Lett A, 1980, 80: 4.

[55] Kenah E, Robins J M. Second look at the epidemic spread[J]. Phys Rev E, 2007, 76(3): 036113.

[56] Karp R,Schindelhauer C,Shenker S,et al. Randomized rumor spreading[Z]. Foundations of Computer Science. Proceedings of 41st Annual Symposium, 2000:565-574.

[57] 許田, 張培培, 姜玉梅, 等. 流行病傳播模型與SARS[J]. 自然雜志, 2004, 26(1): 20-25.

Xu Tian, Zhang Peipei, Jiang Yumei, et al. Epidemic spreading models and SARS[J]. Ziran Zazhi, 2004, 26(1): 20-25.

[58] Saumell-Mendiola A, Serrano M A, Boguna M. Epidemic spreading on interconnected networks[J]. Phys Rev E, 2012, 86(2): 026106.

[59] Rand D A. Correlation equations and pair approximations for spatial ecologies[J]. Advanced Ecological Theory: Principles and Applications, 2009, 12 (34):329-368.

[60] Keeling M J, Rand D A. Spatial correlations and local fluctuations in host parasite systems[J]. From Finite to Infinite Dimensional Dynamical Systems, 2001: 5-57.

[61] Keeling M J.The effects of local spatial structure on epidemiological invasions[J].Proc R SocLond B, 1999, 266 (1421): 859-867.

[62] Morris A J. Representing spatial interactions in simpleecological models[D]. UK: University of Warwick, 1997.

[63] Benoit J, Nunes A, Teloda G M. Pair approximation models for disease spread[J]. Eur Phys J B, 2006, 50 (1/2): 177-181.

[64] Eames K T D, Keeling M J. Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases[J]. Proc Natl Acad Sci, 2002, 99: 13330-13335.

[65] Pugliese E,Castellano C. Heterogeneous pair approximation for voter models on networks[J].EPL, 2009, 88 (5): 58004.

[66] Shaw L B,SchwartzI B. Fluctuating epidemics on adaptive networks[J]. Phys Rev E, 2008, 77 (6): 066101.

[67] Gross T, D’Lima C J D, Blasius B. Epidemic dynamics on an adaptive network[J].Phys Rev Lett, 2006, 96 (20):208701 .

[68] Shaw L B, Schwartz I B. Enhanced vaccine control of epidemics in adaptive networks[J].Phys Rev E,2010, 81 (4): 046120.

[69] Keeling M J, Rand D A, Morris A J. Correlation models for childhood epidemics[J].Proc R Soc B, 1997, 264 (1385): 1149-1156.

[70] CaldarelliG, Pastor-Satorras R, Vespignani A. Structure of cycles and local ordering in complex networks[J].Eur Phys J B,2004, 38 (2): 183-186.

[71] GleissP M, Stadler P F, Wagner A, et al.Relevant cycles in chemical reaction networks[J]. Adv Complex Syst, 2001, 4 (5): 207-226.

[72] Sato K, Matsuda H, Sasaki A. Pathogen invasion and host extinction in lattice structured populations[J]. J Math Biol, 1994, 32 (3): 251-268.

[73] Keeling M J. The ecology and evolution of spatial host-parasite systems[D]. UK: University of Warwick, 1995.

[74] Nakamaru M, Matsuda H, IwasaY. The evolution of cooperation in a lattice-structured population[J].J Theor Biol, 1997, 184 (1): 65-81.

[75] Szabó G,Hauert C. Evolutionary prisoner’s dilemma games with voluntary participation[J].Phys Rev E,2002, 66 (6): 062903.

[76] Perc M,Marhl M. Evolutionary and dynamical coherence resonances in the pair approximated prisoner'sdilemma game[J]. New J Phys, 2006, 8: 142.

[77] Marceau V, No?l P A, Hébert-Dufresne L,et al. Adaptive networks: coevolution of disease and topology [J]. Phys Rev E, 2010, 82(3): 036116.

[79] Durrett R, Gleeson J P, Lloyd A L,et al. Graph fission in an evolving voter model[J]. Proc Natl Acad Sci, 2012, 109: 3682-3687.

[80] Mollison D. Spatial contact models for ecological and epidemic spread[J]. J Roy Stat Soc B, 1977, 39: 283-326.

[81] Grassberger P. On the critical behavior of the general epidemic process and dynamical percolation[J]. Math Biosci, 1982, 63: 157-172.

[82] Newman M E J. Spread of epidemic disease on networks[J]. Phys Rev E, 2002, 66(1): 016128.

[83] Newman M E J. Threshold effects for two pathogens spreading on a network[J]. Phys Rev Lett, 2005： 95(10):108701.

[84] Meyers L A, Newman M E J, Martin M,et al. Applying network theory to epidemics: Control measures for mycoplasma pneumoniae outbreaks[J]. EID, 2003, 92: 204.

[85] Madar N, Kalisky T, Cohen R, et al. Immunization and epidemic dynamics in complex networks[J]. Eur Phys J B, 2004, 38: 269.

[86] No?l P A, Davoudi B, Brunham R C, et al. Time evolution of epidemic disease on finite and infinite networks[J]. Phys Rev E,2009, 79(2): 026101.

[87] Marder M. Dynamics of epidemics on random networks[J]. Phys Rev E, 2007, 75: 066103.

[88] Miller J C. Epidemic size and probability in populations with heterogeneous infectivity and susceptibility[J]. Phys Rev E, 2007, 76(1): 010101.

[90] Newman M E J. Random graphs with clustering[J]. Phys Rev Lett, 2009, 103(5): 058701.

[91] Newman M E J. Competing epidemics on complex networks[J]. Phys Rev E, 2011, 84(3): 036106.

[92] Allard A, No?l P A, DubéL J. Heterogeneous bond percolation on multitype networks with an application to epidemic dynamics[J]. Phys Rev E, 2009, 79(3): 036113.

[93] Funk S, Jansen V A. Interacting epidemics on overlay networks[J]. Phys Rev E, 2010, 81(3): 036118.

[94] Dickison M, Havlin S, Stanley H E. Epidemics on interconnected networks[J]. Phys Rev E, 2012, 85(6), 066109.

[95] Wang Y, Xiao G. Epidemics spreading in interconnected complex networks[J]. Phys Lett A, 2012, 376: 42-43.

[96] Pearl J. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[M]. San Francisco:Morgan Kaufmann,1988.

[97] Mézard M, Parisi M, Virasoro M. Spin Glass Theory and Beyond World Scientific[M]. Singapore:Pergamon Press, 1987.

[98] Mézard M, Parisi G. The Bethe lattice spin glass revisited[J], Eur Phys J B,2001, 20(2):217-233.

[99] Karrer B, Newman M E J. Message passing approach for general epidemic models[J]. Phys Rev E, 2010,82(1):016101.

[100]Bianconi G, Gulbahce N. Algorithm for counting large directed loops[J]. J Phys A Math Theor, 2008,41(22):224008.

[101]Sulc P, Zdeborov L. Belief propagation for graph partitioning[J]. J Phys A MathTheor, 2010,43(28):285003.

[102]Shiraki Y, Kabashima Y. Cavity analysis on the robustness of random networks against targeted attacks: influences of degree-degree correlations[J]. Phys Rev E, 2010, 82(3): 036101.

[103]Reichardt J, Alamino J, Saad D. The interplay of microscopic and mesoscopic structure in complex networks[DB/OL].[2014-06-30]. arXiv.org/abs/1012.4524.

[104]Yoon S, Goltsev A V, Dorogovtsev S N, et al. Belief-propagation algorithm and the ising model on networks with arbitrary distributions of motifs[J]. Phys Rev E,2011,84(4):041144.

[105]Pastor-Satorras R, Vespignani A. Epidemic dynamics in ?nite size scale-free networks[J]. Phys Rev E,2002,65(3): 035108.

[106]Shao J, Buldyrev S V, Braunstein L A, et al. Structure of shells in complex networks[J]. Phys Rev E, 2009, 80(3):036105.

[107]張昊,陳超,王長春. 基于空穴理論的復(fù)雜網(wǎng)絡(luò)流行病傳播控制[J]. 電子科技大學(xué)學(xué)報,2011, 40:4.

Zhang Hao, Chen Chao, Wang Changchun. Epidemic spread and control on complex networks using cavity theory[J]. JUESTC, 2011, 40:4.

[108]Altarelli F, Braunstein A, Dall'Asta L, et al. The spread optimization problem[DB/OL].[2014-06-30]. arXiv/org/abs/1203.1426v1.

[109]Volz E. SIR dynamics in random networks with heterogeneous connectivity[J]. Math Biol, 2008, 56:293-310.

[110]Miller J C. A note on the derivation of epidemic final sizes[J]. Bull Math Biol, 2012, 74:2125-2141.

[111]Miller J C. A note on a paper by Erik Volz: SIR dynamics in random networks[J]. Math Biol, 2011, 62:349-358.

[112]Miller J C, Slim A C, Volz M E. Edge-based compartmental modeling for infectious disease spread part I: an overview[J]. Journal of the Royal Society Interface,2011,9(20):890-906.

[113]Miller J C, Volz M E. Edge-based compartmental modeling for epidemic spread Part II: model selection and hierarchies[J]. Journal of Mathematical Biology,2013,67(4):869-899.

[114]Miller J C, Volz M E. Edge-based compartmental modeling for infectious disease spread Part III: disease and population structure[J]. Eprint Arxiv,2011,9(70):890-906.

[115]Volz E, Frost S D W, Rothenberg R, et al. Epidemiological bridging by injection druguse drives an early HIV epidemic[J]. Epidemics,2010,2(3):155-164.

[116]Volz E, Meyers L A. Epidemic thresholds in dynamic contact networks[J]. J R Soc Interface, 2009, 6: 233-241.

[117]Volz E M, Miller J C, Galvani A, et al. Effects of heterogeneous and clustered contact patterns on infectious disease dynamics[J]. PLoS Computational Biology, 2011, 7(6): e1002042.

[118]Miller J C. Epidemics on networks with large initial conditions or changing structure[J]. Plos One,2014,9(7):e101421.

[119]Valdez L D, Macri P A, Braunstein L A. Intermittent social distancing strategy for epidemic control[J].Phys Rev E, 2012, 85(3): 036108.

[120]Valdez L D, Macri P A, Braunstein L A. Temporal percolation of the susceptible network in an epidemic spreading[J].PLoS ONE,2012,7(9): e44188.

[121]Valdez L D, Macri P A, Braunstein L A. Temporal percolation of a susceptible adaptive network[J]. Physica A, 2013, 392(18), 4172-4180.

[123]Wang Y, Chakrabarti D, Wang C, et al. Epidemic spreading in real networks: an eigenvalue viewpoint[C]//Proceedings of the 22nd Symposium on Reliable Distributed Systems, Florece, Italy,2003: 25-34.

[124]Chung F, Lu L, Vu V. Spectra of random graphs with given expected degrees[J]. Proc Natl Acad Sci. 2003, 100 (11): 6313-6318.

[126]Kephart J O, White S R. Directed-graph epidemiological models of computer viruses[C]//Proceedings of the 1991 IEEE Computer Science Symposium on Research in Security and Privacy, New York, USA, 1991: 343-359.

[127]Chakrabarti D, Leskovec J, Faloutsos C, et al. Information survival threshold in sensor and P2P networks[C]//26th IEEE International Conference on Computer Communications, Anchorage, USA, 2007:1316-1324.

[128]Restrepo J G, Ott E, Hunt B R. Emergence of synchronization in complex networks of interacting dynamical systems[J], Physica D, 2006, 224 (1/2): 114-122.

[129]Feng J, Jirsa V K, Ding M. Synchronization in networks with random interactions: theory and applications[J], Chaos, 2006, 16 (1): 015109.

[130]Cohen R, Erez K, Ben-Avraham D, et al. Resilience of the internet to random breakdowns[J]. Phys Rev Lett, 2000, 85 (21): 4626-4628

[131]Vazquez A, Moreno Y. Resilience to damage of graphs with degree correlations[J]. Phys Rev E, 2003, 67 (1): 015101(R).

[133]Molloy M, Reed B. The size of the giant component of a random graph with a given degree sequence[J]. Combinatorics, Prob Comput, 1998, 7(3):295.

[134]Gómez S, Borge-Holthoefer J, Meloni S,et al. Discrete-time Markov chain approach to contact-based disease spreading in complex networks [J]. EPL, 2010, 89 (3): 38009.

[135]Gibson G J. Markov chain Monte Carlo methods for fitting spatiotemporal stochastic in plant epidemiology[J]. Appl Stat, 1997, 46(2):215-233.

[136]Ganesh A, Massoulié L, Towsley D. The effect of network topology on the spread of epidemics[C]∥ Proceeding of IEEE INFOCOM 2005, 2005: 1455-1466.

[137]O’Neill P D. A tutorial introduction to Bayesian inference for stochastic epidemic models using Markov chain Monte Carlo methods [J]. Math Biosci, 2002, 180(1/2):103-114.

[138]Dorogovtsev S N, Mendes J F F. Evolution of networks[J]. Adv Phys,2002, 51, 1079.

[139]Chakrabarti D,Wang Y, Wang C, et al. Epidemic thresholds in real network[J]. ACM Trans Inf Syst Secur,2008, 10: 1.

[140]Mieghem P V, Omic J, Kooij R. Virus spread in network[C]. IEEE ACMTrans Netw,2009, 17： 1.

[141]Ferreira S C,Castellano C,Pastor-Satorras R. Epidemic thresholds of the susceptible-infected-susceptible model on networks:a comparison of numerical and theoretical results[J]. Phys Rev E,2012, 86(4): 041125.

[142]Cator E, Mieghem P V. Second-order mean-field susceptible-infected-susceptible epidemic threshold[J]. Phys Rev E,2012, 85(5): 056111.

[143]Goltsev A V,Dorogovtsev S N,Oliveira J G, et al. Localization and spreading of diseases in complex networks[J]. Phys Rev Lett,2012, 109(12): 128702.

[144]Lee H K,Shim P S,Noh J D. Epidemic threshold of the susceptible-infected-susceptible model on complex networks[J]. Phys Rev E,2013, 87(6): 062812.

[145]Bogun M, Castellano C, Pastor-Satorras R. Nature of the epidemic threshold for the susceptible-infected-susceptible dynamicsin networks[J]. Phys Rev Lett,2013, 111(6): 068701.

[146]Givan O,Schwartz N,Cygelberg A, et al. Predicting epidemic threshold on complex networks: Limitations of mean-field approaches[J]. J Theor Biol,2011, 288, 21.

[147]Moreno Y,Pastor-Satorras R,Vespignani A. Epidemics outbreaks in complex heterogeneous networks[J]. Eur Phys J B, 2002, 26: 521-529.

[148]Callaway D S, Newman M E J,Strogatz S H,et al.Network robustness and fragility: percolation on random graphs[J]. Phys Rev Lett, 2000, 85, 5468-5471.

[149]Youssef M,Scoglio C. An individual-based approach to SIR epidemics in contact networks[J]. J Theor Biol, 2011, 283.136-144.

[150]Prakash B A, Chakrabarti D, Faloutsos M, et al. Get the flu (or mumps) Check the eigenvalue![DB/OL].[2014-06-30].arxiv.org/abs/10040060.

[151]Zanette D H. Critical behavior of propagation on small-world networks[J].Phys Rev E, 2001, 64(5)： 050901(R).

[152]Ben-Naim E,Krapivsky P L. Size of outbreaks near the epidemic threshold[J]. Phys Rev E, 2004, 69(5): 050901(R).

[153]Ben-Naim E,Krapivsky P L. Scaling behavior of the threshold epidemics[J]. Eur Phys J B, 2012, 85: 1.

[154]Kessler D A, Shnerb N M. Solution of an infection model near threshold[J]. Phys Rev E, 2007, 76(1): 010901(R).

[155]Khalleque A, Sen P. The susceptible-infected-recovered model on a Euclidean network[J]. J Phys A: Math Theor,2013, 46(9): 095007.

[156]Holme P. Extinction times of epidemic outbreaks in networks[J]. PloS ONE,2013, 8: e84429.

[157]Hong H,Ha M,Park H. Finite-size scaling in complex networks[J].Phys Rev Lett, 2007, 98(25): 258701.

[158]Binder K, Heermann D W. Monte carlo simulation in statistical physics,[M]. 5th ed Berlin: Springer-Verlag, 2010.

[159]Ferreira S C,Ferreira R S, Castellano C, et al. Quasistationary simulations of the contact process on quenched networks[J]. Phys Rev E, 2011, 84(6): 066102.

[160]Shu P, Wang W, Tang M, et al, Numerical identification of epidemic thresholds for susceptible-infected-recovered model on finite-size networks [J]. Chaos, 2015, 25 (6), 063104.

[161]Wang W, Liu Q H, Zhong L F, et al, Predicting the epidemic threshold of the susceptible-infected-recovered model [DB/OL]. [2014-06-30]. www.nature.com/articles/srep24676.

[162]Li R, Wang W, Di Z. Effects of human dynamics on epidemic spreading in Cote d'Ivoire[DB/OL]. [2014-06-30]. arxiv.org/abs/160500899v1.

[163]Li R, Tang M, Hui P M. Epidemic spreading on multi-relational networks[J]. Acta Phys Sin, 2013, 62(16): 168903.

[164]Xu E H W, Wang W, Xu C, et al, Suppressed epidemics in multirelational networks [J]. Phys Rev E， 2015, 92 (2): 022812.

[165]Wang W, Tang M, Yang H, et al, Asymmetrically interacting spreading dynamics on complex layered networks [J]. Sci Rep, 2014, 4: 5097.

[166]Liu Q H, Wang W, Tang M, et al, Impacts of complex behavioral responses on asymmetric interacting spreading dynamics in multiplex networks [J]. Food & Nutrition Sciences, 2015,6(6):555-561.

[167]H Yang, M Tang, T Gross, Large epidemic thresholds emerge in heterogeneous networks of heterogeneous nodes [J]. Sci Rep. 2015, 5, 13122.

[168]Wang W, Tang M, Zhang H F, et al, Dynamics of social contagions with memory of nonredundant information [J]. Phys Rev E, 2015, 92(1): 012820.

[169]Wang W, Tang M, Shu P, et al, Dynamics of social contagions with heterogeneous adoption thresholds: crossover phenomena in phase transition [J]. New J Phys, 2016, 18(1): 013029.

[170]Wang W, Shu P, Zhu Y X, et al, Dynamics of social contagions with limited contact capacity [J], Chaos, 2015, 25(10): 103102.

(責(zé)任編輯耿金花)

Review of Threshold Theoretical Analysis About Epidemic Spreading Dynamics on Complex Networks

LI Ruiqi1,2, WANG Wei1, SHU Panpan1,YANG Hui1,PAN Liming1, CUI Aixiang1, TANG Ming1

(1.Web Science Center, University of Electronic Science and Technology of China, Chengdu 611731,China;2.School of Systems Science, Beijing Normal University, Beijing 100875, China)

Abstract:In this Review, we mainly focus on solving the threshold theoretically, and introduce seven common methods, including Mean-Field theory, Pair-wise Approximation, Master equation, Generating Function, Edge Percolation, Cavity method, Edge Classification and Spectral analysis. We also summarize the difference of thresholds between the SIS and SIR model. We are aiming to provide a clear picture for beginners and a good reference for researchers.

Key words:complex networks; epidemic spreading;outbreak threshold; theoretical analysis

文章編號：16723813(2016)01000139;

DOI:10.13306/j.1672-3813.2016.01.001

收稿日期：2014-12-24

基金項目：國家自然科學(xué)基金(11105025,11575041)

作者簡介：李睿琪(1992-)，男，河北張家口人，博士研究生，主要研究方向為復(fù)雜網(wǎng)絡(luò)傳播動力學(xué)、城市演化動力學(xué)、交通分析、人類動力學(xué)、風(fēng)險投資網(wǎng)絡(luò)分析、大數(shù)據(jù)分析與建模。

中圖分類號：N93;N94

文獻標(biāo)識碼：A