高文哲,雒志學

(蘭州交通大學 數理學院, 蘭州 730070)

0 引言

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

1 模型建立

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

(1)

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

令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 基本再生數與平衡點的存在性

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

因此

則定義基本再生數表達式：

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

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

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

(4)

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

(5)

式中：

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

當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判別法證明無病平衡點和邊界平衡點的局部漸近穩(wěn)定性，系統(tǒng)(3)在任意點E(S,V1,V2,I1,I2)處的Jacobi矩陣為：

(6)

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

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

特征方程為：

顯然

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

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

(7)

顯然

特征根p3,p4,p5由如下多項式所確定：

p3+a1p2+a2p+a3=0

式中：

顯然

式中：

利用Hurwitz判別法可知當R1<1

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

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

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

(8)

當t→∞時，I1(t)→0。

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

沿著系統(tǒng)(3)的軌跡對函數V(t)進行求導：

由算數平均數與幾何平均數知

當R2<1時，

式中：

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

由算數平均數與幾何平均數知

由于邊界平衡點E1滿足方程組(4)，經計算可知

5 數值模擬

下面將給出2個數值模擬。

取系統(tǒng)(3)中的參數Λ=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。經計算R1=0.942 5,R2=0.787 5，由定理4知無病平衡E0是全局漸近穩(wěn)定的，其數值模擬如圖1所示。

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

取系統(tǒng)(3)中的參數Λ=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。經計算R1=0.757 6,R2=2.581 8，由定理5知邊界平衡點E1是全局漸近穩(wěn)定的，其數值模擬如圖2所示。

6 結論