具有媒介的多種群炭疽傳播模型的動(dòng)力學(xué)分析

韓夢(mèng)潔, 劉俊利, 劉白茹

(西安工程大學(xué) 理學(xué)院，西安 710048)

0 引 言

炭疽是一種人畜共患的急性傳染病[1]，能夠引起家畜和食草動(dòng)物的突然死亡，造成巨大的經(jīng)濟(jì)損失。世界上多個(gè)國(guó)家都曾爆發(fā)炭疽疫情，在中東、非洲等一些發(fā)展中國(guó)家尤為常見(jiàn)，此外，炭疽也能用于生物戰(zhàn)爭(zhēng)[2]。因此，研究炭疽的傳播動(dòng)態(tài)具有非常重要的意義。

食草動(dòng)物是炭疽最主要的感染源，如牛和馬等。在研究動(dòng)物種群中炭疽的傳播動(dòng)態(tài)時(shí)，一些學(xué)者提出了幾種數(shù)學(xué)模型。1983年，Hahn等以動(dòng)物吸入環(huán)境中的孢子病毒感染炭疽為背景，在動(dòng)物種群中建立了炭疽傳播的微分方程模型[3]。2013年，F(xiàn)riedman等提出了一個(gè)帶有局部擴(kuò)散的微分方程模型，模型既考慮了環(huán)境污染和動(dòng)物死亡，又考慮了動(dòng)物以染病尸體為食的情況[4]。在炭疽爆發(fā)期間，血食性蠅被認(rèn)為是一種潛在的疾病傳播媒介[5]。2013年，Baldacchino等研究了血食性蠅傳播炭疽的方式，指出血食性蠅能夠從嘔吐物和內(nèi)臟中分離出炭疽桿菌，進(jìn)而通過(guò)叮咬傳染易感動(dòng)物[6]。2017年，Mushayabasa等在動(dòng)物種群中建立炭疽傳播模型時(shí)加入了血食性蠅，研究了所建模型的理論結(jié)果與媒介對(duì)炭疽傳播和控制的影響[7]。隨著對(duì)炭疽疾病的深入研究，很多文章表明盡管炭疽是一種食草動(dòng)物疾病，但是接觸炭疽孢子病毒的人和動(dòng)物都有感染炭疽的風(fēng)險(xiǎn)。如1980年，McKendrick研究了當(dāng)血食性蠅存在時(shí)，食草動(dòng)物與人之間的傳播機(jī)制[8]。2010年，F(xiàn)asanella等研究了炭疽在不同物種之間的傳播方式，表明了在炭疽爆發(fā)期間，血食性蠅發(fā)揮著重要的作用[9]。

雖然已經(jīng)有很多文章研究了炭疽在動(dòng)物種群之間的傳播，但是多種群之間炭疽的傳播很少受到關(guān)注。為了尋找控制炭疽傳播的有效措施，本文根據(jù)食草動(dòng)物-血食性蠅-人三個(gè)種群之間的炭疽傳播機(jī)制，建立了一個(gè)具有媒介的多種群炭疽傳播模型。

1 模型的建立

根據(jù)炭疽在食草動(dòng)物-血食性蠅-人之間的傳播機(jī)理建立如下炭疽傳播模型。將模型分為9個(gè)倉(cāng)室，t時(shí)刻易感食草動(dòng)物的數(shù)量為Sa(t)，染病食草動(dòng)物的數(shù)量為Ia(t)，則動(dòng)物的總數(shù)量Na(t)=Sa(t)+Ia(t)。P(t)表示t時(shí)刻環(huán)境中孢子病毒的數(shù)量，C(t)表示t時(shí)刻感染炭疽死亡的動(dòng)物尸體的數(shù)量。假設(shè)t時(shí)刻易感血食性蠅的數(shù)量為SH(t)，染病血食性蠅的數(shù)量為IH(t)，則血食性蠅的總數(shù)量FH(t)=SH(t)+IH(t)。對(duì)于人，考慮易感者(Sh(t))、染病者(Ih(t))和恢復(fù)者(Rh(t))三類(lèi)，且總數(shù)量Nh(t)=Sh(t)+Ih(t)+Rh(t)。為了簡(jiǎn)便，使用FH(t)-IH(t)替代SH(t)，根據(jù)以上假設(shè)，建立如下模型：

