高文哲,雒志學(xué)

(蘭州交通大學(xué) 數(shù)理學(xué)院, 蘭州 730070)

0 引言

長期以來，傳染病都是危害人類健康的最大敵人，人們在努力戰(zhàn)勝它的過程中，已取得了顯著的成果[1-5]。然而隨著疾病的發(fā)展，病毒在傳播過程中可能會發(fā)生變異，從而導(dǎo)致疾病失控。梁桂珍等[6]、渠亞嬌等[7]和李冬梅等[8]研究了具有病毒變異的傳染病模型，討論了疾病流行的閾值和系統(tǒng)各個平衡點(diǎn)的穩(wěn)定性；Cai等[9-10]對具有接種的病毒變異傳染病模型進(jìn)行了穩(wěn)定性分析，但只考慮了接種個體對其中某一階段的病毒完全有效的情況；Baba等[11-12]考慮了對2種病毒分別有阻礙作用的雙接種傳染病問題，卻并未考慮感染這2種病毒的患者之間會發(fā)生轉(zhuǎn)換的情況。因此在病毒變異傳染病問題的研究基礎(chǔ)上，考慮對易感人群進(jìn)行2種疫苗接種的情況，提出一類具有雙接種的病毒變異傳染病模型。

1 模型建立

本節(jié)考慮疫苗接種對預(yù)防傳染病發(fā)生的重要作用，引入一類具有雙接種的病毒變異傳染病模型。假設(shè)第一種疫苗接種者對變異前病毒完全免疫而對變異后病毒有部分抵抗作用，第二種疫苗接種者對變異后病毒完全免疫而對變異前病毒有部分抵抗作用，2類感染者都具有傳染性，且病毒變異前該疾病不足以致命，而病毒變異后該疾病足以致命?；谝陨霞僭O(shè)，建立如下模型：

(1)

式中：S(t),V1(t),V2(t),I1(t),I2(t),R(t)分別表示t時刻易感染者，第一種疫苗接種者，第二種疫苗接種者，病毒變異前感染者，病毒變異后感染者和恢復(fù)者的數(shù)量；Λ表示種群的輸入率；β1,β2分別表示病毒變異前后感染者對易感人群的傳染系數(shù)；d表示種群的自然死亡率；φ1,φ2表示第一種疫苗與第二種疫苗的接種率；k2表示病毒變異前感染者對第二種疫苗接種者的傳染率；k1表示病毒變異后感染者對第一種疫苗接種者的傳染率；γ1和γ2分別表示病毒變異前后感染者的恢復(fù)率；ε表示病毒變異前感染者轉(zhuǎn)換為病毒變異后感染者的速率；δ表示病毒變異后感染者的因病死亡率。根據(jù)模型的生物意義，假設(shè)所有參數(shù)均為正數(shù)。

令N(t)=S(t)+V1(t)+V2(t)+I1(t)+I2(t)+R(t)，由系統(tǒng)(1)可得

由常微分方程的比較原理知

所以

(2)

因系統(tǒng)(1)的前5個方程中不含變量R，所以只需要考慮如下系統(tǒng)：

(3)

式中：λ=d+φ1+φ2,α1=d+γ1+ε,α2=d+γ2+δ，記S=S(t),V1=V1(t),V2=V2(t),I1=I1(t)，I2=I2(t)。

2 基本再生數(shù)與平衡點(diǎn)的存在性

F和V在E0處的Jacobi矩陣為：

因此

則定義基本再生數(shù)表達(dá)式：

R0=ρ(FV-1)=max{R1,R2}

以下將對系統(tǒng)(3)各平衡點(diǎn)的存在性進(jìn)行討論。

證明系統(tǒng)(3)的平衡點(diǎn)滿足如下方程：

(4)

當(dāng)I1=0,I2≠0時，由方程組(4)的前3個式子可得

(5)

式中：

a1=α2k1β2,a2=α2dβ2+α2λk1-β2k1Λ,a3=α2dλ-k1φ1Λ-β2dΛ

