周 瑤,呂貴臣
(重慶理工大學 理學院,重慶 400054)
近年來,登革熱作為一種典型的媒介傳染病,經常在熱帶及亞熱帶國家肆虐,給人們帶來經濟上的損失和生活上的不便[1-2]。登革熱是由登革病毒引起的急性傳染病,主要通過伊蚊叮咬傳播給人類。感染登革病毒之后可導致患者出現不同的臨床表現,輕者表現為隱性感染,重者可能發生致命的登革出血熱(DHF)及登革休克綜合征(DSS)。隨著疾病流行地區的不斷擴大,登革熱已從過去的散發性疾病演變為全球性公共衛生問題,對人類健康造成嚴重威脅,在全球范圍內的發病率不斷增加[3]。如何更好地模擬登革熱病毒的傳播機理并提供有效控制措施去抑制病毒的傳播,已經引起了學者們關注。
最近幾十年來,利用數學模型研究登革熱病毒的傳播特征具有廣闊的應用前景[4]。1998年,Esteva等[5]建立了關于登革熱傳播過程的經典數學模型,隨后他們進一步改進了模型,分別考慮了人類種群變化[6]和病毒可在蚊蟲中垂直傳播[7]的數學模型。此后,基于他們的模型,研究者們對登革熱病毒傳播過程進行更為細致的刻畫,建立了一系列更復雜且切合實際的數學模型[8-10]。根據相關報道,雖然登革熱病毒主要通過媒介傳播到人類身上,但是沒有明顯證據可以排除病毒通過輸血或骨髓移植等方式由感染者直接傳給易感者的可能性[11]。Martcheva等[12]考慮上述因素后建立了一類包含媒介傳染病可在人類中不經媒介而直接進行傳播的數學模型。鑒于在非洲和東南亞的一些疾病較為流行的熱帶地區醫療條件不理想,在某些極端情況下,登革熱病毒可能直接發生有限的人際傳播。因此本文考慮一類人類種群規模變化、病毒可在人群中不經蚊蟲而發生人際傳播且疾病可致患者死亡的數學模型是合適的。
將研究對象聚焦于人和蚊蟲,并對Esteva等[2]的模型進行推廣。根據種群患病和易感狀態,將人類種群分為3類,蚊蟲分為2類。設SH、IH、RH分別為t時刻人類中易感者、感染者和移出者的數量;SV、IV分別表示t時刻易感蚊子和感染蚊子的數量。再令NH=SH+IH+RH和NV=SV+IV分別為t時刻人類和蚊蟲的總數量。考慮到蚊蟲的生命周期極短,其一旦染病,將終生不能恢復,故不考慮蚊子的移出者類。假設在單位時間內患病人類接觸易感人類和易感蚊子的數量有限,從而考慮病毒以標準發生率進行傳播;另外假設患病蚊子傳染易感人類的能力與環境內總人口數成正比,故考慮病毒以雙線性發生率進行傳播。基于以上假設和登革熱病毒的傳播特性,可知該病在各倉室之間的傳播流程見圖1。
圖1中,β1SHIV表示單位時間內易感人類被患病蚊子所感染的人數表示易感人類和患病人類通過輸血等方式傳播病毒,在單位時間內新增的患病人數為單位時間內易感蚊子叮咬患病人類之后受到感染的蚊蟲數量。另外,考慮到模型的實際生物學意義和登革熱病毒發生人際傳播的有限性,對相關參數作出如下假設:α>dH和dH≥β2。
基于傳播流程(圖1),可以寫出登革熱病毒傳播的微分方程模型如下:

圖1 登革熱病毒在人和蚊蟲之間的傳播流程框圖

模型中所涉及參數的具體意義見表1,其中所有參數都是非負的。

表1 模型(1)中相關參數
根據模型(1)前3個方程,可以得到人類種群滿足N′H=S′H+I′H+R′H=(α-μH)NH-dHIH。由后2個方程得到蚊蟲滿足N′V=S′V+I′V=A-μVNV。另外,容易證明。因此,考慮在的區域中對模型進行動力學分析。為便于后續分析,令

結合x+y=1,系統(1)經過上述尺度變換化為如下系統

