張靖文，王浩華，2，3

（1.海南大學(xué) 理學(xué)院，海南 海口 570228；2.海南大學(xué) 海南省工程建模與統(tǒng)計計算重點實驗室，海南 ?？?570228；3.海南大學(xué) 熱帶特色林木花卉遺傳與種質(zhì)創(chuàng)新教育部重點實驗室，海南 ?？?570228）

流行性傳染病的爆發(fā)不僅會導(dǎo)致大規(guī)模的死亡，而且會影響社會穩(wěn)定，嚴(yán)重影響經(jīng)濟和生活，因此，研究流行性傳染病的傳播動力學(xué)模型，在預(yù)測傳染病的動向和阻斷其傳播病源等方面至關(guān)重要［1?6］.乙肝病毒（HBV）通過血液傳播、母嬰傳播、無保護性行為傳播等方式感染肝臟細(xì)胞，誘發(fā)炎癥，從而引起多種疾?。?］.據(jù)世界衛(wèi)生組織報告，全球約有20億人攜帶乙肝病毒，約200萬人患有慢性肝感染，每年超過78萬人因此死亡［8］.

如何有效預(yù)防和控制乙肝的傳播是關(guān)心的熱點問題之一.目前，眾多研究者基于傳染病模型對HBV的傳播動力學(xué)及相關(guān)性質(zhì)進行了深入的研究.Khan等［9］在考慮疫苗接種因素下，對HBV傳播的阻斷進行了定量研究；Zhao等［10］研究表明，有效的信息干預(yù)手段（包括媒體宣傳，病例教育等），可以有效降低感染人數(shù)峰值；Bao等［11］分析了環(huán)境噪聲對隨機HBV感染模型傳播動力學(xué)的影響，并提出了干預(yù)策略；Wang等［12］研究具有細(xì)胞間傳播和免疫應(yīng)答的隨機HBV感染模型，分析了疾病了滅絕和持續(xù)性存在的條件.筆者考察了有效疫苗接種率以及信息干預(yù)下的HBV傳染病模型，給出了新模型平衡點的局部和全局穩(wěn)定性的充分條件及參數(shù)的動力學(xué)特性.此外，為了控制HBV的傳播，添加了信息傳播狀態(tài)變量，再利用最優(yōu)控制方法建立了基于疫苗有效率、信息反應(yīng)強度以及有效治療率的多維控制動力學(xué)模型，在論證最優(yōu)控制存在的基礎(chǔ)上，通過數(shù)值模型驗證的相關(guān)理論結(jié)果.

1 模型建立

為了有效控制疾病的傳播，將建立一個帶有疫苗接種和信息干預(yù)的HBV傳染模型來研究疾病的傳播動力學(xué).

首先，Khan等［9］基于疫苗接種率提出了HBV的傳播動力學(xué)模型

其中，S(t)，I(t)，R(t)分別代表易感人群，感染人群和恢復(fù)人群，λ代表出生率，α代表從易感人群到感染人群的轉(zhuǎn)移概率，μ0代表自然死亡率，μ1代表疾病導(dǎo)致的死亡率，γ1代表感染人群的治愈率，v代表疫苗接種率.Khan等將總?cè)丝趧澐譃?類：易感人群，感染人群和恢復(fù)人群.新出生的孩子都為易感人群，易感人群與感染人群接觸會以一定的比例轉(zhuǎn)移到感染人群中，易感人群接種疫苗會轉(zhuǎn)移到恢復(fù)人群中，3種人群都有固定的自然死亡率，感染人群有疾病導(dǎo)致的死亡以及被成功治療轉(zhuǎn)移到恢復(fù)人群中，恢復(fù)人群都有永久免疫力.

再考慮有效信息干預(yù)，得到如下帶有信息干預(yù)和疫苗接種的HBV傳染病模型（圖1）

圖1 模型（2）示意圖

其中，Z(t)表示信息密度，m表示信息干預(yù)率，w表示信息的反應(yīng)強度，a表示信息的增長率，b表示飽和常量，a0表示信息的自然衰退率.

2 模型平衡點與基本再生數(shù)

計算模型（2）的平衡點以及基本再生數(shù)，是研究其傳播動力學(xué)的基本步驟.基本再生數(shù)R0定義為患病期內(nèi)病人傳染致病的平均人數(shù).根據(jù)文獻［13］中的方法，模型（2）的基本再生數(shù)為

