弱Allee效应下的不育蚊子的两种释放策略

Two Release Strategies of Sterile Mosquitoes under Weak Allee Effect

摘要: 在本文中，我们研究了弱Allee效应下不育蚊子释放技术(SIT)的两种释放策略模型，重点分析了不育蚊子的恒定释放策略和脉冲释放策略的可行性。针对恒定释放策略，我们得到了一个与野生蚊子种群灭绝相对应的边界平衡点的存在性与稳定性结果，并且探讨了代表野生蚊子和不育蚊子共存的内部平衡点的存在性和稳定性。与此同时，我们通过考察脉冲释放策略的平凡周期解存在性与稳定性以及非平凡解的存在性，进一步深入研究了脉冲释放策略的可行性。最后，数值模拟结果为在特定参数条件下存在非平凡周期解提供了有力证据，强烈显示出该系统中可能存在非平凡周期解现象。

Abstract: In this paper, we investigate two models of the Sterile Insect Technique (SIT) under weak Allee effect, and focus on two strategies of constant release and pulse release of sterile mosquitoes. For the constant release strategy, we establish the existence and stability of a boundary equilibrium corresponding to the extinction of wild mosquito populations, as well as an internal equilibrium representing the coexistence of both wild and sterile mosquito types. Correspondingly, we investigate the pulse release strategy by examining the existence and stability of its trivial periodic solution, as well as non-trivial solutions. Our numerical simulations provide compelling evidence for the existence of nontrivial periodic solutions under certain parameter conditions, strongly suggesting the presence of periodic solutions within the system.

1. 引言

蚊子传播的疾病，如疟疾[1]和登革热[2]，对全球公共卫生构成了重大威胁。蚊子是这些疾病的主要媒介，被列为四大害虫之一。目前，蚊子可以分为三类：按蚊、白纹伊蚊和库蚊。它们成为疟疾、登革热和黄热病等疾病传播的危险载体。全球预计有40亿人面临蚊媒疾病的风险，2025年8月底，登革热报告的病例超过1230万，几乎是2023年650万例的两倍[3]。

每年，许多社会资源被用于预防和控制蚊媒疾病的传播。最初的方法是使用化学农药，但其对环境的危害以及蚊子抵抗力的增加促使人们转向基于分子生物学的解决方案。因此，提出了不育昆虫技术(SIT)，即在实验室中对蚊子进行不育处理后释放到野外，在与野生蚊子的交配过程中干扰其正常的繁殖过程，从而控制野生蚊子的数量[4]。近年来，数学家们建立了各种数学模型来研究SIT对种群的影响，以提出更好的控制策略。Li通过微分方程系统描述了野生蚊子与不育蚊子之间的相互作用，建立了一个连续时间阶段结构的数学模型，考虑了模型中三种不同的不育蚊子释放策，研究了模型的动力学行为[5]，Huang [6]提出了两种脉冲释放的数学模型，第一个模型只考虑周期脉冲释放，得到周期解的一些定性结论，第二个数学模型则是通过密切监测野生蚊子的状态来决定脉冲释放，证明了周期解的存在以及一个局部吸引域。Li [7]为了研究不育蚊子周期性脉冲释放对野生蚊子与不育蚊子相互作用动力学的影响，在建模过程中只考虑那些性活跃的不育蚊子，建立了一种新的时滞模型比较不同策略下模型的动力学行为。Allee效应最早是在20世纪30年代以Allee命名的。Allee [8]通过对金鱼的实验发现金鱼越多，在水槽中的存活率越高。所以种群之间的相互作用也对种群数量也会产生很大影响。Allee效应包括强Allee效应和弱Allee效应两种类型。强Allee效应是指当种群密度低于某一阈值时，种群密度呈现负增长；弱Allee效应意味着稀疏的人口密度会降低人口增长率，但不会出现负增长。对于强Allee效应，Luo [9]研究了在强Allee效应下的野生蚊子与不育蚊子在不同释放速率下的相互作用动力学，得到了不同释放速率下的释放阈值，发现该系统在一定条件下是双稳态的。对于弱Allee效应，Cai [10]研究了弱Allee效应下三种不同释放速率下野生蚊子与不育蚊子的相互作用动力学。在文献[10]的基础上，Wang等在文献[11]中分析了野生蚊子和释放的不育蚊子之间相互作用的常微分方程模型下的分支情况。结果表明，该模型经历了一系列的分岔，显示出单稳态、双稳态或三稳态的动态。Kabor’e等[12]提出了一个通过持续和定期释放感染沃尔巴克氏体的雄性蚊子来控制蚊子种群动态的数学模型，表明结合机械控制和幼虫剂处理，定期释放感染沃尔巴克氏体的雄蚊可以有效减少野生蚊虫数量。Li [13]提出对不育蚊子的三种释放策略，无Allee效应的常数释放模型、考虑Allee效应时不育蚊子与野生蚊子成比例释放模型，考虑Allee效应时具有饱和度比例释放模型。这些模型都假设野生蚊和不育蚊的死亡率与密度无关[14]-[16]。为了在野外蚊虫较少的地区找到更优化、更经济的不育蚊子释放策略[17] [18]，我们在模型中加入了Allee效应。由于即使蚊子很难找到配偶，但它们的繁殖率仍然是正的，所以本文考虑的Allee效应是弱Allee效应[19]。我们的模型与已有文献[20] [21]中关于Allee效应对不育蚊子释放技术(SIT)模型中平衡点存在性与稳定性结论表现出一致性。这种一致性不仅增强了我们模型的可靠性，也为进一步探索Allee效应在SIT应用中的重要性提供了理论支持。

