崔玉美，陳姍姍，傅新楚

(上海大學(xué)數(shù)學(xué)系，上海 200444)

0 引言

傳染病歷來是人類健康的大敵，傳染病的防治一直是人類社會面臨的永恒課題。歷史上，如鼠疫、天花、登革熱、非典型性肺炎、肺結(jié)核、艾滋病等疾病多次肆虐全球，嚴(yán)重威脅著人類的生命安全[1-5]。近年來，諸如梅麗莎、灰鴿子、熊貓燒香、機(jī)器狗等計算機(jī)病毒的泛濫，因特網(wǎng)、社交網(wǎng)絡(luò)上各種謠言的傳播給人類社會帶來了新的威脅和挑戰(zhàn)。因此，研究疾病、病毒的傳播和擴(kuò)散機(jī)制以及相應(yīng)的預(yù)防措施是當(dāng)前復(fù)雜系統(tǒng)和傳染病動力學(xué)研究領(lǐng)域的熱點問題。早在1760年，D.Bernouli就用數(shù)學(xué)模型來研究天花的傳播。1911年，Ross為研究瘧疾傳播過程建立的倉室模型為今后傳染病動力學(xué)的研究奠定了基礎(chǔ)。1927年兩位數(shù)學(xué)家Kermack和McKendrick[6]研究黑死病及瘟疫的流行規(guī)律，建立了SIR倉室模型，提出預(yù)測傳染病傳播的“閾值定理”，直到今天仍被廣泛應(yīng)用和不斷完善。

通過微分方程構(gòu)建的傳染病動力學(xué)模型雖然可以刻畫許多疾病的動力學(xué)行為，然而這種基于均勻混合假設(shè)的模型有很大的局限性，它忽略了不同個體傳染能力的差異。不論是疾病、計算機(jī)病毒或謠言等，他們的傳播和擴(kuò)散必然會受到網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的影響。通過將個體抽象為網(wǎng)絡(luò)節(jié)點，個體間的相互作用關(guān)系抽象為節(jié)點之間的連邊，我們能夠更真實地刻畫疾病的傳播過程。1998年Watts和Strogatz[7]提出小世界網(wǎng)絡(luò)及1999年Barabasi 和Albert[8]提出無標(biāo)度網(wǎng)絡(luò)，自此復(fù)雜網(wǎng)絡(luò)理論獲得迅猛發(fā)展。眾多網(wǎng)絡(luò)模型的提出和推導(dǎo)為傳染病動力學(xué)的研究提供了新的方向。其中最經(jīng)典的有3種：一類是2001年由Pastor-Satorras和Vespignani[9-11]提出的平均場模型；一類是2002年由Newman[12]提出的滲流模型以及相關(guān)的分支過程方法和概率生成函數(shù)方法；一類是2003年由Wang[13]等人提出的離散概率模型。此外還有對逼近方法[14]、度有效法[15]等等。文獻(xiàn)[16]中總結(jié)了很多流行病建模的方法，比較全面地分析了這些動力學(xué)解析方法的前提假設(shè)、具體思路、步驟、及其應(yīng)用局限。

在傳染病動力學(xué)模型中，基本再生數(shù)R0是十分重要的參量，它表示在無病平衡狀態(tài)時，引入一個新的染病者，在其平均染病周期內(nèi)所能感染的人數(shù)。R0是決定一種病毒是否流行的標(biāo)志?；驹偕鷶?shù)的計算對疾病的預(yù)防和控制、免疫策略等具有指導(dǎo)性的意義。本文主要介紹了傳染病模型中基本再生數(shù)的導(dǎo)出方法，并且以平均場模型及其相應(yīng)的衍生模型度相關(guān)網(wǎng)絡(luò)模型、對逼近網(wǎng)絡(luò)模型、多菌株網(wǎng)絡(luò)模型、加權(quán)網(wǎng)絡(luò)模型以及時滯微分方程模型傳播模型作為實例，計算了其基本再生數(shù)。

1 基本再生數(shù)的導(dǎo)出方法

我們可以通過以下幾種方法導(dǎo)出基本再生數(shù)：定義法，根據(jù)初始時刻染病者的單調(diào)性導(dǎo)出基本再生數(shù)，根據(jù)正平衡點的存在性導(dǎo)出基本再生數(shù)，根據(jù)無病平衡點的局部穩(wěn)定性即通過在無病平衡點處求解再生矩陣或者雅可比矩陣的方法導(dǎo)出基本再生數(shù)。

相比較而言，根據(jù)定義法導(dǎo)出基本再生數(shù)僅適用于比較簡單的傳染病模型。根據(jù)初始時刻染病者的單調(diào)性導(dǎo)出基本再生數(shù)并不是一個精確地結(jié)果，它并不是判斷疾病爆發(fā)或消亡的充要條件，僅是一個充分條件。根據(jù)無病平衡點的局部穩(wěn)定性，通過計算無病平衡點處的再生矩陣求解基本再生數(shù)的方法需要滿足一系列的條件[17]，通過計算無病平衡點處的雅克比矩陣求解基本再生數(shù)的方法要求所有特征值均具有負(fù)實部。所以本文的第二章中以傳染病模型中比較典型的幾種模型為例子給了出比較精確的導(dǎo)出基本再生數(shù)的方法和具體步驟。

1.1 基本再生數(shù)的定義及性質(zhì)

基本再生數(shù)R0是刻畫傳染病發(fā)病初期的一個重要參量，它表示在易感者人群中，引入一個染病者，該染病者在其平均染病周期內(nèi)所能感染的人數(shù)的期望。當(dāng)R0<1時，染病者在其平均染病期內(nèi)能感染的人數(shù)小于1，那么疾病會自行逐漸消亡。當(dāng)R0>1時，染病者在平均染病期內(nèi)能傳染的人數(shù)大于1，則對于SIS類易感者的補(bǔ)充模型，疾病將一直持續(xù)形成地方?。欢鴮τ赟IR模型，患病者數(shù)量將出現(xiàn)增長，到達(dá)頂峰后再單調(diào)減少而最終趨向于零。

