具有雙接種的病毒變異傳染病模型穩(wěn)定性分析

高文哲,雒志學(xué)

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

0 引言

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

1 模型建立

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

(1)

式中：S(t),V1(t),V2(t),I1(t),I2(t),R(t)分別表示t時(shí)刻易感染者，第一種疫苗接種者，第二種疫苗接種者，病毒變異前感染者，病毒變異后感染者和恢復(fù)者的數(shù)量；Λ表示種群的輸入率；β1,β2分別表示病毒變異前后感染者對(duì)易感人群的傳染系數(shù)；d表示種群的自然死亡率；φ1,φ2表示第一種疫苗與第二種疫苗的接種率；k2表示病毒變異前感染者對(duì)第二種疫苗接種者的傳染率；k1表示病毒變異后感染者對(duì)第一種疫苗接種者的傳染率；γ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個(gè)方程中不含變量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}

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

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

(4)

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

(5)

式中：

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

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

式中：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)定性分析

下面通過(guò)Jacobi矩陣與Hurwitz判別法證明無(wú)病平衡點(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時(shí)，特征根均為負(fù)實(shí)數(shù)，即當(dāng)R0=max{R1,R2}<1時(shí)，無(wú)病平衡點(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é)將通過(guò)構(gòu)造Lyapunov函數(shù)并借助LaSalle不變集原理對(duì)系統(tǒng)(3)的2個(gè)平衡點(diǎn)E0和E1的全局漸近穩(wěn)定性進(jìn)行分析。

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

(8)

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

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

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

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

當(dāng)R2<1時(shí)，

式中：

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

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

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

5 數(shù)值模擬

下面將給出2個(gè)數(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)計(jì)算R1=0.942 5,R2=0.787 5，由定理4知無(wú)病平衡E0是全局漸近穩(wěn)定的，其數(shù)值模擬如圖1所示。

圖1 無(wú)病平衡點(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)計(jì)算R1=0.757 6,R2=2.581 8，由定理5知邊界平衡點(diǎn)E1是全局漸近穩(wěn)定的，其數(shù)值模擬如圖2所示。

6 結(jié)論