本文研究了Allee效应下不育蚊子的常数释放策略和脉冲释放策略。在第2节中，我们提供了模型的建立过程。第3节和第4节中，我们研究了不育蚊子的常数释放和脉冲释放模型，建立了常数释放模型的平衡点并研究了平衡点的存在性和稳定性。对于脉冲释放模型，讨论了其平凡周期解的存在性和全局性质，并利用比较定理研究了非平凡周期解的潜在的存在性。在第5节中，我们使用数值示例来验证我们的结果。最后讨论了我们的结果。

2. 模型的建立

本节中，我们主要叙述野生蚊子和不育蚊子相互作用脉冲模型的建立过程，Li [13]建立了一个无Allee效应的不育蚊子常数释放到野生蚊群的互动模型如下：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v,\\ \frac{\text{d}{S}_g}{\text{d}t}=B(\cdot )-{\mu }_v{S}_g.\end{array}$ (1)

其中，$S_v$ 表示野生蚊子在$t$ 时刻的数量，$S_g$ 表示不育蚊子在$t$ 时刻的数量，$a$ 代表每只蚊子生产的野生后代数量，假设野生蚊和不育蚊的死亡率都和密度有关，${\mu }_v$ 代表蚊子种群与密度无关的死亡率，${\xi }_v$ 代表蚊子种群与密度相关的死亡率，$B(\cdot )$ 代表不育蚊子的释放率。

考虑到野生当野生蚊子数量较小时，在较小的种群中野生蚊子寻找配偶时可能表现出多样性，Li假

设存在野生蚊群此时存在弱Allee效应。使得导致模型(1)中的交配率变化为$\frac{a{S}_v}{1+S_v}$ ，呈现Holling-II型。

