李淑萍，劉花芳

(中北大學(xué) 理學(xué)院，太原 030051)

數(shù)千年以來，結(jié)核病對人類的健康造成了重大的危害，從近年來全球與中國關(guān)鍵數(shù)據(jù)可以看出[1-2]，我國的結(jié)核病情況依然很嚴(yán)重。為了研究結(jié)核病的傳播，人們通過考慮不同的因素建立了結(jié)核病傳播的各種動力學(xué)模型：

Bowong等[3]和Huo等[4]應(yīng)用一般接觸率分別考慮了簡單的SEI模型及含有耐藥菌株的SE1E2I1I2的結(jié)核病模型；考慮到治療對結(jié)核病的作用，宋妮等[5]、Das等[6]、Saputra等[7]和Feng等[8]研究了帶治療(恢復(fù))的SEIT (SEI R)模型；隨著卡介苗的出現(xiàn)，Gao等[9]和李春等[10]探究了帶免疫的SVEIT模型，但輸入全部進(jìn)入了易感者倉室。考慮到輸入人群中除易感人群外，還有已經(jīng)接種疫苗的人群，本文考慮了將按比例輸入與不完全免疫及治療相結(jié)合的結(jié)核病傳播動力學(xué)模型[11]。對于該模型，運(yùn)用穩(wěn)定性理論[12]，對其無病平衡點(diǎn)和正平衡點(diǎn)的局部穩(wěn)定性與全局穩(wěn)定性進(jìn)行分析，并應(yīng)用最優(yōu)控制理論得到最優(yōu)控制解的存在唯一性及其表達(dá)形式；其次，通過數(shù)值模擬驗(yàn)證理論結(jié)果并對參數(shù)的敏感性進(jìn)行分析；最后，總結(jié)本文主要的研究工作。

1 模型的建立

將總?cè)丝贜(t)分為5 類，分別是：易感者S(t)，疫苗接種者V(t)，潛伏者E(t)，染病者I(t)，治療者T(t)，顯然有N(t)=S(t)+V(t)+E(t)+I(t)+T(t)。假設(shè)Λ為總?cè)丝谥械妮斎肼?，q表示輸入人口中進(jìn)入易感者倉室的比例，其余輸入人口均為疫苗接種者，a表示針對易感者的接種率，b表示被接種者免疫力的喪失率，μ表示自然死亡率，d表示因病死亡率，β表示染病者對易感者的傳染率系數(shù)，k表示潛伏者到染病者的轉(zhuǎn)化系數(shù)，α表示治療率?；谝陨霞僭O(shè)，建立如圖1所示的結(jié)核病傳播機(jī)制。

圖1 結(jié)核病傳播機(jī)制示意圖

對應(yīng)的結(jié)核病傳播動力學(xué)模型為

(1)

將式(1)的所有方程相加可得

(2)