結合尺度變換,定義系統(2)的可行域為Σ={(s,i,r,x,y)∈R5+|0≤s+i+r≤1,0≤x+y≤1}。
由于系統(2)中的第1、2、5個方程中均不含r和x,于是考慮分析降維后的系統動力學行為。降維后的3維系統為

考慮到模型的實際生物學意義,在Ω={(s,i,y)∈R3+|0≤y≤1,0≤s,0≤i,s+i≤1}中討論系統(3)的解的性態,易證明Ω是系統(3)的正向不變集。
根據基本再生數R0的下一代矩陣計算方法[13],模型(3)可寫為如下形式:

式中:h=(h1,h2)=(i,y);p=s;Fj(h,p)表示染病倉室j的新患病感染率;Vj(h,p)表示染病倉室j中剩余的轉移項(如出生率、死亡率、治愈率等)。此時

和

由Martcheva的計算公式[13],基本再生數與無病平衡點有關。易知,模型(4)總有唯一無病平衡點P0=(0,0,1)。根據文獻[13]的計算公式,可得

和

進一步,得到

因此

本節討論系統平衡點的存在性。
定理1系統(3)的無病平衡點P0=(1,0,0)始終存在。
證明:根據系統(3),無病平衡點必然滿足以下關系

代入i=0到上面等式中,得到s=1和y=0,因此無病平衡點P0始終存在。
定理2當R0>1時,系統(3)存在唯一的地方病平衡點Pe。
證明:由三維系統的第1、3個方程,地方病平衡點必然滿足以下關系


由于平衡點應被限制在可行域Ω內,滿足上面等式的地方病平衡點也應滿足下面的不等式:

為了分析g(i)的根的存在區間,根據αH的正負情況分別進行討論。
Case 1:αH>0(dH>β2)
前面有假設α>dH(α>αH=dH-β2),式(5)和(6)顯然成立;當0<i<1且時,式(7)成立;當0<i<1且時,式(8)成立。
進一步,為了得到滿足上面4個不等式關系的i,需要找到g(i)在區間(0,i1)上的根,其中


其中:U=β2α(μV+β3i),V=αβ3i-αHβ3i2,W=(λHβ3-αHμV)i。
顯然U、V、W均為正數。
接下來對i1可能存在的3種情況分別計算與討論:
1)當i1=1時





由于(α+γH+dH(1-i))(αμV+V+W)>λHαβ3+U,故
因此,由1)~3)可得,g(i1)<0。
此外,容易得到g(i)取得局部最大值時i為

其中

由上面不等式得到

因此g(i)的最大值在區間[0,i1]右側取得。由于g(i1)<0,當且僅當g(0)>0才能使得g(i)有唯一的1個根i*∈(0,i1)。由于

綜上,只有當R0>1時,g(i)有唯一的正根i*∈(0,i1),顯然原系統(2)此時也只有唯一的正根i∈(0,1),系統正平衡點的唯一性得證。
Case 2:αH=0(dH=β2)
同Case 1做法類似,代入αH=0到g(i)的表達式中,此時有

以及

經過簡單計算,可得到g(i3)=g(1)<0以及。進一步分析,顯然只有當g(0)=αμVM(R0-1)>0時,g(i)在區間(0,i3)范圍才會存在唯一的正根。即當R0>1時,系統存在唯一的正根i∈(0,1),系統正平衡點的唯一性得證。
定理1已經證明系統(3)始終存在1個無病平衡點P0(1,0,0),下面給出它的穩定性。
2.2.1 局部穩定性
定理3當R0<1時,系統(3)的無病平衡點P0局部漸近穩定。
證明:求出系統(3)在P0處的雅可比矩陣為

顯然系統(9)有1個特征值為λ1=-α,設其另2個特征值分別為λ2、λ3,且滿足

R0<1時,有,故

因此λ2+λ3<0顯然成立。此外,再結合,可得λ2·λ3>0。顯然,特征根λ2,λ3都具有負實部。綜上,P0的雅可比矩陣的特征值全部具有負實部,于是P0是局部漸近穩定的。
2.2.2 全局穩定性
定理4當R0<1時,系統(3)的無病平衡點P0全局漸近穩定。
證明:構造Lyapunov函數,令

