一類具有常數(shù)免疫輸入的病毒-免疫應(yīng)答動(dòng)力學(xué)模型分析

李亞男，王曉琴，李建全

(陜西科技大學(xué) 文理學(xué)院，陜西 西安 710021)

0 引言

眾所周知，乙肝病毒(HBV)、丙肝病毒(HCV)和艾滋病病毒(HIV)等在人體內(nèi)可以建立長(zhǎng)久的慢性感染，對(duì)人們的健康甚至生命造成了嚴(yán)重的威脅.病毒在人體內(nèi)的慢性感染程度，通常與個(gè)體的免疫應(yīng)答強(qiáng)度密切相關(guān).目前，在臨床上已基本了解人體免疫系統(tǒng)對(duì)大多數(shù)病毒的作用機(jī)理，但具體的作用方式仍未完全理解.于是，借助數(shù)學(xué)模型來(lái)刻畫免疫系統(tǒng)對(duì)病毒感染作用的具體形式，對(duì)認(rèn)識(shí)和了解以及控制病毒對(duì)人體的感染有著重要的現(xiàn)實(shí)意義，并且已取得一些顯著的成果[1,2].

人體的免疫形式分為體液免疫和細(xì)胞免疫兩類.體液免疫通常是以抗體對(duì)病毒的抑制作用為主，而細(xì)胞免疫則主要以免疫T細(xì)胞對(duì)作為抗原的被感染細(xì)胞進(jìn)行溶解等作用來(lái)控制病毒的產(chǎn)生[3].對(duì)此，用數(shù)學(xué)模型描述病毒的感染過(guò)程時(shí)，至少需要引入未感染細(xì)胞、被感染細(xì)胞和病毒三種變量[2,4,5].若再考慮免疫應(yīng)答的作用，就需要增加數(shù)學(xué)模型中的變量個(gè)數(shù)，這樣模型的結(jié)構(gòu)變得比較復(fù)雜，同時(shí)也增加了模型動(dòng)力學(xué)性態(tài)的分析難度.進(jìn)一步，若是以Holling II 或 III 等非線性函數(shù)來(lái)刻畫相互作用形式時(shí)，模型的動(dòng)力性質(zhì)往往難以得到較全面的理論分析結(jié)果[6,7].為了能清晰地研究病毒和免疫應(yīng)答的直接作用方式，一些不考慮宿主細(xì)胞，而僅考慮病毒和免疫應(yīng)答相互作用的簡(jiǎn)單數(shù)學(xué)模型已經(jīng)提出，這就為后續(xù)的數(shù)學(xué)分析提供了方便[8-12].

在文獻(xiàn)[10,11,12]中，對(duì)于免疫細(xì)胞的輸入，均假設(shè)僅因個(gè)體體內(nèi)存在的病毒刺激免疫系統(tǒng)而產(chǎn)生，這就意味著病毒不存在時(shí)，免疫系統(tǒng)的免疫力消失.同時(shí)以具有飽和形式Holling II函數(shù)描述病毒對(duì)免疫細(xì)胞的刺激產(chǎn)生率.但臨床檢測(cè)已經(jīng)發(fā)現(xiàn)，當(dāng)病毒不存在時(shí)，由免疫系統(tǒng)(如淋巴系統(tǒng))分泌的免疫細(xì)胞會(huì)處于一種動(dòng)態(tài)平衡狀態(tài).因此，本文將在已有文獻(xiàn)所述模型的基礎(chǔ)上，再將不依賴于病毒的免疫細(xì)胞輸入引入模型之中，進(jìn)而分析其動(dòng)力學(xué)性態(tài).

1 模型及其有界性

(1)

對(duì)于模型(1)有

定理1對(duì)t≥0，模型(1)在非負(fù)初始條件下的解(v(t),z(t))非負(fù)且最終一致有界.

證明：由模型(1)的第一個(gè)方程可得，當(dāng)v(0)=0時(shí)，對(duì)所有的t≥0，都有v(t)≡0.

由于v(t)≡0是模型(1)的解，所以根據(jù)自治系統(tǒng)解軌線的唯一性知，當(dāng)v(0)>0時(shí)，對(duì)所有的t≥0，都有v(t)>0.

