一類具有時滯的HIV-1型傳染病模型的穩(wěn)定性與Hopf分岔分析

郭少軍, 寇春海, 申明圓

(東華大學 理學院, 上海 201620)

一類具有時滯的HIV-1型傳染病模型的穩(wěn)定性與Hopf分岔分析

郭少軍, 寇春海, 申明圓

(東華大學 理學院, 上海 201620)

利用特征值理論分析了無病平衡點和地方病平衡點的局部漸近穩(wěn)定性. 同時, 由于時滯τ的出現, Hopf分岔行為在系統(tǒng)中產生, 應用中心流形定理和規(guī)范形理論, 給出系統(tǒng)的分岔方向及分岔周期解穩(wěn)定性計算公式. 最后, 通過數值模擬驗證了理論分析結果.

傳染病模型; 時滯; Hopf分岔; 穩(wěn)定性

近年來, 在研究和分析HIV(艾滋病病毒)感染的動力學性態(tài)過程中, 數學模型的重要作用已被諸多研究所證明[1-3]. CD4+T細胞感染HIV模型的研究有了若干研究成果, 2014年, 文獻[4]建立了一個五維系統(tǒng)并討論了該系統(tǒng)的穩(wěn)定性. 同時, 數學家對免疫時滯的病毒動力學模型也有了一些研究[5-9]. 由于從健康T細胞和病毒結合到能釋放出病毒之間需要一段時間τ, 因此, 本文將在文獻[4]研究的基礎上考慮具有Holling II型感染率且?guī)в袝r滯的HIV模型, 模型見式(1)所示. 式(1)中各參量的意義及取值見表1所示.

(1)

定義該模型的初始條件為

y1(θ)=φ1(θ), y2(θ)=φ2(θ),

y3(θ)=φ3(θ), y4(θ)=φ4(θ),

y5(θ)=φ5(θ)

(2)

表1 各參數的含義和模擬中的取值[4]Table 1 The meaning and the simulation values of the parameters

在生物學的意義下, 假定φi(θ)≥0, θ∈[-τ, 0], i=1, 2, 3, 4, 5, 且式(1)中的所有系數均為正. 不難證明, 具有式(2)所示初始條件的式(1)模型的解總存在且為正.

1 平衡點的穩(wěn)定性和Hopf分岔分析

(3)

1.1 無病平衡點E0的穩(wěn)定性分析

定理1.1.1 當R0<1時, 對任意τ≥0, E0是局部漸近穩(wěn)定的; 當R0>1時, E0是不穩(wěn)定的．

證明: 式(1)模型在E0處的特征方程為

(4)

顯然, ξ1, 2, 3=-d, -b, -h是方程式(4)的3個單重負根．考慮下面的超越方程

(5)

(1) 當τ=0時, 有

(6)

又

因此, 當τ=0時, 方程式(6)有兩個負根, 即E0局部漸近穩(wěn)定.

(2) 若τ>0時, 式(5)為超越方程, 設ξ=±iω為式(5)的純虛根, ω>0, 將ξ代入式(5)得

(7)

分離實部和虛部并消去τ, 且令ω2=x, 有

(8)

1.2 地方病平衡點E2的穩(wěn)定性分析

式(1)模型在平衡點E2處的特征方程為

ξ5+p1ξ4+p2ξ3+p3ξ2+p4ξ+p5+

e-ξ τ(q1ξ3+q2ξ2+q3ξ+q4)=0

(9)

(1) 當τ=0時, 式(9)變形為

ξ5+p1ξ4+(p2+q1)ξ3+(p3+q2)ξ2+

(p4+q3)ξ+p5+q4=0

(10)

由Routh-Hurwiz判據可知, 方程(10)的根均有負實部的充要條件為

其中: a1=p1, a2=p2+q1, a3=p3+q2, a4=p4+q3, a5=p5+q4.

(2) 當τ>0時, 將ξ=iω(ω>0)代入式(9)中, 分離實部和虛部得

(11)

即

(12)

從式(12)中消去τ, 得

ω10+J1ω8+J2ω6+J3ω4+J4ω2+J5=0

(13)

m5+J1m4+J2m3+J3m2+J4m+J5=0

(14)

(15)

在方程式(9)中, 對τ進行微分可得

(16)

將ξ=iω0代入式(16), 有

(17)

上面的分析可以得到下面定理1.2.1的結果.

定理1.2.1 當τ=τ0時, 式(1)模型在地方病平衡點E2處產生Hopf分岔.

2 Hopf分岔的方向及其穩(wěn)定性

本節(jié)主要運用規(guī)范形原理和中心流形原理[10], 給出式(1)模型在分岔點τ=τ0處的準確公式, 這會決定Hopf分岔的方向、穩(wěn)定性和周期解的周期性.

(18)

(19)

這里u(t)=(u1(t), u2(t), u3(t), u4(t), u5(t))T∈C([-τ, 0], R5), Lμ: C→R5定義為: 對φ∈C([-τ, 0], R5),

Lμφ=Q1φ(0)+Q2φ(-τ)

(20)

其中:

且