此时这个具有弱Allee效应的野生蚊和不育蚊的互动模型如下：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{1+S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v,\\ \frac{\text{d}{S}_g}{\text{d}t}=B(\cdot )-{\mu }_v{S}_g.\end{array}$ (2)

在模型(1)和(2)中，Li研究了[13]中的三种不同释放策略。首先是无Allee效应的不育蚊的常数释放策略；随后在存在Allee效应的背景下，分析了不育蚊子和野生蚊子的比例释放策略；最后分析了在Allee效应影响下结合饱和的比例释放策略。

随着现代技术的进步，野生蚊子的数量有所减少，但它们的繁殖率仍然是正的。因此我们认为在本研究中，考虑蚊子种群中存在弱Allee效应是必要的。基于模型(1)和(2)的优势，我们提出了新的策略，常数释放不育蚊子和脉冲释放不育蚊子。如果不育蚊子的释放量是恒定常数$B(\cdot )=c$ ，我们有系统如下：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{1+S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v,\\ \frac{\text{d}{S}_g}{\text{d}t}=c-{\mu }_v{S}_g.\end{array}$

如果不育蚊的释放是脉冲释放，则有如下系统：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{1+S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v\begin{array}{ccc}& & t\in \left(iT,\left(i+1\right)T\right]\end{array},\\ \frac{\text{d}{S}_g}{\text{d}t}=-{\mu }_v{S}_g\begin{array}{cc}\begin{array}{ccc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& & \end{array}& \end{array}\begin{array}{ccc}& & t\in \left(iT,\left(i+1\right)T\right],\end{array}\\ S_g\left(n{T}^{+}\right)=S_g\left(nT\right)+m\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& & & \begin{array}{cc}& t=iT.\end{array}\end{array}\end{array}$

在这两个野生蚊和不育蚊的模型中，假设野生蚊子和不育蚊子的死亡率都是一致的[21]。同样的，$S_v$ 表示野生蚊子在$t$ 时刻的数量，$S_g$ 表示不育蚊子在$t$ 时刻的数量，$a$ 代表每只蚊子生产的野生后代数量，假设野生蚊和不育蚊的死亡率都和密度有关，${\mu }_v$ 代表蚊子种群与密度无关的死亡率，${\xi }_v$ 代表蚊子种群与密度相关的死亡率。$m$ 代表不育蚊子在脉冲阶段的释放率。

3. 常数释放模型

在本节中，我们假设不育蚊子的释放量是常数$B(\cdot )=c$ 。因此，野生蚊子和不育蚊子在弱Allee效应下的常数放模型如下：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{1+S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v,\\ \frac{\text{d}{S}_g}{\text{d}t}=c-{\mu }_v{S}_g.\end{array}$ (3)

我们定义了一个新的变量：$S{S}_g=1+S_g$ 。系统(3)被重写为：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{S_v+S{S}_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v,\\ \frac{\text{d}S{S}_g}{\text{d}t}=b-{\mu }_{v}S{S}_g.\end{array}$ (4)

其中$b=c+{\mu }_v$ 。现在我们将研究系统(4)所描述的常数释放模型。

3.1. 系统(4)平衡点的存在性

本节我们研究系统(4)平衡的存在性，通过分析系统的动力学性质可以确定系统能够达到稳态的条件。

定理3.1：

(1) 系统(4)始终具有一个边界平衡$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ ；

(2) 当$b<c_0\doteq \frac{r^2}{4a{\xi }_v}$ ，系统(4)有两个正平衡点$E_c^{-}=\left(S_v^{-},\frac{b}{{\mu }_v}\right)$ ，$E_c^{+}=\left(S_v^{+},\frac{b}{{\mu }_v}\right)$ 和$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ ；

(3) 当$b>c_0\doteq \frac{r^2}{4a{\xi }_v}$ ，系统(4)有一个正平衡点$E_c=\left(\frac{r}{2a{\xi }_v},\frac{b}{{\mu }_v}\right)$ 和$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 。

证明：我们知道系统(3)和(4)是等价的，所以在接下来的证明中可以通过分析系统(4)来分析系统(3)。

首先，系统(4)的边界平衡为：$E^0:=\left(0,\frac{c+{\mu }_v}{{\mu }_v}\right)$ 。系统(4)的正平衡点满足：

$\frac{a{S}_v}{S_v+\frac{c+{\mu }_v}{{\mu }_v}}\left(1-{\xi }_v{S}_v\right)-{\mu }_v=0.$

假设

$\frac{a{S}_v}{S_v+\frac{c+{\mu }_v}{{\mu }_v}}\left(1-{\xi }_v{S}_v\right)-{\mu }_v:=-\frac{{\mu }_{v}Q\left(S_v\right)}{{\mu }_v\left(S_v+1\right)+c}=0,$ (5)

其中

$Q\left(S_v\right):=a{\xi }_v{S}_v{}^2-r{S}_v+b\begin{array}{c},\end{array}r=a-{\mu }_v.$

由于野生蚊子的出生率总是大于自然死亡率，所以我们假设$r>0$ ，则$Q\left(S_v\right)$ 的解为：

$S_v{}^{\pm }=\frac{r\pm \sqrt{r^2-4a{\xi }_{v}b}}{2a{\xi }_v},$ (6)

因此，系统(4)有两个正平衡点：

$E_c{}^{+}=\left(\frac{r+\sqrt{r^2-4a{\xi }_{v}b}}{2a{\xi }_v},\frac{b}{{\mu }_v}\right)\begin{array}{cc},& E_c{}^{-}=\left(\frac{r-\sqrt{r^2-4a{\xi }_{v}b}}{2a{\xi }_v},\frac{b}{{\mu }_v}\right)\end{array},$

取决于$r^2-4a{\xi }_{v}b$ 。当$r^2-4a{\xi }_{v}b=0$ ，我们可以得到释放量阈值$c_0$ 如下：

$c_0=\frac{r^2}{4a{\xi }_v}.$ (7)

因此，当不育蚊的释放量与释放量阈值之间分别为$b=c_0、b<c_0、b>c_0$ 时分别对应存在一个正平衡点、两个正平衡点、不存在正平衡点。

3.2. 系统(4)边界平衡的稳定性

在本节中，我们分别研究了系统(4)边界平衡的局部和全局渐近稳定性。

3.2.1. 边界平衡的局部稳定性

定理3.2 $E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 对于系统(4)总是局部渐近稳定的。

证明：系统(4)的雅可比矩阵如下：

${\Lambda }_0=\left(\begin{array}{cc}a{S}_v\frac{S{S}_g\left(1-{\xi }_v{S}_v\right)}{{\left(S_v+S{S}_g\right)}^2}-\frac{{\xi }_v{S}_v}{S_v+S{S}_g}-{\mu }_v& \frac{a{S}_v^2\left({\xi }_v{S}_v-1\right)}{S_v+S{S}_g}\\ 0& -{\mu }_v\end{array}\right).$ (8)

将$E_c^0$ 代入(8)中可以得到，

${\Lambda }_0=\left(\begin{array}{cc}-{\mu }_v& 0\\ 0& -{\mu }_v\end{array}\right).$ (9)

显然，在边界平衡处$-{\mu }_v<0$ 。所以$E_c^0$ 总是局部渐近稳定的。

3.2.2. 边界平衡的全局渐近稳定性

在这一节我们将会证明平衡点的全局性质，我们引进一个引理。

引理3.2 [22] 考虑以下系统：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=f\left(x\right),\\ \frac{\text{d}y}{\text{d}t}=g\left(x,y\right),\end{array}$ (10)

其中$\left(x,y\right)\in R^n\times R^m$ ，令$\left(x^{\ast },y^{\ast }\right)$ 是一个平衡点。若$x^{\ast }$ 对于系统$\frac{\text{d}x}{\text{d}t}=f\left(x\right)$ 是全局渐近稳定的；$y^{\ast }$ 对于$\frac{\text{d}y}{\text{d}t}=g\left(x^{\ast },y\right)$ 是全局渐近稳定的，则$\left(x^{\ast },y^{\ast }\right)$ 是渐近稳定的。此外，如果(10)的所有轨迹都是正向有界的，

那么$\left(x^{\ast },y^{\ast }\right)$ 是全局渐近稳定的。

定理3.3 当$b>c_0\doteq \frac{r^2}{4a{\xi }_v}$ ，$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 对系统(4)具有全局渐近稳定的。

证明：我们将系统(4)重写为：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=F\left(S_v,S{S}_g\right),\\ \frac{\text{d}S{S}_g}{\text{d}t}=G\left(S{S}_g\right),\end{array}$ (11)

其中$\left(S_v,S{S}_g\right)\in R\times R$ ，$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 是系统(11)的平衡点。

当$b>c_0$ 时，$S{S}_g^{*}=\frac{b}{{\mu }_v}$ 对于系统$\frac{\text{d}S{S}_g}{\text{d}t}=G\left(S{S}_g\right)$ 是全局渐近稳定的；$S_v^{*}$ 对于系统$\frac{\text{d}{S}_v}{\text{d}t}=F\left(S_v,S{S}_g\right)$ 是

全局渐近稳定的。则$\left(S_v^{*},S{S}_g^{*}\right)$ 对于系统(11)是渐近稳定的。这时(11)的所有轨迹都是向前有界的。由引理3.2，我们有$\left(S_v^{*},S{S}_g^{*}\right)$ 对于系统(11)是全局渐近稳定的。

综上所述，当$b>c_0$ 时，$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 是系统(11)的全局渐近稳定的平衡点。

3.3. 系统(4)的正平衡点的稳定性

定理3.4 $b<c_0\doteq \frac{r^2}{4a{\xi }_v}$ ，系统(4)有两个正平衡点$E_c^{-}=\left(S_v^{-},\frac{b}{{\mu }_v}\right)$ ，$E_c^{+}=\left(S_v^{+},\frac{b}{{\mu }_v}\right)$ 和$E_c^0=\left(0,\frac{b}{{\mu }_v}\right)$ 其中$E_c^{-}$ 是一个不稳定的鞍点，$E_c^{+}$ 是一个局部渐近稳定的节点。

证明：我们将雅可比矩阵(8)改写为：

${\Lambda }_0=\left(\begin{array}{cc}{a}_{11}& \frac{a{S}_v^2\left({\xi }_v{S}_v-1\right)}{S_v+S{S}_g}\\ 0& -{\mu }_v\end{array}\right),$

其中，

$\begin{array}{l}{a}_{11}=\frac{\partial \left(S_v\left(\frac{a{S}_v}{S_v+S{S}_g}\left(1-{\xi }_v{S}_v\right)-{\mu }_v\right)\right)}{\partial \left(S_v\right)}\\ \begin{array}{ccc}& =\left(\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{S_v+S{S}_g}-{\mu }_v\right)+\frac{\partial \left(\frac{a{S}_v}{S_v+S{S}_g}\left(1-{\xi }_v{S}_v\right)-{\mu }_v\right)}{\partial \left(S_v\right)}& \end{array}\end{array}$ (12)

由于，在平衡点处，

$S_v+S{S}_g=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v}.$ (13)

将(13)代入(12)中，(12)中的$\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{S_v+S{S}_g}-{\mu }_v$ 可化简为

$\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{S_v+S{S}_g}-{\mu }_v=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v}}-{\mu }_v\text{=}{\mu }_v-{\mu }_v=0.$

所以

$\begin{array}{l}{a}_{11}=0+\frac{\partial \left(\frac{a{S}_v}{S_v+S{S}_g}\left(1-{\xi }_v{S}_v\right)-{\mu }_v\right)}{\partial \left(S_v\right)}\\ \begin{array}{cc}& =\end{array}a{S}_v\left(\frac{S{S}_g\left(1-{\xi }_v{S}_v\right)}{{\left(S_v+S{S}_g\right)}^2}-\frac{{\xi }_v{S}_v}{S_v+S{S}_g}\right)\\ \begin{array}{cc}& =\end{array}\frac{a{S}_v}{{\left(S_v+S{S}_g\right)}^2}\left(S{S}_g\left(1-{\xi }_v{S}_v\right)-{\xi }_v{S}_v\left(S_v+S{S}_g\right)\right).\end{array}$ (14)

显然$\frac{a{S}_v}{{\left(S_v+S{S}_g\right)}^2}>0$ ，现在只需对$\left(S{S}_g\left(1-{\xi }_v{S}_v\right)-{\xi }_v{S}_v\left(S_v+S{S}_g\right)\right)$ 进行化简得到正负。因此，将 $S{S}_g=\frac{b}{{\mu }_v}$ 代入(14)中得到

$\begin{array}{l}{a}_{11}=\frac{a{S}_v}{S_v+S{S}_g}\left(\frac{b}{{\mu }_v}\left(1-{\xi }_v{S}_v\right)-{\xi }_v{S}_v\left(S_v+\frac{b}{{\mu }_v}\right)\right)\\ \begin{array}{cc}& =\end{array}\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(b-a{\xi }_v{S}_v^2\right).\end{array}$ (15)

由(13)，我们可得

$\begin{array}{l}a{\xi }_v{S}_v^2=a{S}_v-{\mu }_v\left(S_v+S{S}_g\right)\\ \begin{array}{ccc}& & =r\end{array}{S}_v-{\mu }_{v}S{S}_g\\ \begin{array}{ccc}& & =\end{array}{S}_v-{\mu }_v\frac{b}{{\mu }_v}\\ \begin{array}{ccc}& & =r\end{array}{S}_v-b,\end{array}$ (16)

将(15)代入(16)中可得

$a_{11}=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(2b-r{S}_v\right).$ (17)

将$S_v^{\pm }$ 代入(17)中可得

$\begin{array}{l}{a}_{11}=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(2b-r{S}_v\right)\\ \begin{array}{cc}& =\end{array}\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(2b-r\cdot \frac{r\pm \sqrt{r^2-4a{\xi }_{v}b}}{2a{\xi }_v}\right)\\ \begin{array}{cc}& =\end{array}\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)}{2a{\xi }_v{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(4a{\xi }_{v}b-r^2\mp r\sqrt{r^2-4a{\xi }_{v}b}\right)\\ \begin{array}{cc}& =\end{array}\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)\sqrt{r^2-4a{\xi }_{v}b}}{2a{\xi }_v{\mu }_v{\left(S_v+S{S}_g\right)}^2}\left(-\sqrt{r^2-4a{\xi }_{v}b}\mp r^2\right).\end{array}$ (18)

则有

$\begin{array}{l}{a}_{11}\left(S_v^{+}\right)=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)\sqrt{r^2-4a{\xi }_{v}b}}{{\mu }_v{\left(S_v+S{S}_g\right)}^{2}2a{\xi }_v}\left(-\sqrt{r^2-4a{\xi }_{v}b}-r^2\right)<0,\\ a_{11}\left(S_v^{-}\right)=\frac{a{S}_v\left(1-{\xi }_v{S}_v\right)\sqrt{r^2-4a{\xi }_{v}b}}{{\mu }_v{\left(S_v+S{S}_g\right)}^{2}2a{\xi }_v}\left(-\sqrt{r^2-4a{\xi }_{v}b}+r^2\right)>0.\end{array}$

