具有人傳人的登革熱傳染病模型的動(dòng)力學(xué)分析

周 瑤，呂貴臣

（重慶理工大學(xué) 理學(xué)院，重慶 400054）

近年來(lái)，登革熱作為一種典型的媒介傳染病，經(jīng)常在熱帶及亞熱帶國(guó)家肆虐，給人們帶來(lái)經(jīng)濟(jì)上的損失和生活上的不便［1-2］。登革熱是由登革病毒引起的急性傳染病，主要通過(guò)伊蚊叮咬傳播給人類(lèi)。感染登革病毒之后可導(dǎo)致患者出現(xiàn)不同的臨床表現(xiàn)，輕者表現(xiàn)為隱性感染，重者可能發(fā)生致命的登革出血熱（DHF）及登革休克綜合征（DSS）。隨著疾病流行地區(qū)的不斷擴(kuò)大，登革熱已從過(guò)去的散發(fā)性疾病演變?yōu)槿蛐怨残l(wèi)生問(wèn)題，對(duì)人類(lèi)健康造成嚴(yán)重威脅，在全球范圍內(nèi)的發(fā)病率不斷增加［3］。如何更好地模擬登革熱病毒的傳播機(jī)理并提供有效控制措施去抑制病毒的傳播，已經(jīng)引起了學(xué)者們關(guān)注。

最近幾十年來(lái)，利用數(shù)學(xué)模型研究登革熱病毒的傳播特征具有廣闊的應(yīng)用前景［4］。1998年，Esteva等［5］建立了關(guān)于登革熱傳播過(guò)程的經(jīng)典數(shù)學(xué)模型，隨后他們進(jìn)一步改進(jìn)了模型，分別考慮了人類(lèi)種群變化［6］和病毒可在蚊蟲(chóng)中垂直傳播［7］的數(shù)學(xué)模型。此后，基于他們的模型，研究者們對(duì)登革熱病毒傳播過(guò)程進(jìn)行更為細(xì)致的刻畫(huà)，建立了一系列更復(fù)雜且切合實(shí)際的數(shù)學(xué)模型［8-10］。根據(jù)相關(guān)報(bào)道，雖然登革熱病毒主要通過(guò)媒介傳播到人類(lèi)身上，但是沒(méi)有明顯證據(jù)可以排除病毒通過(guò)輸血或骨髓移植等方式由感染者直接傳給易感者的可能性［11］。Martcheva等［12］考慮上述因素后建立了一類(lèi)包含媒介傳染病可在人類(lèi)中不經(jīng)媒介而直接進(jìn)行傳播的數(shù)學(xué)模型。鑒于在非洲和東南亞的一些疾病較為流行的熱帶地區(qū)醫(yī)療條件不理想，在某些極端情況下，登革熱病毒可能直接發(fā)生有限的人際傳播。因此本文考慮一類(lèi)人類(lèi)種群規(guī)模變化、病毒可在人群中不經(jīng)蚊蟲(chóng)而發(fā)生人際傳播且疾病可致患者死亡的數(shù)學(xué)模型是合適的。

1 模型建立

將研究對(duì)象聚焦于人和蚊蟲(chóng)，并對(duì)Esteva等［2］的模型進(jìn)行推廣。根據(jù)種群患病和易感狀態(tài)，將人類(lèi)種群分為3類(lèi)，蚊蟲(chóng)分為2類(lèi)。設(shè)SH、IH、RH分別為t時(shí)刻人類(lèi)中易感者、感染者和移出者的數(shù)量；SV、IV分別表示t時(shí)刻易感蚊子和感染蚊子的數(shù)量。再令NH=SH+IH+RH和NV=SV+IV分別為t時(shí)刻人類(lèi)和蚊蟲(chóng)的總數(shù)量?？紤]到蚊蟲(chóng)的生命周期極短，其一旦染病，將終生不能恢復(fù)，故不考慮蚊子的移出者類(lèi)。假設(shè)在單位時(shí)間內(nèi)患病人類(lèi)接觸易感人類(lèi)和易感蚊子的數(shù)量有限，從而考慮病毒以標(biāo)準(zhǔn)發(fā)生率進(jìn)行傳播；另外假設(shè)患病蚊子傳染易感人類(lèi)的能力與環(huán)境內(nèi)總?cè)丝跀?shù)成正比，故考慮病毒以雙線性發(fā)生率進(jìn)行傳播?；谝陨霞僭O(shè)和登革熱病毒的傳播特性，可知該病在各倉(cāng)室之間的傳播流程見(jiàn)圖1。