令模型（2）所有方程右端為0，可得模型（2）的2個平衡點：

a無病平衡點

b疾病平衡點

其中，

I*為下列方程的正根

即模型（2）的所有解都是非負(fù)的.

令總?cè)丝贜=S+I+R，由模型（2）可知

因此，得到不變集Γ，

3 平衡點分析

探討了模型2個平衡點的局部穩(wěn)定性，分叉和全局穩(wěn)定性，此外還對模型的參數(shù)進行了局部敏感性分析.

3.1平衡點的局部穩(wěn)定性模型（2）的變分矩陣為

定理1

1）若R0＜1，模型（2）的無病平衡點E1是局部漸進穩(wěn)定的；若R0＞1，則E1不穩(wěn)定.

2）若R0＞1，則模型有唯一的疾病平衡點E2，此外，若A1B1＞C1且A1(B1C1-A1D1)＞C21，則E2是局部漸進穩(wěn)定的.

證明

1）計算變分矩陣J在無平衡點E1的特征值為

2）計算變分矩陣J在疾病平衡點E2的特征方程為

3.2分叉

引理1［15］考慮下列帶參數(shù)φ的ODE方程

其中，0是系統(tǒng)的穩(wěn)定點，對任意φ都有f(0，φ)=0.假設(shè)fk是f的第k個分量，并且

其中，y，z分別表示特征值0的右特征向量和左特征向量，上述系統(tǒng)在0附近的局部動力學(xué)完全由a1，b1決定.

1）a1＞0，b1＞0.當(dāng)φ＜0并且|φ|?1時，0是局部漸進穩(wěn)定點，并且存在一個正的不穩(wěn)定的平衡點；當(dāng)0＜φ?1時，0是不穩(wěn)定點，并且存在一個負(fù)的局部漸進穩(wěn)定平衡點.

2）a1＜0，b1＜0.當(dāng)φ＜0并且|φ|?1時，0是不穩(wěn)定點；當(dāng)0＜φ?1時，0是局部漸進穩(wěn)定點，并且存在一個正的不穩(wěn)定平衡點.

3）a1＞0，b1＜0.當(dāng)φ＜0并且|φ|?1時，0是不穩(wěn)定點，并且存在一個負(fù)的局部漸進穩(wěn)定平衡點；當(dāng)0＜φ?1時，0是穩(wěn)定點，并且有一個正的不穩(wěn)點平衡點.

4）a1＜0，b1＞0.當(dāng)φ從負(fù)值變化到正值時，0從穩(wěn)定點變成不穩(wěn)定點.相應(yīng)地一個負(fù)的不穩(wěn)定平衡點變成正的并且局部漸近穩(wěn)定點.

定理2當(dāng)R0=1時，模型（2）有一個前叉分支.

證明令x1=S，x2=I，x3=R，x4=Z，再令φ=α為分支參數(shù).在φ=φ*=α*，由R0=1可以推出α*=模型（2）可以重新寫成

顯然，當(dāng)R0=1時，J x*(α*)有一特征值為0，并且其他特征值都非負(fù)，矩陣J x*(α*)對應(yīng)于特征值0的右特征向量為y=(y1，y2，y3，y4)＇，其中

另外，左特征向量為z=(z1，z2，z3，z4)，其中，z1=0，z2=1，z3=0，z4=0.根據(jù)引理1可知，當(dāng)R0=1時，可以利用a1，b1的值判別分支的方向.模型（4）中f=(f1，f2，f3，f4)在(x*，α*)上用于計算a1，b1的非零二階偏導(dǎo)有根據(jù)引理1可得

顯然a1＜0，b1＞0.根據(jù)引理1，當(dāng)R0=1時，模型（2）有一個前叉分支.

3.3平衡點的全局穩(wěn)定性假設(shè)系統(tǒng)可以寫成如下形式

其中，X∈R3，Y∈R分別代表未感染者和感染者.設(shè)U0=(X0，0)為模型（5）的無病平衡點.

引理2［16］如果R0＜1，且滿足下列2個條件：

2）對?(X，Y)∈Γ，有其中DYG(X0，0)是一個M-矩陣，Γ是前文定義的變量邊界，那么模型（5）的無病平衡點U0=(X0，0)是全局漸進穩(wěn)定的.

定理3若R0＜1，則模型（2）中的無病平衡點E1是全局漸近穩(wěn)定的.

證明根據(jù)引理2，把模型（2）改寫，其中

