具有飽和發(fā)生率的兩毒株時(shí)滯HIV 感染模型的閾值動(dòng)力學(xué)?

陳偉，張龍

(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，新疆 烏魯木齊 830017)

0 引言

HIV,人類免疫缺陷病毒,一直是人類社會(huì)生命健康的主要威脅之一.因此,生物數(shù)學(xué)學(xué)者越來越重視對(duì)病毒動(dòng)力學(xué)的研究[1?2],并且多種HIV 感染模型對(duì)理解病毒感染機(jī)理方面發(fā)揮著重要作用[3?5].在1996 年,Perelson等構(gòu)造了一種模型,其已被廣泛應(yīng)用于模擬HIV 感染患者的血漿病毒載量[6].但是在實(shí)際治療HIV 感染患者的過程中,Ribeiro 和Bonhoeffer 發(fā)現(xiàn)在HAART (高效抗逆轉(zhuǎn)錄病毒療法)中出現(xiàn)了一種不良現(xiàn)象,即在開始治療之前,HIV 毒株的單一突變(或多突變組合)也可能導(dǎo)致耐藥性[7].為了更深入研究耐藥型毒株的演變,Rong 等提出了一類包含敏感型毒株和耐藥型毒株的病毒動(dòng)力學(xué)模型[8]:

其中:T(t),Ts(t),Vs(t),Tr(t) 和Vr(t) 分別是易感細(xì)胞,敏感型感染細(xì)胞,敏感型毒株,耐藥型感染細(xì)胞,耐藥型毒株在時(shí)刻t的濃度.λ和d分別為易感細(xì)胞的補(bǔ)充率和死亡率,ks和kr分別為敏感型感染細(xì)胞和耐藥型感染細(xì)胞的感染率,Ns和Nr分別是敏感型毒株和耐藥型毒株的釋放率,δ和c分別為感染細(xì)胞的死亡率和游離病毒的清除率.參數(shù)u(0 ≤u<1) 是Ts到Tr的轉(zhuǎn)化率,也就是在逆轉(zhuǎn)錄過程中,敏感型毒株發(fā)生突變,變成耐藥型毒株,使得感染細(xì)胞轉(zhuǎn)化的比例(簡(jiǎn)稱SR 轉(zhuǎn)換).目前,已有學(xué)者研究了一些關(guān)于模型(1) 全局動(dòng)力學(xué)包括平衡點(diǎn)全局穩(wěn)定性和模型一致持續(xù)性的有趣結(jié)論[9].

在易感細(xì)胞和游離病毒間的非線性發(fā)生率函數(shù)有多種具體形式,例如,雙線性發(fā)生率,飽和發(fā)生率,B?D發(fā)生率等,值得注意的是,B?D發(fā)生率函數(shù)可以更客觀地描述易感細(xì)胞和游離病毒的動(dòng)態(tài)演化過程[10].B?D發(fā)生率函數(shù)可表示為,(a,b>0),特別的,當(dāng)a=0,b>0 時(shí),該發(fā)生率函數(shù)可退化為飽和發(fā)生率函數(shù).因此Chen 等擴(kuò)展了上述模型,并且建立了全局動(dòng)力學(xué)的閾值準(zhǔn)則[11].時(shí)滯在許多病毒感染模型中得到了廣泛的應(yīng)用,對(duì)研究實(shí)際生物過程起決定性作用[12?14],很容易注意到,相比于離散時(shí)滯,分布潛伏時(shí)滯更適用于實(shí)際情況.受此啟發(fā),建立了一類具有飽和發(fā)生率和分布潛伏時(shí)滯的敏感型和耐藥型毒株HIV 感染模型:

在本文中,ω1,α1,τs,τr都是非負(fù)常數(shù),且均為有限值.當(dāng)ω1=0 或α1=0 時(shí),對(duì)Vs或Vr相應(yīng)的發(fā)生率退化成雙線性發(fā)生率.這里假設(shè)ω1+α1>0,以及hi(η) (i=s,r) 在[0,τi] 上是非負(fù)連續(xù)函數(shù),并且同樣定義

本文首先根據(jù)泛函微分方程理論方法建立了模型解的正性和有界性,得到了平衡點(diǎn)的存在性,進(jìn)一步利用下代矩陣算法定義了基本再生數(shù).其次對(duì)具有SR 轉(zhuǎn)換的模型(2),分別通過線性化方法和構(gòu)造Lyapunov 泛函方法,建立并證明了無病平衡點(diǎn)及耐藥型感染平衡點(diǎn)局部和全局穩(wěn)定性的閾值準(zhǔn)則.然后運(yùn)用動(dòng)力系統(tǒng)持續(xù)性理論,研究了模型(2) 所有正解的一致持續(xù)性.最后,在數(shù)值模擬中,三維相圖驗(yàn)證了滿足地方病平衡點(diǎn)存在的條件,其也是全局漸近穩(wěn)定.

