登革熱傳播模型的最優(yōu)綜合控制研究*

李雅芝，劉利利

（1.黔南民族師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，貴州 都勻 558000；2.黔南州復(fù)雜系統(tǒng)與智能優(yōu)化重點實驗室，貴州 都勻 558000；3.山西大學(xué) 復(fù)雜系統(tǒng)研究所，太原 030006）

引 言

近些年來，登革熱在世界各地迅速流行.據(jù)估計，每年約1 億有癥狀登革熱病例，另有約3 億無癥狀感染病例，負(fù)擔(dān)最大的是亞洲(75%)[1]，這不僅威脅到了人們的身體健康，而且造成了嚴(yán)重的經(jīng)濟(jì)負(fù)擔(dān)，因此需要尋找有效的防控措施以緩解登革熱在世界范圍內(nèi)的傳播.

目前還沒有針對登革熱的特定疫苗或抗病毒治療，故媒介控制是防控登革熱病毒傳播的重要策略.這些策略有：通過環(huán)境治理防止蚊子獲得產(chǎn)卵地、妥善處理固體廢物，去除人為制造的蚊蟲棲息地；對室外儲水容器施用適當(dāng)殺蟲劑；使用個人家庭防護(hù)，如紗窗、長袖衣服、經(jīng)殺蟲劑處理的材料、蚊香和噴霧式殺蟲劑；不育技術(shù)和基因控制等[2].在實際中，往往會同時實施多種控制措施以達(dá)到更好的控制效果.但是任何一種控制措施都要耗費一定的人力、物力和財力，如何在花費盡可能少的前提下達(dá)到最好的控制效果是一個值得研究的問題.

關(guān)于登革熱的最優(yōu)控制已有一些科研工作者進(jìn)行了研究，如Agusto 和Khan[3]考慮了使用殺蟲劑和疫苗的最優(yōu)控制.Prasetyo 等[4]建立了一個SVIR 模型研究最優(yōu)疫苗控制策略.Xue 等[5]建立了一個多菌株的登革熱模型并考慮了最優(yōu)疫苗控制.Fister 等[6]研究了投放絕育蚊子和生活環(huán)境改造對蚊媒傳染病的最優(yōu)控制.Masud 等[7]建立了一個SIR-SI 模型，加入兩種控制措施，研究其在降低蚊子對人的叮咬率和蚊子被人感染的有效性，并考慮了最優(yōu)控制.從現(xiàn)有文獻(xiàn)中發(fā)現(xiàn)，在模型中加入基因控制研究最優(yōu)控制策略的研究還較少，關(guān)于基因控制目前主要是投放攜帶Wolbachia 的蚊子.有實驗發(fā)現(xiàn)Wolbachia 可以抑制登革熱病毒在蚊子體內(nèi)的復(fù)制，所以讓蚊子攜帶Wolbachia 有助于抑制登革熱的傳播.而只要成年雌蚊感染了Wolbachia，其后代也會感染此菌，并且被感染的雄性與正常雌性的后代會因為細(xì)胞質(zhì)不相容(CI)而死亡[8]，這些特性都有助于Wolbachia 在蚊子種群的傳播.相比傳統(tǒng)的控制措施，使用Wolbachia 是一種可持續(xù)且環(huán)境友好的控制方式.

為了研究登革熱的最優(yōu)控制策略，本文將建立登革熱在蚊子和人之間的傳播模型，加入自我保護(hù)、殺蟲劑和Wolbacia 三種控制措施，研究最優(yōu)綜合控制策略.

1 模型建立

關(guān)于登革熱傳播模型的建立，有以下幾點解釋：

(Ⅰ)對于人采用SIS 倉室模型，因為登革熱的潛伏期相對人的壽命是很短的，在此忽略.而登革熱有四種類型，人感染了一種類型的登革熱并恢復(fù)后，仍可能感染其他類型的登革熱，所以SIS 模型是合理的.

(Ⅱ)對于蚊子采用SEI 倉室模型，因為登革熱病毒在蚊子體內(nèi)的潛伏期(約7 d)相對蚊子的整個生命周期(30～100 d)是不能被忽略的.

(Ⅲ)考慮Wolbachia 在蚊子種群中的垂直傳播和CI 影響.

