郭少軍, 寇春海, 申明圓

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

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

郭少軍, 寇春海, 申明圓

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

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

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

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

(1)

定義該模型的初始條件為

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

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

y5(θ)=φ5(θ)

(2)

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

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

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

(3)

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

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

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

(4)

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

(5)

(1) 當(dāng)τ=0時, 有

(6)

又

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

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

(7)

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

(8)

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

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

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

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

(9)

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

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

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

(10)

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

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

(2) 當(dāng)τ>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)中, 對τ進(jìn)行微分可得

(16)

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

(17)

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

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

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

本節(jié)主要運(yùn)用規(guī)范形原理和中心流形原理[10], 給出式(1)模型在分岔點(diǎn)τ=τ0處的準(zhǔn)確公式, 這會決定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的有界變差矩陣函數(shù)η(θ, μ), θ∈[-τ, 0], 使得

(21)

選擇

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

(22)

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

(23)

且

(24)

方程式(18)可以改寫為

(25)

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

(26)

并定義雙線性內(nèi)積映射

(27)

(28)

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

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

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

定義

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)], 比較系數(shù)得

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

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

(39)

由式(35)和(38)可得

(40)

解方程組(40)得

(41)

(42)

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

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

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

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

(43)

由特征值和特征向量的關(guān)系知

(44)

由式(41)～(44)得

M41, M51)T

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

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

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

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

3 數(shù)值模擬

本節(jié)通過數(shù)值仿真來證明本文理論分析結(jié)果的正確性. 式(1)模型中參數(shù)的選取參見表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 無病平衡點(diǎn)E0處的數(shù)值模擬

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

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

3.2 地方病平衡點(diǎn)E2處的數(shù)值模擬

將各參數(shù)代入式(1)模型中, 據(jù)平衡點(diǎn)E2產(chǎn)生的條件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)會產(chǎn)生Hopf分岔現(xiàn)象. 通過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．由計算結(jié)果和定理2.1可知, 此分岔周期解是亞臨界的, 如圖2和3所示. 分岔周期解的周期隨著τ的增大而增大, 如圖4和5所示(圖4為圖2的局部放大圖).

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

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

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

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

3.3 不同抑制率α的數(shù)值模擬

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

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

4 結(jié) 語

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

[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]. 北京: 北京大學(xué)出版社, 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

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項基金資助項目(CUSF-DH-D-2015061)

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

O 175.6

A