一类媒介具有阶段结构的西尼罗河病毒模型的敏感性分析

Sensitivity Analysis of a Model Formulated for West Nile Virus with a Stage-Structured Vector Population

摘要: 本文研究了一类结合鸟类–蚊子种群结构的西尼罗河病毒(WNv)传播模型，用PRCC方法对模型的蚊虫再生数和WNv疾病暴发再生数进行了敏感性分析，对防蚊灭蚊和控制病毒传播提出更好的防控措施。结果表明：控制WNv传播最好办法是降低媒介蚊子的数量，而限制蚊子的最佳策略是从幼虫阶段开始。另外未受感染鸟类的补给率的下降将有益于WNv的传播，因此在WNv暴发期间，利用控制鸟类的数量来抑制其传播是危险的。相反，我们应该适当加大未感染的鸟类的投放。

Abstract: A mathematical model is formulated for the transmission dynamics of West Nile virus (WNv) between a stage-structured mosquito and bird population. Through sensitivity analysis of the basic reproduction number for mosquito and for outbreaks of WNv by the PRCC method, preventive and control of anti-mosquito and WNv disease transmission are discussed. The results show that the best way to control the spread of WNv is to decrease the number of mosquitoes and control mosquitoes at the larval stage. Besides decreasing the recruitment rate of the uninfected birds is beneficial to the prevalence of the WNv. This observation suggests that it is a risk factor for the spread of WNv to control the birds during the period that the WNv prevails. By contraries, we should increase the recruitment rate of uninfected birds.

1. 引言

西尼罗河病毒由感染的蚊虫叮咬而传播，可使人和马等患上致命的脑炎，是一种人兽共患传染病。由于WNv对公共社会的巨大危害，增加对WNv传播动力学的理解是一个重要课题。自1999以来，已经有大量的数学模型来研究WNv的传播动力学 [1] [2] [3] [4]。Bowman等 [3] 提出了一种结合蚊虫–鸟类–人类种群动力学的模型来评估对病毒的控制策略。Blayneh等 [5] 考虑了一类确定性WNv动力学模型，模型中WNv疾病作为人畜共染病，考虑了蚊–鸟–人之间的循环传播。Abdelrazec等 [1] 将传播宿主鸟类分为乌鸦类和非乌鸦类，建立了仓室模型研究宿主物种对蚊子–鸟类WNv传播动力学的影响。Wang等 [6] 研究了因暴雨产生池塘的管理对传播和暴发的WNv的影响。这些模型为了解WNv传播机理提供了有用的信息。

然而，蚊子的生命周期是非常复杂的，经历了完全变态阶段，雄蚊种群可能会显著影响雌蚊种群大小，这可能会影响WNv的传播动力学。此外，成年雌蚊平均存活两到三周，而成年雄蚊通常存活大约一周，它们都比WNv主要寄主鸟类的寿命短得多。这意味着媒介蚊子种群的子系统比宿主鸟类(和人类)的种群系统快得多。因此，作为WNv的媒介，蚊虫种群动态在快速时间尺度上对WNv的传播起着重要的作用。近些年，有学者把蚊虫种群动力学引入其他蚊媒传染病的模型研究中 [7] [8] [9] [10]。Ai等 [8] 考虑了蚊子的四个不同的变质阶段，构建了一个具阶段结构的蚊子种群系统的疟疾模型来研究疟疾的传播动力学。Hughes等 [11] 和Ndii等 [12] 通过研究蚊子种群动力学模型，并在埃及伊蚊中引入Wolbachia病毒来防控登革热流行。

本文建立了一类结合鸟类–蚊子种群结构的西尼罗河病毒(WNv)传播的6维模型，模型结合了蚊子性别结构和2个变态阶段：成虫和水生阶段，分析了模型的动力学性态，通过敏感性分析得到对于基本再生数影响较大的参数，以此确定影响疾病传播的主要因素。

2. 模型建立