(Ⅳ)控制變量u1(t)表示通過使用蚊帳、穿長袖衣服等保護(hù)人類免受蚊蟲叮咬的效果，控制變量u2(t)表示為減少蚊子數(shù)量而進(jìn)行的捕殺工作.

綜合以上幾種因素，可得到如下模型：

其中Nmi=Smi+Emi+Imi,Nmu=Smu+Emu+Imu,Nm=Nmi+Nmu,Nh=Sh+Ih.各狀態(tài)變量和參數(shù)的含義分別見表1和表2，“W”表示“Wolbachia”，“D”表示“登革熱”.由于有Wolbachia 菌株的垂直傳播率和CI 發(fā)生可能性均接近1，故在模型(1)中，將兩種影響的可能性視為1，即只要成年雌蚊感染了Wolbachia，其后代也一定會感染，感染雄蚊和正常雌蚊的后代一定會因為CI 影響而死亡.模型(1)第二個方程中的項表示正常雌性產(chǎn)生正常后代.另外，由于對蚊子的捕殺不可能超過蚊子的出生，否則蚊子種群將會由于捕殺而滅絕，故假設(shè)bm>u2.

表1 模型(1)中各狀態(tài)變量的含義Table 1 The meanings of variables in model (1)

表2 模型(1)中各參數(shù)的含義Table 2 The meanings of parameters in model (1)

定義

關(guān)于模型(1)的適定性，有如下結(jié)果.

定理1模型(1)的解非負(fù)且在D內(nèi)是有界的.進(jìn)一步地，D是模型(1)的正不變集.

證明由于系統(tǒng)(1)的初值非負(fù)，易知模型(1)的解非負(fù)，故只需證明Nm和Nh是有界的.

分別將系統(tǒng)(1)中的前六個方程和后兩個方程相加，可以得到

積分可得

得證.

2 常數(shù)控制分析

令u1(t)=u1,u2(t)=u2，其中u1,u2為常數(shù).下考慮常數(shù)控制對模型(1)基本再生數(shù)的影響.

在系統(tǒng)(1)中令Emi=0,Imi=0,Emu=0,Imu=0,Ih=0，則得到如下子系統(tǒng)：

定義g(0,0)=0，則系統(tǒng)(2)被連續(xù)延拓到 (0,0).分別令f(Smi,Smu)=0,g(Smi,Smu)=0,v(Sh)=0，通過計算模型(2)的平衡態(tài)可知模型(1)最多存在三個無病平衡態(tài) (,,0,0,0,0,0)，分別為

由解對初值的連續(xù)依賴性可知E00是不穩(wěn)定的.記E01和E02對應(yīng)的基本再生數(shù)分別為R01和R02，使用下一代矩陣方法可以計算得到

通過計算系統(tǒng)(1)在E0i(i=1,2)處 的Jacobi 矩陣對應(yīng)特征多項式的特征根，可以得到：當(dāng)R01<1時，E01是局部漸近穩(wěn)定的；當(dāng)R02<1時，E02是局部漸近穩(wěn)定的.

比較R01和R02的表達(dá)式，其區(qū)別在于參數(shù)和 βh，由表2可知，<βh，故有R01

下面分析R0i(i=1,2)對 控制參數(shù)u1和u2的 敏感性，以R01為例.

在此采用Chitnis 等[9]定義的標(biāo)準(zhǔn)正向靈敏度指數(shù)衡量變量P對參數(shù)u的敏感性.結(jié)合R01的表達(dá)式，通過計算可以得到

可以看出，R01與u1和u2均為負(fù)相關(guān)的，即增大常數(shù)控制參數(shù)u1和u2，均會使得R0i(i=1,2)減小，進(jìn)而減少登革熱的感染人數(shù).又? ||/?u1=1/(1?u1)2>0,?||/?u2=bm/2(bm?u2)2>0，所以u1和u2越 大，R01與它們的相關(guān)性越強.