所以，我们可得到$E_c^{-}$ 是一个不稳定的鞍点，$E_c^{+}$ 是一个局部渐近稳定的节点。

4. 脉冲释放模型

在实际应用中，脉冲释放模型能更精准的描述不育蚊子和野生蚊子之间的动力学作用。本节研究了弱Allee效应下不育蚊的脉冲释放模型。野生蚊子与不育蚊子相互作用的动力学模型如下：

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=\frac{a{S}_v}{1+S_v+S_g}\left(1-{\xi }_v{S}_v\right)S_v-{\mu }_v{S}_v\begin{array}{ccc}& & t\in \left(iT,\left(i+1\right)T\right]\end{array},\\ \frac{\text{d}{S}_g}{\text{d}t}=-{\mu }_v{S}_g\begin{array}{cc}\begin{array}{ccc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& & \end{array}& \end{array}\begin{array}{ccc}& & t\in \left(iT,\left(i+1\right)T\right],\end{array}\\ S_g\left(n{T}^{+}\right)=S_g\left(nT\right)+m\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& & & \begin{array}{cc}& t=iT.\end{array}\end{array}\end{array}$ (19)

其中，$T$ 是不育蚊子的释放时间，$m$ 为$t$ 时释放的不育蚊子数。现在我们将对系统(19)所描述的脉冲模型进行分析。

4.1. 系统(19)的周期解的存在性

系统(19)的周期解描述了野生蚊群的规模被控制在一定范围之内。

当$S_v\left(t\right)=0$ 时，系统(19)被简化为

$\{\begin{array}{l}\frac{\text{d}{S}_v}{\text{d}t}=0,\\ \frac{\text{d}{S}_g}{\text{d}t}=-{\mu }_v{S}_g,\\ S_g\left(n{T}^{+}\right)=S_g\left(nT\right)+m\end{array}$ (20)

结合(20)最后两个方程，我们可以得到：$S_g\left(\left(n+1\right)T^{+}\right)=S_g\left(n{T}^{+}\right){\text{e}}^{-{\mu }_{v}T}+m$ 。

这是一个离散系统，当$n\to \infty$ 时，此离散系统满足

$S_g^{\ast }=\frac{m}{1-{\text{e}}^{-{\mu }_{v}T}}.$ (21)

所以系统(19)存在一个与野生蚊灭绝对应的周期解

$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)=\left(0,S_g^{\ast }{\text{e}}^{-{\mu }_v(tmodT)}\right).$ (22)

4.2. 系统(19)的平凡周期解的稳定性

在这一部分中，我们讨论了系统(19)的平凡周期解的稳定性。

4.2.1. 平凡周期解的局部渐近稳定性

定理4.1 平凡周期解$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)=\left(0,S_g^{\ast }{\text{e}}^{-{\mu }_v(tmodT)}\right).$对于系统(19)总是局部渐近稳定的。

