登革热传播模型的最优综合控制研究*

李雅芝，刘利利

（1.黔南民族师范学院 数学与统计学院，贵州 都匀 558000；2.黔南州复杂系统与智能优化重点实验室，贵州 都匀 558000；3.山西大学 复杂系统研究所，太原 030006）

引 言

近些年来，登革热在世界各地迅速流行.据估计，每年约1 亿有症状登革热病例，另有约3 亿无症状感染病例，负担最大的是亚洲(75%)[1]，这不仅威胁到了人们的身体健康，而且造成了严重的经济负担，因此需要寻找有效的防控措施以缓解登革热在世界范围内的传播.

目前还没有针对登革热的特定疫苗或抗病毒治疗，故媒介控制是防控登革热病毒传播的重要策略.这些策略有：通过环境治理防止蚊子获得产卵地、妥善处理固体废物，去除人为制造的蚊虫栖息地；对室外储水容器施用适当杀虫剂；使用个人家庭防护，如纱窗、长袖衣服、经杀虫剂处理的材料、蚊香和喷雾式杀虫剂；不育技术和基因控制等[2].在实际中，往往会同时实施多种控制措施以达到更好的控制效果.但是任何一种控制措施都要耗费一定的人力、物力和财力，如何在花费尽可能少的前提下达到最好的控制效果是一个值得研究的问题.

关于登革热的最优控制已有一些科研工作者进行了研究，如Agusto 和Khan[3]考虑了使用杀虫剂和疫苗的最优控制.Prasetyo 等[4]建立了一个SVIR 模型研究最优疫苗控制策略.Xue 等[5]建立了一个多菌株的登革热模型并考虑了最优疫苗控制.Fister 等[6]研究了投放绝育蚊子和生活环境改造对蚊媒传染病的最优控制.Masud 等[7]建立了一个SIR-SI 模型，加入两种控制措施，研究其在降低蚊子对人的叮咬率和蚊子被人感染的有效性，并考虑了最优控制.从现有文献中发现，在模型中加入基因控制研究最优控制策略的研究还较少，关于基因控制目前主要是投放携带Wolbachia 的蚊子.有实验发现Wolbachia 可以抑制登革热病毒在蚊子体内的复制，所以让蚊子携带Wolbachia 有助于抑制登革热的传播.而只要成年雌蚊感染了Wolbachia，其后代也会感染此菌，并且被感染的雄性与正常雌性的后代会因为细胞质不相容(CI)而死亡[8]，这些特性都有助于Wolbachia 在蚊子种群的传播.相比传统的控制措施，使用Wolbachia 是一种可持续且环境友好的控制方式.

为了研究登革热的最优控制策略，本文将建立登革热在蚊子和人之间的传播模型，加入自我保护、杀虫剂和Wolbacia 三种控制措施，研究最优综合控制策略.

1 模型建立

关于登革热传播模型的建立，有以下几点解释：

(Ⅰ)对于人采用SIS 仓室模型，因为登革热的潜伏期相对人的寿命是很短的，在此忽略.而登革热有四种类型，人感染了一种类型的登革热并恢复后，仍可能感染其他类型的登革热，所以SIS 模型是合理的.

(Ⅱ)对于蚊子采用SEI 仓室模型，因为登革热病毒在蚊子体内的潜伏期(约7 d)相对蚊子的整个生命周期(30～100 d)是不能被忽略的.

(Ⅲ)考虑Wolbachia 在蚊子种群中的垂直传播和CI 影响.

(Ⅳ)控制变量u1(t)表示通过使用蚊帐、穿长袖衣服等保护人类免受蚊虫叮咬的效果，控制变量u2(t)表示为减少蚊子数量而进行的捕杀工作.

综合以上几种因素，可得到如下模型：

