復(fù)雜網(wǎng)絡(luò)上非馬爾科夫易感-感染模型的二階平均場求解

祁 婷,林詔華,馮 秘,唐 明

（1. 華東師范大學(xué) 通信與電子工程學(xué)院, 上海 200241; 2. 華東師范大學(xué) 物理與電子科學(xué)學(xué)院,上海 200241; 3. 香港浸會大學(xué) 物理系, 香港 999077）

0 引 言

疾病傳播是當今復(fù)雜網(wǎng)絡(luò)科學(xué)非常熱門的研究課題之一, 傳統(tǒng)的網(wǎng)絡(luò)動力學(xué)模型是具有無記憶特性的馬爾科夫過程, 即對于網(wǎng)絡(luò)中所有的個體, 疾病傳播過程和恢復(fù)過程都可以看成是典型的泊松過程[1-3]. 在馬爾科夫過程中, 兩個連續(xù)的感染過程或者恢復(fù)過程之間的時間間隔遵循指數(shù)分布, 換言之, 對當前的等待時間并不會影響未來的等待時間. 通過對傳統(tǒng)的SIS(Susceptible-Infected-Susceptible)模型和SIR(Susceptible-Infected-Recovered)模型的平均場方法的研究, 可以發(fā)現(xiàn)馬爾科夫假設(shè)極大地促進了復(fù)雜網(wǎng)絡(luò)上動力學(xué)數(shù)學(xué)理論的發(fā)展[2-3].

然而, 越來越多的實證研究與經(jīng)驗分析表明, 現(xiàn)實生活中的人類活動往往呈現(xiàn)出具有“胖尾”時間分布的非馬爾科夫特性[4-6]. 因此, 經(jīng)典的馬爾科夫框架僅僅是對真實世界的一個近似刻畫.

網(wǎng)絡(luò)中的非馬爾科夫傳播動力學(xué)近幾年來引起了人們廣泛的研究熱情[7-15]. 馬爾科夫理論無法精確地描述人與人之間的交互活動在很早以前就已經(jīng)被注意到了, 同時也發(fā)現(xiàn)感染過程中的異質(zhì)性會阻礙網(wǎng)絡(luò)中疾病的擴散[8]. 消息傳播理論求解了具有任意感染時間分布和恢復(fù)時間分布的SIR疾病傳播模型[9]. 通過研究具有固定恢復(fù)時間分布以及“胖尾”感染時間分布的SIR模型, 發(fā)現(xiàn)了感染過程中的時序異質(zhì)性能夠顯著抑制疾病的傳播[10]. 一個比較重要的研究是, 通過重新定義有效感染率, 能夠發(fā)現(xiàn)網(wǎng)絡(luò)中傳播動力學(xué)的馬爾科夫過程與非馬爾科夫過程在穩(wěn)態(tài)時存在著等價性[11-13]. 兩種經(jīng)典的SIS模型的連邊激活機制被提出后, 該等價性存在的條件被證明為激活邊的產(chǎn)生與激活邊指向的易感節(jié)點之間不存在時序相關(guān)性[13]. 另外, 在非馬爾科夫傳播動力學(xué)方面, SIR模型的點對近似理論以及SIS模型的淬火平均場理論最近也被提了出來[16-18].

求解網(wǎng)絡(luò)傳播動力學(xué)比較常用的方法是平均場理論. 早期的平均場方法假設(shè)了網(wǎng)絡(luò)中的所有節(jié)點都能夠在統(tǒng)計意義上等價[19]. 另外, 為了考慮網(wǎng)絡(luò)度分布的異質(zhì)特性, 異質(zhì)平均場理論將網(wǎng)絡(luò)中所有度相同的節(jié)點都視為等價[1]. 而淬火平均場方法則將網(wǎng)絡(luò)中每個節(jié)點的連接情況都考慮到理論框架中[20-21]. 然而, 以上的理論方法都忽略了傳播過程中節(jié)點與節(jié)點之間的動力學(xué)相關(guān)性, 為了克服這個困難, 點對近似理論深入地探究了兩個鄰接節(jié)點之間的相關(guān)特性[22-25]. 另外, 主方程理論能夠得到更為精確的結(jié)果, 并且在一定條件下能夠退化為平均場理論或者點對近似理論[26].

