脉冲释放沃尔巴克氏体感染雄蚊模型动力学分析

Dynamic Analysis of the Wolbachia-Infected Male Mosquitoes Model by Pulsed Release

摘要: 本文以昆虫不相容技术为背景研究了具有阶段和性别结构的野生与不育蚊子相互作用脉冲模型。其目的是研究采用周期释放不育蚊子的方法控制野生蚊子数量的可行性。首先通过动力学分析得到了系统平凡周期解的存在性，并分别利用Floquet理论和Lyapunov稳定性定理给出了相应的局部稳定性条件和全局稳定性条件，从而验证了脉冲释放不育蚊子可以使野生蚊子灭绝。其次得到了系统非平凡周期解的存在性，并在具体参数下讨论了其局部稳定性，发现系统在一定阈值条件下出现平凡周期解和非平凡周期解共存的双稳现象，这表明脉冲控制也可以在不消灭野生蚊子的情况下控制它们的数量。最后利用数值模拟验证了相关理论结果。

Abstract: In this paper, we study a wild and sterile mosquito interaction impulsive model with stage-structure and sex-structure based on insect incompatibility technology. The aim is to study the feasibility of controlling the number of wild mosquitoes by releasing sterile mosquitoes periodically. Firstly, the existence of trivial periodic solutions is obtained by dynamic analysis, and the corresponding local stability and global stability conditions are proved by Floquet theory and Lyapunov stability theorem respectively. The results show that the pulsed release of sterile mosquitoes could make the wild mosquitoes extinct. Then, the existence of nontrivial periodic solutions is obtained, and its local stability is discussed under specific parameter values. It is found that the system has the bistable phenomenon in which trivial periodic solution and non-trivial periodic solution coexist under certain threshold conditions. This shows that the population of mosquitoes can be controlled by pulse control without eliminating them. Finally, numerical simulation verifies the theoretical results.

1. 引言

多年来蚊媒传染病一直危害着人类的健康，雌性蚊子可以通过叮咬传播细菌，病毒或寄生虫等，使易感者患病。蚊媒传染病每年给多个国家的大量人员带来感染和死亡，严重危害人们的身体健康和经济。蚊媒传染病种类繁多，常见的有流行性乙型脑炎、疟疾、登革热、黄热病等，其特点是发病率很高且难以利用疫苗或其他方法有效克服。目前控制蚊子的数量是最有效的方法 [1] [2]。

为了减少蚊虫数量，人们曾致力于大范围使用杀虫剂，然而这种方法不仅会对环境造成严重的污染，还会使蚊子产生耐药性。不育昆虫释放技术(Sterile Insect Technique, SIT)是一种通过释放辐射不育昆虫，使与其交配的野生昆虫不能产生后代，SIT已经在多种农业和畜牧业的害虫控制中成功应用 [3]。然而该项技术存在一定的缺陷，大量的辐射常常会使不育蚊子的交配竞争力和生存能力降低。

基于沃尔巴克氏体感染的昆虫不相容技术(Incompatible Insect Technique, IIT)弥补了这项缺陷。沃尔巴克氏体(Wolbachia)是一种母系遗传的胞内共生菌，当释放的携菌雄蚊与野生的不携菌雌蚊交配时，会产生胞质不相容(Cytoplasmic Incompatibility, CI)现象，使其产的卵不能孵化，同时沃尔巴克氏体的感染对蚊子的交配竞争力和生存能力影响极小。但由于同时携菌的雌雄蚊子可以产生携菌的后代，这可能会产生携菌蚊子种群完全替代原生蚊子种群的风险。为了避免种群替代的风险，通常会在不育蚊子释放前严格区分雌雄或在不育蚊子的生产中联合使用SIT和IIT [4]。例如2016至2017年，中山大学研究团队在中国登革热主要流行地，广州2个独立岛的居民居住区域，针对白纹伊蚊利用SIT-IIT相结合的技术进行了开放式田间试验。成功使蚊卵数年均下降了94%，成蚊数年均下降了83%~94%，而且连续13周没有发现活卵，连续6周没有监测到成蚊 [5]。

不育蚊子技术完备的同时利用一些生物数学方法，比如建立和解析数学模型，能够帮助生物学家解释害虫防治的很多棘手问题。在进行实践和现场试验的同时，近几十年人们提出了一些动力学模型来研究野生与不育蚊子相互作用下种群的动态和控制 [6] - [17]，这些研究通常有助于了解蚊子种群变化特征，以及不育蚊子的释放对野生蚊子种群的影响，为不育蚊子释放策略的制定提供理论依据。