對L沿著t求全導數,可化簡為

其中E、F、M分別為

顯然滿足E≤0,F<0。因此,當R0≤1時,。易證系統(3)在集合上的最大不變集為單點集 {P }0。根據LaSalle不變集原理,當R0≤1時,P0全局漸近穩定。
為證明地方病平衡點的穩定性,需要用到合作系統相關理論,下面給先出合作系統定義。
定義1考慮一般的n維系統:

其中f∶Rn→Rn是一個連續可微的映射,且有x=(x1,x2,…,xn)∈Rn,其過初值x(0)=x0的解記作x(t,x0)。如果其雅可比矩陣Df(x)的非對角線上元素均為非負,即對所有的i≠j(i,j=1,…,n),當x∈Rn+時,有。則系統(10)稱為合作系統。另外,若其雅可比矩陣Df(x)為不可約矩陣,則稱系統(10)為不可約系統。
根據合作系統定義,可得如下定理:
定理5系統(2)為Σ上的合作不可約系統。
證明:根據s+i+r=1和x+y=1可以得到y=1-x,i=1-s-r。將其代入系統(2)的第1個方程和第4個方程之后,得到系統(2)的雅可比矩陣為

其中主對角線元素用*表示。顯然,上述雅可比矩陣的非對角線元素均為非負,根據合作系統的定義,系統(2)為合作系統。容易證明其雅可比矩陣DF為不可約矩陣,因此系統(2)為不可約系統。
2.3.1 系統(3)的持久性
為證明系統地方病平衡點的全局穩定性,先證明系統的持久性。
定理6當R0>1時,系統(3)在Ω內一致持久,即存在c1、c2、c3>0使得系統(3)的任意解(s,i,y)滿足

證明:根據Freedman在文獻[14]中的定理4.3,令X=R3+,X1=Ω,X1是系統(3)的閉正向不變子集。為便于后續證明,記u(t)=(s,i,y)是系統(3)的解,,和M?={u(0)∈?X1|u(t)∈?X1},={(s,0,0)∈?X1|0≤s≤1}。顯然,系統在邊界上的最大閉不變子集N是單點集P0。欲證是相對開集,可以取X=R,Y=Ω,V={(s,i,y)|-1<s<1,i>0,y>0},從而,即證。
下用反證法證明

假設


易知(11)是(3)后2個方程的極限系統,則t→∞時,i(t)→0,y(t)→0。
由于R0>1,得到



因此

根據極限系統理論,這與t→∞時,s(t)→1,i(t)→0,y(t)→0矛盾,故系統(3)是一致持久的。
2.3.2 全局穩定性
定理7當R0>1時,系統(3)的地方病平衡點Pe全局漸近穩定。
證明:由定理6可知,系統(2)是一致持久的,又因為系統(2)合作不可約,根據Jiang[15]的定理C可知系統(2)的唯一地方病平衡點是全局漸近穩定的,從而得到系統(3)的唯一地方病平衡點是全局漸近穩定的。
運用數值模擬的方法對3維系統平衡點的動力學性質進行分析,模擬中參數取值:α=0.02,A=7 000,μH=0.000 023,dH=0.001 24,μV=0.014 5,γH=0.001 67,β1=0.000 012,β2=0.000 000 1,β3=0.000 534。經過計算得到地方病平衡點為pe=(0.106,0.822,0.029),閾值參數為R0=9.312 36>1,對地方病平衡點的全局穩定性進行模擬。然后分別取4組初值為(0.3,0.2,0.25),(0.35,0.25,0.28),(0.15,0.34,0.17),(0.23,0.17,0.38)并得到系統的解曲線如圖2所示。

圖2 R0>1時,地方病平衡點的全局穩定性曲線
模擬圖顯示當t→+∞時,系統的解曲線趨于地方病平衡點pe,從而對前面的理論證明給出了很好的數值驗證。
以登革熱病毒在人群和蚊蟲之間的傳播機理建立動力學模型,還考慮了人際傳播對模型的影響,得到模型的閾值并分析了平衡點的穩定性。分析表明該傳播模型的無病平衡點和地方病平衡點均表示全局漸近穩定。最后通過數值模擬進一步驗證了地方病平衡點的全局穩定性。