引　　言

近些年来，登革热在世界各地迅速流行。据估计，每年约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影响。

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

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

其中$N_{\rm{mi}}=S_{\rm{mi}}+E_{\rm{mi}}+I_{\rm{mi}}, N_{\rm{mu}}=S_{\rm{mu}}+E_{\rm{mu}}+I_{\rm{mu}}, N_{\rm{m}}=N_{\rm{mi}}+N_{\rm{mu}}, N_{\rm{h}}=S_{\rm{h}}+I_{\rm{h}}$。各状态变量和参数的含义分别见表1和表2，“W”表示“Wolbachia”，“D”表示“登革热”。由于有Wolbachia菌株的垂直传播率和CI发生可能性均接近1，故在模型(1)中，将两种影响的可能性视为1，即只要成年雌蚊感染了Wolbachia，其后代也一定会感染，感染雄蚊和正常雌蚊的后代一定会因为CI影响而死亡。模型(1)第二个方程中的项$b_{\rm{m}}\left(1-\dfrac{N_{\rm{mi}}}{N_{\rm{m}}}\right)N_{\rm{mu}}$表示正常雌性产生正常后代。另外，由于对蚊子的捕杀不可能超过蚊子的出生，否则蚊子种群将会由于捕杀而灭绝，故假设$b_{\rm{m}}>u_{2}$。

表  1  模型(1)中各状态变量的含义

Table  1.  The meanings of variables in model (1)

state variable	biological meaning
$S_{{\rm{mi}}}$	mosquito population infected with W and susceptible to D
$S_{{\rm{mu}}}$	mosquito population not infected with W and susceptible to D
$E_{{\rm{mi}}}$	mosquito population infected with W and in the incubation period of D
$E_{{\rm{mu}}}$	mosquito population not infected with W and in the incubation period of D
$I_{{\rm{mi}}}$	mosquito population infected with W and D
$I_{{\rm{mu}}}$	mosquito population not infected with W but infected with D
$S_{{\rm{h}}}$	human population susceptible to D
$I_{{\rm{h}}}$	human population infected with D

下载: 导出CSV 

| 显示表格

表  2  模型(1)中各参数的含义

Table  2.  The meanings of parameters in model (1)

parameter	biological meanings
$w$	recruitment rate of human beings
$b_{{\rm{m}}}$	birth rate of mosquitoes
$d_{{\rm{h}}}$	natural death rate of human beings
$d_{{\rm{m}}}$	natural death rate of mosquitoes
$\beta_{{\rm{h}}}$	D infection rate of human beings bitten by W-free mosquitoes
$e_{{\rm{m}}}$	rate at which exposed mosquitoes become infectious
$\overline{\beta}_{{\rm{h}}}=q\beta_{{\rm{h}}}$	D infection rate of human beings bitten by W-infected mosquitoes (q<1, which reflects the inhibition effect of W)
$\beta_{{\rm{m}}}$	D infection rate of mosquitoes biting D-infected human beings
$\gamma_{{\rm{h}}}$	recovery rate of infectious human beings

下载: 导出CSV 

| 显示表格

定义

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

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

证明　由于系统(1)的初值非负，易知模型(1)的解非负，故只需证明$N_{\rm{m}}$和$N_{\rm{h}}$是有界的。

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

积分可得

得证。

2.   常数控制分析

令$u_{1}(t)=u_{1}, u_{2}(t)=u_{2}$，其中$u_{1}, u_{2}$为常数。下考虑常数控制对模型(1)基本再生数的影响。

在系统(1)中令$E_{\rm{mi}}=0, I_{\rm{mi}}=0, E_{\rm{mu}}=0, I_{\rm{mu}}=0, I_{\rm{h}}=0$，则得到如下子系统：