文献 [11] 以常微分方程为基础建立了野生与不育蚊子相互作用连续时间模型，模型包含了蚊子种群的阶段性结构和性别结构。本文在此基础上建立脉冲模型，通过对系统的动力学分析，讨论脉冲控制下不育蚊子对野生蚊子数量动态的影响。我们在文章的第二部分，描述了脉冲模型的建立；在文章的第三部分，对模型进行动力学分析得到了系统周期解的存在性及稳定性；在文章的第四部分，对所得理论结果进行了数值模拟。

2. 脉冲模型

本节中，我们主要叙述昆虫不相容技术下野生与不育蚊子相互作用脉冲模型的建立过程。

首先考虑不存在不育蚊子释放的情况。蚊子的生命周期分为四个阶段，包括虫卵，幼虫，虫蛹，成虫。处于前三个阶段的蚊子不产生后代且依赖于水生环境，会因为水生环境的拥挤产生依赖于密度的死亡率，于是将这三个水生阶段的蚊子定义为“蚊子幼虫”，用L表示，而处于最后一个阶段的蚊子会飞离水生环境，不受环境容量的限制，同时雌性成虫会产卵，而雄性成虫不产卵，于是将这个阶段的蚊子定义为“雌性成虫”和“雄性成虫”，分别用F和M表示。则可以建立一个具有阶段和性别结构的蚊子种群模型 [11]