隨著網(wǎng)絡(luò)傳染病動力學(xué)的發(fā)展，物理上常用傳播閾值λc作為刻畫病毒、信息等是否流行的參量，其作用與基本再生數(shù)R0是等價的，因為

(1)

當(dāng)傳播率λ<λc即R0<1時，疾病將會消亡；而當(dāng)傳播率λ>λc即R0>1時，疾病將會流行。因此，R0是判斷疾病是否流行的重要標(biāo)準(zhǔn)。為控制疾病流行，必須盡量減小R0也就是提高傳播閾值λc。

1.2 根據(jù)定義導(dǎo)出R0

對于一個常見的總?cè)丝诤愣ǖ腟IR模型：

(2)

(3)

解上述函數(shù)可得I(α)=I(0)e-(μ+γ)α。令ξ表示一個染病者的感染壽命，可以得到概率函數(shù)

(4)

1.3 根據(jù)無病平衡點的局部穩(wěn)定性導(dǎo)出

在本節(jié)中僅僅對傳統(tǒng)的傳染病模型進(jìn)行一個簡單的分析，最常用的方法就是根據(jù)無病平衡點的局部穩(wěn)定性，采用再生矩陣的方法求出基本再生數(shù)R0。傳統(tǒng)的傳染病模型可以看作均勻網(wǎng)絡(luò)上的傳染病模型，則下文中對于網(wǎng)絡(luò)傳染病模型中基本再生數(shù)的其他求法當(dāng)然也適用于傳統(tǒng)的傳染病模型，在這里就不一一贅述了，在后面會進(jìn)行詳細(xì)的介紹。

設(shè)有

(5)

考慮一個一般的SIR模型：

(6)

其中，S，I,R分別為群體中易感者、染病者及恢復(fù)(移出)者的個數(shù)。這里假定新出生的人口數(shù)為A，自然死亡率為d，恢復(fù)率為γ，傳染率系數(shù)為β。很容易求得該模型的無病平衡點為(0,A/d,0)，而且F=βA/d，V=d+γ，則

(7)

上邊用再生矩陣方法求解了最簡單而且最常見的SIR模型的基本再生數(shù)，下面利用同樣的方法求解推廣的SI模型的基本再生數(shù)。

(8)

其中，S,I分別代表群體中易感者、染病者的個數(shù)在總?cè)丝谥兴嫉谋戎兀@然S+I=1。在該模型中根據(jù)不同的感染階段將I分成m組，易感者被染病者感染進(jìn)入第1階段I1，然后I1以一定的感染概率進(jìn)入第2階段I2……，一直持續(xù)到第m階段Im，且Im不會再繼續(xù)傳染下去。這里設(shè)vi為從i階段向i+1階段發(fā)展的平均概率，di為染病者在第i階段的死亡率，b為個體的出生率和易感者的自然死亡率，βk為處于第k階段的個體的傳染率系數(shù)。顯然，該模型的無病平衡點為(0,…,0,1)，而且可以求得

(9)

(10)

令V-1=(αij)，則

通過計算再生矩陣的譜半徑可以求得基本再生數(shù)

(11)

1.4 網(wǎng)絡(luò)傳染病模型中由初始時刻染病者單調(diào)性導(dǎo)出R0

對于一般的SIR網(wǎng)絡(luò)傳染病動力學(xué)模型：

(12)

(13)

(14)

1.5 網(wǎng)絡(luò)傳染病模型中由正平衡點的存在性導(dǎo)出R0

對于易感者有補(bǔ)充的模型，當(dāng)R0>1時，即一個染病者在平均染病期內(nèi)能傳染的人數(shù)大于1時，染病者數(shù)量逐漸增加，疾病一致持續(xù)形成地方病平衡點。則可以根據(jù)地方病平衡點的存在性判定疾病是否流行導(dǎo)出基本再生數(shù)。這方面工作主要來自文獻(xiàn)[18]。對于一般的SIS網(wǎng)絡(luò)傳染病動力學(xué)模型：

(15)

(16)

顯然Θ=0始終是該自治方程的一個平凡解。下面導(dǎo)出自治方程(16)存在正解0<Θ<1的條件。

記

(17)

經(jīng)計算有：

(18)

(19)

則F(Θ)為凹函數(shù)，因F(0)=0，且

則F(Θ)=0在(0,1)上有唯一正解的充分必要條件是在Θ=0時有

(20)

從而得到基本再生數(shù)：

(21)

當(dāng)R0>1時，系統(tǒng)(5)存在唯一的正平衡點。

1.6 網(wǎng)絡(luò)傳染病模型中由無病平衡點的穩(wěn)定性導(dǎo)出R0

根據(jù)初始時刻染病者的單調(diào)性，或者正平衡點的存在性導(dǎo)出的R0，只是對基本再生數(shù)的一種相對粗略的估計。比如根據(jù)正平衡點的存在性導(dǎo)出的R0，條件R0>1只能保證方程(15)存在正平衡點，并沒有證明該正平衡點的穩(wěn)定性。也就是說當(dāng)R0>1時，疾病不一定會持續(xù)形成地方病。同理，R0<1只說明方程(15)不存在正平衡點，疾病不一定會逐漸消亡。為導(dǎo)出基本再生數(shù)的相對精確的表達(dá)式，可以應(yīng)用微分方程穩(wěn)定性理論的相關(guān)方法[19]。R0作為區(qū)分疾病是否流行的參考量，當(dāng)R0<1時，少量感染者進(jìn)入系統(tǒng)后疾病會逐漸消亡，即無病平衡點(零解)是局部漸近穩(wěn)定的。根據(jù)系統(tǒng)局部漸近穩(wěn)定性的條件即可導(dǎo)出R0的表達(dá)式。例如對于1.5中的式(15)，

設(shè)網(wǎng)絡(luò)的最大度為n，則該方程在無病平衡點Ik=0處的Jacobian矩陣為