定义$g(0,0)=0$，则系统(2)被连续延拓到$(0,0)$。分别令$f(S_{\rm{mi}},S_{\rm{mu}})=0, g(S_{\rm{mi}},S_{\rm{mu}})=0, v(S_{\rm{h}})=0$，通过计算模型(2)的平衡态可知模型(1)最多存在三个无病平衡态 $(S^{*}_{\rm{mi}},S^{*}_{\rm{mu}},0,0,0,0,S^{*}_{\rm{h}},0)$，分别为

由解对初值的连续依赖性可知$E_{00}$是不稳定的。记$E_{01}$和$E_{02}$对应的基本再生数分别为$R_{01}$和$R_{02}$，使用下一代矩阵方法可以计算得到

通过计算系统(1)在$E_{0i}(i=1,2)$处的Jacobi矩阵对应特征多项式的特征根，可以得到：当$R_{01}<1$时，$E_{01}$是局部渐近稳定的；当$R_{02}<1$时，$E_{02}$是局部渐近稳定的。

比较$R_{01}$和$R_{02}$的表达式，其区别在于参数$\overline{\beta}_{\rm{h}}$和$\beta_{\rm{h}}$，由表2可知，$\overline{\beta}_{\rm{h}}<\beta_{\rm{h}}$，故有$R_{01}<R_{02}$。所以Wobachia有助于降低登革热的传播力度和规模。

下面分析$R_{0i}(i=1,2)$对控制参数$u_{1}$和$u_{2}$的敏感性，以$R_{01}$为例。

在此采用Chitnis等^[9]定义的标准正向灵敏度指数$\varUpsilon^{P}_{u}:=\dfrac{\partial P}{\partial u}\times\dfrac{u}{P}$衡量变量P对参数u的敏感性。结合$R_{01}$的表达式，通过计算可以得到

可以看出，$R_{01}$与$u_{1}$和$u_{2}$均为负相关的，即增大常数控制参数$u_{1}$和$u_{2}$，均会使得$R_{0i}(i=1,2)$减小，进而减少登革热的感染人数。又$\partial|\varUpsilon^{R_{01}}_{u_{1}}|/\partial u_{1}=1/(1-u_{1})^{2}>0 , \partial|\varUpsilon^{R_{01}}_{u_{2}}|/\partial u_{2}=b_{\rm{m}}/2(b_{\rm{m}}-u_{2})^{2}>0$，所以$u_{1}$和$u_{2}$越大，$R_{01}$与它们的相关性越强。

图1为当$R_{0i}(i=1,2)>1$时，控制参数$u_{1}$和$u_{2}$对登革热感染人数的影响。参数取值为：$b_{{\rm{m}}}=0.8, d_{{\rm{m}}}=4{\rm{E}}-$$2, \beta_{{\rm{m}}}=0.001\;5, e_{{\rm{m}}}=0.05, w=2, d_{{\rm{h}}}=4{\rm{E}}-5, \gamma_{{\rm{h}}}=0.002, \beta_{{\rm{h}}}=0.000\;1, q_{3}=0.9$。从图中可以看出：控制参数取值越大，登革热感染人数越少，且病例数呈显著上升趋势所需的时间越长。

图  1  常数控制对登革热感染人数的影响

Figure  1.  Effects of constant control on the number of dengue infections

下载: 全尺寸图片 幻灯片

3.   最优控制分析

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

首先建立如下目标泛函：

其中$(u_{1}(t),u_{2}(t))\in K$，这里的K为控制集，定义为

A表示由携带登革热病毒的病人造成的人均损失，$B_{1}$和$B_{2}$分别表示人的自我保护和灭蚊工作相关的成本，T表示控制的终端时间。假设时间区间为$[0,T]$，我们需要寻找最优的$u_{1}(t)$和$u_{2}(t)$使得目标泛函$J(u_{1}(t),u_{2}(t))$达到最小。

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

定理2　对于模型(1)，存在最优控制$u_{1}^{*}(t), u_{2}^{*}(t)$，使得