其中Nmi=Smi+Emi+Imi,Nmu=Smu+Emu+Imu,Nm=Nmi+Nmu,Nh=Sh+Ih.各状态变量和参数的含义分别见表1和表2，“W”表示“Wolbachia”，“D”表示“登革热”.由于有Wolbachia 菌株的垂直传播率和CI 发生可能性均接近1，故在模型(1)中，将两种影响的可能性视为1，即只要成年雌蚊感染了Wolbachia，其后代也一定会感染，感染雄蚊和正常雌蚊的后代一定会因为CI 影响而死亡.模型(1)第二个方程中的项表示正常雌性产生正常后代.另外，由于对蚊子的捕杀不可能超过蚊子的出生，否则蚊子种群将会由于捕杀而灭绝，故假设bm>u2.

表1 模型(1)中各状态变量的含义Table 1 The meanings of variables in model (1)

表2 模型(1)中各参数的含义Table 2 The meanings of parameters in model (1)

定义

关于模型(1)的适定性，有如下结果.

定理1模型(1)的解非负且在D内是有界的.进一步地，D是模型(1)的正不变集.

证明由于系统(1)的初值非负，易知模型(1)的解非负，故只需证明Nm和Nh是有界的.

分别将系统(1)中的前六个方程和后两个方程相加，可以得到

积分可得

得证.

2 常数控制分析

令u1(t)=u1,u2(t)=u2，其中u1,u2为常数.下考虑常数控制对模型(1)基本再生数的影响.

在系统(1)中令Emi=0,Imi=0,Emu=0,Imu=0,Ih=0，则得到如下子系统：

定义g(0,0)=0，则系统(2)被连续延拓到 (0,0).分别令f(Smi,Smu)=0,g(Smi,Smu)=0,v(Sh)=0，通过计算模型(2)的平衡态可知模型(1)最多存在三个无病平衡态 (,,0,0,0,0,0)，分别为

由解对初值的连续依赖性可知E00是不稳定的.记E01和E02对应的基本再生数分别为R01和R02，使用下一代矩阵方法可以计算得到

通过计算系统(1)在E0i(i=1,2)处 的Jacobi 矩阵对应特征多项式的特征根，可以得到：当R01<1时，E01是局部渐近稳定的；当R02<1时，E02是局部渐近稳定的.

比较R01和R02的表达式，其区别在于参数和 βh，由表2可知，<βh，故有R01

下面分析R0i(i=1,2)对 控制参数u1和u2的 敏感性，以R01为例.

在此采用Chitnis 等[9]定义的标准正向灵敏度指数衡量变量P对参数u的敏感性.结合R01的表达式，通过计算可以得到

可以看出，R01与u1和u2均为负相关的，即增大常数控制参数u1和u2，均会使得R0i(i=1,2)减小，进而减少登革热的感染人数.又∂ ||/∂u1=1/(1−u1)2>0,∂||/∂u2=bm/2(bm−u2)2>0，所以u1和u2越 大，R01与它们的相关性越强.

图1为当R0i(i=1,2)>1时，控制参数u1和u2对 登革热感染人数的影响.参数取值为：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.从图中可以看出：控制参数取值越大，登革热感染人数越少，且病例数呈显著上升趋势所需的时间越长.

图1 常数控制对登革热感染人数的影响Fig.1 Effects of constant control on the number of dengue infections

3 最优控制分析

由于Wolbachia 主要是通过抑制登革热病毒在蚊子体内的复制，进而达到抑制登革热传播的目的，对于蚊子种群属于内部控制，控制效果与最初的投放量有关，故在此仅研究自我保护和杀虫剂两种外部控制措施的最优综合控制策略.

首先建立如下目标泛函：

其中(u1(t),u2(t))∈K，这里的K为控制集，定义为A表示由携带登革热病毒的病人造成的人均损失，B1和B2分别表示人的自我保护和灭蚊工作相关的成本，T表示控制的终端时间.假设时间区间为 [0,T]，我们需要寻找最优的u1(t)和u2(t)使 得目标泛函J(u1(t),u2(t))达到最小.

为了研究系统(1)在目标泛函(3)下的最优控制策略，现在给出最优控制的存在性结论如下.

定理2对于模型(1)，存在最优控制(t),(t)，使得

证明易知目标函数J(u1(t),u2(t))中 的被积函数关于u1(t)和u2(t)是凸函数.由定理1 知系统(1)的解是有界的，故系统关于变量Smi(t),Smu(t),Emi(t),Emu(t),Imi(t),Imu(t),Sh(t)，Ih(t)满足Lipschitz 性质，进而存在最优的二元组((t),(t)).得证.