证明：首先，我们为系统(19)定义一个小扰动振幅：${\tilde{S}}_v=S_v-{\overline{S}}_v\begin{array}{cc},& {\tilde{S}}_g=S_g-{\overline{S}}_g\end{array}.$

我们有${\tilde{S}}_g\left(n{T}^{+}\right)=\left(S_g\left(nT\right)+b\right)-\left({\tilde{S}}_g\left(nT\right)+b\right)=S_g\left(nT\right).$ 此时(19)中的脉冲项消失。

系统(19)可以被重写为：

$\{\begin{array}{l}\frac{\text{d}{\tilde{S}}_v}{\text{d}t}=\frac{a{S}_v}{1+{\tilde{S}}_v+{\overline{S}}_g+{\tilde{S}}_g}\left(1-{\xi }_v{\tilde{S}}_v\right){\tilde{S}}_v-{\mu }_v{\tilde{S}}_v,\\ \frac{\text{d}{\tilde{S}}_g}{\text{d}t}=-{\mu }_v{\tilde{S}}_g,\end{array}$ (23)

将(23)重写为：

$\left(\begin{array}{l}{\tilde{S}}_v\left(t\right)\\ {\tilde{S}}_g\left(t\right)\end{array}\right)=\Phi \left(t\right)\left(\begin{array}{l}{\tilde{S}}_v\left(0\right)\\ {\tilde{S}}_g\left(0\right)\end{array}\right).$ (24)

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

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

其中

${\Lambda }_1\left(t\right)=\left(\begin{array}{cc}\frac{2a{\tilde{S}}_v\left(1-{\xi }_v{\tilde{S}}_v\right)-a{\xi }_v{S}_v}{1+{\tilde{S}}_v+{\overline{S}}_g+{\tilde{S}}_g}-\frac{a{\tilde{S}}_v^2\left(1-{\xi }_v{\tilde{S}}_v\right)}{{\left(1+{\tilde{S}}_v+{\overline{S}}_g+{\tilde{S}}_g\right)}^2}-{\mu }_v& \frac{a{\tilde{S}}_v^2\left(1-{\xi }_v{\tilde{S}}_v\right)}{{\left(1+{\tilde{S}}_v+{\overline{S}}_g+{\tilde{S}}_g\right)}^2}\\ 0& -{\mu }_v\end{array}\right).$

当$S_v\left(t\right)=0$ 时，我们有

${\Lambda }_1\left(t\right)=\left(\begin{array}{cc}-{\mu }_v& 0\\ 0& -{\mu }_v\end{array}\right),$ (26)

且

$\left(\begin{array}{l}{\tilde{S}}_v\left(T^{+}\right)\\ {\tilde{S}}_g\left(T^{+}\right)\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{l}{\tilde{S}}_v\left(nT\right)\\ {\tilde{S}}_g\left(nT\right)\end{array}\right).$ (27)

单值矩阵为

$\Gamma \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\Phi \left(T\right)=\Phi \left(T\right),$ (28)

可知$\Phi \left(T\right)=\Phi \left(0\right)\exp \left({\displaystyle {\int }_0^{\text{T}}tr{\Lambda }_1\left(t\right)\text{d}t}\right)$ ，$\Phi (0)$ 是单位矩阵。假设$\exp \left({\displaystyle {\int }_0^{\text{T}}tr{\Lambda }_1\left(t\right)\text{d}t}\right)$ 的特征值为${\tau }_1,{\tau }_2.$ 则有：

${\tau }_1={\text{e}}^{-{\displaystyle {\int }_0^{\text{T}}{\mu }_v\text{d}t}}<1\begin{array}{cc},& {\tau }_2={\text{e}}^{-{\displaystyle {\int }_0^{\text{T}}{\mu }_v\text{d}t}}<1\end{array}.$

根据Floquet理论[23]可知，系统(19)的平凡周期解总是局部渐近稳定的。

4.2.2. 平凡周期解的全局渐近稳定性

定理4.2 平凡周期解$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)=\left(0,S_g^{\ast }{\text{e}}^{-{\mu }_v(tmodT)}\right).$对于系统(19)是全局渐近稳定的。

