廖 夢，王興玉，王連文,施 欣

(湖北民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，湖北 恩施 445000)

結(jié)核病是因感染結(jié)核分枝桿菌而引起的一種慢性傳染病，主要由患者咳嗽、咳痰等方式向空氣中釋放細(xì)菌而廣泛傳播.在新型冠狀病毒肺炎(COVID-19)大流行之前，結(jié)核病是單一傳染病導(dǎo)致死亡的主要原因，排在人類免疫缺陷病毒(HIV) /艾滋病(AIDS)之前[1].2019年，約有1 000萬人患上結(jié)核病，并有140萬人死亡[1].更糟糕的是，到2020年，全世界約有25%的人口是潛伏性結(jié)核病感染者(latent tuberculosis infection，LTBI)[1].

1 SVLEIT結(jié)核病模型的建立

個體感染結(jié)核桿菌后會經(jīng)過幾個月、幾年甚至幾十年的潛伏期，長短因人而異.短期潛伏性結(jié)核病感染者稱為快LTBI，長期潛伏性結(jié)核病感染者稱為慢LTBI，前者發(fā)展成活動性結(jié)核病的風(fēng)險(約95%)遠(yuǎn)高于后者(約5%)[2].這對隱藏在人群中的數(shù)量龐大的潛伏性結(jié)核病感染者的及時發(fā)現(xiàn)和檢測工作帶來巨大挑戰(zhàn)，也是結(jié)核病難以快速消除的一個重要原因.慶幸的是，結(jié)核病是可以預(yù)防和治愈的[1,3].事實(shí)上，在過去的一個世紀(jì)里，卡介苗(Bacille Calmette-Guérin，BCG)已被廣泛接種[1,4]，全世界絕大多數(shù)國家建議新生兒出生后就立即接種該疫苗，對于未接種的嬰幼兒和學(xué)齡兒童予以補(bǔ)種，可用于預(yù)防嬰兒和10歲以下兒童的嚴(yán)重結(jié)核病，如結(jié)核性腦膜炎和粟粒性結(jié)核病.然而，卡介苗對人群不能提供完全有效地保護(hù)，而且疫苗效率不完美，平均為80%[5].治療雖然是控制結(jié)核病流行的重要措施之一，但是對大多數(shù)結(jié)核病患者而言治療效果往往是不完全的，即治療成功后患者不再具有傳染性，但體內(nèi)仍然攜帶結(jié)核分枝桿菌[6].考慮到治療環(huán)境等因素的影響，處于治療階段的患者可能具有傳染性[7-9].此外未被治療的患者有一定的自愈能力，自愈后可能再次轉(zhuǎn)化為潛伏者[9-10].基于文獻(xiàn)[10]，本文建立了一類具有快慢潛伏感染、不完全接種、不完全治療和自愈的SVLEIT結(jié)核病模型：

(1)

圖1 SVLEIT模型流程Fig.1 The flow diagram for SVLEIT model

其中S、V、L、E、I和T分別表示易感者、接種者、慢LTBI、快LTBI、感染者和治療者.模型(1)的流程和參數(shù)的具體含義如圖1和表1所示，模型中的所有參數(shù)都是非負(fù)的.

表1 參數(shù)的生物意義Tab.1 Biological description of parameters

許多學(xué)者通過構(gòu)造Lyapunov函數(shù)的方法來解決模型平衡點(diǎn)的全局穩(wěn)定性問題，其中李健全等[11-12]應(yīng)用代數(shù)方法來確定Volterra型Lyapunov函數(shù)的組合系數(shù)，高大鵬等[13]運(yùn)用該代數(shù)方法解決了文獻(xiàn)[7]所建立的SVLIT結(jié)核病模型地方病平衡點(diǎn)的全局穩(wěn)定性.本文的貢獻(xiàn)在于如下兩方面:① 由于模型(1)的維數(shù)高，運(yùn)用傳統(tǒng)的Routh-Hurwitz判據(jù)難以有效研究其平衡點(diǎn)局部穩(wěn)定性，因此將充分利用控制再生數(shù)和平衡點(diǎn)方程，通過反證法來解決其局部穩(wěn)定性；② 由于模型(1)具有多個非線性項(xiàng)，并考慮了不完全接種和不完全治療，應(yīng)用代數(shù)方法[11-12]來確定Volterra型Lyapunov函數(shù)的組合系數(shù)計(jì)算過于繁瑣，將呈現(xiàn)一個簡單的方法來確定這些組合系數(shù).