(22)

通過無病平衡點穩(wěn)定性導(dǎo)出R0是目前最精確的計算傳播閾值的方法。但對于相對復(fù)雜的微分方程系統(tǒng)，求解無病平衡點處Jacobian矩陣的所有特征值計算量較大，2002年Van den Driessche和Willboordse[17]提出了再生矩陣方法，通過計算再生矩陣的譜半徑求得無病平衡點的漸近穩(wěn)定性，從而導(dǎo)出基本再生數(shù)。

對于異質(zhì)網(wǎng)絡(luò)模型，記x=(x1,x2,…,xn)t，其中xi>0(i=1,2,…,n)表示每一組個體的數(shù)量或者比例。為方便計算將x分成兩部分，前m項xi,…,xm表示感染組個體，剩下的xm+1,…,xn為非感染組個體。感染組與非感染組個體的劃分應(yīng)根據(jù)流行病學(xué)的生物意義，非單純通過數(shù)學(xué)方程劃分。設(shè)有如下微分方程：

(23)

R0=ρ(FV-1)

例如，對于帶有出生死亡的異質(zhì)網(wǎng)絡(luò)SIR模型：

(24)

利用再生矩陣的譜半徑計算其基本再生數(shù)：

胃癌是嚴(yán)重威脅我國居民健康的重大疾病，其發(fā)病率居惡性腫瘤第 2 位，死亡率居第 3 位[1]。辨證論治是中醫(yī)臨床診療的核心，但目前國內(nèi)對胃癌的證候尚無統(tǒng)一的認(rèn)知，也未形成確能提高療效的共識。因此，基于中醫(yī)自身的實踐性與經(jīng)驗性，在臨床中提出新的學(xué)術(shù)觀點并在實踐中進(jìn)行驗證，是提高胃癌中醫(yī)診治水平的必由之路。海軍軍醫(yī)大學(xué)（第二軍醫(yī)大學(xué)）長征醫(yī)院中醫(yī)科作為國家教育部中西醫(yī)結(jié)合臨床重點學(xué)科，在名老中醫(yī)魏品康教授帶領(lǐng)下，幾十年來致力于胃癌的中西醫(yī)結(jié)合防治研究，在實踐中提出并構(gòu)建了較為系統(tǒng)的胃癌痰證理論。

第2步：求方程(24)對應(yīng)的無病平衡點：

第3步：在無病平衡點E0處求解：

第4步：計算譜半徑ρ(FV-1)得，系統(tǒng)(24)的基本再生數(shù)為

(25)

1.7 網(wǎng)絡(luò)傳染病模型中由數(shù)值估算閾值

對于網(wǎng)絡(luò)傳染病模型的不同理論分析方法，由于分析角度不同，可能會導(dǎo)出不同的基本再生數(shù)或者疾病的傳播閾值，甚至有些傳染病模型不能夠通過理論分析求出流行閾值,正確且具體的表達(dá)形式。而數(shù)值模擬能夠更加真實地刻畫疾病傳播的規(guī)律，得到更加真實有效的數(shù)據(jù)。因此數(shù)值模擬是理論分析的一個有效估計和補(bǔ)充，對于求解傳染病模型的傳播閾值起著非常重要的指導(dǎo)作用。2012年Silvio C. Ferreira等[20]通過對有限網(wǎng)絡(luò)上SIS模型進(jìn)行了大量的數(shù)值分析，提出可以通過敏感性峰值的方法有效估計傳染病的流行閾值，并與異質(zhì)平均場理論(HMF)、淬火平均場理論(QMF)得到的結(jié)果進(jìn)行了比較分析。2014年P(guān)anpan Shu等[21]對其進(jìn)行了補(bǔ)充說明，針對有限網(wǎng)絡(luò)上SIR模型提出了估計傳染病流行閾值的新方法——變異性峰值方法。2015年Can Liu等[22]將變異性峰值的方法，應(yīng)用到基于人們行為反應(yīng)的傳染病模型，估計出合理的流行閾值，并且指出了利用平均場方法建模的局限性。本節(jié)主要介紹一下參考文獻(xiàn)[20],[21]中提出的兩種典型的估計傳染病流行閾值的方法。

敏感性峰值方法：文獻(xiàn)[20]中，將敏感性χ定義為

(26)

敏感性峰值方法是預(yù)測SIS模型流行病閾值的有效方法，但是由于SIR模型中λ>λc時爆發(fā)大小的雙峰分布導(dǎo)致敏感性峰值方法用于預(yù)測SIR模型流行病閾值時，往往高于SIR模型實際的流行病閾值。所以參考文獻(xiàn)[21]中提出了預(yù)測流行病閾值的新方法——變異性峰值方法。

變異性峰值方法：變異性

(27)

文獻(xiàn)[22]通過比較在隨機(jī)規(guī)則網(wǎng)絡(luò)(RRN)SIS模型、SIR模型中應(yīng)用變異性峰值方法估計的流行閾值與異質(zhì)平均場理論(HMF)得到的流行閾值，發(fā)現(xiàn)二者是一致的，證明了變異性峰值方法估計SIS模型的流行閾值的準(zhǔn)確性。通過研究無標(biāo)度網(wǎng)絡(luò)和實際網(wǎng)絡(luò)Facebook等SIR模型的模擬閾值進(jìn)一步表明，變異性峰值方法估計流行閾值比敏感性峰值方法準(zhǔn)確。

2 幾類典型網(wǎng)絡(luò)傳播模型的基本再生數(shù)的計算

2.1 度相關(guān)網(wǎng)絡(luò)模型