(1)

式中：假設(shè)Λa為食草動(dòng)物的輸入率；Λh為易感人群的輸入率；bH(FH)為血食性蠅的出生率函數(shù)；di(i=a,H,h)為食草動(dòng)物、血食性蠅和人的自然死亡率；δa為染病動(dòng)物的因病死亡率；δh為感染炭疽的人的因病死亡率，ηC是染病動(dòng)物尸體釋放孢子病毒的速率；孢子病毒對(duì)易感動(dòng)物的傳染率為βPa。

由于染病動(dòng)物的糞便中攜帶孢子病毒，因此ηa表示染病動(dòng)物釋放孢子病毒的速率。假設(shè)孢子病毒的衰減率為k，染病動(dòng)物尸體的腐爛率為μ，則1/k為每克孢子病毒存活的平均時(shí)間，1/μ是染病尸體完全腐爛的平均時(shí)間。血食性蠅叮咬染病動(dòng)物后攜帶孢子病毒，之后在叮咬動(dòng)物的時(shí)候可以傳染食草動(dòng)物，則βaH表示染病動(dòng)物對(duì)血食性蠅的傳染率，βHa表示染病的血食性蠅對(duì)易感動(dòng)物的傳染率。人不僅能夠通過(guò)吸入孢子病毒感染炭疽，而且能夠通過(guò)接觸染病動(dòng)物皮毛等制品或食用染病動(dòng)物感染炭疽，則βPh、βah和βCh分別表示孢子病毒對(duì)易感人群的傳染率、染病動(dòng)物對(duì)易感人群的傳染率和染病動(dòng)物尸體對(duì)易感人群的傳染率。此外，大部分的炭疽感染者都能恢復(fù)[10]，則ρI表示染病者的恢復(fù)率，恢復(fù)者的免疫喪失率為ρR。假設(shè)δa≥0，其他參數(shù)為正常數(shù)。

下面證明系統(tǒng)(1)的非負(fù)有界性。

定理1對(duì)于系統(tǒng)(1)的任意初值(Sa(0),Ia(0),P(0),C(0),FH(0),IH(0),Sh(0),Ih(0),Rh(0))∈Ω，當(dāng)t≥0時(shí)系統(tǒng)(1)存在唯一的非負(fù)有界解，其中

此外，系統(tǒng)(1)的正向不變集為：

證明由微分方程基本理論可知，對(duì)任意的初值(Sa(0),Ia(0),P(0),C(0),FH(0),IH(0),Sh(0),Ih(0),Rh(0))∈Ω，系統(tǒng)(1)在其最大存在區(qū)間[0,T)，T≤∞上存在唯一的非負(fù)解。

將系統(tǒng)(1)中的Sa(t)和Ia(t)相加可得：

(2)

則對(duì)任意的t∈[0,T)，Na(t)有界。同理可得，對(duì)任意的t∈[0,T)，Nh(t)有界。由于

(3)

根據(jù)比較定理可知，P(t)和C(t)在t∈[0,T)上有界。

血食性蠅增長(zhǎng)所滿(mǎn)足的方程：

(4)

在文獻(xiàn)[11]中曾經(jīng)提到出生函數(shù)bH(FH)的假設(shè)如下：令血食性蠅的出生率函數(shù)bH(FH)=B(FH)FH，對(duì)FH∈(0,+∞)，假設(shè)B(FH)滿(mǎn)足以下3個(gè)性質(zhì)：

注：在生物學(xué)文獻(xiàn)中滿(mǎn)足性質(zhì)(A1)～(A3)的出生率函數(shù)B(FH)的例子有：

(ii)B2(FH)=be-aFH，a>0，b>dH

