按比例輸入的結核病模型的動力學研究

李淑萍，劉花芳

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

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

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

1 模型的建立

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

圖1 結核病傳播機制示意圖

對應的結核病傳播動力學模型為

(1)

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

(2)

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

對于系統(1)，有以下結論。

命題1Ω在R5內為系統(1)的正向不變集。

證明定義函數

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

則

(3)

由于系統(1)的前4個方程均與T無關，故系統(1)可導出下面的等價系統

(4)

在下述討論中只考慮等價系統(4)。

2 平衡點和基本再生數

通過簡單計算可得，系統(4)的無病平衡點為P0=(S0，V0，0，0)，其中

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

為了得到系統(4)的正平衡點P*=(S*，V*，E*，I*)，令

(5)

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

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

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

(6)

顯然，式(6)有2個根

(ⅰ)I0=0；對應無病平衡點;

當R0>1時，系統(4)存在唯一的正平衡點P*=(S*，V*，E*，I*)，其中S*，V*，E*，I*的表達式如下

(7)

綜上，有下面的定理。

定理1 當R0≤1時，系統(4)僅有一個無病平衡點P0；當R0>1時，系統(4)有兩個平衡點，分別是無病平衡點P0和正平衡點P*。

3 平衡點的穩定性

3.1 無病平衡點的穩定性

定理2 當R0<1時，系統(4)的無病平衡點P0在Ω內全局漸近穩定；當R0>1時，無病平衡點不穩定。

證明為了證明P0的全局穩定性，構造如下的Lyapunov函數

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

它沿著系統(4)的解的導數為

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

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

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

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

下面證明R0>1時，P0不穩定。系統(4)在P0處的Jacobian矩陣為

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

(8)

則

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

化簡可得

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

3.2 正平衡點的穩定性

定理3 當R0>1時，系統(4)唯一的正平衡點P*局部漸近穩定。

證明系統(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 當R0>1時，系統(4)唯一的正平衡點P*全局漸近穩定。

(9)

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

構造如下的Lyapunov函數

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

E*(z-1+lnz)+

則

F(x，y，z，u)

當bV*<(a+μ)S*時，上式可化為

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

當bV*>(a+μ)S*時，上式可化為

4 最優控制策略

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

(10)

目標集X為

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

控制集U為

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

0

考慮目標函數

(11)

(12)

由文獻[15-16]可得系統最優控制解的存在性，接下來根據Pontryagin最大值原理，給出最優控制解的表達形式。

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

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

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

① 控制方程滿足

② 協態方程滿足

③ 極小值條件為

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

④ 極小值條件對于?t∈[0，T]成立，構造哈密爾頓函數

λ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

由文獻[15-17]可得以下定理。

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

橫截條件為

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

最優控制解的表達形式為

5 數值模擬

5.1 平衡點的穩定性

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

圖2 當R0<1時，P0全局漸近穩定曲線

在系統(4)中選取另外一組參數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，通過計算可得R0=2.763 8>1。此時，系統(4)存在惟一的正平衡點P*，由定理可得P*是全局漸近穩定的，且在第25天(A點)時染病者數量達到峰值，如圖3所示。

圖3 當R0>1時，P*全局漸近穩定曲線

5.2 敏感性分析

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

圖4 R0與參數的相關性直方圖

圖5 Imax與參數的相關性直方圖

圖6 峰值時間與參數的相關性直方圖

5.3 最優控制策略

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

圖7 控制變量u(t)與t的函數關系

圖8 施加控制前后染病者的數量

6 結論

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