证明：我们将系统(23)重写为

$\{\begin{array}{l}\frac{\text{d}{\tilde{S}}_v}{\text{d}t}=F\left({\tilde{S}}_v,{\tilde{S}}_g\right),\\ \frac{\text{d}{\tilde{S}}_g}{\text{d}t}=G\left({\tilde{S}}_g\right).\end{array}$ (29)

事实上，0是系统(29)中的第一个式子的全局渐近稳定的，现我们只需要证明0对于系统(29)的第二项是全局渐近稳定的。

通过计算可知$S_g{}^{\ast }{\text{e}}^{-{\mu }_{v}T}\le S_g{}^{\ast }{\text{e}}^{-{\mu }_v(t\begin{array}{c}\end{array}mood\begin{array}{c}\end{array}T)}\le S_g{}^{\ast }.$

考虑以下两个比较系统

$\frac{\text{d}{\tilde{S}}_{v_{\min }}\left(t\right)}{\text{d}t}=F\left({\tilde{S}}_{v_{\min }}\left(t\right),S_g^{\ast }\right),$ (30)

$\frac{\text{d}{\tilde{S}}_{v_{\max }}\left(t\right)}{\text{d}t}=F\left({\tilde{S}}_{v\max }\left(t\right),S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}\right).$ (31)

其中${\tilde{S}}_v{}_{{}_{\min }}\left(t\right)\doteq \inf {\tilde{S}}_v,{\tilde{S}}_v{}_{{}_{\max }}\left(t\right)\doteq \sup {\tilde{S}}_v$ ，其中${\tilde{S}}_v{}_{{}_{\min }}\left(0\right)={\tilde{S}}_v\left(0\right)={\tilde{S}}_v{}_{{}_{\max }}\left(0\right)={\tilde{S}}_v\left(0\right)$ ，利用比较定理我们有

$0\le {\tilde{S}}_v{}_{{}_{\min }}\left(t\right)\le {\tilde{S}}_v\left(t\right)\le {\tilde{S}}_v{}_{{}_{\max }}\left(t\right).$

我们已知，当$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v}$ ，$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v{\text{e}}^{-{\mu }_{v}t}}$ 时，系统(30) (31)存在全局渐近稳定的平凡周期解。利用比较定理可知，当$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v}$ 时，0对于系统(23)是全局渐近稳定的。由引理3.2可知，当$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v}$ 时，$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)=\left(0,S_g^{\ast }{\text{e}}^{-{\mu }_v(tmodT)}\right)$是对于系统(23)是全局渐近稳定的。

4.3. 非平凡解的存在性