2 平衡點(diǎn)分析

平衡點(diǎn)是研究系統(tǒng)(1)動(dòng)力學(xué)性質(zhì)的基礎(chǔ)，本節(jié)將討論系統(tǒng)(1)的平衡點(diǎn)的存在性。

當(dāng)血食性蠅不存在時(shí)，當(dāng)且僅當(dāng)Rp>1時(shí)，系統(tǒng)(1)存在地方病平衡點(diǎn)E2=(Sa2,Ia2,P2,C2,0,0,Sh2,Ih2,Rh2)，其中

下面計(jì)算系統(tǒng)(1)的地方病平衡點(diǎn)。

記m1=βPaP+βHaIH，m2=βaHIa，m3=βPhP+βChC+βahIa，此時(shí)系統(tǒng)(1)的正平衡點(diǎn)的分量為：

把上面的表達(dá)式代入到m1和m2的表達(dá)式中，得：

(5)

(6)

將式(6)代入式(5)中，得：

(7)

下面僅考慮m1≠0的情況。式(7)兩邊同時(shí)除以m1，可得：

(8)

3 平衡點(diǎn)的穩(wěn)定性和疾病的持久性

無(wú)病平衡點(diǎn)表示疾病的消亡，地方病平衡點(diǎn)表示疾病最終會(huì)長(zhǎng)期存在而成為一種地方病。研究平衡點(diǎn)的穩(wěn)定性能夠?yàn)榭刂铺烤业膫鞑ヌ峁├碚撘罁?jù)。本節(jié)將從理論上討論系統(tǒng)(1)的平衡點(diǎn)的穩(wěn)定性。

(1)如果R0<1成立，則平衡點(diǎn)E3是局部漸近穩(wěn)定的；

(2)當(dāng)R0>1時(shí)，正平衡點(diǎn)E*存在。若δa=0成立，則平衡點(diǎn)E*是局部漸近穩(wěn)定的。

證明對(duì)于平衡點(diǎn)E3，它的特征方程為：

(9)

其中I4×4是四階單位矩陣，且

由于地方病平衡點(diǎn)E*的穩(wěn)定性不易證明，因此在證明其穩(wěn)定性時(shí)只考慮δa=0的情況。

對(duì)于平衡點(diǎn)E*，當(dāng)δa=0時(shí)，其雅克比矩陣的形式為：

其中

由分塊矩陣的性質(zhì)可知，J(E*)的特征值由J11和J22的特征值組成。J11的特征方程為:

(10)

其中I4×4是四階單位矩陣，且

-M2的順序主子式為

此時(shí)-M2的順序主子式均大于0，又因?yàn)?M2的主對(duì)角線(xiàn)元素均大于零，非主對(duì)角線(xiàn)元素非正，則-M2是M-矩陣。由文獻(xiàn)[13]可知，M2的所有特征根均具有負(fù)實(shí)部。由于方程(10)的另外兩個(gè)特征根均小于0，則矩陣J11的所有特征根均有負(fù)實(shí)部。容易證明-J22是M-矩陣。因此J22的所有特征根均具有負(fù)實(shí)部,從而矩陣J(E*)的所有特征根均具有負(fù)實(shí)部，則E*是局部漸近穩(wěn)定的。定理證畢。

局部穩(wěn)定性只能描述當(dāng)初始值位于平衡點(diǎn)附近時(shí)系統(tǒng)的動(dòng)力學(xué)行為，而全局穩(wěn)定性則可以將系統(tǒng)的初始值擴(kuò)大到更大范圍來(lái)研究解的長(zhǎng)期行為。下面將分析系統(tǒng)(1)的平衡點(diǎn)的全局穩(wěn)定性。

(1)當(dāng)R0<1時(shí)，如果FH(0)>0，則E3在Ω內(nèi)部是全局漸近穩(wěn)定的;

(2)當(dāng)正平衡點(diǎn)E*存在時(shí)，對(duì)于δa=0，若FH(0)>0成立，則E*在Ω內(nèi)部是全局漸近穩(wěn)定的。