1 預(yù)備知識(shí)

定義τ=max{τs,τr} 以及={(x1,x2,x3,x4,x5) :xi≥0,i=1,2,3,4,5}.令C=C([?τ,0],R5) 表示從[?τ,0] 到R5的全體連續(xù)函數(shù)所構(gòu)成的Banach 空間,并且對(duì)任意?∈C,有‖?‖=max?τ≤θ≤0|?(θ)|.在這里,?(θ)=(?1(θ),?2(θ),?3(θ),?4(θ),?5(θ)) 并且定義是C的正錐.則模型(2) 的初始條件可表示為

這里θ∈[?τ,0] 以及?=(?1,?2,?3,?4,?5)∈C+.

首先,關(guān)于模型(2) 解的非負(fù)性和有界性,有如下結(jié)論.

定理1任意初始函數(shù)?=(?1,?2,?3,?4,?5)∈C+,模型(2) 具有初始條件(3) 的解(T(t),Ts(t),Vs(t),Tr(t),Vr(t)) 在t∈[0,∞) 上有定義,非負(fù)且最終有界.

證明根據(jù)泛函微分方程基本理論[15],容易得知模型(2) 有唯一的解(T(t),Ts(t),Vs(t),Tr(t),Vr(t)) 滿足初始條件(3),并且在t∈[0,t0) 上有定義,其中t0<∞或t0=∞.為了得到解(T(t),Ts(t),Vs(t),Tr(t),Vr(t)) 的非負(fù)性,根據(jù)解對(duì)初始函數(shù)的連續(xù)依賴性,僅需要證明當(dāng)?i(0)>0 (i=1,2,3,4,5) 時(shí),對(duì)所有t∈[0,t0) 有(T(t),Ts(t),Vs(t),Tr(t),Vr(t))>0.由模型(2) 的第一個(gè)方程有:

2 模型（2）的穩(wěn)定性

運(yùn)用線性化方法和構(gòu)造適當(dāng)?shù)睦钛牌罩Z夫泛函方法,考慮平衡點(diǎn)的局部和全局穩(wěn)定性,定義函數(shù)g(x)=x?1?lnx.令是(2) 的任意平衡點(diǎn),則可得如下線性系統(tǒng)

首先,關(guān)于無病平衡點(diǎn)E0的穩(wěn)定性,有如下定理.

定理3(i) 如果(1?u)Rs<1 且Rr<1,則無病平衡點(diǎn)E0局部漸近穩(wěn)定;

(ii) 如果Rs≤1 且Rr≤1,則E0全局漸近穩(wěn)定;

(iii) 如果(1?u)Rs>1 或Rr>1,則E0不穩(wěn)定.

證明由(7) 可知在平衡點(diǎn)E0處的特征方程為

顯然(8) 式的一個(gè)根為λ0=?d<0.剩余的根滿足下列方程

矛盾.因此如果(1?u)Rs<1 且Rr<1,則方程(8) 的所有根均具有負(fù)實(shí)部.故無病平衡點(diǎn)E0局部漸近穩(wěn)定.

如果(1?u)Rs>1 或Rr>1,直接計(jì)算有f1(0)=δc(1?(1?u)Rs)<0,f1(+∞)=+∞或f2(0)=δc(1?Rr)<0,f2(+∞)=+∞.因此,特征方程(8) 至少有一個(gè)根具有正實(shí)部.故E0不穩(wěn)定.

為了證明E0的全局穩(wěn)定性,定義如下李雅普諾夫函數(shù)H1(t)

注2在定理4 中,我們僅得到了當(dāng)Rr>max{1,Rs+(Rs?1)} 時(shí),Er全局漸近穩(wěn)定.然而,結(jié)合定理4 的結(jié)論(i),一個(gè)有趣的問題是:當(dāng)Rr>max{1,(1?u)Rs+((1?u)Rs?1)} 時(shí),Er是否全局漸近穩(wěn)定.