本文在前人工作的基礎(chǔ)上, 提出了關(guān)于非馬爾科夫 SI 模型的二階平均場方法; 根據(jù)SI模型中易感態(tài)節(jié)點能轉(zhuǎn)換為感染態(tài)節(jié)點而感染態(tài)節(jié)點無法轉(zhuǎn)變?yōu)橐赘袘B(tài)節(jié)點的特點, 提出了閑置邊的概念, 該方法能夠求解SI模型的非馬爾科夫爆發(fā)過程, 同時也能夠預(yù)測每個節(jié)點被感染的平均時刻. 本文分別在真實網(wǎng)絡(luò)和人工網(wǎng)絡(luò)上進行了實驗?zāi)M, 并驗證了以上研究成果的正確性和準確性.

1 非馬爾科夫SI模型二階平均場理論

1.1 非馬爾科夫SI模型

在復(fù)雜網(wǎng)絡(luò)的SI傳播模型中, 節(jié)點存在著兩種狀態(tài)—易感態(tài)(S態(tài))與感染態(tài)(I態(tài)), 其中I態(tài)節(jié)點有一定的概率將與其相連的S態(tài)節(jié)點感染成I態(tài), I態(tài)節(jié)點卻無法恢復(fù)為S態(tài)節(jié)點. 為了研究方便, 我們將無向網(wǎng)絡(luò)中的無向邊看作是兩條方向相反的有向邊, 對于有向邊而言, 疾病的傳播具有方向性, 只能是箭尾處的I態(tài)節(jié)點才能感染箭頭處的S態(tài)節(jié)點. 在非馬爾科夫SI模型中, 1條有向邊存在著4種狀態(tài)—〈SS〉態(tài)、〈SI〉態(tài)、〈IS〉態(tài)和〈II〉態(tài), 分別表示該有向邊的箭頭節(jié)點為S態(tài)并且箭尾節(jié)點為S態(tài)、箭頭節(jié)點為S態(tài)并且箭尾節(jié)點為I態(tài)、箭頭節(jié)點為I態(tài)并且箭尾節(jié)點為S態(tài)、箭頭節(jié)點為I態(tài)并且箭尾節(jié)點為I態(tài).

如果1條有向邊的箭頭節(jié)點為易感態(tài), 箭尾節(jié)點為感染態(tài)時, 即〈SI〉態(tài), 這條邊就被稱為激活邊,其他狀態(tài)的邊則被稱為非激活邊. 同時將當前時刻與該有向邊成為激活邊的時刻的時間差定義為激活邊的狀態(tài)年齡, 用K表示. 疾病的感染的傳播滿足特定的時間分布—感染時間分布ψinf(κ), 其物理意義在于, 1 條激活邊傳播疾病時的狀態(tài)年齡K處在區(qū)間 [κ,κ+dκ)之內(nèi)概率為ψinf(κ)dκ. 同時, 我們用ωinf(κ)表示激活邊的感染率, 其物理意義為條激活邊的狀態(tài)年齡K如果處在區(qū)間 [κ,κ+dκ) 之內(nèi),則該激活邊傳播疾病將其指向的易感態(tài)節(jié)點感染的概率為ωinf(κ)dκ. 另外,ωinf(κ)dκ還可以理解為是, 激活邊在狀態(tài)年齡的區(qū)間 [0,κ)內(nèi)沒有傳播疾病的前提下, 在狀態(tài)年齡區(qū)間 [κ,κ+dκ) 傳播疾病的概率. 則ψinf(κ) 與ωinf(κ) 的關(guān)系可以表示為

1.2 閑置邊的定義

復(fù)雜網(wǎng)絡(luò)可以由鄰接矩陣A表示, 即

鄰接矩陣A中的元素aij=1表示復(fù)雜網(wǎng)絡(luò)中存在著有向邊i←j(即由節(jié)點j指向節(jié)點i的有向邊),aij=0 表示復(fù)雜網(wǎng)絡(luò)中不存在有向邊i←j. 值得注意的是, 在SI模型中, 某些有向邊在整個爆發(fā)過程中始終不會成為激活邊, 即始終不參與疾病的傳播. 如圖1所示, 紅色節(jié)點為種子節(jié)點, 黃色節(jié)點為非種子節(jié)點, 藍色邊為無向邊, 灰色和綠色的兩條邊是將無向邊看作的兩條有向邊(所有無向邊在理論計算時都應(yīng)看作是兩條方向相反的有向邊, 但圖中只畫出了一條邊被看作兩條有向邊). 對于疾病爆發(fā)而言, 疾病只能從節(jié)點j通過連邊i←j傳向節(jié)點i, 當節(jié)點i被感染后, 由于SI模型中的I態(tài)節(jié)點無法轉(zhuǎn)變?yōu)镾態(tài)節(jié)點, 故連邊j←i將永遠不可能變?yōu)榧せ钸厖⑴c傳播, 我們稱這樣的連邊為閑置邊, 見圖1, 其中灰色的邊即為閑置邊.

