異質網絡上SIS傳染病模型兩種逼近方法的比較

羅曉峰，李 星

(1.中北大學 理學院， 太原 030051； 2.山西劍橋國際學校， 山西 晉中 030600)

傳染病始終伴隨和影響著人類的生產生活，如新冠肺炎、SARS和艾滋病等[1]。由于傳染病的不可實驗性，傳染病動力學模型是研究傳染病在人群中爆發和流行的規律，進而預防和控制傳染病的重要工具之一[2-3]。基于均勻混合的傳染病動力學模型研究已有百年的歷史，但其忽略了個體行為和個體間的接觸異質性對傳染病傳播的影響。隨著近十幾年網絡科學的發展，網絡傳染病動力學模型的出現很好地克服了這一缺陷[4]。

網絡是把人群內的個體作為節點，個體間的接觸作為邊，若一個節點的度為k，則該節點有k條邊[5]。網絡的度分布pk表示網絡中隨機選一個節點的度為k的概率。

目前，網絡傳染病動力學模型主要是基于節點、度、邊和滲流理論建模[6]。由于人群間的接觸對傳染病傳播和控制有著實質性的影響，它隨時間的變化影響個體隨時間的變化，因此，考慮接觸邊和個體的耦合演化更能揭示傳染病在人群中的傳播規律。對(邊)模型是一類基于邊建模的動力學模型，涵蓋了不同狀態的節點和它們之間形成的對(或邊，也稱二元組)[7-8]。然而，由于低元組的演化離不開高元組，模型無法封閉，因此在單節點和二元組的基礎上應用對逼近方法給出三元組的近似公式封閉方程，這樣得到的封閉模型被稱為對逼近模型[9-10]。近年來也被推廣到離散時間的自適應網絡研究[11]、局部重連的影響研究[12]、聚類網絡的對模型閾值研究[13]，以及疫苗接種博弈[14]和疾病的實際傳播研究中，如新冠肺炎[15]。

不同的逼近方法與網絡結構和鄰居狀態滿足的分布有關，因此所得對逼近模型反映的疾病傳播情況和網絡拓撲屬性(如網絡平均度、聚類系數等)也不同[4]。Morris[16]給出了服從泊松分布和多項分布的逼近公式，前者無網絡拓撲參數，后者含網絡平均度。House等[17]研究了由不同度的節點構成對(邊)的網絡逼近問題，給出了異質對逼近方法。Simon等基于線性假設提出了一種超緊對逼近方法。考慮了更多的網絡拓撲參數，包括節點的度和網絡的一階矩(即平均度)， 2階矩和3階矩[18]。不同的逼近公式捕捉的網絡拓撲信息不同，且所得模型的復雜度也不同。精度高、含更多的網絡拓撲參數，所得模型維數低是評價對逼近公式優劣的關鍵。

因此，本文主要基于SIS傳染病模型比較異質網絡上Keeling的異質對逼近與Kiss的超緊對逼近方法的精度。首先推導出2種逼近下的SIS對逼近模型，進而分析基本再生數和無病平衡點的穩定性，最后通過隨機模擬和誤差分析比較2種逼近方法，給出結論。

1 構建異質網絡對逼近模型

下面基于異質對逼近和超緊對逼近2種方法，推導相應的SIS傳染病模型。

1.1 異質SIS對逼近模型(H-PW)

1.1.1異質對逼近公式

對于SIS對逼近傳染病模型，設個體總數即網絡節點總數為N且保持不變，個體按狀態分為易感者S和染病者I，染病者恢復后變為易感者。 一個染病者和一個易感者每次接觸傳染的概率為τ，染病者的恢復率為γ。設網絡最大度為K

(1)

其中，[SkI]表示網絡中度為k的易感者Sk與染病者I構成的所有二元組的數量。 根據反卷積逼近[1]有

(2)

進而由方程式(1)可得

(3)

對式(3)中的變量[SI]，利用主方程得

(4)

網絡中不同狀態節點的連邊有如下恒等關系

(5)

式中：[SI]+[SS]表示狀態為S得節點發出的總邊數即[S·]；n1表示網絡的平均度。因此，方程(4)中[II]為

因方程(4)中含三元組變量[SSI]，[ISI]，模型不能封閉，為了封閉模型，下面給出異質網絡的三元組逼近，即異質對逼近公式。

1.1.2H-PW模型

(6)

式中：A∈{S，I},式(2)可得異質對逼近公式為

(7)

(8)

式中：k，j=1，2，…，K。

1.2 K+1和三維異質SIS超緊對逼近模型(HS-PW和HSL-PW)

以上基于異質對逼近建立的模型(8)維數高且僅含網絡平均度，既不利于數學分析又不能反映網絡的異質性。為了使模型更易分析同時包含更多的網絡拓撲參數，下面推導基于線性假設的超緊對逼近公式，并給出低維模型。

1.2.1超緊對逼近公式

對于給定度分布pk的配置網絡，文獻[18]假設sk/pk與k呈線性關系，即

(9)