$\{\begin{array}{l}\frac{\text{d}L\left(t\right)}{\text{d}t}=bF-\beta L-d_{1}L-\kappa L^2,\\ \frac{\text{d}F\left(t\right)}{\text{d}t}=\eta \beta L-d_{2}F,\\ \frac{\text{d}M\left(t\right)}{\text{d}t}=\left(1-\eta \right)\beta L-d_{3}M,\end{array}$ (2.1)

其中b为雌性成虫产卵率； $\beta$ 为蚊子幼虫成熟率； $d_1$ 为蚊子幼虫自然死亡率； $\kappa L$ 表示依赖于密度线性变化的蚊子幼虫死亡率； $\eta$ 为雌性比例； $d_2$ 为雌性成虫自然死亡率； $d_3$ 为雄性成虫自然死亡率。

接着考虑存在不育蚊子释放的情况。假设 $M_w\left(t\right)$ 为不育蚊子在t时刻的数量，这部分蚊子在IIT作用下不存在自然出生率，它们的出生率即不育蚊子的释放率。于是不育蚊子的种群数量变化可表达如下：

$\frac{\text{d}{M}_w}{\text{d}t}=\theta (\cdot )-d_4{M}_w,$ (2.2)

其中 $\theta (\cdot )$ 为不育蚊子释放率； $d_4$ 为不育蚊子自然死亡率。

此时压倒性数量的不育雄蚊被释放进环境中，可能导致蚊子种群中雌雄比例悬殊，进而为雄蚊寻找配偶带来困难，因此我们在模型中引入Allee效应参数 $\gamma$。同时考虑沃尔巴克氏体引导CI概率q，雌性成虫的产卵率可以表示为两部分，一部分为不携菌蚊子之间的产卵率 $bFM/\left(M+M_w+\gamma \right)$，另一部分为CI失败产生的产率 $bF\left(1-q\right)M_w/\left(M+M_w+\gamma \right)$。则在不育蚊子作用下野生蚊子出生率可表达为

$bF\left(\frac{M+\left(1-q\right)M_w}{M+M_w+\gamma }\right).$ (2.3)

大多数野生与不育蚊子相互作用模型基于完全连续的动力系统，在描述不育蚊子释放过程时通常使用一个连续的增长率，比如最常见的释放方式经常假设释放率为常数，既 $\theta (\cdot )=\theta$。这相当于假设不育蚊子的释放是不间断的。然而，现实中蚊子的释放过程是难以做到连续的。不育蚊子的释放通常是分多次执行的，例如2017年美国弗雷斯诺县曾在20周中每周释放100万只不育蚊子 [18]。不育蚊子的每次释放表现为数量的一次瞬间变化，这会使不育蚊子数量出现脉冲式的变化。因此，相比于连续系统，脉冲系统更能合理的描述整个释放过程。不育蚊子脉冲式释放受到了人们的广泛关注 [19] [20] [21] [22]。于是我们假设不育蚊子释放每隔固定时间T执行一次，且每次释放量相同。则不育蚊子的种群数量变化可利用脉冲微分方程表达

$\{\begin{array}{l}\frac{\text{d}{M}_w}{\text{d}t}=-d_4{M}_w,\\ M_w\left(n{T}^{+}\right)=M_w\left(nT\right)+\theta .\end{array}$ (2.4)

根据式(2.1)，(2.3)，(2.4)我们给出昆虫不相容技术下不育雄蚊与野生蚊子相互作用脉冲模型

$\{\begin{array}{l}\frac{\text{d}L\left(t\right)}{\text{d}t}=bF\left(\frac{M+\left(1-q\right)M_w}{M+M_w+\gamma }\right)-\beta L-d_{1}L-\kappa L^2,\\ \frac{\text{d}F\left(t\right)}{\text{d}t}=\eta \beta L-d_{2}F,\\ \frac{\text{d}M\left(t\right)}{\text{d}t}=\left(1-\eta \right)\beta L-d_{3}M,\\ \frac{\text{d}{M}_w\left(t\right)}{\text{d}t}=-d_4{M}_w,\\ M_w\left(n{T}^{+}\right)=M_w\left(nT\right)+\theta \end{array}$ (2.5)

其中T为不育蚊子释放周期； $\theta$ 为不育蚊子每次释放量。

3. 动力学分析

3.1. 平凡周期解及其稳定性

3.1.1. 平凡周期解存在性

首先讨论当野生蚊子被全部消灭时系统所对应平凡周期解的存在性。

当 $L\left(t\right)=0$， $F\left(t\right)=0$， $M\left(t\right)=0$ 时系统(2.5)可简化成

$\{\begin{array}{l}\frac{\text{d}L\left(t\right)}{\text{d}t}=0,\\ \frac{\text{d}F\left(t\right)}{\text{d}t}=0,\\ \frac{\text{d}M\left(t\right)}{\text{d}t}=0,\\ \frac{\text{d}{M}_w\left(t\right)}{\text{d}t}=-d_4{M}_w,\\ M_w\left(n{T}^{+}\right)=M_w\left(nT\right)+\theta .\end{array}$ (3.1)

由系统(3.1)后两式得到

$M_w\left[\left(n+1\right)T^{+}\right]=M_w\left(n{T}^{+}\right){\text{e}}^{-d_4\left(t\mathrm{mod}T\right)},$

这个离散系统满足当 $n\to \infty$ 时

$M_w^{\ast }=\frac{\theta }{1-{\text{e}}^{-d_{4}T}},$

则系统(2.5)存在对应于野生蚊子灭绝的平凡周期解

$\left(\bar{L}\left(t\right),\bar{F}\left(t\right),\bar{M}\left(t\right),{\bar{M}}_w\left(t\right)\right)=\left(0,0,0,M_w^{\ast }{\text{e}}^{-d_4\left(t\mathrm{mod}T\right)}\right).$ (3.2)

3.1.2. 平凡周期解局部稳定性

我们在这部分讨论系统(2.5)平凡周期解的局部稳定性。

取系统(2.5)的一个小幅度扰动 $\tilde{L}=L-\bar{L}$， $\tilde{F}=F-\bar{F}$， $\tilde{M}=M-\bar{M}$， ${\tilde{M}}_w=M_w-{\bar{M}}_w$，此时脉冲项消失：

${\tilde{M}}_w\left(n{T}^{+}\right)=\left[M_w\left(nT\right)+\theta \right]-\left[{\bar{M}}_w\left(nT\right)+\theta \right]={\tilde{M}}_w\left(nT\right).$

系统变成

$\{\begin{array}{l}\frac{\text{d}\tilde{L}\left(t\right)}{\text{d}t}=b\tilde{F}\left(\frac{\tilde{M}+\left(1-q\right)\left({\tilde{M}}_w+{\bar{M}}_w\right)}{\tilde{M}+\left({\tilde{M}}_w+{\bar{M}}_w\right)+\gamma }\right)-\beta \tilde{L}-d_1\tilde{L}-\kappa {\tilde{L}}^2,\\ \frac{\text{d}\tilde{F}\left(t\right)}{\text{d}t}=\eta \beta \tilde{L}-d_2\tilde{F},\\ \frac{\text{d}\tilde{M}\left(t\right)}{\text{d}t}=\left(1-\eta \right)\beta \tilde{L}-d_3\tilde{M},\\ \frac{\text{d}{\tilde{M}}_w\left(t\right)}{\text{d}t}=-d_4{\tilde{M}}_w.\end{array}$ (3.3)

将上式记作

$\left(\begin{array}{c}\tilde{L}\left(t\right)\\ \tilde{F}\left(t\right)\\ \tilde{M}\left(t\right)\\ {\tilde{M}}_w\left(t\right)\end{array}\right)=\Phi \left(t\right)\left(\begin{array}{c}\tilde{L}\left(0\right)\\ \tilde{F}\left(0\right)\\ \tilde{M}\left(0\right)\\ {\tilde{M}}_w\left(0\right)\end{array}\right),$

这里 $\Phi \left(t\right)$ 满足

$\frac{\text{d}\Phi \left(t\right)}{\text{d}t}=\Lambda \left(t\right)\Phi \left(t\right),$

其中

$\Lambda \left(t\right)=\left(\begin{array}{cccc}-\beta -d_1-2\kappa \tilde{L}& b\left(\frac{\tilde{M}+\left(1-q\right)\left({\tilde{M}}_w+{\bar{M}}_w\right)}{\tilde{M}+\left({\tilde{M}}_w+{\bar{M}}_w\right)+\gamma }\right)& bF\left(\frac{\gamma +q\left({\tilde{M}}_w+{\bar{M}}_w\right)}{{\left[\tilde{M}+\left({\tilde{M}}_w+{\bar{M}}_w\right)+\gamma \right]}^2}\right)& bF\left(\frac{\left(1-q\right)\gamma -q\tilde{M}}{{\left[\tilde{M}+\left({\tilde{M}}_w+{\bar{M}}_w\right)+\gamma \right]}^2}\right)\\ \eta \beta & -d_2& 0& 0\\ \left(1-\eta \right)\beta & 0& -d_3& 0\\ 0& 0& 0& -d_4\end{array}\right).$

当 $\tilde{L}=0$， $\tilde{F}=0$， $\tilde{F}=0$ 时

$\Lambda \left(t\right)=\left(\begin{array}{cccc}-\beta -d_1& b\left(\frac{\left(1-q\right){\bar{M}}_w}{\bar{M}{}_w+\gamma }\right)& 0& 0\\ \eta \beta & -d_2& 0& 0\\ \left(1-\eta \right)\beta & 0& -d_3& 0\\ 0& 0& 0& -d_4\end{array}\right),$

并且

$\left(\begin{array}{c}\tilde{L}\left(n{T}^{+}\right)\\ \tilde{F}\left(n{T}^{+}\right)\\ \tilde{M}\left(n{T}^{+}\right)\\ {\tilde{M}}_w\left(n{T}^{+}\right)\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}\tilde{L}\left(nT\right)\\ \tilde{F}\left(nT\right)\\ \tilde{M}\left(nT\right)\\ {\tilde{M}}_w\left(nT\right)\end{array}\right),$

单值矩阵 [23] [24] [25] 为

$\Gamma =\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\Phi \left(T\right)=\Phi \left(T\right).$

可知 $\Phi \left(T\right)=\Phi \left(0\right)\exp \left({\displaystyle {\int }_0^{T}tr\Lambda \left(t\right)\text{d}t}\right)$， $\Phi \left(0\right)$ 是单位矩阵。记 $\exp \left({\displaystyle {\int }_0^{T}tr\Lambda \left(t\right)\text{d}t}\right)$ 的特征值为 ${\lambda }_1$， ${\lambda }_2$， ${\lambda }_3$， ${\lambda }_4$，相应的特征多项式为

$\left(\lambda +d_3\right)\left(\lambda +d_4\right)\left({\lambda }^2+A\lambda +B\right)=0,$

其中

$\begin{array}{l}A=\beta +d_1+d_2>0,\\ B=\frac{d_2\gamma \left(\beta +d_1\right)\left[\gamma -\left(r_0-1\right){\bar{M}}_w\right]}{{\bar{M}}_w+\gamma }.\end{array}$

这里

$r_0=\frac{b\eta \beta \left(1-q\right)}{d_2\left(\beta +d_1\right)}.$ (3.4)

显然 ${\lambda }_1=-d_3<0$， ${\lambda }_2=-d_4<0$，下面我们讨论 ${\lambda }_3$， ${\lambda }_4$ 的符号。

当 $r_0<1$ 时 $B>0$，由一元二次方程韦达定理可知 ${\lambda }_3+{\lambda }_4=-A<0$， ${\lambda }_3{\lambda }_4=B>0$，则 ${\lambda }_3<0$， ${\lambda }_4<0$。当 $r_0>1$ 时，我们定义一个阈值

$M_w^0=\frac{\gamma }{r_0-1}.$ (3.5)

${\bar{M}}_{w\max }<M_w^0$ 时 $B>0$，可知 ${\lambda }_3<0$， ${\lambda }_4<0$， ${\bar{M}}_{w\min }>M_w^0$ 时 $B<0$，可知 ${\lambda }_3{\lambda }_4<0$。则由Floquet理论，当 $r_0<1$ 或当 $r_0>1$ 且满足阈值条件 ${\bar{M}}_{w\max }<M_w^0$ 时，系统(2.5)的平凡周期解是局部渐近稳定的。而当 $r_0>1$ 且 ${\bar{M}}_{w\min }>M_w^0$ 时，系统(2.5)的平凡周期解是不稳定的。

定理3.1系统(2.5)平凡周期解 $\left(0,0,0,M_w^{\ast }{\text{e}}^{-d_4\left(t\mathrm{mod}T\right)}\right)$ 的局部稳定性结论如下：

1) 如果 $r_0<1$，则总是局部渐近稳定的；

2) 如果 $r_0>1$，则在 ${\bar{M}}_{w\max }<M_w^0$ 时是局部渐近稳定的，而在 ${\bar{M}}_{w\min }>M_w^0$ 时是不稳定的。

3.1.3. 平凡周期解全局稳定性

定理3.2当 $\frac{b\eta \beta }{d_2\left(\beta +d_1\right)}\le 1$ 时，系统(2.5)的平凡周期解 $\left(0,0,0,M_w^{\ast }{\text{e}}^{-d_4\left(t\mathrm{mod}T\right)}\right)$ 是全局渐近稳定的。

证明：构造李雅普诺夫函数 $V=d_{2}L+bF$，则

$\begin{array}{c}{V}^{\prime }=d_2\left[bF\left(\frac{M+\left(1-q\right)M_w}{M+M_w+\gamma }\right)-\beta L-d_{1}L-\kappa L^2\right]+b\left(\eta \beta L-d_{2}F\right)\\ \le \left(b\eta \beta -d_2\left(\beta +d_1\right)\right)L\\ =d_2\left(\beta +d_1\right)\left(r_1-1\right)L.\end{array}$

其中

$r_1=\frac{b\eta \beta }{d_2\left(\beta +d_1\right)}.$ (3.6)

可知当 $r_1\le 1$ 时， $V^{\prime }\left(t\right)\le 0$，且仅当 $L=0$， $F=0$， $M=0$ 时 $V^{\prime }\left(t\right)=0$。所以 $r_1\le 1$ 时，系统(2.5)的平凡周期解是全局渐近稳定的。

3.2. 非平凡周期解及其稳定性

3.2.1. 非平凡周期解存在性

在实践中有时我们可以不完全消灭蚊子种群，而是将它们的数量控制在一定范围之内。基于这个目的，我们继续研究系统(2.5)对应于野生蚊子部分存活的非平凡周期解。下面继续利用系统(3.3)的正平衡点讨论系统(2.5)非平凡周期解的存在性。系统(3.3)在正平衡点满足

$\begin{array}{l}b\tilde{F}\left(\frac{\tilde{M}+\left(1-q\right)\left({\tilde{M}}_w+{\bar{M}}_w\right)}{\tilde{M}+\left({\tilde{M}}_w+{\bar{M}}_w\right)+\gamma }\right)-\beta \tilde{L}-d_1\tilde{L}-\kappa {\tilde{L}}^2=0,\\ \eta \beta \tilde{L}-d_2\tilde{F}=0,\\ \left(1-\eta \right)\beta \tilde{L}-d_3\tilde{M}=0,\\ {\tilde{M}}_w=0.\end{array}$ (3.7)

由式(3.6)前两项得出关于 $\tilde{J}$ 的方程，记为

$\frac{\left(A_1{\tilde{L}}^2+B_1\tilde{L}+C_1\right)\tilde{L}}{d_2\left(\left(1-\eta \right)\beta \tilde{L}+d_3{\bar{M}}_w+d_m\gamma \right)}=0,$ (3.8)

其中

$\begin{array}{l}{A}_1=\kappa d_2\left(1-\eta \right)\beta >0,\\ B_1=\kappa d_2{d}_3\left({\bar{M}}_w+\gamma \right)-\beta d_2\left(\beta +d_1\right)\left(1-\eta \right)\left(r_1-1\right),\\ C_1=d_3{d}_2\left(\beta +d_1\right)\left(\gamma -\left(r_0-1\right){\bar{M}}_w\right).\end{array}$

假设 ${\bar{M}}_{w\max }<M_w^0$， $M_w^0$ 由式(3.5)给出，则 $A_1>0$， $C_1>0$，此时系统(3.3)正平衡点的存在性取决于 $B_1$ 的符号。若 $r_1\le 1$，则 $B_1>0$，系统(3.3)不存在正平衡点，对应于系统(2.5)不存在正的非平凡周期解。若 $r_1>1$，则 $B_1$ 的符号取决于 ${\bar{M}}_w$，由此我们定义一个系统(3.3)正平衡点存在性阈值

$M_w^1=\frac{d_2\beta \left(\beta +d_1\right)\left(1-\eta \right)\left(r_1-1\right)-\kappa d_2{d}_3\gamma }{\kappa d_2{d}_3},$ (3.9)

如果 ${\bar{M}}_{w\max }<M_w^1$，则 $B_1<0$，系统(3.3)存在0，1或2个正平衡点，对应于系统(2.5)存在0，1或2个正的非平凡周期解，如果 ${\bar{M}}_{w\max }>M_w^1$，则 $B_1>0$，系统(3.3)不存在正平衡点，对应于系统(2.5)不存在正的非平凡周期解。

下面我们假设 ${\bar{M}}_{w\max }<M_w^0$ 且 ${\bar{M}}_{w\max }<M_w^1$，即系统(2.5)存在0，1或2个正的非平凡周期解。式(3.7)正根的个数取决于 $\Delta =B_1^2-4A_1{C}_1$，若 $\Delta >0$ 式(3.7)存在两个不相同的正根，若 $\Delta <0$ 式(3.7)不存在正根，而由于 $\Delta$ 的值依赖于周期变化的 ${\bar{M}}_w$，所以 $\Delta$ 的值不能恒等于0。由此定义一个正平凡周期解存在个数阈值

$R_0^c=\frac{B_1{}^2}{4A_1{C}_1}.$ (3.10)

如果 $R_0^c>1$，则系统(2.5)存在两个正的非平凡周期解，如果 $R_0^c<1$ 则系统(2.5)不存在正的非平凡周期解。

定理3.2假设 ${\bar{M}}_{w\max }<M_w^0$，则系统(2.5)的正非平凡周期解结论总结如下

其中 $r_1$ 由式(3.6)给出， $M_w^1$ 由式(3.9)给出， $R_0^c$ 由式(3.10)给出。

3.2.2. 非平凡周期解的进一步讨论

下面假设 ${\bar{M}}_{w\max }<M_w^0$， ${\bar{M}}_{w\max }<M_w^1$ 且 $R_0^c>1$，利用Floquet理论讨论非平凡周期解的稳定性。此时系统(3.3)存在两个正平衡点

$\begin{array}{l}{E}_1=\left(\frac{-B_1-\sqrt{\Delta }}{2A_1},\frac{\eta \beta {\tilde{L}}_1}{d_2},\frac{\left(1-\eta \right)\beta {\tilde{L}}_1}{d_3},{\bar{M}}_w\right),\\ E_2=\left(\frac{-B_1+\sqrt{\Delta }}{2A_1},\frac{\eta \beta {\tilde{L}}_2}{d_2},\frac{\left(1-\eta \right)\beta {\tilde{L}}_2}{d_3},{\bar{M}}_w\right),\end{array}$

对应系统(2.5)存在两个正的非平凡周期解 $P_i=\left({\bar{L}}_i,{\bar{F}}_i,{\bar{M}}_i,{\bar{M}}_w\right)$， $i=1,2$。

类似定理3.1的证明，记 $\exp \left({\displaystyle {\int }_0^T\Lambda \left(t\right)\text{d}t}\right)$ 的特征值为 ${\lambda }_{11}$， ${\lambda }_{12}$， ${\lambda }_{13}$， ${\lambda }_{14}$。扰动系统(3.3)存在两个正平衡点 $E_i=\left({\bar{L}}_i\left(t\right),{\bar{F}}_i\left(t\right),{\bar{M}}_i\left(t\right),{\bar{M}}_w\left(t\right)\right)$， $i=1,2$ 并且在正平衡点满足 ${\tilde{M}}_w=0$。此时

${\Lambda }_1\left(t\right)=\left(\begin{array}{cccc}-\beta -d_1-2\kappa {\bar{L}}_i& b\left(\frac{{\bar{M}}_i+\left(1-q\right)}{{\bar{M}}_i+{\bar{M}}_w+\gamma }\right)& b{\bar{F}}_i(\frac{\gamma +q+{\bar{M}}_w}{{\left({\bar{M}}_i+{\bar{M}}_w+\gamma \right)}^2})& b{\bar{F}}_i\left(\frac{\left(1-q\right)\gamma -q{\bar{M}}_i}{{\left({\bar{M}}_i+{\bar{M}}_w+\gamma \right)}^2}\right)\\ \eta \beta & -d_2& 0& 0\\ \left(1-\eta \right)\beta & 0& -d_3& 0\\ 0& 0& 0& -d_4\end{array}\right),$

记 ${\Lambda }_1\left(t\right)$ 的特征值为 ${\bar{\lambda }}_{11}$， ${\bar{\lambda }}_{12}$， ${\bar{\lambda }}_{13}$， ${\bar{\lambda }}_{14}$，其特征多项式为

$\left({\bar{\lambda }}_1-d_4\right)\left({\bar{\lambda }}_1{}^3+A_1{\bar{\lambda }}_1{}^2+B_1{\bar{\lambda }}_1+C_1\right)=0,$ (3.11)

其中

$\begin{array}{l}{A}_2=d_1+d_2+d_3+\beta +2\kappa {\tilde{L}}_i>0,\\ B_2=\left(\beta d_2+d_1{d}_2+d_2{d}_3+\beta d_3+d_1{d}_3+2\kappa {\tilde{L}}_i{d}_2+2\kappa {\tilde{L}}_i{d}_3\right)-\frac{b\eta \beta \left({\tilde{M}}_i+\left(1-q\right){\bar{M}}_w\right)}{{\tilde{M}}_i+r+{\bar{M}}_w}-\frac{b\eta \beta \tilde{F}\left(r+q{\bar{M}}_w\right)}{{\left({\tilde{M}}_i+r+{\bar{M}}_w\right)}^2},\\ C_2=\beta d_2{d}_3+d_1{d}_2{d}_3+2\kappa {\tilde{L}}_i{d}_2{d}_3-\frac{b\eta \beta d_3\left({\tilde{M}}_i+\left(1-q\right){\bar{M}}_w\right)}{{\tilde{M}}_i+r+{\bar{M}}_w}-\frac{b\eta \beta {\tilde{F}}_i{d}_2\left(r+q{\bar{M}}_w\right)}{{\left({\tilde{M}}_i+{\bar{M}}_w+r\right)}^2},\end{array}$

则 $\exp \left({\displaystyle {\int }_0^{T}tr{\Lambda }_1\left(t\right)\text{d}t}\right)$ 的特征值为

${\lambda }_{11}={\text{e}}^{{\displaystyle {\int }_0^T{\bar{\lambda }}_{11}\text{d}t}},{\lambda }_{12}={\text{e}}^{{\displaystyle {\int }_0^T{\bar{\lambda }}_{12}\text{d}t}},{\lambda }_{13}={\text{e}}^{{\displaystyle {\int }_0^T{\bar{\lambda }}_{13}\text{d}t}},{\lambda }_{14}={\text{e}}^{{\displaystyle {\int }_0^T{\bar{\lambda }}_{14}\text{d}t}}.$

易知 ${\lambda }_{11}<1$，而 ${\lambda }_{12}$， ${\lambda }_{13}$， ${\lambda }_{14}$ 的符号取决于 $B_2\left({\tilde{L}}_i\right)$ 和 $C_2\left({\tilde{L}}_i\right)$。比如，当 $B_2{\left({\tilde{L}}_i\right)}_{\min }>0$ 且 $C_2{\left({\tilde{L}}_i\right)}_{\min }>0$ 时，根据一元三次方程根与系数的关系可知式(3.11)的根都有负实部，进而得到 ${\lambda }_{12}$， ${\lambda }_{13}$， ${\lambda }_{14}$ 均小于1，此时对应的非平凡周期解则是稳定的。由于计算过于复杂，无法给出具体稳定性条件。我们将在数值模拟部分赋予具体参数，举例说明非平凡周期解的稳定性。

4. 数值模拟

前面我们给出了昆虫不相容技术下不育雄蚊与野生蚊子相互作用脉冲模型周期解存在性及稳定性分析，给出了平凡周期解局部稳定性及全局稳定性条件，以及非平凡周期解存在性阈值 $M_w^1$ 和存在个数阈值 $R_0^c$。

下面我们在具体参数下举例说明非平凡周期解 $P_i=\left({\bar{L}}_i\left(t\right),{\bar{F}}_i\left(t\right),{\bar{M}}_i\left(t\right),{\bar{M}}_w\left(t\right)\right)$， $i=1,2$ 的稳定性，并数值模拟理论结果。参数取值为

$b=0.8,q=0.9,\beta =\frac{1}{6},\kappa =0.001,d_1=\frac{1}{12},d_2=\frac{1}{13},d_3=\frac{1}{14},d_4=\frac{1}{14}.$

首先假设 $\theta =3$， $T=14$，在此参数条件下 $r_1=3.4667>1$， ${\bar{M}}_{w\max }<M_w^1=719.3444$， $R_{0\min }^c=15.3621>1$，系统(2.5)存在两个正的非平凡周期解。

此时，经计算 $B_2{\left({\tilde{L}}_1\right)}_{\min }>0$， $C_2{\left({\tilde{L}}_1\right)}_{\min }>0$，则 ${\bar{\lambda }}_{12}<0$， ${\bar{\lambda }}_{13}<0$， ${\bar{\lambda }}_{14}<0$，进而得到 ${\lambda }_{12}<1$， ${\lambda }_{13}<1$， ${\lambda }_{14}<1$，由Floquet理论可知平凡周期解 $P_2=\left({\bar{L}}_2\left(t\right),{\bar{F}}_2\left(t\right),{\bar{M}}_2\left(t\right),{\bar{M}}_w\left(t\right)\right)$ 是局部渐近稳定的。分析发现在该参数条件下系统(3.3)同时存在0平衡点 $E_0$ 和正平衡点 $E_1$， $E_2$，且 $E_0$， $E_2$ 是局部渐近稳定的，而 $E_1$ 是不稳定的，对应于系统(2.5)同时存在平凡周期解 $P_0$，和正的非平凡周期解 $P_1$， $P_2$，且 $P_0$， $P_2$ 是局部渐近稳定的，而 $P_1$ 是不稳定的，此时系统(2.5)出现脉冲双稳现象。

如图1的(a)~(d)分别为当 $\theta =3$， $T=14$ 时蚊子幼虫 $L\left(t\right)$，雌性成虫 $F\left(t\right)$，雄性成虫 $M\left(t\right)$ 和不育蚊子 $M_w\left(t\right)$ 数量随时间变化曲线。可见 $P_2=\left({\bar{L}}_2\left(t\right),{\bar{F}}_2\left(t\right),{\bar{M}}_2\left(t\right),{\bar{M}}_w\left(t\right)\right)$ 为局部渐近稳定的非平凡周期解。如图2(a)~(c)我们设置其他参数不变，取 $\theta$ 作为分支参数模拟系统(2.5)产生的分支现象。可见系统(2.5)在该参数条件下出现脉冲双稳现象。如图3(a)，设置 $\theta =30$， $T=1$，此时 ${\bar{M}}_{w\max }=435.1786<M_w^1$， $R_0^c<1$，由定理3.3，系统(2.5)不存在正的非平凡周期解。

接下来我们重新设置参数为

$b=\frac{1}{3},q=0.9,\beta =\frac{1}{6},\kappa =0.001,d_1=\frac{1}{7},d_2=\frac{1}{10},d_3=\frac{1}{11},d_4=\frac{1}{11}.$

如图3(b)，此时 $r_1=0.8974<1$，由定理3.3系统(2.5)的平凡周期解是全局渐近稳定的。

(a) (b) (c) (d)

Figure 1. (a)~(c) and (d) show the time-varying curves of mosquito larvae, female adults, male adults and sterile mosquitoes when T = 14, θ = 3

图1. 图(a)~(d) 分别为当T = 14, θ = 3时蚊子幼虫，雌性成虫，雄性成虫和不育蚊子数量随时间变化曲线

(a) (b) (c)

Figure 2. In figures (a), (b) and (c), the blue line represents the stable periodic solution P₀, P₂ and the black dotted line represents the unstable periodic solution P₁

图2. 图(a)~(c)蓝色线条表示稳定的周期解P₀和P₂，黑色虚线表示不稳定的周期解P₁

(a) (b)

Figure 3. The parameter condition of figure (a) makes ${\bar{M}}_{w\max }<M_w^1$, $R_0^c<1$, and there is no positive nontrivial periodic solution. The parameter condition of figure (b) makes $r_1<1$, and the nontrivial periodic solution of the system is globally asymptotically stable

图3. 图(a)的参数条件使 ${\bar{M}}_{w\max }<M_w^1$， $R_0^c<1$，此时系统不存在正的非平凡周期解。图(b)的参数条件使 $r_1<1$，此时系统的平凡周期解是全局渐近稳定的

5. 结论

本文在文献 [11] 基础上建立了昆虫不相容技术下不育雄蚊与野生蚊子相互作用脉冲模型。首先发现系统总是存在对应于蚊子种群灭绝的平凡周期解，并分别利用Floquet理论和Lyapunov稳定性定理给出了相应的局部稳定性条件和全局稳定性条件，从而验证了脉冲释放不育蚊子可以使野生蚊子灭绝。其次发现在满足一定的阈值条件时系统可以存在两个对应于蚊子种群部分存活的非平凡周期解，并在具体参数值下讨论了它们的局部稳定性，研究表明有且只有一个非平凡周期解是局部渐近稳定的，此时系统出现平凡周期解和非平凡周期解共存的双稳现象。这表明我们通过控制不育蚊子的释放量和释放周期既可以完全消灭野生蚊子，也可以在不消灭野生蚊子的情况下控制它们的数量。

基金项目

2020年吉林省教育厅“十三五”科学技术项目(JJKH20200726KJ)。