实际上，完全根除野生蚊子是十分困难的，但控制野生蚊子种群的总数[24] [25]也具有一定的意义。

由系统(22)可知${\overline{S}}_{g{}_{\max }}\doteq S_g^{\ast }\begin{array}{cc},& {\overline{S}}_{g{}_{\min }}\doteq \end{array}{S}_g^{\ast }{\text{e}}^{-{\mu }_{v}T}.$

其最大值$S_g^{\ast }$ 和最小值$S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}$ 取决于释放周期$T$ 和释放量$m$ 。将${\overline{S}}_{g_{\min }}=S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}$代入(23)，我们有

$\{\begin{array}{l}\frac{\text{d}{\tilde{S}}_v}{\text{d}t}=\frac{a{S}_v}{1+{\tilde{S}}_v+S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}+{\tilde{S}}_g}\left(1-{\xi }_v{\tilde{S}}_v\right){\tilde{S}}_v-{\mu }_v{\tilde{S}}_v,\\ \frac{\text{d}{\tilde{S}}_g}{\text{d}t}=-{\mu }_v{\tilde{S}}_g,\end{array}$ (32)

由系统(32)可得到

$\{\begin{array}{l}\frac{a{S}_v}{1+{\tilde{S}}_v+S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}+{\tilde{S}}_g}\left(1-{\xi }_v{\tilde{S}}_v\right){\tilde{S}}_v={\mu }_v{\tilde{S}}_v,\\ {\tilde{S}}_g=0.\end{array}$

我们假设$\Omega \left(S_v\right):=\frac{a{\tilde{S}}_v\left(1-{\xi }_v{\tilde{S}}_v\right)}{1+{\tilde{S}}_v+S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}+{\tilde{S}}_g}-{\mu }_v=0,$$\Omega \left(S_v\right)$ 的两根为

$S_v^{\pm }:=\frac{r\pm \sqrt{r^2-4a{\xi }_v{\mu }_v\left(S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}+1\right)}}{2a{\xi }_v{\mu }_v}.$ (33)

从方程(33)中我们可以看到当$r^2-4a{\xi }_v{\mu }_v\left(S_g^{\ast }{\text{e}}^{-{\mu }_{v}T}+1\right)>0,$ 系统(32)有两个正解，此时

$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v{\text{e}}^{-{\mu }_{v}t}}.$

类似的，我们将${\overline{S}}_{g{}_{\max }}=S_g^{\ast }$ 代入(23)可得

$\{\begin{array}{l}\frac{\text{d}{\tilde{S}}_v}{\text{d}t}=\frac{a{S}_v}{1+{\tilde{S}}_v+S_g^{\ast }+{\tilde{S}}_g}\left(1-{\xi }_v{\tilde{S}}_v\right){\tilde{S}}_v-{\mu }_v{\tilde{S}}_v,\\ \frac{\text{d}{\tilde{S}}_g}{\text{d}t}=-{\mu }_v{\tilde{S}}_g,\end{array}$ (34)

可知(34)有两个正平衡点

$S_v^{\pm }:=\frac{r\pm \sqrt{r^2-4a{\xi }_v{\mu }_v\left(S_g^{\ast }+1\right)}}{2a{\xi }_v{\mu }_v},$ (35)

此时$r^2-4a{\xi }_v{\mu }_v\left(S_g^{\ast }+1\right)>0,$ 所以有$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v}.$

我们得出结论，当${\tilde{S}}_g$ 达到其最大值和最小值时，相应的系统(32)和系统(34)都存在两个正解，它们对应于系统(19)的两个非平凡周期解的存在性。因此，根据比较定理$S_g{}^{\ast }{\text{e}}^{-{\mu }_{v}T}\le S_g{}^{\ast }{\text{e}}^{-{\mu }_v(t\begin{array}{c}\end{array}mood\begin{array}{c}\end{array}T)}\le S_g{}^{\ast }.$

当$S_g=S_g^{\ast }{\text{e}}^{-{\mu }_v(tmodT)}$ 时，系统(19)存在两个非平凡解。此时，系统(19)的动力学行为与连续系统相似，表现出双稳状态，其中平凡周期解和其中一个非平凡周期解都是局部稳定的。

数值模拟结果支持了这一结论，如图3(a)中的蓝色曲线所示，它显示出周期性波动的特征，虽然我们不能推导出精确的非平凡周期解，但我们利用这些数值模拟的结果可以验证非平凡周期解的潜在存在性。图中箭头所示的区域强调了周期性行为，反映了系统在指定参数下的复杂的动力学特征。

5. 数值模拟

本文利用二维齐次模型，分析了野生蚊子与不育蚊子相互作用的常数模型和脉冲模型中周期解的存在性与稳定性。通过数值模拟，证明了两种释放策略在特定参数值下模型动力学特征的影响。

参数取值为$a=10,{\xi }_v=0.3,{\mu }_{\text{v}}=0.5,T=14,b=7.$

恒定释放系统中的释放阈值$c_0=7.521$ 。当释放量$b=7<c_0$ 时，系统中不存在正平衡。所有的解都接近于边界平衡$E_0$ 。如图1所示。

当$b>c_0$ 时，系统存在两个正平衡点$E_c^{-}$ 和$E_c^{+}$ 。其中，$E_c^{-}$ ，为鞍点，$E_c^{+}$ 为局部渐近稳定节点，如图2所示。

(a) (b)

