韓夢潔, 劉俊利, 劉白茹

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

0 引 言

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

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

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

1 模型的建立

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

(1)

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

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

下面證明系統(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)

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

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

(4)

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

注：在生物學文獻中滿足性質(A1)～(A3)的出生率函數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)定性能夠為控制炭疽的傳播提供理論依據。本節(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 控制措施

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

4.1 彈性指數

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

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

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

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

表1 基本再生數R0的彈性指數值

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

4.2 目標再生數

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

(15)

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

(16)

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

5 結 論

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