(11)

考慮系統(tǒng)(11)的輔助系統(tǒng)：

(12)

(13)

f1(ρυ1,ρυ2,ρυ3,ρυ4)>ρf1(υ1,υ2,υ3,υ4)，f2(ρυ1,ρυ2,ρυ3,ρυ4)=ρf2(υ1,υ2,υ3,υ4)

f3(ρυ1,ρυ2,ρυ3,ρυ4)=ρf3(υ1,υ2,υ3,υ4)，f4(ρυ1,ρυ2,ρυ3,ρυ4)>ρf4(υ1,υ2,υ3,υ4)

(14)

定理4中討論了δa=0時(shí)地方病平衡點(diǎn)的全局穩(wěn)定性，對(duì)于δa≥0的一般情況，正平衡點(diǎn)的全局穩(wěn)定性不易得到，下面的定理給出了R0>1時(shí)疾病的持久性。

定理5若R0>1，則存在ε>0，使得對(duì)于系統(tǒng)(1)中具有初值條件Ia(0)>0，P(0)>0，C(0)>0，IH(0)>0，Ih(0)>0的解(Sa(t),Ia(t),P(t),C(t),FH(t),IH(t),Sh(t),Ih(t),Rh(t))滿(mǎn)足

由于定理5的證明方法與文獻(xiàn)[15]類(lèi)似，因此這里不再贅述。

4 控制措施

在傳染病建模中，依據(jù)所建模型的理論結(jié)果與實(shí)際生活中傳染病的病例數(shù)據(jù)研究疾病的控制措施是極其重要的。結(jié)合文獻(xiàn)[16]中給出的參數(shù)值，利用彈性指數(shù)和目標(biāo)再生數(shù)兩種方法尋找食草動(dòng)物-血食性蠅-人之間炭疽傳播的有效控制措施。

4.1 彈性指數(shù)

基本再生數(shù)R0在炭疽控制方面起著非常重要的作用。當(dāng)R0>1時(shí)炭疽流行成為一種地方病持續(xù)存在，R0<1時(shí)炭疽消亡，因此，通過(guò)研究R0尋找控制炭疽傳播的有效措施極其重要。在許多流行病學(xué)研究中，通常使用彈性指數(shù)研究基本再生數(shù)R0對(duì)參數(shù)的敏感性。此時(shí)，變量為R0，根據(jù)文獻(xiàn)[16]的計(jì)算方法可得彈性指數(shù)為：

顯然，上述所有的彈性指數(shù)都有一個(gè)確定的符號(hào)，即基本再生數(shù)R0隨著參數(shù)βPa、ηa、ηC、βaH和βHa的增加而增加；隨著參數(shù)k、μ、dH、da和δa的增加而降低。

為了從數(shù)值上確定系統(tǒng)(1)的彈性指數(shù)，下面以土耳其的人炭疽病例數(shù)為例，對(duì)系統(tǒng)(1)的相關(guān)參數(shù)值進(jìn)行估計(jì)。2004年土耳其接種疫苗的動(dòng)物數(shù)量為60 387[17]，感染炭疽的動(dòng)物數(shù)量為376[17]，其中337只動(dòng)物病死[17]，動(dòng)物出生率為1.216 7/年[16]，由于染病動(dòng)物與死亡動(dòng)物的數(shù)值均是在已接種動(dòng)物中監(jiān)測(cè)得到的，則Sa(0)=60 387，Ia(0)=376，C(0)=337，從而食草動(dòng)物的常數(shù)輸入為Λa=Sa(0)×1.216 7=73 472.862 9，染病動(dòng)物的因病死亡率為δa=C(0)/Ia(0)=0.896 3/年。假設(shè)血食性蠅出生率函數(shù)為Ricker型，即B(FH)=be-aFH。血食性蠅一生產(chǎn)卵量約500粒，生命周期約30天，則b=6 083.33/年，dH=12.166 7/年。據(jù)統(tǒng)計(jì)，土耳其2004年總?cè)丝跒? 115萬(wàn)人[18]，出生率為19.9‰[18]，人均壽命為70歲，則Λh=7.115×107×19.9‰=1 358 965，dh=0.014 3/年。由于染病者恢復(fù)后不會(huì)永久免疫，則ρR=1/年。由文獻(xiàn)[10]和文獻(xiàn)[18]知，參數(shù)δh=0.023 5/年，ρI=0.976 5/年，βph=0.038/年，βah=0.952/年，βCh=0.5/年。其余參數(shù)值見(jiàn)表1，根據(jù)世界衛(wèi)生組織發(fā)布的土耳其1996～2004年人炭疽病例數(shù)[18]進(jìn)行擬合，對(duì)參數(shù)ηa、βaH和βHa進(jìn)行估計(jì)，具體參數(shù)值如表1所示。