关于动力系统的最优控制方法很多，动力系统类型不同，求解最优控制的方法也不同[10-11].模型(1)为一个高维ODE 系统，可以采用Pontryagin 极值原理[12]寻找具体的最优控制((t),(t)).

首先构造Hamilton 函数H如下：

其中λi(t)(i=1,2,···,8)称为协态变量，满足下面的伴随方程：

横截条件为

利用Hamilton 函数，我们将状态方程(1)、目标泛函(3)和伴随方程(4)联立起来构成一个最优控制系统.由最优条件可得

综上，关于最优控制的具体形式有如下结果.

定理3在控制集K上使得J(u1(t),u2(t))取 得最小值的最优控制((t),(t))为

4 数值模拟

本部分采用四阶Runge-Kutta 前推回代法计算最优控制和状态值[13].记状态变量向量为

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

协态变量向量为

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

控制变量向量为

u=(u1,u2).

其基本算法步骤如下：

步1 在时间范围[0,T]内对u进行初步估计.

步2 使用初始条件和估计的u值，依据系统(1)关于时间向前求解x.

步3 使用横截条件(5)和u，x的值, 依据系统(4)关于时间向后求解λ.

步4 通过在最优控制的表征(6)中输入新的x和λ 值来更新u.

步5 检查收敛性.重复前面的步骤，直到当前的状态值、伴随值和控制值充分收敛.

取参数值为：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.本组参数值使得无控制系统的基本再生数为232.5204，也即登革热会在人群中流行.控制的终端时间T=30.初值取为：Smi(0)=1 000,Smu=1 000,Emi=100,Emu=100,Imi=20,Imu=20,Sh=100,Ih=5.

控制策略对登革热感染人数的影响见图2，权重取值为：A=1，B1=1，B2=1.可以看出，通过采取控制措施，在一个月内被感染人数有了显著降低，且最优综合控制策略的效果最好，蚊子捕杀的控制效果优于人的自我保护的控制效果.这充分说明直接对传播媒介本身采取例如捕杀等控制措施是非常有必要的.

图2 控制策略对登革热感染人数的影响Fig.2 Effects of the control strategy on the number of dengue infections

最优综合控制与单一控制的比较见图3，权重取值为：A=1,B1=1,B2=1.由图3(a)可知，对于人的自我保护，如果同时还对蚊子进行捕杀，在前2 d，人的自我保护在控制下界，之后大约需要在上界持续12 d 左右，随后逐渐下降到控制下界；如果只有人的自我保护，没有捕杀蚊子，可以发现几乎在整个模拟周期内人的自我保护都要维持在上界.由图3(b)可知，对于蚊子的捕杀，如果同时还有人的自我保护，对蚊子的捕杀最开始就要在控制上界，而且要持续24 d 左右，之后再缓慢下降到控制下界；如果只捕杀蚊子，没有人的自我保护，对蚊子的捕杀要在控制上界持续27 d 左右.总之，相比综合控制策略，单一控制策略在上界的时间要长.

图3 最优综合控制与单一控制的比较：(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

不同权重对控制变量的影响见图4.B1和B2取值越大表示对应控制的成本越高，A取值越大表示感染登革热造成的损失越大.从图4(a)、(b)可以发现，随着控制成本增大，两种控制措施的控制力度都在减小.从图4(c)、(d)可以发现，随着疾病造成的损失增大，人的自我保护力度要增大到最大.

图4 不同权重对控制变量的影响：(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 结 语

本文建立了一个SEI-SIS 模型描述登革热在人和蚊子之间的传播，加入了Wolbachia、人的自我保护和蚊子捕杀三种控制措施，研究最优综合控制策略，并采用数值模拟探讨控制策略对登革热感染人数的影响，对比综合控制策略与单一控制策略，对比权重对控制策略的影响.在实际中，登革热的控制方法有多种，可以根据实际情况，对不同的控制方法进行组合来研究最优综合控制策略.