圖1為當(dāng)R0i(i=1,2)>1時，控制參數(shù)u1和u2對 登革熱感染人數(shù)的影響.參數(shù)取值為：bm=0.8,dm=4E?2,βm=0.001 5,em=0.05,w=2,dh=4E?5,γh=0.002,βh=0.000 1,q3=0.9.從圖中可以看出：控制參數(shù)取值越大，登革熱感染人數(shù)越少，且病例數(shù)呈顯著上升趨勢所需的時間越長.

圖1 常數(shù)控制對登革熱感染人數(shù)的影響Fig.1 Effects of constant control on the number of dengue infections

3 最優(yōu)控制分析

由于Wolbachia 主要是通過抑制登革熱病毒在蚊子體內(nèi)的復(fù)制，進(jìn)而達(dá)到抑制登革熱傳播的目的，對于蚊子種群屬于內(nèi)部控制，控制效果與最初的投放量有關(guān)，故在此僅研究自我保護(hù)和殺蟲劑兩種外部控制措施的最優(yōu)綜合控制策略.

首先建立如下目標(biāo)泛函：

其中(u1(t),u2(t))∈K，這里的K為控制集，定義為A表示由攜帶登革熱病毒的病人造成的人均損失，B1和B2分別表示人的自我保護(hù)和滅蚊工作相關(guān)的成本，T表示控制的終端時間.假設(shè)時間區(qū)間為 [0,T]，我們需要尋找最優(yōu)的u1(t)和u2(t)使 得目標(biāo)泛函J(u1(t),u2(t))達(dá)到最小.

為了研究系統(tǒng)(1)在目標(biāo)泛函(3)下的最優(yōu)控制策略，現(xiàn)在給出最優(yōu)控制的存在性結(jié)論如下.

定理2對于模型(1)，存在最優(yōu)控制(t),(t)，使得

證明易知目標(biāo)函數(shù)J(u1(t),u2(t))中 的被積函數(shù)關(guān)于u1(t)和u2(t)是凸函數(shù).由定理1 知系統(tǒng)(1)的解是有界的，故系統(tǒng)關(guān)于變量Smi(t),Smu(t),Emi(t),Emu(t),Imi(t),Imu(t),Sh(t)，Ih(t)滿足Lipschitz 性質(zhì)，進(jìn)而存在最優(yōu)的二元組((t),(t)).得證.

關(guān)于動力系統(tǒng)的最優(yōu)控制方法很多，動力系統(tǒng)類型不同，求解最優(yōu)控制的方法也不同[10-11].模型(1)為一個高維ODE 系統(tǒng)，可以采用Pontryagin 極值原理[12]尋找具體的最優(yōu)控制((t),(t)).

首先構(gòu)造Hamilton 函數(shù)H如下：

其中λi(t)(i=1,2,···,8)稱為協(xié)態(tài)變量，滿足下面的伴隨方程：

橫截條件為

利用Hamilton 函數(shù)，我們將狀態(tài)方程(1)、目標(biāo)泛函(3)和伴隨方程(4)聯(lián)立起來構(gòu)成一個最優(yōu)控制系統(tǒng).由最優(yōu)條件可得

綜上，關(guān)于最優(yōu)控制的具體形式有如下結(jié)果.

定理3在控制集K上使得J(u1(t),u2(t))取 得最小值的最優(yōu)控制((t),(t))為

4 數(shù)值模擬

本部分采用四階Runge-Kutta 前推回代法計算最優(yōu)控制和狀態(tài)值[13].記狀態(tài)變量向量為

x=(Smi,Smu,Emi,Emu,Imi,Imu,Sh,Ih),

協(xié)態(tài)變量向量為

λ=(λ1,λ2,λ3,λ4,λ5,λ6,λ7,λ8),

控制變量向量為

u=(u1,u2).

其基本算法步驟如下：

步1 在時間范圍[0,T]內(nèi)對u進(jìn)行初步估計.

步2 使用初始條件和估計的u值，依據(jù)系統(tǒng)(1)關(guān)于時間向前求解x.

步3 使用橫截條件(5)和u，x的值, 依據(jù)系統(tǒng)(4)關(guān)于時間向后求解λ.

步4 通過在最優(yōu)控制的表征(6)中輸入新的x和λ 值來更新u.

步5 檢查收斂性.重復(fù)前面的步驟，直到當(dāng)前的狀態(tài)值、伴隨值和控制值充分收斂.

