異質(zhì)網(wǎng)絡(luò)上SIS傳染病模型兩種逼近方法的比較

羅曉峰，李 星

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

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

網(wǎng)絡(luò)是把人群內(nèi)的個體作為節(jié)點，個體間的接觸作為邊，若一個節(jié)點的度為k，則該節(jié)點有k條邊[5]。網(wǎng)絡(luò)的度分布pk表示網(wǎng)絡(luò)中隨機選一個節(jié)點的度為k的概率。

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

不同的逼近方法與網(wǎng)絡(luò)結(jié)構(gòu)和鄰居狀態(tài)滿足的分布有關(guān)，因此所得對逼近模型反映的疾病傳播情況和網(wǎng)絡(luò)拓?fù)鋵傩?如網(wǎng)絡(luò)平均度、聚類系數(shù)等)也不同[4]。Morris[16]給出了服從泊松分布和多項分布的逼近公式，前者無網(wǎng)絡(luò)拓?fù)鋮?shù)，后者含網(wǎng)絡(luò)平均度。House等[17]研究了由不同度的節(jié)點構(gòu)成對(邊)的網(wǎng)絡(luò)逼近問題，給出了異質(zhì)對逼近方法。Simon等基于線性假設(shè)提出了一種超緊對逼近方法。考慮了更多的網(wǎng)絡(luò)拓?fù)鋮?shù)，包括節(jié)點的度和網(wǎng)絡(luò)的一階矩(即平均度)， 2階矩和3階矩[18]。不同的逼近公式捕捉的網(wǎng)絡(luò)拓?fù)湫畔⒉煌宜媚Ｐ偷膹?fù)雜度也不同。精度高、含更多的網(wǎng)絡(luò)拓?fù)鋮?shù)，所得模型維數(shù)低是評價對逼近公式優(yōu)劣的關(guān)鍵。

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

1 構(gòu)建異質(zhì)網(wǎng)絡(luò)對逼近模型

下面基于異質(zhì)對逼近和超緊對逼近2種方法，推導(dǎo)相應(yīng)的SIS傳染病模型。

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

1.1.1異質(zhì)對逼近公式

對于SIS對逼近傳染病模型，設(shè)個體總數(shù)即網(wǎng)絡(luò)節(jié)點總數(shù)為N且保持不變，個體按狀態(tài)分為易感者S和染病者I，染病者恢復(fù)后變?yōu)橐赘姓摺?一個染病者和一個易感者每次接觸傳染的概率為τ，染病者的恢復(fù)率為γ。設(shè)網(wǎng)絡(luò)最大度為K

(1)

其中，[SkI]表示網(wǎng)絡(luò)中度為k的易感者Sk與染病者I構(gòu)成的所有二元組的數(shù)量。 根據(jù)反卷積逼近[1]有

(2)

進而由方程式(1)可得

(3)

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

(4)

網(wǎng)絡(luò)中不同狀態(tài)節(jié)點的連邊有如下恒等關(guān)系

(5)

式中：[SI]+[SS]表示狀態(tài)為S得節(jié)點發(fā)出的總邊數(shù)即[S·]；n1表示網(wǎng)絡(luò)的平均度。因此，方程(4)中[II]為

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

1.1.2H-PW模型

(6)

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

(7)

(8)

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

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

以上基于異質(zhì)對逼近建立的模型(8)維數(shù)高且僅含網(wǎng)絡(luò)平均度，既不利于數(shù)學(xué)分析又不能反映網(wǎng)絡(luò)的異質(zhì)性。為了使模型更易分析同時包含更多的網(wǎng)絡(luò)拓?fù)鋮?shù)，下面推導(dǎo)基于線性假設(shè)的超緊對逼近公式，并給出低維模型。

1.2.1超緊對逼近公式

對于給定度分布pk的配置網(wǎng)絡(luò)，文獻[18]假設(shè)sk/pk與k呈線性關(guān)系，即

(9)

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

即

(10)

(11)

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

(12)

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

(13)

其中，k，j=1，2，…，K。顯然，模型(13)相比異質(zhì)SIS對逼近模型(8)維數(shù)沒變，但包含了網(wǎng)絡(luò)的度分布的一階矩、二階矩和三階矩等網(wǎng)絡(luò)拓?fù)鋮?shù)。

(14)

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

(15)

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

2 基本再生數(shù)

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

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

顯然，該方程有K-1個負(fù)實根-γ，當(dāng)且僅當(dāng)

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

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

(16)

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

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

定理1當(dāng)R0<1時，系統(tǒng)(8)，(13)和(15)的無病平衡點是局部漸近穩(wěn)定的。

綜上，系統(tǒng)H-PW，HSH-PW和HSL-PW的基本再生數(shù)相同，三系統(tǒng)在基本再生數(shù)上無區(qū)別，但三維異質(zhì)SIS超緊對逼近模型HSL-PW維數(shù)低且包含了更多的網(wǎng)絡(luò)拓?fù)鋮?shù)，代替高維系統(tǒng)既易進行數(shù)學(xué)分析，又便于研究網(wǎng)絡(luò)拓?fù)洚愘|(zhì)性的影響。

3 模擬與誤差分析

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

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

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

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

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

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

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

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

圖3 當(dāng)R0>1時，[SI]的數(shù)值解與隨機均值的相對誤差曲線

4 結(jié)論

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

2) 數(shù)值上，三維異質(zhì)SIS超緊對逼近模型HSL-PW精度明顯高于其他2個高維模型，即超緊對逼近方法導(dǎo)出的模型維數(shù)低，精度高，且包含了更多的網(wǎng)絡(luò)拓?fù)鋮?shù)。所得結(jié)果豐富了網(wǎng)絡(luò)傳染病傳播動力學(xué)的研究成果。