由Riesz表示引理[11]知, 存在一個5×5的有界變差矩陣函數η(θ, μ), θ∈[-τ, 0], 使得

(21)

選擇

η(θ, μ)=Q1δ(θ)+Q2δ(θ+τ)

(22)

其中: δ(θ)是狄拉克函數. 因為φ∈C([-τ, 0], R5), 定義

(23)

且

(24)

方程式(18)可以改寫為

(25)

對于ψ∈C([-τ, 0], R5),A的伴隨算子A*定義為

(26)

并定義雙線性內積映射

(27)

(28)

由于A*的特征向量q*的表達式過長, 這里從略.

為了保證〈q*, q〉=1, 由式(27)可得

采用與文獻[10]中同樣的符號, 首先計算出中心流形C0在μ=0處的坐標.

定義

ut(θ)-2Re{Z(t)q(t)}

(29)

在中心流形C0上有

(30)

這里

(31)

(32)

(33)

(34)

(35)

(36)

另外, 在C0上有

(37)

(38)

ut(θ)=(u1(t+θ), u2(t+θ), u3(t+θ),

u4(t+θ), u5(t+θ))

式(34)兩邊同乘以[k16+k15φ1(-τ)][k16+k15φ1(0)], 比較系數得

為了得到g21, 需要計算當θ∈[-τ, 0)時W20(θ)和W11(θ)的表達式. 從式(35)可知

將上式與式(36)比較系數得

(39)

由式(35)和(38)可得

(40)

解方程組(40)得

(41)

(42)

將方程式(36)兩邊同乘以[k16+k15φ1(-τ)][k16+k15φ1(0)], 與式(42)比較系數可得

2(M11, M21, M31, M41, M51)T

(M13, M23, M33, M43, M53)T

由A的定義和式(35)和(38), 得

(43)

由特征值和特征向量的關系知

(44)

由式(41)～(44)得

M41, M51)T

從式(41)計算可得g21, 經過簡單的計算可得:

定理2.1 (1) Hopf分岔的方向可以由μ2確定, 如果μ2<0(μ2>0), 則Hopf分岔是亞臨界(超臨界)的; 并且當τ<τ0(τ>τ0)時, 分岔周期解存在.

(2) β2確定分岔周期解的穩(wěn)定性. 若β2>0(β2<0), 那么分岔周期解是軌道不穩(wěn)定(穩(wěn)定)的.

(3) T2決定分岔周期解的周期性．若T2>0(T2<0), 那么分岔周期解的周期增大(減小).

3 數值模擬

本節(jié)通過數值仿真來證明本文理論分析結果的正確性. 式(1)模型中參數的選取參見表1: y1=0.5, y2=1.5, y3=0.1, y4=1.5, y5=1.5, d=0.1, p=1, β=0.007 5, a=1.5, u=8, c=0.1, k=25, q=0.1, b=0.2, h=0.1, α=0.031 25.

3.1 無病平衡點E0處的數值模擬

將各參數代入式(1)模型中, 據定理1.1.1的結論, 要求R0<1, 此時得到λ<6.400 0, 選取λ=λ1=4且τ=4進行數值模擬, 經計算可得E0=(40, 0, 0, 0, 0), 此時系統(tǒng)是漸近穩(wěn)定的, 如圖1所示.

圖1 τ=4, λ=4, 式(1)模型中各分量趨于平衡點的情況Fig.1 Each amount converge to their equilibrium with the parametic values τ=4, λ=4 in model (1)

3.2 地方病平衡點E2處的數值模擬

將各參數代入式(1)模型中, 據平衡點E2產生的條件R0>R1, 計算可得λ>7.511 1, 選取λ=λ2=10, 計算可得E2=(87.785 7, 2.222 2, 6.944 4, 0.223 4, 0.049 6). 此時由于時滯τ的影響, 系統(tǒng)會產生Hopf分岔現象. 通過Matlab計算可得: ω0≈0.145 5, τ0≈7.256 1, C1(0)≈0.379 2-6.485 7i, β2=0.758 4, μ2≈-5.712 8, T2≈3.657 1．由計算結果和定理2.1可知, 此分岔周期解是亞臨界的, 如圖2和3所示. 分岔周期解的周期隨著τ的增大而增大, 如圖4和5所示(圖4為圖2的局部放大圖).

圖2 τ=6, λ=10, 式(1)模型中各分量趨于平衡點的情況Fig.2 Each amount converge to their equilibrium with the parametic values τ=6, λ=10 in model (1)

圖3 τ=8, λ=10, 式(1)模型中各分量趨于平衡點的情況Fig.3 Each amount converge to their equilibrium with the parametic values τ=8, λ=10 in model (1)

圖4 τ=6, λ=10, 式(1)模型中各分量趨于平衡點的情況Fig.4 Each amount converge to their equilibrium with the parametic values τ=6, λ=10 in model (1)

圖5 τ=1, λ=10, 式(1)模型中各分量趨于平衡點的情況Fig.5 Each amount converge to their equilibrium with the parametic values τ=1, λ=10 in model (1)

3.3 不同抑制率α的數值模擬