在研究網(wǎng)絡(luò)上的動力學(xué)行為時，網(wǎng)絡(luò)的結(jié)構(gòu)對節(jié)點的動力學(xué)行為具有重要作用。對于現(xiàn)實網(wǎng)絡(luò)，還有很多情況需要考慮。真實世界網(wǎng)絡(luò)中節(jié)點之間的連接選擇并不均等，而是存在明顯的偏好。大量實驗研究表明[23]，蛋白質(zhì)交互網(wǎng)絡(luò)和神經(jīng)網(wǎng)絡(luò)等生物網(wǎng)絡(luò)以及互聯(lián)網(wǎng)等技術(shù)網(wǎng)絡(luò)是異配的，也就是說網(wǎng)絡(luò)中度大的節(jié)點更傾向于和度小的節(jié)點相連。而包括科研人員合作網(wǎng)或者電影演員合作網(wǎng)，微博，F(xiàn)acebook等在內(nèi)的許多社交網(wǎng)絡(luò)則多呈現(xiàn)較為明顯的同配特性。從社會學(xué)與心理學(xué)角度分析，精英人士往往更傾向于和不低于自己地位的人交往，從而形成了節(jié)點度的正相關(guān)性。2002年，Marian Boguna和Romualdo Pastor-Satorras[24]研究了度相關(guān)網(wǎng)絡(luò)動力學(xué)模型，并給出了基本再生數(shù)的計算方法，本節(jié)主要介紹文獻(xiàn)[24]，[25]的相關(guān)工作。

kp(k′|k)p(k)=k′p(k|k′)p(k′)≡〈k〉p(k,k′)

(28)

(29)

則Θk表示一條連接度為k的節(jié)點的邊的另一端指向染病者的概率。由此可以建立度相關(guān)異質(zhì)網(wǎng)絡(luò)SIS平均場模型：

(30)

其無病平衡點Ik=0，k=1,2,…,n處的Jacobian矩陣L={lkk′}，其中

lkk′=-γδkk′+λkp(k′|k)

(31)

(32)

其中,Λm是網(wǎng)絡(luò)的連接矩陣C={ckk′}，ckk′=kp(k′|k)的最大特征值。由此可見，度相關(guān)網(wǎng)絡(luò)的基本再生數(shù)不僅與網(wǎng)絡(luò)尺寸有關(guān)，還受到網(wǎng)絡(luò)的相關(guān)性的影響，即條件概率p(k′|k)的表達(dá)式?jīng)Q定著R0的取值。

2.2 對逼近網(wǎng)絡(luò)模型

在考慮傳染性疾病在人群中的傳播時，疾病的傳染只有在染病者(I)和易感者(S)接觸(即有邊相連)時才可能進(jìn)行。比如性傳播疾病，在忽略同性性行為的假設(shè)下，只有異性之間發(fā)生性交(即網(wǎng)絡(luò)中兩點相連形成邊)時，疾病才有可能傳播?；诠?jié)點建立的平均場模型不能充分刻畫疾病或者信息的傳播必須通過網(wǎng)絡(luò)連邊進(jìn)行這一基本特征。因此科學(xué)家建立了以邊為研究對象的對逼近(Pair approximation)模型。記[S]、[I]、[SI]分別代表易感者、染病者及二者構(gòu)成的二元組的數(shù)量。如果是完全圖，顯然有[SI]=[S][I]。但如果不是完全圖，則上述微分方程不能封閉，需要給出[SI]隨時間變化的動力學(xué)方程，這又涉及到[SS]、[II]以及三元組[SSI]等隨時間的變化，最終將得到一組無限維的微分方程組。為了對模型進(jìn)行動力學(xué)分析，就需要在誤差允許的范圍內(nèi)對方程組進(jìn)行封閉，這個過程稱為“矩封閉”。對逼近是最簡單的一種“矩封閉模型”，是關(guān)于二元組層面上的封閉，得到的是二元組變化的微分方程組。在矩封閉過程中，必須考慮網(wǎng)絡(luò)的結(jié)構(gòu)特性。在靜態(tài)網(wǎng)絡(luò)中，對于規(guī)則網(wǎng)絡(luò)和節(jié)點的度波動不太大的隨機(jī)網(wǎng)絡(luò)的研究，已有大量研究成果。

1995年Altmann[26]首次建立了研究性病傳播的對關(guān)系SIR模型，并利用Markov過程計算疾病傳播的基本再生數(shù)。隨后Keeling和Rand[27]等建立了SEIR對逼近模型來研究麻疹的傳播。2001年Filipe[28]提出了正方形逼近方法(SA方法)、對角線逼近方法(DA方法)、混合成對團(tuán)簇逼近方法(HPA方法)以及Sato團(tuán)簇逼近方法。2003年Thomson[29]建立了有死亡的SIS網(wǎng)絡(luò)傳染病模型，提出了雙邊緣的對逼近方法(PEA方法)。2007年，Trapman[30]給出了均勻網(wǎng)絡(luò)結(jié)構(gòu)下的對封閉模型理論推導(dǎo)，定義了該網(wǎng)絡(luò)下的有效再生數(shù)及其計算方法。下面將結(jié)合[30],[31]在均勻網(wǎng)絡(luò)下建立一般的對逼近模型并求解。

因疾病的傳播只有染病者(I)和易感者(S)接觸才可能進(jìn)行，據(jù)此可建立總節(jié)點數(shù)N不變的SIS對逼近傳染病模型：

(33)

其中,[S]，[I]，[SI]分別代表易感者，染病者及二者構(gòu)成二元組的數(shù)量。二元組[SI]的存在使方程(33)不能封閉，需要給出[SI]隨時間變化的動力學(xué)方程。因為節(jié)點狀態(tài)的變化受到與之相連的節(jié)點的影響，則根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)可以用三元組來表示二元組[SI]的變化：

(34)

方程右端出現(xiàn)三元組使方程不能封閉，在此做近似截斷如下：在均勻網(wǎng)絡(luò)如ER隨機(jī)網(wǎng)絡(luò)中，如果節(jié)點的染病者鄰居滿足多項分布，則三元組[ABC]可由二元組逼近，近似公式為

(35)

其中,n表示網(wǎng)絡(luò)的最大度。則對(34)中的三元組[SSI]與[ISS]可以近似為

(36)

(37)