圖1中，β1SHIV表示單位時(shí)間內(nèi)易感人類(lèi)被患病蚊子所感染的人數(shù)表示易感人類(lèi)和患病人類(lèi)通過(guò)輸血等方式傳播病毒，在單位時(shí)間內(nèi)新增的患病人數(shù)為單位時(shí)間內(nèi)易感蚊子叮咬患病人類(lèi)之后受到感染的蚊蟲(chóng)數(shù)量。另外，考慮到模型的實(shí)際生物學(xué)意義和登革熱病毒發(fā)生人際傳播的有限性，對(duì)相關(guān)參數(shù)作出如下假設(shè)：α＞dH和dH≥β2。

基于傳播流程（圖1），可以寫(xiě)出登革熱病毒傳播的微分方程模型如下：

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

模型中所涉及參數(shù)的具體意義見(jiàn)表1，其中所有參數(shù)都是非負(fù)的。

表1 模型（1）中相關(guān)參數(shù)

根據(jù)模型（1）前3個(gè)方程，可以得到人類(lèi)種群滿足N′H=S′H+I′H+R′H=（α-μH）NH-dHIH。由后2個(gè)方程得到蚊蟲(chóng)滿足N′V=S′V+I′V=A-μVNV。另外，容易證明。因此，考慮在的區(qū)域中對(duì)模型進(jìn)行動(dòng)力學(xué)分析。為便于后續(xù)分析，令

結(jié)合x(chóng)+y=1，系統(tǒng)（1）經(jīng)過(guò)上述尺度變換化為如下系統(tǒng)

結(jié)合尺度變換，定義系統(tǒng)（2）的可行域?yàn)棣?｛（s，i，r，x，y）∈R5+|0≤s+i+r≤1，0≤x+y≤1｝。

2 平衡點(diǎn)的存在性和穩(wěn)定性

由于系統(tǒng)（2）中的第1、2、5個(gè)方程中均不含r和x，于是考慮分析降維后的系統(tǒng)動(dòng)力學(xué)行為。降維后的3維系統(tǒng)為

考慮到模型的實(shí)際生物學(xué)意義，在Ω=｛（s，i，y）∈R3+|0≤y≤1，0≤s，0≤i，s+i≤1｝中討論系統(tǒng)（3）的解的性態(tài)，易證明Ω是系統(tǒng)（3）的正向不變集。

根據(jù)基本再生數(shù)R0的下一代矩陣計(jì)算方法［13］，模型（3）可寫(xiě)為如下形式：

式中：h=（h1，h2）=（i，y）；p=s；Fj（h，p）表示染病倉(cāng)室j的新患病感染率；Vj（h，p）表示染病倉(cāng)室j中剩余的轉(zhuǎn)移項(xiàng)（如出生率、死亡率、治愈率等）。此時(shí)

和

由Martcheva的計(jì)算公式［13］，基本再生數(shù)與無(wú)病平衡點(diǎn)有關(guān)。易知，模型（4）總有唯一無(wú)病平衡點(diǎn)P0=（0，0，1）。根據(jù)文獻(xiàn)［13］的計(jì)算公式，可得

和

進(jìn)一步，得到

因此

2.1 平衡點(diǎn)的存在性

本節(jié)討論系統(tǒng)平衡點(diǎn)的存在性。

定理1系統(tǒng)（3）的無(wú)病平衡點(diǎn)P0=（1，0，0）始終存在。

證明：根據(jù)系統(tǒng)（3），無(wú)病平衡點(diǎn)必然滿足以下關(guān)系

代入i=0到上面等式中，得到s=1和y=0，因此無(wú)病平衡點(diǎn)P0始終存在。

定理2當(dāng)R0＞1時(shí)，系統(tǒng)（3）存在唯一的地方病平衡點(diǎn)Pe。

證明：由三維系統(tǒng)的第1、3個(gè)方程，地方病平衡點(diǎn)必然滿足以下關(guān)系

由于平衡點(diǎn)應(yīng)被限制在可行域Ω內(nèi)，滿足上面等式的地方病平衡點(diǎn)也應(yīng)滿足下面的不等式：

為了分析g（i）的根的存在區(qū)間，根據(jù)αH的正負(fù)情況分別進(jìn)行討論。

Case 1：αH＞0（dH＞β2）

前面有假設(shè)α＞dH（α＞αH=dH-β2），式（5）和（6）顯然成立；當(dāng)0＜i＜1且時(shí)，式（7）成立；當(dāng)0＜i＜1且時(shí)，式（8）成立。