证明　易知目标函数$J(u_{1}(t),u_{2}(t))$中的被积函数关于$u_{1}(t)$和$u_{2}(t)$是凸函数。由定理1知系统(1)的解是有界的，故系统关于变量$S_{\rm{mi}}(t), S_{\rm{mu}}(t), E_{\rm{mi}}(t), E_{\rm{mu}}(t), I_{\rm{mi}}(t), I_{\rm{mu}}(t), S_{\rm{h}}(t)$，$I_{\rm{h}}(t)$满足Lipschitz性质，进而存在最优的二元组$(u_{1}^{*}(t), u_{2}^{*}(t))$。得证。

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

首先构造Hamilton函数H如下：

其中$\lambda_{i}(t)(i=1,2,\cdots,8)$称为协态变量，满足下面的伴随方程：

横截条件为

利用Hamilton函数，我们将状态方程(1)、目标泛函(3)和伴随方程(4)联立起来构成一个最优控制系统。由最优条件$\dfrac{\partial H}{\partial u_{1}}=0$和$\dfrac{\partial H}{\partial u_{2}}=0$可得

其中，$S_{\rm{mi}}^{*}(t),S_{\rm{mu}}^{*}(t),E_{\rm{mi}}^{*}(t),E_{\rm{mu}}^{*}(t),I_{\rm{mi}}^{*}(t),I_{\rm{mu}}^{*}(t)$，$S_{\rm{h}}^{*}(t),I_{\rm{h}}^{*}(t)$分别是相应状态变量的最优解；$\lambda_{i}^{*}(t)(i=1,2,\cdots,8)$是系统(4)在条件(5)下的解。

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

定理3　在控制集K上使得$J(u_{1}(t),u_{2}(t))$取得最小值的最优控制$(u_{1}^{*}(t),u_{2}^{*}(t))$为

4.   数 值 模 拟

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

协态变量向量为

控制变量向量为

其基本算法步骤如下：

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

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

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

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

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

取参数值为：$b_{\rm{m}}=0.8, d_{\rm{m}}=4{\rm{E}}-5, \beta_{\rm{m}}=0.001\;5, e_{\rm{m}}=0.05, w=2, d_{\rm{h}}=4{\rm{E}}-5, \gamma_{\rm{h}}=0.002, \beta_{\rm{h}}=0.001, q_{3}=0.9$。本组参数值使得无控制系统的基本再生数为232.5204，也即登革热会在人群中流行。控制的终端时间$T=30$。初值取为：$S_{\rm{mi}}(0)=1\;000, S_{\rm{mu}}=1\;000, E_{\rm{mi}}=100, E_{\rm{mu}}=100, I_{\rm{mi}}=20, I_{\rm{mu}}=20, S_{\rm{h}}=100, I_{\rm{h}}=5$。

图  2  控制策略对登革热感染人数的影响

Figure  2.  Effects of the control strategy on the number of dengue infections

下载: 全尺寸图片 幻灯片

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

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

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

图  3  最优综合控制与单一控制的比较：(a) 两种控制策略下$u_{1}$的变化图；(b) 两种控制策略下$u_{2}$的变化图

Figure  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  不同权重对控制变量的影响：(a) 改变$B_{1}$对$u_1$的影响，$A=1, B_{2}=1$；(b) 改变$B_{2}$对$u_2$的影响，$A=1, B_{1}=1$；(c) A=1时，$u_{1}, u_{2}$的变化图，$B_{1}=1, B_{2}=1$；(d) $A=5$时，$u_{1}, u_{2}$的变化图，$B_{1}=1, B_{2}=1$

4.  The influences of different weights on control variables: (a) the influence of $B_{1}$ on $u_{1}$, $A=1, B_{2}=1$; (b) the influence of $B_{2}$ on $u_{2}$, $A=1, B_{1}=1$; (c) the variation diagram of $u_{1}\; {\rm{and}} \;u_{2},$ $A=1, B_{1}=1, B_{2}=1$; (d) the variation diagram of $u_{1}\; {\rm{and}} \;u_{2},$ $A=5$, $B_{1}=1, B_{2}=1$

下载: 全尺寸图片 幻灯片

5.   结　　语

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