圖1 閑置邊示意圖Fig. 1 Schematic diagram of idle edges

但是在理論計算中, 這種閑置邊也會參與計算, 而在實驗?zāi)M中這種邊卻不會參與疾病的傳播,因此閑置邊的存在會使計算誤差增大, 故有必要將連邊j←i對應(yīng)的矩陣元素aji置為0. 我們將閑置邊對應(yīng)的矩陣元素值置為0后得到的矩陣記為, 即

為了得到, 我們需要判斷鄰接矩陣A中每個為1的元素所對應(yīng)的邊是否是閑置邊, 然后將滿足條件的元素置為0.

判定一條邊j←i是否為閑置邊的算法如下.

(1) 如果節(jié)點j為種子節(jié)點, 則j←i一定為閑置邊, 反之, 則進行步驟(2);

(2) 將節(jié)點i標記, 如果存在節(jié)點k(kj)使得連邊i←k存在, 則將節(jié)點k標記;

(3) 得到被標記的節(jié)點集合M, 如果存在節(jié)點l(l∈M)與節(jié)點k(k∈/M)使得連邊l←k存在, 并且連邊l←k不為連邊i←j, 則將節(jié)點k標記;

(4) 不斷重復(fù)步驟(3), 直到不再有新的節(jié)點被標記, 則進行步驟(5);

(5) 如果種子節(jié)點沒有被標記, 則連邊j←i是閑置邊, 反之則不是.

通過以上算法我們可以判斷網(wǎng)絡(luò)中哪些邊為閑置邊, 進而可以得到矩陣.

而此方法的基本思想是, 如果某條邊j←i為閑置邊, 那么說明種子節(jié)點不可能通過連邊j←i將節(jié)點j感染; 同樣, 如果將疾病反向傳播, 節(jié)點i在連邊j←i不存在的條件下也無法將疾病傳染給種子節(jié)點.

1.3 二階平均場方法

對于節(jié)點而言, 我們定義Ii(t)和Si(t)分別為節(jié)點i在時刻t處于感染態(tài)和易感態(tài)的概率. 對于有向邊, 我們定義 〈SS〉i←j(t) 、〈IS〉i←j(t) 、〈II〉i←j(t)分別為連邊i←j在時刻t處于 〈SS〉態(tài)、〈IS〉態(tài)、〈II〉態(tài)的概率, 另外, 定義 〈SI〉i←j(κ;t)為連邊i←j在時刻t處于〈SI〉態(tài)并且其狀態(tài)年齡為κ的概率密度.

對于 〈SI〉 態(tài)的具有任意狀態(tài)年齡的有向邊, 只要它傳播了疾病, 就能夠?qū)⒓^節(jié)點感染為S態(tài), 因此, 能夠得到

被感染的S態(tài)節(jié)點會轉(zhuǎn)化為I態(tài)節(jié)點, 因此可以得到關(guān)于Ii(t) 的表達式

對于〈SS〉態(tài)連邊而言, 其箭頭節(jié)點被感染會轉(zhuǎn)化為〈IS〉態(tài)連邊, 其箭尾節(jié)點被感染會轉(zhuǎn)化為〈SI〉態(tài)連邊, 則 〈SS〉i←j(t) 的表達式可以表示為

其中, 式子右端第一項表示 〈SS〉態(tài)連邊箭頭節(jié)點被感染, 第二項表示箭尾節(jié)點被感染.

對于一個狀態(tài)年齡為κ的〈SI〉態(tài)連邊, 其箭頭節(jié)點的感染方式有兩種: 一種是箭頭節(jié)點被除箭尾節(jié)點以外的其他節(jié)點感染; 另一種是箭頭節(jié)點被箭尾節(jié)點感染. 如果其箭頭節(jié)點在時刻t沒有被感染,則在t+dt, 它的狀態(tài)年齡將變?yōu)棣?dκ, 其中 dκ=dt. 因此可以得到

進而可以得到

由于 〈SS〉態(tài)連邊箭尾節(jié)點被感染會轉(zhuǎn)變?yōu)闋顟B(tài)年齡為0的 〈SI〉 態(tài)連邊, 因此可以得到