將近似等式(36)(38)代入等式(34)封閉方程得到均勻網(wǎng)絡(luò)下對逼近SIS模型：

(38)

其中,2[SI]+[II]+[SS]=nN,[S]+[I]=N。下面應(yīng)用再生矩陣譜半徑的方法計算R0。

首先求得方程(38)的無病平衡點：

E0=([S]*,[I]*,[SS]*,[SI]*,[II]*)=(N,0,nN,0,0)

(39)

當(dāng)R0<1時，無病平衡點漸近穩(wěn)定，疾病趨于消亡；而當(dāng)R0>1時，平衡點不穩(wěn)定。

2.3 多菌株網(wǎng)絡(luò)模型

疾病傳播過程中，對于不同的疾病，病菌的傳染力和疾病的流行機(jī)理不同。傳染病往往是多種病原體或多種菌株相互作用的結(jié)果。比如引起艾滋病(AIDS)的HIV病毒，有多種互不相同的菌株比如HIV-1，HIV-2等。以下主要根據(jù)文獻(xiàn)[32]，[33]介紹兩種簡單的異質(zhì)網(wǎng)絡(luò)多菌株傳染病模型及其基本再生數(shù)的計算。

2.3.1 不同感染力的多菌株SIS模型

在疾病的傳播過程中，不同的人群具有相異的感染力。把具有很強(qiáng)傳染力的人群稱為核心人群。本小節(jié)主要介紹S.Nobuaki與A.Kazuyuki在文獻(xiàn)[32]中的工作。在易感者與染病者每次接觸過程中，易感者Sk一部分以概率1-p轉(zhuǎn)化為具有較弱傳染力λ1的一般感染者Uk，其余以概率p轉(zhuǎn)化為具有較強(qiáng)傳染力λ2的核感染者Vk。假設(shè)網(wǎng)絡(luò)是度不相關(guān)的，在此情形下建立模型：

(40)

此處用正平衡點存在性導(dǎo)出基本再生數(shù)。首先令方程(40)右邊等于零，得到模型的正平衡點：

(41)

(42)

因為F(0)=0，

(43)

2.3.2 競爭病毒的SIS模型

在多菌株異質(zhì)網(wǎng)絡(luò)傳染病模型中，相互作用的兩個病毒之間主要存在3種關(guān)系：1)Competing pathogens，表示兩種病毒相互競爭同一個寄主，不能共存；2)Super-infection，表示在同一個寄主中，一種病毒可以戰(zhàn)勝另一種病毒取而代之；3)Co-infection，表示兩病毒可以共存于同一個寄主中。以下介紹文獻(xiàn)[33]中關(guān)于異質(zhì)網(wǎng)絡(luò)競爭病毒的SIS模型及其基本再生數(shù)的計算。

(44)

(45)

其中,

(J1)kk′=-γ1δkk′+β1kp(k′)/〈k〉

(46)

(J2)kk′=-γ2δkk′+β2kp(k′)/〈k〉

(47)

(48)

(49)

(50)

2.4 帶有權(quán)重的無標(biāo)度網(wǎng)絡(luò)模型

現(xiàn)實世界的真實網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)一般是異質(zhì)的，大多具有無標(biāo)度特性，而這種不均勻性不僅表現(xiàn)為連邊的數(shù)量不同，同時也體現(xiàn)在邊上權(quán)重的差異。以接觸網(wǎng)絡(luò)為例，接觸網(wǎng)絡(luò)上連邊的權(quán)重，往往表示連接節(jié)點之間的親密程度或接觸時間的長度，從而可以用來衡量疾病傳播幾率的大小，權(quán)重越大，連接的兩個節(jié)點越親密，則被感染的概率越大[34-37]。同時，在很多傳染病模型中，一定時間內(nèi)，感染節(jié)點能夠接觸到的鄰居數(shù)(傳染力)為節(jié)點度。但實際情況下，感染節(jié)點的傳染力應(yīng)該是與節(jié)點度有關(guān)的函數(shù)。文章[38]基于以上兩點介紹了無標(biāo)度網(wǎng)絡(luò)上帶有權(quán)重的SIR模型，下面主要介紹該文中對基本再生數(shù)的計算。

(51)

(52)

(53)

顯然，為了得到基本再生數(shù)需要計算出Φ(t)的值。根據(jù)(52)及(53)得，

(54)

(55)

(56)

(57)

2.5 時滯微分方程網(wǎng)絡(luò)模型

常微分方程是研究傳染病傳播機(jī)制和動力學(xué)行為的一個重要手段，但這是一種理想化的方法?，F(xiàn)實生活中，如流感、SARS病毒、艾滋病、計算機(jī)病毒等，雖然一個易感者與一個染病者接觸會以一定概率感染疾病，但此時病毒或病菌只是在其體內(nèi)發(fā)展。經(jīng)過一段時間(即潛伏期)后，才可能具有感染力或表現(xiàn)出癥狀。另外，染病者患病后也不可能立即恢復(fù)成易感者或移除者，這需要一段時間(患病期)的治療或自我恢復(fù)。針對于此，學(xué)者們在模型中加入時滯來刻畫潛伏期或恢復(fù)期對疾病傳播的影響。對于經(jīng)典的傳染病模型，這方面的研究成果很多[39]。由于網(wǎng)絡(luò)模型都是高維微分方程，再將時滯引入模型會使模型變得更加復(fù)雜。因此網(wǎng)絡(luò)中利用時滯模型研究傳染病的相關(guān)文章很少。文獻(xiàn)[40]與文獻(xiàn)[41]分別用離散和連續(xù)時滯微分方程討論了異質(zhì)網(wǎng)絡(luò)上傳染病的動力學(xué)行為。這兩篇文獻(xiàn)所用的求基本再生數(shù)的方法相似，我們以文獻(xiàn)[40]為例介紹網(wǎng)絡(luò)時滯微分方程的基本再生數(shù)的計算方法。