Ω={(S，V，E，I，T)|S≥0，V≥0，E≥0，

對于系統(tǒng)(1)，有以下結(jié)論。

命題1Ω在R5內(nèi)為系統(tǒng)(1)的正向不變集。

證明定義函數(shù)

N(t)=S(t)+V(t)+E(t)+I(t)+T(t)

則

(3)

由于系統(tǒng)(1)的前4個(gè)方程均與T無關(guān)，故系統(tǒng)(1)可導(dǎo)出下面的等價(jià)系統(tǒng)

(4)

在下述討論中只考慮等價(jià)系統(tǒng)(4)。

2 平衡點(diǎn)和基本再生數(shù)

通過簡單計(jì)算可得，系統(tǒng)(4)的無病平衡點(diǎn)為P0=(S0，V0，0，0)，其中

利用下一代矩陣法[13]，可得系統(tǒng)(4)的基本再生數(shù)為

為了得到系統(tǒng)(4)的正平衡點(diǎn)P*=(S*，V*，E*，I*)，令

(5)

由式(5)的前2個(gè)方程可得

由式(5)的最后一個(gè)方程可得

將S*與E*的表達(dá)式代入式(5)的第3個(gè)方程可得

(6)

顯然，式(6)有2個(gè)根

(ⅰ)I0=0；對應(yīng)無病平衡點(diǎn);

當(dāng)R0>1時(shí)，系統(tǒng)(4)存在唯一的正平衡點(diǎn)P*=(S*，V*，E*，I*)，其中S*，V*，E*，I*的表達(dá)式如下

(7)

綜上，有下面的定理。

定理1 當(dāng)R0≤1時(shí)，系統(tǒng)(4)僅有一個(gè)無病平衡點(diǎn)P0；當(dāng)R0>1時(shí)，系統(tǒng)(4)有兩個(gè)平衡點(diǎn)，分別是無病平衡點(diǎn)P0和正平衡點(diǎn)P*。

3 平衡點(diǎn)的穩(wěn)定性

3.1 無病平衡點(diǎn)的穩(wěn)定性

定理2 當(dāng)R0<1時(shí)，系統(tǒng)(4)的無病平衡點(diǎn)P0在Ω內(nèi)全局漸近穩(wěn)定；當(dāng)R0>1時(shí)，無病平衡點(diǎn)不穩(wěn)定。

證明為了證明P0的全局穩(wěn)定性，構(gòu)造如下的Lyapunov函數(shù)

V(t)=kE+(k+μ)I

它沿著系統(tǒng)(4)的解的導(dǎo)數(shù)為

k(βSI-(k+μ)E)+

(k+μ)(kE-(d+μ+α)I)=

kβ(S-S0)I+(k+μ)·

(d+μ+α)I(R0-1)

下面證明R0>1時(shí)，P0不穩(wěn)定。系統(tǒng)(4)在P0處的Jacobian矩陣為

矩陣J(P0)的特征方程為

(8)

則

(α+k+d+2μ)2>4kβS0+(α-k+d)2

化簡可得

即R0<1，與定理中R0>1矛盾，故假設(shè)錯(cuò)誤，所以λ1>0，從而式(8)有一個(gè)正實(shí)根，因此R0>1時(shí)，無病平衡點(diǎn)不穩(wěn)定。

3.2 正平衡點(diǎn)的穩(wěn)定性

定理3 當(dāng)R0>1時(shí)，系統(tǒng)(4)唯一的正平衡點(diǎn)P*局部漸近穩(wěn)定。

證明系統(tǒng)(4)在P*處的Jacobian矩陣為

特征方程為

λ4+a1λ3+a2λ2+a3λ+a4=0

其中

a1=a+μ+b+μ+βI*+d+μ+α+k+μ>0

a2=(a+μ)(b+μ)+βI*(b+μ)-

ab+(d+μ+α+k+μ)(a+μ+b+μ+βI*)>0

a3=(d+μ+α+k+μ)[(a+μ)(b+μ)+

βI*(b+μ)-a*b]+

(d+μ+α)(k+μ)βI*>0

a4=(d+μ+α)(k+μ)(b+μ)βI*>0

定理4 當(dāng)R0>1時(shí)，系統(tǒng)(4)唯一的正平衡點(diǎn)P*全局漸近穩(wěn)定。

(9)

利用式(9)，系統(tǒng)(4)可被改寫為

構(gòu)造如下的Lyapunov函數(shù)

V(t)=S*(x-1+lnx)+

E*(z-1+lnz)+

則

F(x，y，z，u)

當(dāng)bV*<(a+μ)S*時(shí)，上式可化為

F(x，y，z，u)=((a+μ)S*-bV*)·

當(dāng)bV*>(a+μ)S*時(shí)，上式可化為

4 最優(yōu)控制策略

最優(yōu)控制方法有助于找到經(jīng)濟(jì)有效的控制各種疾病的策略，從而達(dá)到盡可能減少患病人數(shù)與相應(yīng)戰(zhàn)略成本的目的。因此，我們可以通過接種疫苗，提高治療率及隔離等策略以不同的成本實(shí)現(xiàn)對結(jié)核病的控制。令X=(S，V，E，I，T)，定義U={ui|i=1，2，3}，其中u1表示可以提高個(gè)體免疫力的疫苗接種策略，u2表示可以提高染病者治療率的控制策略，u3表示減少易感者與感染者接觸的隔離策略。由于醫(yī)療技術(shù)與成本的控制，每一種控制策略ui都是有上界uimax的。具有控制策略的模型由以下非線性微分方程組給出：

(10)

目標(biāo)集X為

X={X(·)∈W1，1([0，T];R5)|滿足式(1)與(10)}

控制集U為

U={U(·)∈L1([0，T];R3)|

0

考慮目標(biāo)函數(shù)

(11)

(12)

由文獻(xiàn)[15-16]可得系統(tǒng)最優(yōu)控制解的存在性，接下來根據(jù)Pontryagin最大值原理，給出最優(yōu)控制解的表達(dá)形式。

若U(·)∈X為在最終固定時(shí)間T上的最優(yōu)解，則存在非平凡，全連續(xù)的映射λ∶[0，T]→R5。

λ(t)=(λ1(t)，λ2(t)，λ3(t)，λ4(t)，λ5(t))

稱為協(xié)態(tài)向量，即存在不為零的協(xié)態(tài)向量λ(t)，使得

① 控制方程滿足

② 協(xié)態(tài)方程滿足

③ 極小值條件為

H(X◇(t)，U◇(t)，λ◇(t))=

④ 極小值條件對于?t∈[0，T]成立，構(gòu)造哈密爾頓函數(shù)

λ1(qΛ+bV-(1-u3)βSI-(u1+μ)S)+

λ2((1-q)Λ-(b+μ)V+u1S)+

λ3((1-u3)βSI-(k+μ)E)+

λ4(kE-(d+μ+u2)I)+λ5(u2I-μT)

⑤ 橫截條件為

λi(t)=0，i=1，…，5

由文獻(xiàn)[15-17]可得以下定理。

定理5 如果系統(tǒng)(10)存在最優(yōu)控制U(t)及相應(yīng)的最優(yōu)解(S◇(·)，V◇(·)，E◇(·)，I◇(·)，T◇(·))，則存在協(xié)態(tài)向量λi，i=1，…，5使得

橫截條件為

λi(T)=0，i=1，…，5

最優(yōu)控制解的表達(dá)形式為

5 數(shù)值模擬

5.1 平衡點(diǎn)的穩(wěn)定性

在系統(tǒng)(4)中選取參數(shù)q=0.01，Λ=100，b=0.03，β=0.025，a=0.35，μ=0.1，k=0.032，d=0.12，α=0.3，通過計(jì)算可得R0=0.752 7<1，此時(shí)，系統(tǒng)(4)存在惟一的無病平衡點(diǎn)P0，由定理(2)可得P0是全局漸近穩(wěn)定的，如圖2所示。

圖2 當(dāng)R0<1時(shí)，P0全局漸近穩(wěn)定曲線

在系統(tǒng)(4)中選取另外一組參數(shù)q=0.15，Λ=1 500，b=0.45，a=0.55，β=0.000 375，k=0.48，μ=0.096，d=0.002 6，α=0.65，通過計(jì)算可得R0=2.763 8>1。此時(shí)，系統(tǒng)(4)存在惟一的正平衡點(diǎn)P*，由定理可得P*是全局漸近穩(wěn)定的，且在第25天(A點(diǎn))時(shí)染病者數(shù)量達(dá)到峰值，如圖3所示。

圖3 當(dāng)R0>1時(shí)，P*全局漸近穩(wěn)定曲線

5.2 敏感性分析

敏感性分析常用來確定模型對參數(shù)值的穩(wěn)定性。圖4表明Λ對R0的影響很小，β、b、k、q均為R0的正相關(guān)變量，μ、a、d、α均為R0的負(fù)相關(guān)變量。顯然，在所有變量中，β對R0的影響最大，即β是R0中更重要的因素。從圖5和圖6可以看出：β也是Imax與患病人數(shù)到達(dá)峰值時(shí)間的關(guān)鍵參數(shù)。換言之，敏感性分析告訴我們預(yù)防勝于治療，在控制肺結(jié)核傳播方面，加強(qiáng)預(yù)防的努力顯得尤為重要。

圖4 R0與參數(shù)的相關(guān)性直方圖

圖5 Imax與參數(shù)的相關(guān)性直方圖

圖6 峰值時(shí)間與參數(shù)的相關(guān)性直方圖

5.3 最優(yōu)控制策略

為了驗(yàn)證最優(yōu)控制的有效性，利用前后掃描法與4階龍格庫塔公式相結(jié)合的方法對其進(jìn)行了數(shù)值模擬?？紤]到醫(yī)療技術(shù)與成本的因素，每種控制策略都有上限，分別取控制變量的最大值如下：u1=0.6，u2=0.8，u3=0.7。由圖7可知，隨著時(shí)間的增加，各種控制策略的投入可以適當(dāng)?shù)販p少，從而減少成本；圖8更為直接地反映了無論是R0>1還是R0<1，采取最優(yōu)控制方案后，患病者的數(shù)量都急劇減少，顯然最優(yōu)策略的應(yīng)用可以使結(jié)核病得到有效的控制。

圖7 控制變量u(t)與t的函數(shù)關(guān)系

圖8 施加控制前后染病者的數(shù)量

6 結(jié)論

研究了一種按比例輸入的結(jié)核病模型的全局動力學(xué)行為，通過構(gòu)造Lyapunov函數(shù)，利用LaSalle不變集原理的方法，研究了系統(tǒng)的無病平衡點(diǎn)和正平衡點(diǎn)的穩(wěn)定性。當(dāng)R0<1 時(shí)，無病平衡點(diǎn)是全局漸近穩(wěn)定的，此時(shí)疾病將會逐漸消失；當(dāng)R0>1時(shí)，正平衡點(diǎn)是全局漸近穩(wěn)定的，在這種情況下，結(jié)核病將一直存在。通過對基本再生數(shù)，最高患病人數(shù)及患病人數(shù)到達(dá)峰值時(shí)間作敏感性分析得到傳染率系數(shù)β是一個(gè)關(guān)鍵參數(shù)，對于結(jié)核病的預(yù)防及控制有重要意義。因此，在人們的日常生活中，可以通過提高個(gè)人的預(yù)防意識來降低β的數(shù)值，從而達(dá)到有效預(yù)防和控制結(jié)核病的目的。