進(jìn)一步，為了得到滿足上面4個(gè)不等式關(guān)系的i，需要找到g（i）在區(qū)間（0，i1）上的根，其中

其中：U=β2α（μV+β3i），V=αβ3i-αHβ3i2，W=（λHβ3-αHμV）i。

顯然U、V、W均為正數(shù)。

接下來(lái)對(duì)i1可能存在的3種情況分別計(jì)算與討論：

1）當(dāng)i1=1時(shí)

由于（α+γH+dH（1-i））（αμV+V+W）＞λHαβ3+U，故

因此，由1）～3）可得，g（i1）＜0。

此外，容易得到g（i）取得局部最大值時(shí)i為

其中

由上面不等式得到

因此g（i）的最大值在區(qū)間［0，i1］右側(cè)取得。由于g（i1）＜0，當(dāng)且僅當(dāng)g（0）＞0才能使得g（i）有唯一的1個(gè)根i*∈（0，i1）。由于

綜上，只有當(dāng)R0＞1時(shí)，g（i）有唯一的正根i*∈（0，i1），顯然原系統(tǒng)（2）此時(shí)也只有唯一的正根i∈（0，1），系統(tǒng)正平衡點(diǎn)的唯一性得證。

Case 2：αH=0（dH=β2）

同Case 1做法類(lèi)似，代入αH=0到g（i）的表達(dá)式中，此時(shí)有

以及

經(jīng)過(guò)簡(jiǎn)單計(jì)算，可得到g（i3）=g（1）＜0以及。進(jìn)一步分析，顯然只有當(dāng)g（0）=αμVM（R0-1）＞0時(shí)，g（i）在區(qū)間（0，i3）范圍才會(huì)存在唯一的正根。即當(dāng)R0＞1時(shí)，系統(tǒng)存在唯一的正根i∈（0，1），系統(tǒng)正平衡點(diǎn)的唯一性得證。

2.2 無(wú)病平衡點(diǎn)的穩(wěn)定性

定理1已經(jīng)證明系統(tǒng)（3）始終存在1個(gè)無(wú)病平衡點(diǎn)P0（1，0，0），下面給出它的穩(wěn)定性。

2.2.1 局部穩(wěn)定性

定理3當(dāng)R0＜1時(shí)，系統(tǒng)（3）的無(wú)病平衡點(diǎn)P0局部漸近穩(wěn)定。

證明：求出系統(tǒng)（3）在P0處的雅可比矩陣為

顯然系統(tǒng)（9）有1個(gè)特征值為λ1=-α，設(shè)其另2個(gè)特征值分別為λ2、λ3，且滿足

R0＜1時(shí)，有，故

因此λ2+λ3＜0顯然成立。此外，再結(jié)合，可得λ2·λ3＞0。顯然，特征根λ2，λ3都具有負(fù)實(shí)部。綜上，P0的雅可比矩陣的特征值全部具有負(fù)實(shí)部，于是P0是局部漸近穩(wěn)定的。

2.2.2 全局穩(wěn)定性

定理4當(dāng)R0＜1時(shí)，系統(tǒng)（3）的無(wú)病平衡點(diǎn)P0全局漸近穩(wěn)定。

證明：構(gòu)造Lyapunov函數(shù)，令

對(duì)L沿著t求全導(dǎo)數(shù)，可化簡(jiǎn)為

其中E、F、M分別為

顯然滿足E≤0，F(xiàn)＜0。因此，當(dāng)R0≤1時(shí)，。易證系統(tǒng)（3）在集合上的最大不變集為單點(diǎn)集 {P }0。根據(jù)LaSalle不變集原理，當(dāng)R0≤1時(shí)，P0全局漸近穩(wěn)定。

2.3 地方病平衡點(diǎn)的穩(wěn)定性

為證明地方病平衡點(diǎn)的穩(wěn)定性，需要用到合作系統(tǒng)相關(guān)理論，下面給先出合作系統(tǒng)定義。

定義1考慮一般的n維系統(tǒng)：

其中f∶Rn→Rn是一個(gè)連續(xù)可微的映射，且有x=（x1，x2，…，xn）∈Rn，其過(guò)初值x（0）=x0的解記作x（t，x0）。如果其雅可比矩陣Df（x）的非對(duì)角線上元素均為非負(fù)，即對(duì)所有的i≠j（i，j=1，…，n），當(dāng)x∈Rn+時(shí)，有。則系統(tǒng)（10）稱(chēng)為合作系統(tǒng)。另外，若其雅可比矩陣Df（x）為不可約矩陣，則稱(chēng)系統(tǒng)（10）為不可約系統(tǒng)。