滿足G(X，0)=0，X=(S，R，Z)＇，Y=I.

對于模型（2）中的疾病平衡點E2，可以構(gòu)建一個Lyapunov函數(shù)來判別全局穩(wěn)定性.

定理4當(dāng)R0＞1時，模型（2）中的疾病平衡點E2是全局漸進穩(wěn)定的.

3.4局部敏感性分析模型（2）中的參數(shù)具有不確定性，會影響基本再生數(shù)R0的大小，局部敏感性分析揭示了基本再生數(shù)R0與模型中參數(shù)之間的關(guān)系.局部敏感性指標(biāo)定義為［17］

其中，Π為模型中的參數(shù).該分析提供了參數(shù)值與基本再生數(shù)之間的變化，為疾病的控制和傳播提供了依據(jù).利用式（6）和下列的參數(shù)值

得到下列敏感度指標(biāo)

由上式可知，參數(shù)α與R0成正比，參數(shù)α的增加或減少都會使R0增加或者減少.另一方面，參數(shù)v，μ1，γ1與R0成反比，參數(shù)v，μ1，γ1的增加或減少都會使R0減少或者增加.因此，在局部敏感性分析的基礎(chǔ)上，可以建立一個最優(yōu)控制系統(tǒng).

4 最優(yōu)控制

為了探究如何有效控制HBV的傳播，通過增加控制變量改進模型（2）.在模型（2）中添加3個控制變量：u1(t)表示疫苗的有效率，u2(t)表示信息反應(yīng)強度大小，u3(t)表示對感染人群治療的有效率.定義如下控制集

其中，T是控制策略結(jié)束的最終時間.方便起見，用u1，u2，u3表示3個控制變量.用k1I(t)表示感染人群的權(quán)重，k2，k3，k4分別表示疫苗的有效率，信息反應(yīng)強度大小，治療的有效率的權(quán)重，k2u12，k3u22，k4u32分別描述了疫苗接種，信息傳播和治療的相關(guān)成本.為了降低感染人群的數(shù)量，并使疫苗接種，信息傳播和治療成本最低.因此，控制問題可以描述為

受約束于

其中，k1，k2，k3，k4都是正常數(shù).式（7）表示為了控制疾病傳播的總代價，因此，尋找最優(yōu)控制u*=(u1*，u2*，u3*)，即使得J(u*)=minJ(u1，u2，u3).

定理5對于最優(yōu)控制問題（7）和（8），在U中一定存在最優(yōu)控制u*=(u1*，u2*，u3*)，得J(u*)=minJ(u1，u2，u3).

證明為了證明最優(yōu)控制的存在性，利用文獻［18］中的方法.根據(jù)定義，控制集U是凸閉集.由于控制變量和狀態(tài)變量都是非負(fù)的值，所以在最小值的問題上，目標(biāo)函數(shù)（7）在u1，u2，u3中的必要凸性是滿足的.因此，可以看出最優(yōu)系統(tǒng)是有界的，確保了最優(yōu)控制存在的緊密性.目標(biāo)函數(shù)（7）中的被積函數(shù)在控制集U上是凸的.

證畢.

為了找到最優(yōu)控制問題（7）和（8）的最優(yōu)解，定義如下Hamilton函數(shù)

其中，λ1，λ2，λ3，λ4是協(xié)態(tài)變量，根據(jù)Pontryagin最大值原則，得到如下定理.

定理6如果u1*，u2*，u3*是最優(yōu)控制問題（7）和（8）的最優(yōu)控制變量，且S*，I*，R*，Z*是相應(yīng)的最優(yōu)狀態(tài)變量，則存在協(xié)態(tài)變量λ1，λ2，λ3，λ4滿足下列協(xié)態(tài)方程

且邊界條件為λ1(T)=0，λ2(T)=0，λ3(T)=0，λ4(T)=0.

此外相應(yīng)的最優(yōu)控制變量u1*，u2*，u3*為

證明如果u1*，u2*，u3*是最優(yōu)控制問題（7）和（8）的最優(yōu)控制變量，且S*，I*，R*，Z*是相應(yīng)的最優(yōu)狀態(tài)變量.根據(jù)Pontryagin最大值原則，那么存在協(xié)態(tài)變量λ1，λ2，λ3，λ4滿足下列協(xié)態(tài)方程

根據(jù)定義的控制集U，有

可以改寫成式（9），（10）和（11）的形式.

證畢.