設(shè)網(wǎng)絡(luò)節(jié)點數(shù)為保持不變，節(jié)點的最大度為。建立帶有時滯的SIS微分方程模型，當(dāng)時有：

(58)

(59)

對等式(59)兩邊積分可得：

(60)

初始階段即0≤t≤τ時，系統(tǒng)(59)或(60)服從一般的SI微分方程：

(61)

(62)

(63)

顯然Θ*=0是方程(63)的解，即為該時滯系統(tǒng)的無病平衡點。

則f(Θ*)是(0,1)上的凸函數(shù)，在區(qū)間(0,1)上有正解的充要條件是：

由此導(dǎo)出網(wǎng)絡(luò)上時滯微分方程的基本再生數(shù)為

(64)

特別地，考慮年齡結(jié)構(gòu)的網(wǎng)絡(luò)傳播模型和時滯微分方程模型形式是一致的，在導(dǎo)出具有年齡結(jié)構(gòu)的微分方程模型對應(yīng)的基本再生數(shù)時，基本的方法是首先將偏微分方程化成常微分方程，然后通過基本再生數(shù)的定義或根據(jù)無病平衡點的穩(wěn)定性來導(dǎo)出年齡結(jié)構(gòu)模型的基本再生數(shù)。

2.6 有因病死亡的SIR脈沖預(yù)防接種模型

該模型一般根據(jù)無病周期解的存在性和局部漸近穩(wěn)定來求解基本再生數(shù)。由于對網(wǎng)絡(luò)上脈沖接種模型進(jìn)行分析存在一定的困難,這方面的文獻(xiàn)資料還不夠完善，所以這里考慮當(dāng)預(yù)防接種在離散的時間點(t=k,k=0,1,2…)以周期為1進(jìn)行時，可以用以下脈沖系統(tǒng)來修正傳統(tǒng)的SIR模型。

(65)

(66)

下面討論無病1周期解的局部漸近穩(wěn)定性。

引理1[42]設(shè)有周期系統(tǒng)

(67)

其中，f∈C[R×Rn,Rn]，φk∈C[Rn,Rn]，f(t+1,y)=f(t,y),?t∈R；φk+1=φk；且系統(tǒng)(65)關(guān)于其周期解y(t)的線性近似系統(tǒng)為

(68)

并設(shè)Φ(t)是(68)中方程的一個基解矩陣，即滿足

若矩陣M=B1Φ(1)的一切特征根的絕對值均小于1，則系統(tǒng)(68)的零解，也即系統(tǒng)(67)的周期解y(t)局部漸近穩(wěn)定。

下面做變換x(t)=S(t)-S*(t),y(t)=I(t)-0,z(t)=R(t)-R*(t)，可以得到原系統(tǒng)的關(guān)于1周期解(S*(t),0,R*(t))的線性近似系統(tǒng)具有以下形式：

(69)

其基解矩陣為

(70)

其中,

由方程組(66)可得

(71)

應(yīng)用引理1，

(72)

可以求出M的3個特征根分別為

λ1=(1-p)e-b,λ2=φ22(1),λ2=e-b

由于00,所以0<λ1<1,0<λ3<1,

2.7 動態(tài)免疫網(wǎng)絡(luò)模型

在上一節(jié)中提到對于網(wǎng)絡(luò)上的脈沖接種模型的分析是具有一定困難的，目前還沒有特別固定的模型和方法供我們使用。在本節(jié)中以文獻(xiàn)[43]為例，考慮一個動態(tài)免疫網(wǎng)絡(luò)模型的基本再生數(shù)的估計方法，利用近似函數(shù)根據(jù)無病平衡點的穩(wěn)定性導(dǎo)出基本再生數(shù)R0，它為我們估計網(wǎng)絡(luò)上脈沖接種模型的傳播閾值提供了思路。

設(shè)網(wǎng)絡(luò)節(jié)點總數(shù)為N，而且它可以抽象為一個由點集V和邊集E構(gòu)成的圖G=(V,E)。A表示G的鄰接矩陣，即當(dāng)G中節(jié)點i和節(jié)點j有邊相連時，aij=1否則aij=0。

一個節(jié)點如果與一個或者多個節(jié)點相連時，應(yīng)該先被免疫。否則，它很有可能被鄰居節(jié)點所感染，這樣的節(jié)點稱為高危節(jié)點。直覺上，我們可以通過控制高危節(jié)點來控制疾病的傳播。這一模型中主要考慮某一集合Ω中的高危節(jié)點，顯然?Ω?V。

文中通過馬爾可夫鏈的方法構(gòu)建了如下傳染病模型:

(73)

(74)

其中,β為傳染率系數(shù)，δ為免疫概率，Ν(i)為節(jié)點i的鄰域，即與節(jié)點i有邊相連的點的集合。

首先考慮模型的平衡點。由于模型方程的第2個等式和第3個等式都是Ω中節(jié)點的方程，所以我們先考慮這兩個等式。令這兩個等式右端項等于0，且qi,t=qi，pi,t=pi可得pi=0。

下面我們考慮無病平衡點pi=0，i∈V/Ω的局部穩(wěn)定性。利用當(dāng)a<<1，b<<1時(1-a)(1-b)≈1-a-b這一近似，可以把模型簡化成以下形式：

(75)

要使無病平衡點pi=0，i∈V/Ω局部漸進(jìn)穩(wěn)定，則

其中ΛmaxA表示矩陣A的最大特征值。該模型的基本再生數(shù)

3 時變的基本再生數(shù)