2 模型的局部穩(wěn)定性

對于模型(1)，集合

是模型(1)的正向不變集，其中S0=q1μΛ/(μ+ν),V0=q2Λ+q1νΛ/(μ+ν).易知模型(1)總存在一個無病平衡點(diǎn)ρ0=(S0,V0,0,0,0,0).記α1=μ+ν,α2=μ,α3=μ+η,α4=μ+θ,α5=μ+d+γ,α6=μ+δ,x0=λ(p1S0+V0).利用下一代矩陣法[14]計(jì)算基本再生數(shù)為

(2)

為方便，記x*=λ(I*+κT*)，直接計(jì)算可知地方病平衡點(diǎn)ρ*(S*,V*,E*,L*,I*,T*)滿足下面等式：

定理1當(dāng)且僅當(dāng)Rc>1，模型存在唯一的地方病平衡點(diǎn)ρ*(S*,V*,E*,L*,I*,T*).

定理2由定理1，如下等式明顯成立,

(3)

定理3如果Rc<1，ρ0是局部漸近穩(wěn)定的；而Rc>1，ρ0則變得不穩(wěn)定.

證明模型(1)在ρ0處的特征方程為

(a+α1)(a+α2)(a+α3)(a+α4)(a+α5)(a+α6)=

(a+α1)(a+α2)(a+α6+κk3γ)[x0θη+p2S0λθ(a+α3)]+(a+α1)(a+α2)×

[k1γθη(a+α6)+k2γθ(a+α3)(a+α6)+l1k3γδθη+l2k3γδθ(a+α3)].

(4)

顯然，a1=-α1<0,a2=-α2<0，其他特征值滿足等式

(a+α3)(a+α4)(a+α5)(a+α6)=[(a+α6+κk3γ)x0θη+p2S0λθ(a+α3)]+

k1γθη(a+α6)+k2γθ(a+α3)(a+α6)+l1k3γδθη+l2k3γδθ(a+α3).

(5)

(6)

定理4如果唯一的地方病平衡點(diǎn)ρ*存在，則它是局部漸近穩(wěn)定的.

證明ρ*處的特征方程為

A1A2A3A4A5A6+νλx*S*θη(A6+κk3γ)=

θλS*(a+α1)A2(p1η+p2A3)(A6+κk3γ)+λθηV*A1(a+μ)(A6+κk3γ)+

k1γθηA1A2A6+k2γθA1A2A3A6+l1k3γδθηA1A2+l2k3γδθA1A2A3,

(7)

其中A1=a+x*+α1,A2=a+x*+α1,Ai=a+αi,i=3,4,5,6.

3 模型的全局穩(wěn)定性

定理5(i) 如果Rc<1，ρ0在Ω內(nèi)是全局漸近穩(wěn)定的；(ii) 當(dāng)Rc>1，ρ*在Ω*={φ∈Ω|I>0}內(nèi)是全局漸近穩(wěn)定的.

證明：(i) 首先證明ρ*的全局穩(wěn)定性.為此選擇一個合適的Lyapunov函數(shù)

其中Volterra型的函數(shù)ψ(u)=u-1-lnu對于u>0，滿足ψ(u)≥0，當(dāng)且僅當(dāng)u=1時ψ(u)=0.常數(shù)ci>0,i=1,2,…,6待定.注意到平衡點(diǎn)滿足下面等式：

q2μΛ=-νS*+

c2λ(I+κT)V*-c2λ(I+κT)V+

c3p1λ(I+κT)S+c3

c3

c3

合并上式中λ(I+κT)S，λ(I+κT)V和和的同類項(xiàng)，可得

c2λ(I+κT)V-

c3

其中,

o1=-c3p1λ(I*+κT*)S*-c3λ(I*+κT*)V*-c3k1γI*-c3l1δT*+c4ηL*,

o2=-c4p2λ(I*+κT*)S*-c4k2γI*-c4l2δT*-c4ηL*+c5θE*，

o3=c1λS*I*+c2λV*I*+c3k1γI*+c4k2γI*-c5θE*+c6k3γI*，

