具有媒介的多種群炭疽傳播模型的動力學分析

韓夢潔, 劉俊利, 劉白茹

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

0 引 言

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

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

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

1 模型的建立

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

(1)

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

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

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

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

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

證明由微分方程基本理論可知，對任意的初值(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≤∞上存在唯一的非負解。

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

(2)

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

(3)

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

血食性蠅增長所滿足的方程：

(4)

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

注：在生物學文獻中滿足性質(A1)～(A3)的出生率函數(shù)B(FH)的例子有：

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

2 平衡點分析

平衡點是研究系統(tǒng)(1)動力學性質的基礎，本節(jié)將討論系統(tǒng)(1)的平衡點的存在性。

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

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

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

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

(5)

(6)

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

(7)

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

(8)

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

無病平衡點表示疾病的消亡，地方病平衡點表示疾病最終會長期存在而成為一種地方病。研究平衡點的穩(wěn)定性能夠為控制炭疽的傳播提供理論依據(jù)。本節(jié)將從理論上討論系統(tǒng)(1)的平衡點的穩(wěn)定性。

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

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

證明對于平衡點E3，它的特征方程為：

(9)

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

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

對于平衡點E*，當δa=0時，其雅克比矩陣的形式為：

其中

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

(10)

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

-M2的順序主子式為

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

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

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

(2)當正平衡點E*存在時，對于δa=0，若FH(0)>0成立，則E*在Ω內部是全局漸近穩(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時地方病平衡點的全局穩(wěn)定性，對于δa≥0的一般情況，正平衡點的全局穩(wěn)定性不易得到，下面的定理給出了R0>1時疾病的持久性。

定理5若R0>1，則存在ε>0，使得對于系統(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))滿足

由于定理5的證明方法與文獻[15]類似，因此這里不再贅述。

4 控制措施

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

4.1 彈性指數(shù)

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

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

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

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

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

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

4.2 目標再生數(shù)

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

(15)

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

(16)

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

5 結 論

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