綜上所述，當(dāng)t≥0時(shí)，模型(1)在非負(fù)初始條件下的解(v(t),z(t))是非負(fù)的.

(2)

因此模型(1)在非負(fù)初始條件下的解是最終有界的.

定理1證畢.

因此，只需在Γ上分析模型(1)的動(dòng)力學(xué)性態(tài).

2 平衡點(diǎn)的存在性和局部漸近穩(wěn)定性

模型(1)的平衡點(diǎn)P(v,z)坐標(biāo)滿足方程組：

(3)

當(dāng)v≠0時(shí)，由方程組(3)中的第一個(gè)方程有

(4)

為了保證v的正性，rd>pλ就是模型(1)正平衡點(diǎn)存在的必要條件.

進(jìn)一步，將式(4)代入方程組(3)中的第二個(gè)方程，并整理可得

(5)

綜上所述，關(guān)于模型(1)在正不變集Γ上的平衡點(diǎn)存在性有如下定理.

下面討論平衡點(diǎn)P0和P*的局部漸近穩(wěn)定性.

模型(1)在P0點(diǎn)的Jacobi矩陣為:

矩陣J(P0)的兩個(gè)特征根分別為

所以當(dāng)rdpλ時(shí)P0是不穩(wěn)定的.

當(dāng)rd=pλ時(shí)，矩陣J(P0)的兩個(gè)特征根為λ1=0和λ2<0.因此P0是模型(1)的高階平衡點(diǎn)，同時(shí)在P0附近滿足局部中心流形存在的條件.下面用中心流形來(lái)判斷P0在正不變集Γ上的穩(wěn)定性.

(6)

(7)

其中：

o(|u,y|3).

(8)

系統(tǒng)(7)變?yōu)?/p>

(9)

其中：

g2(w,x)=a1w2+a2wx+a3w3+a4w2x+

o(|w,x|3)

其中a1、a2、a3、a4分別為：

根據(jù)中心流形定理[13]，假設(shè)系統(tǒng)(9)在P0的局部中心流形為

x=h(w)=bw2+o(w2)

(10)

其中：b為待定常數(shù).

對(duì)式(10)兩邊沿著系統(tǒng)(9)的解關(guān)于t求導(dǎo)得

-dx+a1w2+a2wx+a3w3+a4w2x+o(|w,x|3)

(11)

再將x=h(w)=bw2+o(w2)代入式(11)得

(a1-bd)w2+(a2b+a3)w3+a4bw4+o(w4)

(12)

于是系統(tǒng)(9)在P0的局部中心流形為

(13)

進(jìn)一步將式(13)代入系統(tǒng)(9)中的第一方程，得中心流形上系統(tǒng)(9)的退化方程為

(14)

注意到在式(14)右端w的最低次數(shù)是2，且w2的系數(shù)小于零，并結(jié)合式(8)以及平移變換易知當(dāng)rd=pλ時(shí)，邊界平衡點(diǎn)P0在正不變集Γ上是局部漸近穩(wěn)定的.

模型(1)在正平衡點(diǎn)P*的Jacobi矩陣為

借助v*和z*滿足的方程組(3)可將矩陣J(P*)化簡(jiǎn)為

顯然,detJ(P*)>0，trJ(P*)<0.因此正平衡點(diǎn)P*只要存在就是局部漸近穩(wěn)定的.

于是關(guān)于模型(1)平衡點(diǎn)的局部穩(wěn)定性有如下定理.

定理3當(dāng)rd≤pλ時(shí)，模型(1)的邊界平衡點(diǎn)P0在正不變集Γ上是局部漸近穩(wěn)定的，當(dāng)rd>pλ時(shí)P0是不穩(wěn)定的；模型(1)的正平衡點(diǎn)P*只要存在，在正不變集Γ內(nèi)就是局部漸近穩(wěn)定的.

3 平衡點(diǎn)的全局漸近穩(wěn)定性

關(guān)于模型(1)平衡點(diǎn)的全局穩(wěn)定性有如下定理.

定理4當(dāng)rd≤pλ時(shí)，模型(1)的邊界平衡點(diǎn)P0在正不變集Γ上是全局漸近穩(wěn)定；正平衡點(diǎn)P*只要存在，在正不變集Γ內(nèi)就是全局漸近穩(wěn)定的.