以上我們在確定各類參數(shù)以后求解的基本再生數(shù)都是一個固定的值，而傳染病的傳播是一個動態(tài)的過程，因此基本再生數(shù)可能會出現(xiàn)時變的情況。文獻(xiàn)[44]表明，不同時期，基本再生數(shù)是會發(fā)生變化的。這一節(jié)我們就根據(jù)文獻(xiàn)[44]描述的埃博拉病毒的傳播特點建立合理的數(shù)學(xué)模型，通過計算基本再生矩陣的譜半徑求出該模型的基本再生數(shù)。然后根據(jù)該模型中基本再生數(shù)的特點，在社區(qū)和衛(wèi)生保健環(huán)境中通過調(diào)整基本傳染率、診斷率和提高感染控制措施，來說明控制干預(yù)措施在埃博拉病毒爆發(fā)期間的影響。我們將人口分為5類：易感個體(S)、暴露個體(E)、感染個體(I)、住院個體(H)、康復(fù)或者死亡后移除個體(R)。具體模型如式(76)。

(76)

其中,N(t)為t時刻的人口總數(shù)，β(t)為傳染率系數(shù)，l(t)量化了與社區(qū)中的染病者相比，住院患者的相對傳染率。0≤l(t)<1反映了醫(yī)院隔離措施的有效性，將埃博拉病毒的傳播概率降到社區(qū)以下。感染個體以依賴于時間的診斷率γa(t)住院或者以概率γI康復(fù)而不住院，住院個體概率γR康復(fù)。

(77)

R0=Rcomm+Rhosp

因為傳染率系數(shù)β(t)、相對傳染率l(t)、診斷率γa(t)均是隨時間變化的，所以在埃博拉病毒傳播期間，基本再生數(shù)也是隨時間不斷變化的。因此在埃博拉病毒傳播期間，埃博拉病毒爆發(fā)區(qū)域可以通過隔離或者增加防護(hù)裝備、提高醫(yī)療水平等措施，來減少疾病爆發(fā)的可能性。這也說明了埃博拉病毒在發(fā)達(dá)國家不太可能爆發(fā)，而在衛(wèi)生系統(tǒng)極差或者是缺乏的國家比如非洲國家爆發(fā)的可能性相對較大?？梢姰?dāng)?shù)厣鐣?jīng)濟(jì)和社會文化條件是疾病傳播的關(guān)鍵性因素。

4 結(jié)論

隨著1998年小世界網(wǎng)絡(luò)及1999年無標(biāo)度網(wǎng)絡(luò)在真實系統(tǒng)中的發(fā)現(xiàn)和建立，復(fù)雜網(wǎng)絡(luò)上的動力學(xué)行為研究是最近科學(xué)討論的一個熱點。采用復(fù)雜網(wǎng)絡(luò)理論與流行病學(xué)相結(jié)合的方法已成為傳染病動力學(xué)建模的主要趨勢。本文總結(jié)整理了傳染病動力學(xué)中基本再生數(shù)的幾種主要導(dǎo)出方法。即通過基本再生數(shù)的定義，初始時刻染病者的單調(diào)性，正平衡點的存在性，無病平衡點的穩(wěn)定性、數(shù)值模擬導(dǎo)出基本再生數(shù)。作為示例，本文分別介紹了度相關(guān)模型，對逼近模型，多菌株模型，加權(quán)網(wǎng)絡(luò)模型和時滯微分方程模型5種網(wǎng)絡(luò)傳播模型，并導(dǎo)出了基本再生數(shù)的精確表達(dá)式。對于網(wǎng)絡(luò)傳播模型，基本再生數(shù)不僅與傳染參數(shù)(如傳染率，恢復(fù)率)有關(guān)，還與網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)有關(guān)，而且基本再生數(shù)也有能出現(xiàn)時變的情況，因為疾病的爆發(fā)過程是一個動態(tài)過程，模型參數(shù)是隨時間在改變的，所以不同時期，基本再生數(shù)是會發(fā)生變化的。

現(xiàn)有的關(guān)于基本再生數(shù)的計算方法并不精確，從而使由本文給出的幾種方法導(dǎo)出的基本再生數(shù)R0不能同時保證當(dāng)R0<1時疾病消亡，而當(dāng)R0>1時疾病暴發(fā)。這需要進(jìn)一步分析無病平衡點和正平衡點的全局動力學(xué)行為。這將成為今后網(wǎng)絡(luò)傳播動力學(xué)研究中的重點和難點。

由于疾病傳播方式和網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的多樣性，在本文中我們沒有辦法完全列舉出所有傳染病模型的基本再生數(shù)的導(dǎo)出過程。針對于特殊的傳染病模型可能使用幾種方法結(jié)合考慮才能到較為準(zhǔn)確的基本再生數(shù)，甚至只能估計出一個基本再生數(shù)。但本文所總結(jié)的上述方法是主要的，且能夠廣泛應(yīng)用于基本再生數(shù)的導(dǎo)出，相信這些方法能夠給讀者提供有益的參考。

感謝課題組的同仁們，特別是孫孟鋒的熱心幫助和大力支持！

[1]Feng Z L, Jorge X, Velasco-Hernández. Competitive exclusion in a vector-host model for the dengue fever[J]. Journal of Mathematical Biology, 1997, 35: 523-544.

[2]Chan-Yeung M, Xu R H. SARS: epidemiology[J]. Respirology, 2003, 8: 9-14.

[3]Small M, Tse C K, Walker D M. Super-spreaders and the rate of transmission of the SARS virus[J]. Physica D, 2006, 215(2): 146-158.

[4]Castillo-Chavez C, Feng Z L. To treat or not to treat: the case of tuberculosis[J]. Journal of Mathematical Biology, 1997, 35: 629-656.

[5]Castillo-Chavez C, Huang W, Li J. Competitive exclusion in gonorrhea models and other sexually-transmitted diseases[J]. SIAM Journal on Applied Mathematics, 1996, 56: 494-508 .

[6]Kermack W O, Mckendrick A G. Contributions to the mathematical theory of epidemics[J]. Proceedings of the Royal Society A, 1927,115: 700-721.

[7]Watts D J, Strogatz S H. Collective dynamics of small-world networks[J]. Nature, 1998, 393: 440-442.

[8]Barabasi A L, Albert R. Emergence of scaling in random networks[J]. Science, 1999, 286(5439): 509-512.

[9]Pastor-Satorras R, Vespignani A. Epidemic spreading in scale-free networks[J]. Physical Review Letters, 2000, 86: 3200-3203.