5 數(shù)值模擬

使用Runge-Kutta方法，參考Kumar［15］和Khan［18］等人的參數(shù)數(shù)據(jù)，利用Matlab軟件進行數(shù)值模擬.

模型（2）給定初始值S(0)=100，I(0)=40，R(0)=20，Z(0)=10以及參數(shù)［15，18］：λ=0.4，α=0.005，μ0=0.03，μ1=0.002，γ1=0.05，v=0.02，m=0.017，w=0.2，a=0.01，b=1，a0=0.06，可以得到R0=0.487 8＜1，圖2a驗證了定理3.再使用同樣的初始值以及參數(shù)［15，18］：λ=0.4，α=0.02，μ0=0.03，μ1=0.002，γ1=0.05，v=0.02，m=0.017，w=0.2，a=0.01，b=1，a0=0.06，可得R0=1.951 2＞1，圖2b驗證了定理4.

圖2 模型（2）中疾病平穩(wěn)和滅絕性分析

給定初始值S(0)=100，I(0)=20，R(0)=20，Z(0)=10以及下列的參數(shù)［15，18］：λ=0.4，α=0.005，μ0=0.03，μ1=0.002，γ1=0.05，v=0.02，m=0.017，w=0.2，a=0.01，b=1，a0=0.06，k1=1 000，k2=0.45，k3=0.55，k4=0.34.利用前推回代算法［15］解決控制問題（7）和（8），再對比沒有施加控制時的ODE模型（2），易感人群，感染人群和恢復(fù)人群的數(shù)量對比分別表示在圖3a，b和c中.

圖3a結(jié)果表明，三維綜合變量的引入，不僅大幅度地降低易感人群的數(shù)量，而且t=4為最低值，由于感染人群數(shù)量趨于零，疾病已經(jīng)滅絕，不存在人與人之間的傳播，所以易感人群不用再接種疫苗，又因為出生率大于死亡率，導(dǎo)致易感人群數(shù)量會慢慢增加.從圖3b可以看到，感染人數(shù)在控制變量的干預(yù)下幾乎完全改變了趨勢，不僅避免了峰值的出現(xiàn)，而且t=4時感染人數(shù)降為零，即三維綜合控制策略不僅能大幅度地降低感人數(shù)，避免峰值的出現(xiàn)，而且能有效的縮短病源的阻斷時間.圖3c則表明恢復(fù)人群將在短時間內(nèi)達到峰值.

圖3 α=0.005時控制對比圖，

利用相同的方法，把接觸傳染概率系數(shù)提高6倍，即令α=0.03，初始值和其他參數(shù)都不變，同樣對比

沒有施加控制時的ODE模型（2），易感人群，感染人群，恢復(fù)人群的數(shù)量對比分別表示在圖4a，b和c中.

圖4 α=0.03時控制對比圖，

從圖4可以看出，在擴大疾病傳染概率的情況下，本文模型也能較好地降低易感人群的數(shù)量以及加快疾病的滅絕.

圖5為最優(yōu)控制問題（7）和（8）中最優(yōu)控制變量的時間序列圖.從圖5中可以看出，在一開始時達到最大值，在一段時間后急速下降，最后下降到0，表明經(jīng)過開始階段的嚴(yán)厲措施后，隨著時間的延長，控制強度可以慢慢降低趨于常態(tài)化，且u3*衰減最慢，即實際HBV傳染源的阻斷需要對感染人群的治療投入較大，而疫苗有效率次之，信息反應(yīng)強度最少.

圖5 最優(yōu)控制變量時間序列圖

6 小結(jié)

研究基于信息干預(yù)下HBV傳染病模型，綜合考察了信息干預(yù)、疫苗接種以及有效治療率對模型動力學(xué)性質(zhì)的影響.計算了新模型的基本再生數(shù)和模型的平衡點，根據(jù)線性化和構(gòu)造Lyapunov函數(shù)的方法得到2個平衡點全局漸進穩(wěn)定的條件.數(shù)值驗證了2個平衡點的局部和全局漸進穩(wěn)定性，并給出了其分支性質(zhì).通過參數(shù)敏感性分析，引入控制變量，探究最優(yōu)控制問題（7）和（8），驗證了最優(yōu)控制的存在性并且根據(jù)Pontryagin最大值原則得出了最優(yōu)控制變量的表達式.結(jié)果表明，加入控制后能加速疾病的滅絕，有效地控制了疾病的傳播.