(15)

定義Lyapunov函數(shù)

則L2沿著系統(tǒng)(15)解的全導(dǎo)數(shù)為

對(duì)于模型(1)的正平衡點(diǎn)P*(v*,z*)，由方程組(3)有

于是模型(1)可表示為

(16)

F(v)=(v*+ε)ln(v+ε)-v*lnv

易知，v=v*是函數(shù)F(v)在(0,∞)上的唯一極值點(diǎn)，并且F(v)在v=v*處取到最小值.同樣，z=z*是函數(shù)G(z)在(0,∞)上的唯一極值點(diǎn)，并且G(z)在z=z*處取到最小值.

直接計(jì)算可得，函數(shù)L3沿著系統(tǒng)(16)解的全導(dǎo)數(shù)為

定理4證畢.

4 數(shù)值分析與討論

前述定理2、3、4，從平衡點(diǎn)的存在性、局部穩(wěn)定性和全局穩(wěn)定性在理論上說(shuō)明了模型(1)的全局動(dòng)力學(xué)性態(tài).為進(jìn)一步從數(shù)值上進(jìn)行驗(yàn)證，下面將對(duì)rdpλ三種情形予以數(shù)值模擬以顯示理論結(jié)果的正確性.

由圖1(a)、(b)、(c)可以看出，對(duì)于情形rdpλ，模型的可行平衡點(diǎn)有單重?zé)o感染平衡點(diǎn)P0(0,4)和感染平衡點(diǎn)P*(0.1,5.7).這里P0是不穩(wěn)定的，P*在可行域內(nèi)是全局漸近穩(wěn)定的.

(a)情形rd

(b)情形rd=pλ時(shí)的軌線圖(k=2、p=1、λ=2、c=2.5、ε=1.5、d=0.5、r=4)

(c)情形rd>pλ時(shí)的軌線圖(k=2、p=1、λ=2、c=2.5、ε=1.5、d=0.5、r=6)圖1 rdpλ三種情形的軌線圖

下面進(jìn)一步考慮在免疫作用下，感染持續(xù)并處于平衡狀態(tài)時(shí)病毒載量對(duì)各免疫因素p、λ和c的依賴性.

當(dāng)免疫系統(tǒng)起作用時(shí)，若模型(1)的正平衡點(diǎn)存在，則感染平衡點(diǎn)P*(v*,z*)中病毒的分量v*滿足方程f(v)=0.此方程也可寫為

(17)

(a)v*對(duì)p的依賴(r=6,k=2,λ=2,c=2.5,ε=1.5,d=0.5)

(b)v*對(duì)λ的依賴(r=6,k=2,p=1,c=2.5,ε=1.5,d=0.5)

(c)v*對(duì)c的依賴(r=6,k=2,p=1,λ=2,ε=1.5,d=0.5) 圖2 病毒載量平衡值對(duì)參數(shù)p、λ和c的依賴性

應(yīng)當(dāng)注意，根據(jù)定理2、3、4，決定感染能否滅絕的閾值條件不依賴于參數(shù)c，這就意味著增大免疫細(xì)胞的刺激產(chǎn)生率只會(huì)降低病毒載量，不會(huì)使得感染滅絕.相反，增大免疫細(xì)胞對(duì)病毒的殺傷率系數(shù)或免疫細(xì)胞的輸入率都可清除感染.

5 結(jié)論

本文考慮了一類具有常數(shù)免疫輸入的病毒-免疫應(yīng)答動(dòng)力學(xué)模型，通過(guò)對(duì)模型進(jìn)行分析和研究得到模型的邊界平衡點(diǎn)和正平衡點(diǎn)，并對(duì)模型的邊界平衡點(diǎn)和正平衡點(diǎn)的局部穩(wěn)定性和全局穩(wěn)定性進(jìn)行了詳細(xì)的分析和研究，發(fā)現(xiàn)免疫存在時(shí)的病毒內(nèi)稟增長(zhǎng)率決定模型的動(dòng)力學(xué)性態(tài).最后針對(duì)各免疫因素對(duì)病毒持續(xù)的影響進(jìn)行了討論.