依據(jù)上述參數(shù)值，計(jì)算R0的所有參數(shù)的彈性指數(shù)，具體值如表1所示。

表1 基本再生數(shù)R0的彈性指數(shù)值

由表1可知，參數(shù)對(duì)炭疽的爆發(fā)或消亡有不同程度的影響，其中參數(shù)da、δa、dH、βaH和βHa對(duì)R0的影響最顯著，因此及時(shí)對(duì)動(dòng)物接種疫苗、降低食草動(dòng)物的患病率和因病死亡率、盡可能消除蒼蠅的繁殖地點(diǎn)、對(duì)蒼蠅使用殺蟲(chóng)劑、減少染病蒼蠅的數(shù)量是控制炭疽傳播最有效的措施。此外，雖然上述理論分析中沒(méi)有涉及人群對(duì)炭疽傳播的影響，但是對(duì)人群普及教育、不食用染病動(dòng)物等措施能夠降低易感人群感染炭疽的風(fēng)險(xiǎn)，對(duì)炭疽傳播也有一定的抑制作用。

4.2 目標(biāo)再生數(shù)

為了量化根除炭疽而需要控制的特定種群的百分比，依據(jù)文獻(xiàn)[19]描述的方法進(jìn)行計(jì)算。首先計(jì)算食草動(dòng)物進(jìn)行疫苗接種的目標(biāo)再生數(shù)。假設(shè)公共衛(wèi)生部門(mén)決定對(duì)食草動(dòng)物進(jìn)行疫苗接種，此時(shí)目標(biāo)集S={(1,1),(1,2),(1,3),(1,4)}，則目標(biāo)再生數(shù)為：

(15)

血食性蠅是傳播炭疽的一個(gè)重要媒介，在叮咬染病動(dòng)物之后攜帶孢子病毒，進(jìn)而傳染食草動(dòng)物，因此消滅血食性蠅的數(shù)量、降低血食性蠅的傳染率是控制炭疽的一個(gè)有效方法。此時(shí)目標(biāo)集S={(4,1)}，根據(jù)文獻(xiàn)[19]的計(jì)算方法可得目標(biāo)再生數(shù)為：

(16)

根據(jù)式(16)可得，如果Rp<1，則通過(guò)消滅1-1/T41比例的血食性蠅能夠根除炭疽。依據(jù)表1得T41=6.874 6，因此，為了使R0小于1，需要消滅1-1/T41=85.45%的血食性蠅。

5 結(jié) 論

本文研究了具有媒介的食草動(dòng)物和人的炭疽模型，分析了模型的動(dòng)力學(xué)行為。根據(jù)土耳其1996～2004年人炭疽病例數(shù)估計(jì)了模型參數(shù)值，并研究了基本再生數(shù)關(guān)于參數(shù)的敏感性，同時(shí)，依據(jù)目標(biāo)再生數(shù)將控制措施數(shù)值化。結(jié)果表明，對(duì)超過(guò)85.05%的動(dòng)物接種疫苗或消滅85.45%以上的血食性蠅都能夠根除炭疽。