周 瑤，呂貴臣

（重慶理工大學 理學院，重慶 400054）

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

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

1 模型建立

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

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

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

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

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

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

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

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

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

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

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

考慮到模型的實際生物學意義，在Ω=｛（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的下一代矩陣計算方法［13］，模型（3）可寫為如下形式：

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

和

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

和

進一步，得到

因此

2.1 平衡點的存在性

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

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

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

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

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

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

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

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

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

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

進一步，為了得到滿足上面4個不等式關(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ù)。

接下來對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）的最大值在區(qū)間［0，i1］右側(cè)取得。由于g（i1）＜0，當且僅當g（0）＞0才能使得g（i）有唯一的1個根i*∈（0，i1）。由于

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

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

同Case 1做法類似，代入αH=0到g（i）的表達式中，此時有

以及

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

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

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

2.2.1 局部穩(wěn)定性

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

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

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

R0＜1時，有，故

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

2.2.2 全局穩(wěn)定性

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

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

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

其中E、F、M分別為

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

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

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

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

其中f∶Rn→Rn是一個連續(xù)可微的映射，且有x=（x1，x2，…，xn）∈Rn，其過初值x（0）=x0的解記作x（t，x0）。如果其雅可比矩陣Df（x）的非對角線上元素均為非負，即對所有的i≠j（i，j=1，…，n），當x∈Rn+時，有。則系統(tǒng)（10）稱為合作系統(tǒng)。另外，若其雅可比矩陣Df（x）為不可約矩陣，則稱系統(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個方程和第4個方程之后，得到系統(tǒng)（2）的雅可比矩陣為

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

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

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

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

證明：根據(jù)Freedman在文獻［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是單點集P0。欲證是相對開集，可以取X=R，Y=Ω，V=｛（s，i，y）|-1＜s＜1，i＞0，y＞0｝，從而，即證。

下用反證法證明

假設(shè)

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

由于R0＞1，得到

因此

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

2.3.2 全局穩(wěn)定性

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

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

3 數(shù)值模擬

運用數(shù)值模擬的方法對3維系統(tǒng)平衡點的動力學性質(zhì)進行分析，模擬中參數(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)過計算得到地方病平衡點為pe=（0.106，0.822，0.029），閾值參數(shù)為R0=9.312 36＞1，對地方病平衡點的全局穩(wě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時，地方病平衡點的全局穩(wěn)定性曲線

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

4 結(jié)論

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