仍取文獻[4]中的各參數值, 其中λ=10, τ=8. 抑制率α分別取0.02, 0.03, 0.05時, 顯示抑制率與已受HIV病毒感染的細胞濃度之間的關系, 如圖4所示.

圖6 λ=10, τ=8, 已受HIV病毒感染的細胞濃度在不同抑制率影響因素下的變化情況Fig.6 The concentration of infected cells by HIV under different inhibition rates with the parametic values λ=10, τ=8

4 結 語

本文研究了一類具有時滯且伴有非線性發(fā)生率的HIV-1型傳染病模型， 分析了無病平衡點E0的局部漸近穩(wěn)定性和不穩(wěn)定性， 同時對地方病平衡點E2也進行了穩(wěn)定性分析. 在R0>R1的條件下, 當τ=0時, 得到了地方病平衡點E2的局部漸近穩(wěn)定的充分條件; 當τ>0且τ=τ0時, Hopf分岔現象出現在了地方病平衡點E2處, 即疾病會出現周期性爆發(fā).本文運用了中心流形定理和規(guī)范形理論，對式(1)模型所出現的Hopf分岔方向和分岔周期解的穩(wěn)定性以及分岔的周期進行了研究分析. 最終得出, 在式(1)模型中參數滿足一定條件時, 式(1)模型出現了亞臨界的Hopf分岔現象; 隨著時滯量τ的增加, 分岔周期解的周期也在增加. 理論分析的正確性被進一步的數值模擬結果所證明. 為了驗證抑制率α對疾病的影響效果, 通過一定的數值模擬得出, 適當地增大抑制率α, 會使得已受HIV感染的細胞濃度的最小值明顯下降, 這對疾病的治療時機的切入是有利的, 因此式(1)模型中引入抑制率α是必要的, 同時對生物學的研究具有重要意義. 但是, 在平衡點E1處進行平衡性分析過程中出現了異常復雜的數據, 這有待進一步研究論證.

[1] PERELSON A S, NEUMANN A U, MAEKOWITZ M. HIV-1 dynamics in divo: Virion clearance rate, infected cell life-span and viral generation time [J]. Science, 1996, 271(5255): 15821586.

[2] HO D D, NEUMANN A U, PERELSON A S. Rapid turnover of plasma virions and CD4+lymphocytes in HIV-1 infection[J]. Nature, 1995, 373: 123126.

[3] PERELSON A S, NELSON P W. Mathematical analysis of HIV-1 dynamicsin vivo[J]. Society for Industrial and Applied Mathematics, 1999, 41: 344.

[4] YU P, HUANG J N, JIANG J. Dynamics of an HIV-1 infection model with cell mediated immunity[J]. Communications in Nonlinear Science and Numerical Simulation, 2014, 19: 38273844.

[5] NOWAK M A, LLOYD A L, VADQUEZ G M, et al. Viral dynamics of primary viremia and antitroviral therapy in simian immunode ficiency virus infection[J]. Journal of Virology, 1997, 71(10): 75187525.

[6] SHI X, ZHOU X, SONG X. Dynamical behavior of a delay virus dynamics model with CTL immune response[J]. Nonlinear Analysis, 2010, 11: 17951809.

[7] XIE Q, HUANG D, ZHANG S. Analysis of a viral infection model with delayed immune response[J]. Applied Mathematical Modelling, 2010, 34: 23882395.

[8] WANG W D. Global behavior of an SEIRS epidemic model with time delays[J]. Applied Mathematics Letters, 2002, 15: 423428.

[9] KADDAR A. Stability analysis in a delayed SIR epidemic model with a saturated incidence rate[J]. Nonlinear Analysis Model Control, 2010, 15: 299306.

[10] HASSARD B D, KAZAARINOFF N D, WAN Y H. Theory and applications of Hopf bifurcation[M]. London: Cambridge University Press, 1981.

[11] 張恭慶, 郭懋正. 泛函分析[M]. 北京: 北京大學出版社, 1990: 3839.

Stability and Hopf Bifurcation Analysis of a Delay HIV-1 Epidemic Model

GUOShao-jun,KOUChun-hai,SHENMing-yuan

(College of Science, Donghua University, Shanghai 201620, China)

By analyzing the corresponding characteristic equation, the local stability of disease-free equilibrium and endemic equilibrium is discussed. The bifurcation property is obtained as the delay pass through a critical value. Applying the center manifold and normal form theory, some local bifurcation results are obtained and the formulas for determining the bifurcation direction and stability of the bifurcation periodic solution are derived. Finally, numerical simulations are presented to illustrate the theoretical analysis.

epidemic model; delay; Hopf bifurcation; stability

16710444 (2016)06-0906-10

2015-07-27

中央高校基本科研業(yè)務費專項基金資助項目(CUSF-DH-D-2015061)

郭少軍(1988—),男, 山西保德人, 碩士研究生, 研究方向為常微分方程. E-mail：15021986290@163.com 寇春海(聯(lián)系人),男,教授,E-mail：kouchunhai@dhu.edu.cn

O 175.6

A