Figure 1. When $b=7<c_0$ , there exists no positive equilibrium. In figure (a), all solutions approach boundary equilibrium $E_0$ . In figure (b), we show the attraction region of equilibria. The coordinates of points in (b) correspond to the corresponding initial values, while the shape of points indicates whether the initial values converge to $E_0$ . “+” represents convergence to $E_0$ . We can see that $E_0$ is globally asymptotically stable when $b<c_0$

图1. 当$b=7<c_0$ 时，不存在正平衡。在图(a)中，所有的解都接近边界平衡$E_0$ 。。在图(b)中，我们展示了平衡点的吸引域。(b)中的点的坐标对应于对应的初始值，而点的形状表示初始值是否收敛于$E_0$ 。“+”表示收敛到$E_0$ 。我们可以看到，当$b<c_0$ 时，$E_0$ 是全局渐近稳定的

(a) (b)

Figure 2. When $b>c_0$ , there exits positive equilibria $E_c^{-}$ and $E_c^{+}$ . In figure (a), red dot represent $E_c^0$ , purple dot represent $E_c^{+}$ and green dot represent unstable point $E_c^{-}$ . All solutions approach either $E_c^0$ or $E_c^{+}$ depending on their initial values. In figure (b), we show the attraction region of equilibria. The coordinates points in (b) correspond to the corresponding initial values, while the shape of points indicates whether the initial values approach to $E_c^0$ or $E_c^{+}$ . “+” represents the solution approach to $E_c^0$ , “o” represents the solution approach to $E_c^{+}$ . We can see that $E_c^0$ and $E_c^{+}$ are locally asymptotically stable in our given range

图2. 当$b>c_0$ 时，存在正平衡点$E_c^{-}$ 和$E_c^{+}$ 。图(a)中红点代表$E_c^0$ ，紫色点代表$E_c^{+}$ ，绿点代表不稳定点$E_c^{-}$ 。所有的解决方案都可以根据它们的初始值来接近$E_c^0$ 或$E_c^{+}$ 。图(b)展示了平衡点的吸引域。(b)中的坐标点对应于对应的初始值，而点的形状表示初始值是接近$E_c^0$ 还是$E_c^{+}$ 。“+”表示收敛到$E_c^0$ 的解，“o”表示收敛到$E_c^{+}$ 。我们可以看到$E_c^0$ 和$E_c^{+}$ 在我们给定的范围内是局部渐近稳定的

现对脉冲释放模型进行数值模拟，取参数$a=10,{\xi }_v=0.3,{\mu }_v=0.5,T=14,{\overline{SS}}_g=17$。

当$m<\frac{\left(r-4a{\xi }_v{\mu }_v\right)\left(1-{\text{e}}^{-{\mu }_{v}t}\right)}{4a{\xi }_v{\mu }_v}$ 有两个正的非平凡解。${\tilde{E}}^{+}$ 是局部渐近稳定的，而${\tilde{E}}^{-}$ 是不稳定的。系统中发生脉冲双稳态，如图3所示。

(a) (b)

Figure 3. When $m=14$ , there exits two positive nontrivial periodic solutions ${\tilde{E}}^{+}$ and ${\tilde{E}}^{-}$ . In figure (a), solutions approach either $\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ or ${\tilde{E}}^{+}$ depending on their initial values. In (a), the system exhibits distinct periodic fluctuations. The region indicated by the arrow displays periodic characteristics, hinting at the possible existence of nontrivial periodic solutions. In figure (b), we show the attraction region of equilibria. The coordinates of points in (b) correspond to the corresponding initial values, while the shape of points indicates whether the initial values approach to $\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ or ${\tilde{E}}^{+}$ . “+” represents the solution approach to $\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ , “o” represents the solution approach to ${\tilde{E}}^{+}$ . We can see that $\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ and ${\tilde{E}}^{+}$ are locally asymptotically stable in our given range

图3. 当$m=14$ 时，存在两个正的非平凡周期解${\tilde{E}}^{+}$ 和${\tilde{E}}^{-}$ 。在图(a)中，解可以接近$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ 或${\tilde{E}}^{+}$ ，取决于它们的初始值。在(a)中，系统表现出明显的周期性波动，箭头所示的区域显示了周期特征，暗示了可能存在非平凡周期解。图(b)展示了平衡点的吸引区域，(b)中的点的坐标对应于对应的初始值，而点的形状表示初始值是接近于$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ 还是${\tilde{E}}^{+}$ 。“+”表示趋近于$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ 的解；“o”表示趋近于${\tilde{E}}^{+}$ 的解。我们可以看到$\left({\overline{S}}_v\left(t\right),{\overline{S}}_g\left(t\right)\right)$ 和${\tilde{E}}^{+}$ 在我们给定的范围内是局部渐近稳定的

6. 结论

在这项文章，我们考虑了弱Allee效应可能存在于在野生蚊和不育蚊之间的互动。我们不仅对蚊子的常数释放模型进行了详细的动力学分析还研究了脉冲释放模型中平凡周期解的存在性和稳定性，并利用比较定理探讨了非平凡周期解的潜在存在性。我们的数值模拟为非平凡周期解的存在提供了强有力的证据。虽然解析方法未得到精确的非平凡周期解，但数值结果强烈地表明它们存在于系统中。这为进一步的研究奠定了坚实的基础，强调需要加强数值技术或加入额外的分析方法，以更深入地理解系统的动力学，特别是在存在Allee效应的背景下的动力学。