對于 〈IS〉態(tài)連邊i←j而言, 其箭尾節(jié)點被感染會轉(zhuǎn)化為 〈SI〉 態(tài)連邊, 感染方式有兩種: 一種是箭尾節(jié)點被除箭頭節(jié)點以外的其他節(jié)點感染; 另一種是箭頭節(jié)點通過連邊j←i將箭尾節(jié)點感染, 則

其中, 等號右端的第一項表示 〈SS〉態(tài)連邊轉(zhuǎn)變?yōu)?〈IS〉態(tài)連邊, 第二項和第三項表示 〈IS〉 態(tài)連邊的箭尾節(jié)點被感染.

由于 〈SI〉態(tài)連邊的箭尾節(jié)點被感染和 〈IS〉態(tài)連邊的箭頭節(jié)點被感染都會導(dǎo)致連邊轉(zhuǎn)化為 〈II〉 態(tài),因此可以得到

其中, 等號右端的第一項表示 〈SI〉態(tài)連邊轉(zhuǎn)變?yōu)?〈II〉態(tài)連邊, 第二項表示 〈IS〉 態(tài)連邊的箭尾節(jié)點被除箭頭節(jié)點以外的其他節(jié)點感染, 第三項表示箭頭節(jié)點通過連邊j←i將箭尾節(jié)點感染.

1.4 平均感染時刻的預(yù)測

Ii(t)還有另外一層物理意義, 即節(jié)點i的被感染時刻小于t的概率. 根據(jù)概率論理論, 因為在 dt時間內(nèi)Ii(t)的增量 dIi(t)即代表了節(jié)點i從S態(tài)變?yōu)镮態(tài)的概率, 則可以表示為節(jié)點i的被感染時刻等于t的概率密度. 設(shè)節(jié)點i被感染的時間為Ti, 則其期望值E(Ti) 為

通過E(Ti)的計算, 就能夠計算出節(jié)點i被感染的平均時刻, 進而就能夠知道在復(fù)雜網(wǎng)絡(luò)中, 對于擁有特定感染時間分布與特定種子節(jié)點的傳播過程而言, 哪些節(jié)點能夠迅速被感染, 哪些節(jié)點被感染則需要很長的時間.

值得注意的是, 公式(12)只適用于那些狀態(tài)只能單向變換的模型, 如SI模型、SIR模型; 但對SIS模型而言, S態(tài)可以轉(zhuǎn)變?yōu)镮態(tài), I態(tài)能夠恢復(fù)為S態(tài), 則公式(12)不再適用, 因為對于SIS模型,S態(tài)變?yōu)镮態(tài)這個過程可以進行無窮次, 而對于SI模型和SIR模型在整個爆發(fā)過程中只會經(jīng)歷一次.

2 模擬驗證

為了驗證理論的正確性, 我們設(shè)置感染時間分布ψinf(κ) 為韋布爾分布, 即

設(shè)置βI=0.5,αI分別取 0.5, 1, 2, 4.αI和βI的下標I表示感染過程, 對于韋布爾分布,αI值越大,該分布的異質(zhì)性越弱, 圖象的分布越集中; 當αI=1 的時候, 該韋布爾分布為經(jīng)典的指數(shù)分布. 對于4個值不同的αI, 其分布如圖2所示, 隨著αI值越來越大, 其圖象的分布也越來越集中.

圖2 韋布爾分布Fig. 2 Weibull distribution