o4=c1λκS*T*+c2λκV*T*+c3l1δT*+c4l2δT*-c6k3γI*.

通過比較系數(shù)發(fā)現(xiàn):

(8)

通過(i)中類似論證，ρ0在Ω內(nèi)是全局漸近穩(wěn)定的.

4 數(shù)值模擬

本節(jié)通過數(shù)值模擬來檢驗(yàn)前面理論結(jié)果.取定λ=2.411×10-9[16],Λ=1.240 1×109[17],κ=0.024 2[17],μ=1/76.34[18],q1=0.05[19],q2=1-q1,ν=0.508 3[17],=0.35[20],p1=0.90[10],p2=1-p1,d=0.3[21],l2=0.2[18],η=0.5[17],θ=0.016 9[17],γ=0.623 5[17].結(jié)合生物學(xué)意義，可以選取k1=0.3,δ=0.6,k2=0.13,k3=0.57,δ=0.6,l1=0.8,l1=1-l2,使得控制再生數(shù)Rc=0.979 1<1.0和ρ*的全局漸近穩(wěn)定性,如圖2所示.由圖2(a)可知,隨著t→+∞,I→0,T→0,因此無病平衡點(diǎn)ρ0是全局漸近穩(wěn)定的.如果選取λ=3.134 3×10-9而保持其余參數(shù)不變,使得Rc=1.272 9>1.由圖2(b)可知，隨著t→+∞,I→I*=3.21×106,T→T*=1.86×106,因此地方病平衡點(diǎn)ρ*全局漸近穩(wěn)定.

(a) Rc=0.979 1<1時ρ0全局漸近穩(wěn)定 (b) Rc=1.272 9>1時ρ*的全局漸近穩(wěn)定圖2 ρ0和ρ*的全局漸近穩(wěn)定性Fig.2 Globally asymptotically stability of ρ0 and ρ*

接下來，通過對新生兒的接種比例q2、疾病傳染率λ、活動性結(jié)核病例治療比例k3進(jìn)行10%的改變來進(jìn)行參數(shù)敏感性分析，結(jié)果如圖3所示.由圖3可知：疾病傳染率λ是最敏感參數(shù)(見圖3(b))，活動性結(jié)核病例治療比例k3為第2敏感參數(shù)(見圖3(c))，新生兒的接種比例p2為第3敏感參數(shù)(見圖3(a))，但鑒于接種是控制成本最小的措施，上述3種措施均能有效減小活動性結(jié)核病例數(shù)量，應(yīng)當(dāng)同時實(shí)施.

(a) 參數(shù)q2的敏感性分析 (b) 參數(shù)λ的敏感性分析 (c) 參數(shù)k3的敏感性分析圖3 對參數(shù)q2、λ、k3進(jìn)行敏感性分析Fig.3 Sensitivity analysis for the paramenters q2,λ,k3are conducted

5 結(jié)語

綜合文獻(xiàn)[7]和[10]研究的結(jié)核病傳播模型反映的快慢潛伏感染、不完全接種、不完全治療和自愈等結(jié)核病傳播機(jī)理，進(jìn)一步建立一類具有多個非線性項(xiàng)和高維的SVLEIT結(jié)核病模型,通過反證法得到各平衡點(diǎn)的局部穩(wěn)定性，并運(yùn)用一個簡單的方法確定了Volterra型Lyapunov函數(shù)的組合系數(shù)，分析發(fā)現(xiàn)該模型呈現(xiàn)出全局閾值動力學(xué)，即該模型的基本再生數(shù)決定了結(jié)核病滅絕或持久.通常，傳統(tǒng)的代數(shù)方法所確定的Volterra型Lyapunov函數(shù)的組合系數(shù)不總是唯一的，提出的新方法能唯一確定Lyapunov函數(shù)的組合系數(shù)，有效簡化了分析計(jì)算，在對高維的非線性傳染病模型的全局穩(wěn)定性分析中顯得更為突出.此外，數(shù)值模擬探究表明參數(shù)敏感性揭示減少疾病傳染率λ、提高活動性結(jié)核病例治療比例k3和新生兒的接種比例q2均能有效減少活動性結(jié)核病例，這些措施可以為結(jié)核病的防治提供理論參考.