[10] Pastor-Satorras R, Vespignani A. Epidemic dynamics in finite size scale-free networks[J]. Physical Review E, 2002, 65(3): 035108.

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

[12] Newman M E J. Spread of epidemic disease on networks[J]. Physical Review E, 2002, 66(1): 016128.

[13] Wang Y, Chakrabarti D, Wang C, et al. Epidemic spreading in real networks: an eigenvalue viewpoint[J]. International Symposium on Reliable Distributed Systems, 2003, 10(13)：25-43.

[14] Eames K T D. Modelling disease spread through random and regular contacts in clustered populations[J]. Theoretical Population Biology, 2008, 73(1): 104-111.

[15] Lindquist J, Ma J, Driessche P V D, et al. Effective degree network disease models[J]. Journal of Mathematical Biology, 2011, 62: 43-164.

[16] 李睿琪, 王偉, 舒盼盼, 等. 復(fù)雜網(wǎng)絡(luò)上流行病傳播動力學(xué)的爆發(fā)閾值解析綜述[J]. 復(fù)雜系統(tǒng)與復(fù)雜性科學(xué), 2016, 13(1):1-39.

Li Ruiqi, Wang Wei, Shu Panpan, et al. Review of threshold theoretical analysis about epidemic spreading dynamics on complex networks[J]. Complex Systems and Complexity Science, 2016, 13(1):1-39.

[17] Van d D P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[J]. Mathematical Biosciences, 2002, 180(1/2):29-48.

[18] Pastor-Satorras R, Vespignani A. Immunization of complex networks[J]. Physical Review E, 2002, 65(3): 036104.

[19] 馬知恩, 周義倉. 常微分方程定性與穩(wěn)定性方法[M]. 北京: 科學(xué)出版社, 2001.

[20] 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]. Physical Review E, 2012, 86(4): 041125.

[21] Shu P, Wang W, Tang M, et al. Simulated identification of epidemic threshold on finite-size networks[J]. Chaos, 2015, 544(10):167.

[22] Liu C, Xie J R, Chen H S, et al. Interplay between the local information based behavioral responses and the epidemic spreading in complex networks[J]. Chaos, 2015, 25(10): 103111.

[23] Newman M E J. Assortative mixing in networks[J]. Physical Review Letters, 2002, 89(20): 208701.

[24] Boguna M, Pastor-Satorras R. Epidemic spreading in correlated complex networks[J]. Physical Review E, 2002, 66(4): 047104.

[25] Boguna M, Pastor-Satorras R, Vespignani A. Absence of epidemic threshold in scale-free networks with degree correlations[J], Physical Review Letters, 2003, 90(2): 028701.

[26] Altmann M. Susceptible-Infected-Removed epidemic models with dynamic partnerships[J]. Journal of Mathematical Biology, 1995, 33: 661-675.

[27] Keeling M J, Rand D A, Morris A J. Correlation models for childhood epidemics[J]. Proc Biol Sci, 1997, 264(1385): 1149-1156.

[28] Filipe J A N, Gibson G J. Comparing approximations to spatio-temporal models for epidemics with local spread[J]. Bulletin of Mathematical Biology, 2001, 63: 603-624.

[29] Thomson N A, Ellner S P. Pari-edge approximation for heterogeneous lattice population models[J]. Theoretical Population Biology, 2003, 64: 271-280.

[30] Trapman P. Reproduction numbers for epidemics on networks using pair approximation[J]. Mathematical Bioscience, 2007, 210 (2): 464-489.

[31] 靳禎, 孫桂全, 劉茂省. 網(wǎng)絡(luò)傳染病動力學(xué)建模分析[M]. 北京: 科學(xué)出版社, 2014.

[32] Sugimine N, Aihara K. Stability of an equilibrium state in a multi-infectious type SIS model on a truncated network[J]. Artificial Life Robotics, 2007,11: 157-161.

[33] Wu Q C, Fu X C, Yang M. Epidemic thresholds in a heterogenous population with competing strains[J]. Chinese Physics B, 2011, 20: 046401.

[34] Read J M, Eames K T D, Edmunds W J. Dynamic social networks and the implicatins for the spread of infectious disease[J]. Journal of the Royal Society Interface, 2008, 5(26): 1001-1007.

[35] Barrat A, Barthelemy M, Vespignani A. Weighted evolving networks: coupling topology and weight dynamics[J]. Physical Review Letters, 2004,92(22): 228701-228704.

[36] Boccaletti S, Latorab V, Morenod Y, et al. Complex networks: structure and dynamics[J]. Physics Reports, 2006, 424: 175-308.

[37] Newman M E J. Analysis of weighted networks[J]. Physical Review E, 2004,70(5): 056131

[38] Roshani F, Naimi Y. Effects of degree-biased transmission rate and nonlinear infectivity on rumor spreading in complex social networks[J]. Physical Review E, 2012, 85(3): 036109.

[39] Ma Z E, Li J. Dynamical modeling and analysis of epidemics[J].International Association of Geodesy Symposia, 2009, 106(B11):498.

[40] Xia C Y, Wang Z, Sanz J, et al. Effects of delayed recovery and nonuniform transmission on the spreading of diseases in complex networks[J]. Physica A, 2013, 392: 1577-1585.

[41] Liu Q M, Deng C S, Sun M C. The analysis of an epidemic model with time delay on scale-free networks[J]. Physica A, 2014, 410: 79-87.

[42] Heesterbeek J A P, Roberts M G. Threshold quantities for helminth infections[J]. Journal of Mathematical Biology, 1995, 33(4): 415-434.

[43] Wu Q, Fu X, Jin Z, et al. Influence of dynamic immunization on epidemic spreading in networks[J]. Physica A, 2015, 419:566-574.

[44] Chowell G, Nishiura H. Transmission dynamics and control of Ebola virus disease (EVD): a review[J]. BMC Medicine, 2014, 12(1): 196.