取參數(shù)值為：bm=0.8,dm=4E?5,βm=0.001 5,em=0.05,w=2,dh=4E?5,γh=0.002,βh=0.001,q3=0.9.本組參數(shù)值使得無控制系統(tǒng)的基本再生數(shù)為232.5204，也即登革熱會在人群中流行.控制的終端時間T=30.初值取為：Smi(0)=1 000,Smu=1 000,Emi=100,Emu=100,Imi=20,Imu=20,Sh=100,Ih=5.

控制策略對登革熱感染人數(shù)的影響見圖2，權(quán)重取值為：A=1，B1=1，B2=1.可以看出，通過采取控制措施，在一個月內(nèi)被感染人數(shù)有了顯著降低，且最優(yōu)綜合控制策略的效果最好，蚊子捕殺的控制效果優(yōu)于人的自我保護(hù)的控制效果.這充分說明直接對傳播媒介本身采取例如捕殺等控制措施是非常有必要的.

圖2 控制策略對登革熱感染人數(shù)的影響Fig.2 Effects of the control strategy on the number of dengue infections

最優(yōu)綜合控制與單一控制的比較見圖3，權(quán)重取值為：A=1,B1=1,B2=1.由圖3(a)可知，對于人的自我保護(hù)，如果同時還對蚊子進(jìn)行捕殺，在前2 d，人的自我保護(hù)在控制下界，之后大約需要在上界持續(xù)12 d 左右，隨后逐漸下降到控制下界；如果只有人的自我保護(hù)，沒有捕殺蚊子，可以發(fā)現(xiàn)幾乎在整個模擬周期內(nèi)人的自我保護(hù)都要維持在上界.由圖3(b)可知，對于蚊子的捕殺，如果同時還有人的自我保護(hù)，對蚊子的捕殺最開始就要在控制上界，而且要持續(xù)24 d 左右，之后再緩慢下降到控制下界；如果只捕殺蚊子，沒有人的自我保護(hù)，對蚊子的捕殺要在控制上界持續(xù)27 d 左右.總之，相比綜合控制策略，單一控制策略在上界的時間要長.

圖3 最優(yōu)綜合控制與單一控制的比較：(a)兩種控制策略下u 1的 變化圖；(b)兩種控制策略下u 2的變化圖Fig.3 Comparison between the optimal integrated control and the single control: (a)the variation diagram of u 1 under 2 control strategies;(b)the variation diagram of u 2 under 2 control strategies

不同權(quán)重對控制變量的影響見圖4.B1和B2取值越大表示對應(yīng)控制的成本越高，A取值越大表示感染登革熱造成的損失越大.從圖4(a)、(b)可以發(fā)現(xiàn)，隨著控制成本增大，兩種控制措施的控制力度都在減小.從圖4(c)、(d)可以發(fā)現(xiàn)，隨著疾病造成的損失增大，人的自我保護(hù)力度要增大到最大.

圖4 不同權(quán)重對控制變量的影響：(a)改變 B1對u 1的 影響，A =1,B2=1； (b)改變 B2對u 2的 影響，A =1,B1=1；(c)A=1 時，u 1,u2 的變化圖，B 1=1,B2=1； (d)A =5 時，u 1,u2 的 變化圖，B1=1,B2=1Fig.4 The influences of different weights on control variables: (a)the influence of B1 on u 1, A =1,B2=1; (b)the influence of B2 on u 2 , A =1,B1=1;(c)the variation diagram of u 1 and u2, A =1,B1=1,B2=1; (d)the variation diagram of u 1 and u2, A =5 , B1=1,B2=1

5 結(jié) 語

本文建立了一個SEI-SIS 模型描述登革熱在人和蚊子之間的傳播，加入了Wolbachia、人的自我保護(hù)和蚊子捕殺三種控制措施，研究最優(yōu)綜合控制策略，并采用數(shù)值模擬探討控制策略對登革熱感染人數(shù)的影響，對比綜合控制策略與單一控制策略，對比權(quán)重對控制策略的影響.在實際中，登革熱的控制方法有多種，可以根據(jù)實際情況，對不同的控制方法進(jìn)行組合來研究最優(yōu)綜合控制策略.