同時我們分別測試了在BA(Barabási-Albert)無標度網(wǎng)絡(luò)、ER(Erd?s-Rényi)隨機網(wǎng)絡(luò)以及真實網(wǎng)絡(luò)中SI疾病爆發(fā)的理論計算與實驗?zāi)M, 詳見圖3. 圖3中, 橫坐標t表示時間, 縱坐標I(t) 表示整個網(wǎng)絡(luò)在時間t時的感染密度. 圖3a) 為BA網(wǎng)絡(luò)上的理論與模擬結(jié)果, 節(jié)點數(shù)N=10 000, 平均度 〈k〉=4 ;圖3b) 為ER網(wǎng)絡(luò)上的理論與模擬結(jié)果, 節(jié)點數(shù)N=10 000, 平均度 〈k〉=4 ; 圖3c) 為真實網(wǎng)絡(luò)上的結(jié)果. 為了后面能夠在單個圖象上全部展示所有節(jié)點被感染的平均時刻, 我們選取了1個節(jié)點數(shù)只有34的真實網(wǎng)絡(luò), 其平均度 〈k〉≈4.59. 該網(wǎng)絡(luò)為著名的Zachary karate club網(wǎng)絡(luò), 且曾多次被相關(guān)領(lǐng)域的研究者所使用, 其數(shù)據(jù)是在1977年由 Zachary[27]所收集, 數(shù)據(jù)中的節(jié)點代表俱樂部里的成員, 無向連邊表示相連的兩個成員之間存在關(guān)系. 對于所有的實驗?zāi)M, 我們都是實驗10 000次進行平均然后得到結(jié)果. 圖3中的空心菱形、空心圓形、空心三角形、空心倒三角形的散點圖分別代表αI取0.5, 1,2, 4的實驗?zāi)M結(jié)果; 實心菱形、實心圓形、實心三角形、實心倒三角形的圖象分別代表αI取0.5, 1, 2,4的理論計算結(jié)果, 可根據(jù)公式(4)–公式(11)得到. 可以發(fā)現(xiàn), 該理論在4個不同的參數(shù)下, 以及3個不同的網(wǎng)絡(luò)中都可以很好地預(yù)測模擬結(jié)果, 其中隨著αI值的減小, 網(wǎng)絡(luò)爆發(fā)的速率逐步增大, 其主要的原因在于, 在爆發(fā)的初始階段, 較小的αI值會導(dǎo)致更大的感染率, 使得在初始階段其爆發(fā)速率增加,更多的被感染節(jié)點進而引發(fā)更快的傳播速率.

圖3 SI模型理論與模擬對比圖Fig. 3 Comparison between the theory and simulation of the SI model

另外, 我們也測試了在真實網(wǎng)絡(luò)中每個節(jié)點被感染時間的平均值, 并與理論計算值進行了比較,結(jié)果見圖4. 圖4a)、圖4b)、圖4c)與圖4d)分別代表了αI為0.5, 1, 2, 4時的理論與模擬的圖象, 其中橫坐標代表節(jié)點i的編號, 縱坐標代表節(jié)點被感染的平均時間E(Ti), 方形代表理論計算的圖象, 三角形代表實驗?zāi)M的圖象, 理論結(jié)果是根據(jù)公式(4)–公式(12)得到. 可以發(fā)現(xiàn), 我們的理論能夠很好地預(yù)測每個節(jié)點被感染時刻的平均值. 我們發(fā)現(xiàn), 隨著αI的增大, 節(jié)點被感染的平均時刻在整體上逐漸增大, 并且部分節(jié)點在αI比較小的時候, 被感染的平均時刻參差不齊, 卻隨著感染時間分布的異質(zhì)性的減小, 這些節(jié)點被感染的平均時刻趨于相同. 如編號為2到9的節(jié)點, 在αI=0.5 時, 這些節(jié)點被感染的平均時刻各不相同, 但是當αI=4.0, 其對應(yīng)的數(shù)值基本相同. 節(jié)點被感染的平均時刻E(Ti) 的引入, 也讓我們更加清楚地認識到在SI模型的疾病爆發(fā)過程中, 不同異質(zhì)性的感染時間分布下, 疾病從種子節(jié)點到特定目標節(jié)點的時間是呈現(xiàn)各種各樣的非線性性. 在現(xiàn)實生活中, 如果我們掌握了疾病或者消息的傳播網(wǎng)絡(luò)結(jié)構(gòu), 以及對應(yīng)的感染時間分布, 我們就能夠確定出哪些節(jié)點更容易被快速感染.

圖4 真實網(wǎng)絡(luò)中不同參數(shù)下各個節(jié)點被感染平均時刻Fig. 4 Average time for a single node to be infected in real networks with different parameters

3 主要結(jié)論

本文首先給出了閑置邊的概念, 然后通過閑置邊概念的引入提出了一種能夠求解非馬爾科夫SI模型的二階平均場方法, 該方法能夠很好地預(yù)測復(fù)雜網(wǎng)絡(luò)上非馬爾科夫SI模型的爆發(fā)過程, 并且發(fā)現(xiàn)隨著感染時間分布異質(zhì)性的變強, 爆發(fā)的速率也會隨之增加. 通過求解SI模型的時間演化過程,該理論還能夠得到每個節(jié)點被疾病感染的平均時刻, 并且理論計算值也能夠非常準確地預(yù)測實驗?zāi)M結(jié)果, 即我們能夠利用該理論對于特定的感染時間分布, 在特定的網(wǎng)絡(luò)上, 預(yù)測哪些節(jié)點能夠更快地被感染, 哪些節(jié)點被感染需要花費更多的時間, 為控制疾病或謠言在網(wǎng)絡(luò)中的傳播提供理論支撐.