根據(jù)合作系統(tǒng)定義，可得如下定理：

定理5系統(tǒng)（2）為Σ上的合作不可約系統(tǒng)。

證明：根據(jù)s+i+r=1和x+y=1可以得到y(tǒng)=1-x，i=1-s-r。將其代入系統(tǒng)（2）的第1個(gè)方程和第4個(gè)方程之后，得到系統(tǒng)（2）的雅可比矩陣為

其中主對(duì)角線元素用*表示。顯然，上述雅可比矩陣的非對(duì)角線元素均為非負(fù)，根據(jù)合作系統(tǒng)的定義，系統(tǒng)（2）為合作系統(tǒng)。容易證明其雅可比矩陣DF為不可約矩陣，因此系統(tǒng)（2）為不可約系統(tǒng)。

2.3.1 系統(tǒng)（3）的持久性

為證明系統(tǒng)地方病平衡點(diǎn)的全局穩(wěn)定性，先證明系統(tǒng)的持久性。

定理6當(dāng)R0＞1時(shí)，系統(tǒng)（3）在Ω內(nèi)一致持久，即存在c1、c2、c3＞0使得系統(tǒng)（3）的任意解（s，i，y）滿足

證明：根據(jù)Freedman在文獻(xiàn)［14］中的定理4.3，令X=R3+，X1=Ω，X1是系統(tǒng)（3）的閉正向不變子集。為便于后續(xù)證明，記u（t）=（s，i，y）是系統(tǒng)（3）的解，，和M?=｛u（0）∈?X1|u（t）∈?X1｝，=｛（s，0，0）∈?X1|0≤s≤1｝。顯然，系統(tǒng)在邊界上的最大閉不變子集N是單點(diǎn)集P0。欲證是相對(duì)開(kāi)集，可以取X=R，Y=Ω，V=｛（s，i，y）|-1＜s＜1，i＞0，y＞0｝，從而，即證。

下用反證法證明

假設(shè)

易知（11）是（3）后2個(gè)方程的極限系統(tǒng)，則t→∞時(shí)，i（t）→0，y（t）→0。

由于R0＞1，得到

因此

根據(jù)極限系統(tǒng)理論，這與t→∞時(shí)，s（t）→1，i（t）→0，y（t）→0矛盾，故系統(tǒng)（3）是一致持久的。

2.3.2 全局穩(wěn)定性

定理7當(dāng)R0＞1時(shí)，系統(tǒng)（3）的地方病平衡點(diǎn)Pe全局漸近穩(wěn)定。

證明：由定理6可知，系統(tǒng)（2）是一致持久的，又因?yàn)橄到y(tǒng)（2）合作不可約，根據(jù)Jiang［15］的定理C可知系統(tǒng)（2）的唯一地方病平衡點(diǎn)是全局漸近穩(wěn)定的，從而得到系統(tǒng)（3）的唯一地方病平衡點(diǎn)是全局漸近穩(wěn)定的。

3 數(shù)值模擬

運(yùn)用數(shù)值模擬的方法對(duì)3維系統(tǒng)平衡點(diǎn)的動(dòng)力學(xué)性質(zhì)進(jìn)行分析，模擬中參數(shù)取值：α=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。經(jīng)過(guò)計(jì)算得到地方病平衡點(diǎn)為pe=（0.106，0.822，0.029），閾值參數(shù)為R0=9.312 36＞1，對(duì)地方病平衡點(diǎn)的全局穩(wěn)定性進(jìn)行模擬。然后分別取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）并得到系統(tǒng)的解曲線如圖2所示。

圖2 R0＞1時(shí)，地方病平衡點(diǎn)的全局穩(wěn)定性曲線

模擬圖顯示當(dāng)t→+∞時(shí)，系統(tǒng)的解曲線趨于地方病平衡點(diǎn)pe，從而對(duì)前面的理論證明給出了很好的數(shù)值驗(yàn)證。

4 結(jié)論

以登革熱病毒在人群和蚊蟲(chóng)之間的傳播機(jī)理建立動(dòng)力學(xué)模型，還考慮了人際傳播對(duì)模型的影響，得到模型的閾值并分析了平衡點(diǎn)的穩(wěn)定性。分析表明該傳播模型的無(wú)病平衡點(diǎn)和地方病平衡點(diǎn)均表示全局漸近穩(wěn)定。最后通過(guò)數(shù)值模擬進(jìn)一步驗(yàn)證了地方病平衡點(diǎn)的全局穩(wěn)定性。