當(dāng)I1I2≠0時，由方程組(4)的前3個式子可得

式中：a1=α1β1k2,b1=α1β2k2,d1=α1β1d+α1k2λ-β2k2Λ,f1=α1β2d,g1=α1λd-β1Λd-k2φ2Λ,a2=-α2β2k1,b2=β1εk1-α2β1k1,c2=β1εk1,d2=β1dε,e2=β2Λk1-α2β2d-α2λk1,f2=ελk1+β2dε-α2β1d,g2=εdλ,h=β2Λd+k1φ1Λ-α2dλ。

3 局部穩(wěn)定性分析

下面通過Jacobi矩陣與Hurwitz判別法證明無病平衡點(diǎn)和邊界平衡點(diǎn)的局部漸近穩(wěn)定性，系統(tǒng)(3)在任意點(diǎn)E(S,V1,V2,I1,I2)處的Jacobi矩陣為：

(6)

式中：M=-β1I1-β2I2-λ,N=β1S+k2V2-α1,Q=β2S+k1V1-α2。

證明由式(6)知，系統(tǒng)(3)在點(diǎn)E0處的Jacobi矩陣為：

特征方程為：

顯然

故當(dāng)R1<1且R2<1時，特征根均為負(fù)實(shí)數(shù)，即當(dāng)R0=max{R1,R2}<1時，無病平衡點(diǎn)E0局部漸近穩(wěn)定。

證明系統(tǒng)(3)在點(diǎn)E1處的特征方程為：

(7)

顯然

特征根p3,p4,p5由如下多項(xiàng)式所確定：

p3+a1p2+a2p+a3=0

式中：

顯然

式中：

利用Hurwitz判別法可知當(dāng)R1<1

4 全局穩(wěn)定性分析

本節(jié)將通過構(gòu)造Lyapunov函數(shù)并借助LaSalle不變集原理對系統(tǒng)(3)的2個平衡點(diǎn)E0和E1的全局漸近穩(wěn)定性進(jìn)行分析。

證明由于R0=max{R1,R2}<1意味著R1<1，所以由方程組(3)的4個等式知

(8)

當(dāng)t→∞時，I1(t)→0。

故在I1(t)=0的平面上，構(gòu)造如下Lyapunov函數(shù)：

沿著系統(tǒng)(3)的軌跡對函數(shù)V(t)進(jìn)行求導(dǎo)：

由算數(shù)平均數(shù)與幾何平均數(shù)知

當(dāng)R2<1時，

式中：

對函數(shù)V(t)沿著系統(tǒng)(3)的軌跡求導(dǎo)：

由算數(shù)平均數(shù)與幾何平均數(shù)知

由于邊界平衡點(diǎn)E1滿足方程組(4)，經(jīng)計算可知

5 數(shù)值模擬

下面將給出2個數(shù)值模擬。

取系統(tǒng)(3)中的參數(shù)Λ=1,d=0.1,β1=0.15,β2=0.18,γ1=0.48,γ2=0.45,φ1=0.14,φ2=0.3,ε=0.54,δ=0.4,k1=0.16,k2=0.14,S(0)=3,V1(0)=2,V2(0)=1,I1(0)=2,I2(0)=2。經(jīng)計算R1=0.942 5,R2=0.787 5，由定理4知無病平衡E0是全局漸近穩(wěn)定的，其數(shù)值模擬如圖1所示。

圖1 無病平衡點(diǎn)E0的全局漸近穩(wěn)定性

取系統(tǒng)(3)中的參數(shù)Λ=1,d=0.1,β1=0.15,β2=0.22,γ1=0.35,γ2=0.25,φ1=0.5,φ2=0.4,ε=0.54,δ=0.2,k1=0.24,k2=0.15,S(0)=3,V1(0)=2,V2(0)=1,I1(0)=2,I2(0)=2。經(jīng)計算R1=0.757 6,R2=2.581 8，由定理5知邊界平衡點(diǎn)E1是全局漸近穩(wěn)定的，其數(shù)值模擬如圖2所示。

6 結(jié)論