由于自然界中WNv主要在鸟类和雌性成虫蚊子之间 [10]，这里只考虑WNv在鸟类和蚊子之间的传播动力学。将蚊子的生命周期分为成虫阶段和水生阶段两个阶段，其中将卵、幼虫和蛹三个阶段划分为一个水生阶段， $A_M\left(t\right)$ 表示水生阶段蚊子的种群数量。成虫阶段的蚊子分为雄性和雌性两个仓室，数量设为 $M_M\left(t\right)$ 和 $F_M\left(t\right)$。因只有雌性成虫传播WNv病毒，因此将雌性种群分成两类：易感 $S_M\left(t\right)$ 和感染 $I_M\left(t\right)$。设鸟类种群数量为 $N_B\left(t\right)$，将鸟类分为易感 $B_S\left(t\right)$ 和感染 $B_I\left(t\right)$。

下面给出模型一些合理的假设。假定鸟类补给率是常数 ${\Pi }_B$，自然死亡率为 ${\mu }_B$，感染WNv疾病致死率为 $d_B$。设蚊子种群的变化主要来自出生和死亡。此外，假设雌蚊与雄蚊交配是随机的，产生的后代数量用函数表述 $\frac{\left(S_M+I_M\right)M_M}{S_M+I_M+M_M}$。值得注意的是，雌蚊在积水的表面上成群地产卵，因此水生蚊虫种群将受到环境承载力 的限制，这里将密度依赖效应应用到模型中，表示为 $1-\frac{A_M}{K}$。从上述讨论中，选取函数 ${\rho }_M\frac{\left(S_M+I_M\right)M_M}{S_M+I_M+M_M}\left(1-\frac{A_M}{K}\right)$ 作为水生蚊子种群的增长率， ${\rho }_M$ 表示雌性成虫的繁殖率。易感雌性成虫和鸟类感染WNv病毒的比率分别为 $b_M\beta S_M\frac{B_I}{B_S+B_I}$ 和 $b_M{\beta }_{MB}\frac{B_S}{B_S+B_I}{I}_M$。详见图1仓室转移图以及表1参数定义。文中提到的参数都是正的常数。

Figure 1. State transfer diagram of WNv model

图1. WNv模型的仓室转移图

由上分析可以建立结合鸟类–蚊子种群结构的6维的WNv传播模型：

