復雜網絡上非馬爾科夫易感-感染模型的二階平均場求解

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

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

0 引 言

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

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

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

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

本文在前人工作的基礎上, 提出了關于非馬爾科夫 SI 模型的二階平均場方法; 根據SI模型中易感態節點能轉換為感染態節點而感染態節點無法轉變為易感態節點的特點, 提出了閑置邊的概念, 該方法能夠求解SI模型的非馬爾科夫爆發過程, 同時也能夠預測每個節點被感染的平均時刻. 本文分別在真實網絡和人工網絡上進行了實驗模擬, 并驗證了以上研究成果的正確性和準確性.

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

1.1 非馬爾科夫SI模型

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

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

1.2 閑置邊的定義

復雜網絡可以由鄰接矩陣A表示, 即

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

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

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

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

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

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

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

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

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

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

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

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

1.3 二階平均場方法

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

對于 〈SI〉 態的具有任意狀態年齡的有向邊, 只要它傳播了疾病, 就能夠將箭頭節點感染為S態, 因此, 能夠得到

被感染的S態節點會轉化為I態節點, 因此可以得到關于Ii(t) 的表達式

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

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

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

進而可以得到

由于 〈SS〉態連邊箭尾節點被感染會轉變為狀態年齡為0的 〈SI〉 態連邊, 因此可以得到

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

其中, 等號右端的第一項表示 〈SS〉態連邊轉變為 〈IS〉態連邊, 第二項和第三項表示 〈IS〉 態連邊的箭尾節點被感染.

由于 〈SI〉態連邊的箭尾節點被感染和 〈IS〉態連邊的箭頭節點被感染都會導致連邊轉化為 〈II〉 態,因此可以得到

其中, 等號右端的第一項表示 〈SI〉態連邊轉變為 〈II〉態連邊, 第二項表示 〈IS〉 態連邊的箭尾節點被除箭頭節點以外的其他節點感染, 第三項表示箭頭節點通過連邊j←i將箭尾節點感染.

1.4 平均感染時刻的預測

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

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

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

2 模擬驗證

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

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

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

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

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

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

圖4 真實網絡中不同參數下各個節點被感染平均時刻Fig. 4 Average time for a single node to be infected in real networks with different parameters

3 主要結論

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