注3很遺憾沒能建立模型(2)正平衡點(diǎn)Ec局部和全局穩(wěn)定性的閾值條件.原因是很難分析在Ec處的特征方程,且很難構(gòu)造適當(dāng)?shù)睦钛牌罩Z夫函數(shù).故在下節(jié)中,當(dāng)正平衡點(diǎn)Ec存在時(shí),建立了模型(2) 的一致持續(xù)性.

3 一致持續(xù)性

定義

令ω(?) 是解u(t,?) 的ω-極限集.則從以下兩種情形給出證明.

情形(1):(1 ?u)Rs>1 ≥Rr.由定理2 可知,模型(2) 僅有兩個(gè)平衡點(diǎn)E0和Ec.令M0={E0},顯然有對(duì)任意的?∈M?,有u(t,?)∈?X,故當(dāng)t≥0 時(shí),Ts(t,?)≡0 或Vs(t,?)≡0 或Tr(t,?)≡0 或Vr(t,?)≡0.如果Ts(t,?)≡0,則由模型(2) 的第二個(gè)方程可知Vs(t,?)≡0.由此,模型(2) 可退化為下列形式

若?4+?5=0,則由模型(16) 可以得到Tr(t,?)=Vr(t,?)≡0,故(16) 進(jìn)一步退化為

這表明模型(2) 同樣是一致持續(xù)的.

注4一個(gè)有趣的問題是:當(dāng)滿足定理5 的條件時(shí),正平衡點(diǎn)Ec是否全局漸近穩(wěn)定.

4 數(shù)值模擬

本節(jié)給出兩個(gè)數(shù)值算例,驗(yàn)證注4 中提出的問題.對(duì)所有η∈[?τi,0],取核函數(shù)hi(η)=e?0.1η,(i=s,r).在模型(2) 中選取參數(shù)d=0.1,u=0.6,δ=c=1,其余參數(shù)將在具體算例中給出.

例1在模型(2)中,選取初始函數(shù)為(T(θ),Ts(θ),Vs(θ),Tr(θ),Vr(θ))=(0.5+0.2sin(θ)+0.005k,0.06+0.1sin(θ)+0.005k,0.02+0.1sin(θ)+0.004k,0.1+0.01sin(θ)+0.000 2k,0.01+0.04sin(θ)+0.01k),k=1,2,···,20,θ∈[?3,0],并且參數(shù)λ=0.6,Ns=0.4,Nr=0.05,τs=τr=3,ω1=3.5,α1=4.5.通過計(jì)算可得Rs≈6.220 3,(1?u)Rs≈2.488 1 和Rr≈0.777 5.因此,(1?u)Rs>1>Rr,由圖1 和圖2 可知正平衡點(diǎn)Ec=(2.747 2,0.099 5,0.039 8,0.225 8,0.011 3)全局漸近穩(wěn)定,這表明注4 中的結(jié)論正確.

例2在模型(2)中,選取初始函數(shù)為(T(θ),Ts(θ),Vs(θ),Tr(θ),Vr(θ))=(0.8+0.4sin(θ)+0.01k,0.05+0.01sin(θ)+0.02k,0.01+0.01sin(θ)+0.02k,0.3+0.01sin(θ)+0.01k,0.03+0.04sin(θ)+0.02k),k=1,2,···,20,θ∈[?1.5,0],并且參數(shù)λ=0.8,Ns=0.9,Nr=0.4,τs=τr=1.5,ω1=1.5,α1=3.通過計(jì)算可得到Rs≈10.028 9,Rr≈4.457 3 和(1?u)Rs+((1?u)Rs?1)≈6.902 7.于是,(1?u)Rs+((1?u)Rs?1),由圖3 和圖4 可知正平衡點(diǎn)Ec=(2.204 2,0.078,0.070 2,0.501 6,0.200 6) 全局漸近穩(wěn)定,這表明注4 中的結(jié)論正確.

5 結(jié)論

本文研究了一類具有飽和發(fā)生率的兩毒株時(shí)滯HIV 感染模型的閾值動(dòng)力學(xué),主要討論了飽和發(fā)生率和分布潛伏時(shí)滯對(duì)HIV 感染模型閾值動(dòng)力學(xué)的影響.通過構(gòu)造李雅普諾夫泛函和運(yùn)用動(dòng)力系統(tǒng)持續(xù)性理論,得到了一些關(guān)于平衡點(diǎn)全局穩(wěn)定性和模型一致持續(xù)性全局動(dòng)力學(xué)的有趣結(jié)論.今后的工作將集中在具有反應(yīng)擴(kuò)散HIV感染模型的分支分析和敏感性分析,可以更準(zhǔn)確地描述傳染病傳播的實(shí)際情況.