$\{\begin{array}{l}\frac{\text{d}{A}_M}{\text{d}t}={\rho }_M\frac{\left(S_M+I_M\right)M_M}{S_M+I_M+M_M}\left(1-\frac{A_M}{K}\right)-{\mu }_A{A}_M-{\gamma }_M{A}_M\hfill \\ \frac{\text{d}{M}_M}{\text{d}t}={\epsilon }_M{\gamma }_M{A}_M-{\mu }_M{M}_M\hfill \\ \frac{\text{d}{S}_M}{\text{d}t}=\left(1-{\epsilon }_M\right){\gamma }_M{A}_M-{\mu }_F{S}_M-b_M\beta S_M\frac{B_I}{B_S+B_I}\hfill \\ \frac{\text{d}{I}_M}{\text{d}t}=b_M\beta S_M\frac{B_I}{B_S+B_I}-{\mu }_F{I}_M\hfill \\ \frac{\text{d}{B}_S}{\text{d}t}={\Pi }_B-{\mu }_B{B}_S-b_M{\beta }_{MB}\frac{B_S}{B_S+B_I}{I}_M\hfill \\ \frac{\text{d}{B}_I}{\text{d}t}=b_M{\beta }_{MB}\frac{B_S}{B_S+B_I}{I}_M-{\mu }_B{B}_I-d_B{B}_I\hfill \end{array}$ (1)

本文给出系统(1)的动力学分析，得到2个重要阈值：蚊虫基本再生数和WNv染病基本再生数，并用PRCC方法对再生数进行了敏感性分析，分析了各参数对疾病传播的影响程度，更好的了解疾病暴发的形成机理，从而提出更好的防控措施。

Table 1. Parameter definition of WNv model

表1. WNv模型的参数定义

3. 基本再生数

记 $R_M$ 和 $R_0$ 分别表示为蚊虫再生数和疾病WNv暴发的基本再生数，由3.1以及3.2可得

$R_M=\frac{{\rho }_M\left(1-{\epsilon }_M\right){\epsilon }_M{\gamma }_M}{\left({\mu }_A+{\gamma }_M\right)\left({\mu }_M\left(1-{\epsilon }_M\right)+{\epsilon }_M{\mu }_F\right)},R_0=\sqrt{\frac{b_M^2\beta {\beta }_{MB}{\mu }_B\left(1-{\epsilon }_M\right){\gamma }_{M}K}{{\mu }_F^2\left({\mu }_B+d_B\right){\Pi }_B}\left(1-\frac{1}{R_M}\right)}$

3.1. 蚊虫的基本再生数

定理1. 若 $R_M\le 1$，则系统只有一个无蚊虫的平衡点 $E_{00}\left(0,0,0,0,\frac{{\Pi }_B}{{\mu }_B},0\right)$。

根据平衡点的定义，系统正平衡点满足下面代数方程组

$\{\begin{array}{l}{\rho }_M\frac{F_M{M}_M}{F_M+M_M}\left(1-\frac{A_M}{K}\right)-{\mu }_A{A}_M-{\gamma }_M{A}_M=0\hfill \\ {\epsilon }_M{\gamma }_M{A}_M-{\mu }_M{M}_M=0\hfill \\ \left(1-{\epsilon }_M\right){\gamma }_M{A}_M-{\mu }_F{S}_M-b_M\beta S_M\frac{B_I}{N_B}=0\hfill \\ b_M\beta S_M\frac{B_I}{N_B}-{\mu }_F{I}_M=0\hfill \\ {\Pi }_B-{\mu }_B{B}_S-b_M{\beta }_{MB}\frac{B_S}{N_B}{I}_M=0\hfill \\ b_M{\beta }_{MB}\frac{B_S}{N_B}{I}_M-{\mu }_B{B}_I-d_B{B}_I=0\hfill \end{array}$ (2)

这里 $F_M=S_M+I_M,N_B=B_S+B_I$，将方程组(2)中第五个方程和第六个方程相加可得

$N_B=\frac{{\Pi }_B}{{\mu }_B}-\frac{d_B}{{\mu }_B}{B}_I,B_S=\frac{{\Pi }_B}{{\mu }_B}-\left(1+\frac{d_B}{{\mu }_B}\right)B_I$

将上式带入方程组(2)中第六个方程，可以得到

$I_M=\frac{\left({\Pi }_B-d_B{B}_I\right)B_I}{b_M{\beta }_{MB}\left(\frac{{\Pi }_B}{{\mu }_B+d_B}-B_I\right)}$

从(2)的第二、第三和第四个方程，分别得到

$M_M=\frac{{\epsilon }_M{\gamma }_M}{{\mu }_M}{A}_M,F_M=\frac{\left(1-{\epsilon }_M\right){\gamma }_M}{{\mu }_F}{A}_M$

将表达式带入(2)中第一和第四个方程，合并同类项，得到系统(1)的非负解需满足

$\{\begin{array}{l}{A}_M\left(1-\frac{1}{R_M}-\frac{A_M}{K}\right)=0\hfill \\ \begin{array}{l}{B}_I[\frac{d_B\left({\mu }_F{d}_B-b_M\beta {\mu }_B\right)}{b_M{\beta }_{MB}{\mu }_B}{B}_I^2+\left(b_M\beta \frac{\left(1-{\epsilon }_M\right){\gamma }_M}{{\mu }_F}{A}_M-\frac{{\Pi }_B\left(2{\mu }_F{d}_B-b_M\beta {\mu }_B\right)}{b_M{\beta }_{MB}{\mu }_B}\right)B_I\\ +\frac{{\Pi }_B^2{\mu }_F}{b_M{\beta }_{MB}{\mu }_B}\left(1-\frac{b_M^2\beta {\beta }_{MB}\left(1-{\epsilon }_M\right){\gamma }_M{A}_M}{{\Pi }_B{\mu }_F^2\left({\mu }_B+d_B\right)}\right)]=0\end{array}\hfill \end{array}$ (3)

若 $R_M\le 1$，从方程组(3)中第一个方程计算可得，只有唯一非负解 $A_M=0$，系统(1)存在唯一无蚊虫的平衡点 $E_{00}\left(0,0,0,0,\frac{{\Pi }_B}{{\mu }_B},0\right)$。若 $R_M\le 1$，蚊子灭绝，若 $R_M>1$，蚊虫能够存活。因此参数 $R_M$ 可以解释为蚊虫存在的基本再生数。

3.2. WNv疾病暴发的基本再生数

这里假定 $R_M>1$，计算可得系统(1)存在一个无病平衡点 $E_0\left(A_M^{*},M_M^{*},S_M^{*},0,\frac{{\Pi }_B}{{\mu }_B},0\right)$。这里

$A_M^{*}=K\left(1-\frac{1}{R_M}\right),M_M^{*}=\frac{{\gamma }_M{\epsilon }_{M}K}{{\mu }_M}\left(1-\frac{1}{R_M}\right),S_M^{*}=\frac{{\gamma }_M\left(1-{\epsilon }_M\right)K}{{\mu }_F}\left(1-\frac{1}{R_M}\right)$

根据Driessche和Watmoug [13] 的定义，可知系统(1)的新生感染项F和剩余转移项V关于各感染仓室的偏导在 $E_0$ 处取值后的雅克比矩阵可分别表示

$F=\left(\begin{array}{cc}0& \frac{b_M\beta {\mu }_B\left(1-{\epsilon }_M\right){\gamma }_{M}K}{{\mu }_F{\Pi }_B}\left(1-\frac{1}{R_M}\right)\\ b_M{\beta }_{MB}& 0\end{array}\right),V=\left(\begin{array}{cc}{\mu }_F& 0\\ 0& {\mu }_B+d_B\end{array}\right)$

现在用 $R_0$ 表示系统(1)疾病基本再生数，根据Driessche和Watmough中再生数的计算方法可得

$R_0=\rho \left(F{V}^{-1}\right)=\sqrt{\frac{b_M^2\beta {\beta }_{MB}{\mu }_B\left(1-{\epsilon }_M\right){\gamma }_{M}K}{{\mu }_F^2\left({\mu }_B+d_B\right){\Pi }_B}\left(1-\frac{1}{R_M}\right)}$

其中 $\rho \left(F{V}^{-1}\right)$ 表示矩阵 $F{V}^{-1}$ 的谱半径。由文献 [13] 的定理2可得到下面结论。

定理2. 当 $R_0<1$ 时，系统(1)的无病平衡点 $E_0$ 是局部渐近稳定的；当 $R_0>1$ 时，系统(1)的无病平衡点 $E_0$ 是不稳定的。

染病的雌性蚊虫和鸟类的种群数量分别为 $\frac{{\gamma }_M\left(1-{\epsilon }_M\right)K}{{\mu }_F}\left(1-\frac{1}{R_M}\right)$ 和 $\frac{{\Pi }_B}{{\mu }_B}$。在上述条件下，蚊虫叮咬染病的鸟类的频率定义为

${\hat{b}}_M:=\frac{\frac{b_M{\gamma }_M\left(1-{\epsilon }_M\right)K}{{\mu }_F}\left(1-\frac{1}{R_M}\right)}{\frac{{\Pi }_B}{{\mu }_B}}$

染病鸟类的平均寿命为 ${\hat{\tau }}_B:=\frac{1}{{\mu }_B+d_B}$，一个染病的鸟类成功传染给雌蚊的平均数量可以表示为

${\hat{I}}_M:={\hat{b}}_M{\hat{\tau }}_B\beta =\frac{b_M\beta {\mu }_B\left(1-{\epsilon }_M\right){\gamma }_{M}K}{{\mu }_F\left({\mu }_B+d_B\right){\Pi }_B}\left(1-\frac{1}{R_M}\right)$

相类似地，由一个染病的雌性蚊虫成功传染给鸟类的平均数量可以表示为 ${\hat{B}}_I:=b_M{\beta }_M\frac{1}{{\mu }_F}$。因此，当易感蚊虫和鸟类的数量分别为 $\frac{{\gamma }_M\left(1-{\epsilon }_M\right)K}{{\mu }_F}\left(1-\frac{1}{R_M}\right)$ 和 $\frac{{\Pi }_B}{{\mu }_B}$ 时，由一个染病鸟类成功传染给下一代的总的数量可以定义为

${\hat{R}}_0:={\hat{I}}_M{\hat{B}}_I=\frac{b_M^2\beta {\beta }_{MB}{\mu }_B\left(1-{\epsilon }_M\right){\gamma }_{M}K}{{\mu }_F^2\left({\mu }_B+d_B\right){\Pi }_B}\left(1-\frac{1}{R_M}\right)={\left(R_0\right)}^2$

根据基本再生数定义，参数 $R_0$ 可以解释为WNv传播的基本再生数，具有生物学意义。

4. 敏感性分析

敏感性分析指通过一系列的测试，使用不同的输入参数来观察预测参数值的变化如何影响系统的动态行为，识别输入参数与预测参数之间的关系 [14] [15] [16] [17]。LHS/PRCC法是一种全局性的参数敏感性分析法，可以快速有效识别敏感参数。偏秩相关系数(PRCC)绝对值越大相关性越强，正值代表正相关，负值代表负相关。 $R_M$ 和 $R_0$ 分别表示为蚊虫和WNv疾病暴发的基本再生数，决定了蚊虫大量暴发和疾病是否灭绝的阈值，通过采用敏感性分析可得到对于基本再生数影响较大的参数，从而提出更好的防控措施。

首先，设参数均服从某区间的均匀分布： ${\rho }_M\in \left(1,1.5\right)$， ${\epsilon }_M\in \left(0.4,0.6\right)$， ${\gamma }_M\in \left(1/8,1/6\right)$， ${\mu }_A\in \left(1/10,1/7.78\right)$， ${\mu }_M\in \left(1/16,1/7\right)$， ${\mu }_F\in \left(1/12,1/9\right)$， $b_M\in \left(0.05,0.3\right)$， $\beta \in \left(0.1,0.3\right)$， ${\beta }_{MB}\in \left(0.6,0.9\right)$， ${\mu }_B\in \left(1/1000,1/800\right)$， $d_B\in \left(1/5000,1/10000\right)$， ${\Pi }_B\in \left(80,100\right)$， $K\in \left(4000,5000\right)$，参数范围大部分来源于文献 [3] [12]。选取1000个样本，利用Mathematica，给出 $R_M$ 和 $R_0$ 的统计分布直方图，如图所示图2，图3。从图2(b)可得到：蚊虫基本再生数大约85%均大于1，均值和标准差分别为1.76和0.3210；从图3(b)可得到：WNv疾病基本再生数 $R_0$ 大约85%均小于1，均值和标准差分别为0.683和0.340。

(a) (b)

Figure 2. (a) Sensitivity Analysis of the basic reproduction number for mosquito $R_M$ ; (b) Uncertainty analysis of the basic reproduction number for mosquito $R_M$

图2. (a) 蚊虫基本再生数 $R_M$ 的敏感性分析；(b) 蚊虫基本再生数 $R_M$ 的不确定性分析

从 $R_M$ 表达式可以看出，蚊虫基本再生数 $R_M$ 与参数 ${\rho }_M,{\epsilon }_M,{\gamma }_M,{\mu }_M,{\mu }_A,{\mu }_F$ 相关。通过模拟得到 $R_M$ 的敏感性分析，如图2(a)所示，可见 $R_M$ 与参数 ${\rho }_M,{\gamma }_M$ 正相关， $R_M$ 与参数 ${\mu }_M,{\mu }_A,{\mu }_F$ 负相关，参数 $R_M$ 与 ${\epsilon }_M$ 几乎不相关。交配后雌蚊成虫产生水生蚊子的繁殖率 ${\rho }_M$ 、从水生蚊子到成虫蚊子的孵化率 ${\gamma }_M$ 与 $R_M$ 高度正相关，其相关系数分别为+0.95558，+0.71480；而成虫雄蚊的自然死亡率 ${\mu }_M$ 、成虫雌蚊子的自然死亡率 ${\mu }_F$ 、水生蚊子死亡率 ${\mu }_A$ 与 $R_M$ 高度负相关，其相关系数分别为−0.95791，−0.74305，−0.64901。

故降低水生蚊子到成蚊孵化率 ${\gamma }_M$，采用技术降低雌蚊的繁殖率 ${\rho }_M$，加大措施增加水生蚊子和成蚊的死亡率，都可降低蚊虫基本再生数，并且从相关系数大小可以看出，限制蚊子的最佳策略是从幼虫阶段开始。故这与现在采取的消杀蚊虫措施相吻合，即及时清除家里积水、污水以彻底消灭幼蚊的生存空间，这些都是最好、最实际的消杀方法。

(a) (b)

Figure 3. (a) Sensitivity Analysis of the basic reproduction number for outbreaks of WNv $R_0$ ; (b) Uncertainty analysis of the basic reproduction number for outbreaks of WNv $R_0$

图3. (a) 西尼罗河病毒暴发基本再生数 $R_0$ 的敏感性分析；(b) 西尼罗河病毒暴发基本再生数 $R_0$ 的不确定性分析

从 $R_0$ 的表达式可看出，WNv疾病基本再生数 $R_0$ 与参数 ${\rho }_M,{\epsilon }_M,{\gamma }_M,{\mu }_M,{\mu }_A,{\mu }_F,b_M,\beta ,{\beta }_{MB},{\mu }_B,d_B,{\Pi }_B,K$ 相关。通过模拟得到 $R_0$ 的敏感性分析，如图3(a)所示，可见 $R_0$ 与参数 $b_M,\beta ,{\rho }_M,{\gamma }_M,{\beta }_{MB},K,{\mu }_B$ 正相关，与参数 ${\mu }_F,{\mu }_M,d_B,{\epsilon }_M,{\mu }_A,{\Pi }_B$ 负相关。蚊子对鸟类的叮咬率 $b_M$ 、从感染的鸟类到易感蚊子的染病概率 $\beta$ 、交配后雌蚊成虫产生水生蚊子的繁殖率 ${\rho }_M$ 、从水生蚊子到成虫蚊子的孵化率 ${\gamma }_M$ 与 $R_0$ 高度正相关，其相关系数分别为：0.97365，0.79140，0.63339，0.52936。成虫雌蚊子的自然死亡率 ${\mu }_F$ 、成虫雄蚊子的自然死亡率 ${\mu }_M$ 、鸟类因病致死率 $d_B$ 与 $R_0$ 高度负相关，其相关系数分别为：−0.70285，−0.59145，−0.51941。从染病蚊子到易感鸟类的染病概率 ${\beta }_{MB}$ 、成虫转化为雄蚊的比率 ${\epsilon }_M$ 、蚊虫环境容纳量K、水生蚊子死亡率 ${\mu }_A$ 、鸟类种群补给率 ${\Pi }_B$ 、鸟类自然死亡率 ${\mu }_B$ 与 $R_0$ 弱相关，相关系数分别为：0.45959，−0.36667，0.2541，−0.24373，−0.20504，0.07714。

基本再生数是流行病学中的重要参数。一般情况下，希望可采取一些措施使基本再生数降低到1以下来进行疾病预防和控制。从上述分析可得，当蚊子对鸟类的叮咬率 $b_M$ 、从感染的鸟类到易感蚊子的染病概率 $\beta$ 、交配后雌蚊成虫产生水生蚊子的繁殖率 ${\rho }_M$ 、从水生蚊子到成虫蚊子的孵化率 ${\gamma }_M$ 越小，成虫雌蚊子的自然死亡率 ${\mu }_F$ 、成虫雄蚊子的自然死亡率 ${\mu }_M$ 、蚊虫环境容纳量K越大， $R_0$ 越小。

5. 结论

本文建立了一类结合鸟类–蚊子种群结构的西尼罗河病毒(WNv)传播模型，分析了模型的动力学性态，通过LHS/PRCC方法，分别对WNv传播模型的蚊虫再生数和WNv疾病暴发再生数进行了敏感性分析，得出各参数对模型性态的敏感程度，对防蚊灭蚊和控制病毒传播措施提供了理论支撑。由上述分析可得，传染病暴发后控制WNv传播最好办法是降低媒介蚊子的数量，而限制蚊子的最佳策略是从幼虫阶段开始，即保持环境清洁，消除积水和脏水措施，消灭水生蚊虫孳生地，将有效地抑制病毒的传播。注意到未受感染鸟类的补给率的下降将有益于WNv的传播，这一结果显示在WNv暴发期间时，利用控制鸟类的数量来抑制其传播是危险的。相反，应该适当加大未感染的鸟类的投放。

基金项目

浙江省教育厅一般科研项目Y201840260，绍兴文理学院校级科研计划项目2018LG1008。