式中：A、B與時間相關。本文將用不同的推導方法給出超緊對逼近公式，由式(9)可得

即

(10)

(11)

將上式代入式(7)，可得基于線性近似的超緊對逼近公式[18]

(12)

1.2.2HSH-PW模型和HSL-PW模型

(13)

其中，k，j=1，2，…，K。顯然，模型(13)相比異質SIS對逼近模型(8)維數沒變，但包含了網絡的度分布的一階矩、二階矩和三階矩等網絡拓撲參數。

(14)

結合方程(4)，利用超緊對逼近式(2)和恒等式N=[S]+[I]可得如下三維異質SIS超緊對逼近傳染病模型(HSL-PW)，

(15)

至此，建立了K+1維異質SIS對逼近模型，K+1維異質SIS超緊對逼近模型和三維異質SIS超緊對逼近模型。下面通過疾病的基本再生數以及模擬和誤差分析比較3個模型的準確性，進而比較2種逼近方法。

2 基本再生數

令系統(8)的右邊等于零，得系統的無病平衡點為E0(N1，N2，…，NK，0)。該系統在無病平衡點E0處的雅可比矩陣為

對應的特征方程為|λE-J1|=0，E為單位矩陣，即

顯然，該方程有K-1個負實根-γ，當且僅當

所有的特征根有負實部時，無病平衡點E0是局部漸近穩定的。f(λ)的判別式滿足

因此，不妨假設f(λ)的特征根分別為λ1，λ2，則根據圓盤定理有

(16)

同理，可求得模型(13)和(15)的基本再生數也為R0。

此外，上面的計算也證明了3種模型無病平衡點的局部穩定性，定理如下：

定理1當R0<1時，系統(8)，(13)和(15)的無病平衡點是局部漸近穩定的。

綜上，系統H-PW，HSH-PW和HSL-PW的基本再生數相同，三系統在基本再生數上無區別，但三維異質SIS超緊對逼近模型HSL-PW維數低且包含了更多的網絡拓撲參數，代替高維系統既易進行數學分析，又便于研究網絡拓撲異質性的影響。

3 模擬與誤差分析

因傳染病在人群中的不可實驗性，隨機模擬傳染病在人群中的傳播情況進而驗證模型是最常見的方法。這部分選取度異質網絡的代表BA無標度網絡[20]，對3種模型H-PW，HSH-PW和HSL-PW進行數值模擬和隨機模擬(10次)比較它們的準確性，并通過誤差分析驗證了超緊對逼近方法在包含更多的網絡拓撲參數，在不失精度的情形下，可將高維模型轉化為低維模型來研究傳染病的傳播。

首先，構造一個節點數為N=5 000，度分布服從冪律分布pk～k-r(r=2.5)，網絡最大度和最小度分別為kmin=10，kmax=100的配置網絡。網絡的總邊數為n1N=9.96×104，其他參數為a=-101.04，b=34.18，n1=19.92，n2=587.23，n3=2.55×104。

當取參數τ=0.01，γ=0.3且R0=0.95<1時， 隨機選擇2 500個節點染病，得到[Sk](k=15)和[SI]的時間序列，如圖1所示，其中實線和虛線代表數值模擬的結果，灰色區域代表了100次隨機模擬的結果(圓圈為均值)。

圖1 當R0=0.95<1時，[Sk](k=15)、[SI]的時間序列曲線

當取參數τ=0.03，γ=0.1且R0=8.55>1時，隨機選擇9個節點染病，得到[Sk](k=15)和[SI]的時間序列(圖2)(其中實線和虛線代表數值模擬的結果，灰色區域代表了100次隨機模擬的結果，圓圈為均值)。

圖2 當R0=8.55>1時，[Sk](k=15)、[SI]的時間序列曲線

從圖1、2(BA無標度網)可以看出：3個模型無論基本再生數大于1還是小于1，度為k的節點 [Sk]和[SI]的數值解(實線和虛線)均穿越隨機(灰色區域)模擬的中心區域，趨勢均吻合得很好，充分驗證了模型H-PW，HSH-PW 和HSL-PW的合理性，但從3個模型的數值結果對比來看，雖然整體趨勢以及最終染病規模吻合很好，但是有一定的誤差。

圖3進一步顯示：當R0>1時，各系統[SI]的數值解與隨機模擬均值的相對誤差(error=|y-x|/x)隨時間變化，可以看出三維異質SIS超緊對逼近模型HSL-PW精度明顯高于其他2個高維模型，驗證了超緊對逼近方法不僅將高維模型轉化為了低維模型且準確性更高，發展和豐富了網絡傳染病傳播的動力學研究成果。

圖3 當R0>1時，[SI]的數值解與隨機均值的相對誤差曲線

4 結論

1) 理論上，2種逼近下三類模型的基本再生數相同，無病平衡點局部穩定。

2) 數值上，三維異質SIS超緊對逼近模型HSL-PW精度明顯高于其他2個高維模型，即超緊對逼近方法導出的模型維數低，精度高，且包含了更多的網絡拓撲參數。所得結果豐富了網絡傳染病傳播動力學的研究成果。