一类Wolbachia氏菌在蚊群传播的数学模型的动力学研究

Dynamics of a Mathematical Model for the Propagation of Wolbachia Bacteria in Mosquito Populations

摘要: 蚊媒传染病是指由蚊子传播的传染病，主要包括登革热，疟疾等。一种新型的控制和预防蚊媒传染病的方式是利用内共生菌Wolbachia (沃尔巴克氏)来阻断病源体的传播。本文建立了一个新的时滞微分方程模型研究Wolbachia氏菌在蚊群中的传播。首先证明了模型解的非负性和有界性，给出了平衡点的存在条件；在各种不同的参数条件下，得到了平衡点的稳定性态；分析了Hopf分支存在的充分条件；利用中心流形定理和规范型理论给出了确定Hopf分支周期解方向和稳定性的计算公式。最后用数值模拟验证了理论结果。

Abstract: Mosquito-borne diseases are infectious diseases transmitted by mosquitoes, including dengue fever and malaria. A new way to control and prevent mosquito-borne infectious diseases is to use the endosymbiotic bacterium Wolbachia to interrupt the transmission of the pathogen. In this paper, a new delay differential equation model was established to study the transmission of Wolbachia in mosquito population. First of all, the nonnegative and boundedness of solutions of the model are proved, and the existence condition of the equilibrium points is given. The asymptotic stability of the equilibrium is obtained in different situations. Sufficient conditions are analyzed which ensures that Hopf bifurcation occurs. The computational formulas for direction and stability of Hopf bifurcation are given by applying the center manifold theorem and norm form theory. Finally, numerical simulations are provided to demonstrate the theoretical results.

1. 引言

蚊媒传染病指的是由蚊子来传播的传染病，也叫虫媒传染病。蚊媒传染病主要包括登革热、寨卡、疟疾等。近年来，蚊媒传染病在全球时有爆发，其病例数也呈现增长趋势。其中，登革热是发展速度最快的蚊媒传染病，在过去的50年间，其发病的数量迅速增长，目前全世界有将近40%的人面临着被感染的风险。登革热是登革病毒经伊蚊传播引起的一种急性传染病，主要流行于热带和亚热带国家和地区。近年来，我国曾多次出现过登革热的局部暴发，例如在2014年，我国广东省就爆发了很严重的登革热疫情。目前，能够有效减少登革热传播的唯一途径是防蚊灭蚊，比如通过大量喷洒杀虫剂减少蚊子数量。但是这种方法不仅会带来生态问题，还可能导致蚊子对杀虫剂逐渐产生抗药性。考虑到这些情况，防蚊灭蚊的措施亟待改进，应该用更长远、更环保的方法来防止蚊媒传染病的传播。

Wolbachia (沃尔巴克氏)是一类广泛分布于节肢动物和无脊椎动物生殖组织内的内共生细菌，它可以引发宿主产生多种生殖异常反应，如最常见的细胞质不兼容(CI)，等等。因此，有一类新兴的控制蚊媒传染病的方法，就是使自然界中的伊蚊种群感染内共生菌Wolbachia (沃尔巴克氏)。一系列的实验证明这样可以阻断蚊子传播登革热等疾病 [1] [2] [3]。Wolbachia可以通过CI机制成功入侵自然界中的伊蚊种群，使得在伊蚊种群中，由于CI机制的作用，野外雌蚊与感染了Wolbachia的雄蚊交配所产的卵无法正常发育，但对感染Wolbachia的雌蚊无影响，并且其后代大多会感染Wolbachia，因此Wolbachia在入侵自然界的蚊群时拥有传播的优势。而且，基于Wolbachia的蚊媒控制方法是安全的。因此，利用Wolbachia阻断登革热的传播在应用上具有重要价值。早在2005年，生物学家奚志勇就在实验室里成功地使Wolbachia在埃及伊蚊中稳定传播 [3]。随着Wolbachia蚊媒控制技术的发展，建立数学模型研究Wolbachia传播动力学具有了重要意义。通过数学模型分析研究Wolbachia传播动力学，可以为感染Wolbachia的蚊子释放提供很好的策略和方案。

1990年以来，Hoffmann和Turelli建立了世代不重叠的离散模型并做了解释与分析 [4] [5]。2014年，Zheng等建立了如下时滞微分方程来研究Wolbachia传播动力学 [6]：

$\{\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=b_1\left[x\left(t-\tau \right)+2R_F\left(t-\tau \right)\right]-{\delta }_{1}x\left(t\right)\left[x\left(t\right)+y\left(t\right)+R_F\left(t\right)+R_M\left(t\right)\right],\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\frac{b_2{y}^2\left(t-\tau \right)}{x\left(t-\tau \right)+y\left(t-\tau \right)+2R_M\left(t-\tau \right)}-{\delta }_{2}y\left(t\right)\left[x\left(t\right)+y\left(t\right)+R_F\left(t\right)+R_M\left(t\right)\right],\\ \frac{\text{d}{R}_F\left(t\right)}{\text{d}t}=-{\delta }_1{R}_F\left(t\right)\left[x\left(t\right)+y\left(t\right)+R_F\left(t\right)+R_M\left(t\right)\right],\\ \frac{\text{d}{R}_M\left(t\right)}{\text{d}t}=-{\delta }_1{R}_M\left(t\right)\left[x\left(t\right)+y\left(t\right)+R_F\left(t\right)+R_M\left(t\right)\right].\end{array}$ (1)

注意到，模型中出生函数为线性函数。一般说来，随着种群规模的增加，其成年个体的出生率和未成年个体的幸存率将下降，这是由于种群内部的竞争引起的，而蚊群的内部竞争主要发生在未成年阶段，所以单位时间内成年个体不再是线性增长。另外，基于成年个体对资源的竞争假设，模型的死亡函数为二次函数。但是大量的事实表明，成蚊之间的竞争几乎可以忽略。基于以上考虑，本文建立一个新的时滞微分方程模型，出生和幸存函数为Ricker型函数，而死亡函数用线性函数。

具体来说，我们建立以下的时滞微分方程模型

$\{\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=p_1{\text{e}}^{-{\delta }_1\tau }x\left(t-\tau \right){\text{e}}^{-q_1\left[x(t-\tau )+y(t-\tau )\right]}-{\delta }_{1}x\left(t\right)\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=p_2{\text{e}}^{-{\delta }_2\tau }\frac{y^{\text{2}}\left(t-\tau \right)}{x\left(t-\tau \right)+y\left(t-\tau \right)}\cdot {\text{e}}^{-q_2\left[x(t-\tau )+y(t-\tau )\right]}-{\delta }_{2}y\left(t\right)\end{array}$ (2)

初始条件

$x\left(\theta \right)={\phi }_1\left(\theta \right)\ge 0,y\left(\theta \right)={\phi }_2\left(\theta \right)\ge 0,\theta \in \left[-\tau ,0\right],{\phi }_1\left(0\right)>0,{\phi }_2\left(0\right)>0.$ (3)

其中 $x\left(t\right)$ 表示感染Wolbachia的蚊群规模， $y\left(t\right)$ 表示野外蚊群规模， $p_1$， $p_{\text{2}}$ 分别为感染蚊群和野外蚊群

的出生率， ${\delta }_{\text{1}}$， ${\delta }_{\text{2}}$ 分别为感染蚊群和野外蚊群的死亡率， $\frac{\text{1}}{q_1}$， $\frac{\text{1}}{q_{\text{2}}}$ 分别为感染蚊群和野外蚊群以最大繁

殖率繁殖的规模， $\tau$ 是蚊子从交配到后代出现的平均等待时间。

关于Wolbachia入侵动力学及其它传染病模型可参见 [7] - [13] 以及其中的参考文献。例如Keeling等 [8] 建立的常微分方程模型，Huang等 [9] [10] 建立了反应扩散方程及其衍生的常微分方程模型，Hu等 [11] 建立的随机微分方程模型，Li等 [12] 建立的离散竞争模型，Huang等 [13] 建立的时滞微分方程。

本文第二节主要分析了模型(2)的动力学性质，包括解的正性、有界性，平衡点的存在性和稳定性和Hopf分支的存在性。第三节确定了Hopf分支方向和周期解的稳定性计算公式。第四节，通过数值模拟验证了理论的结论。

2. 模型(2)的平衡点的存在性、稳定性与Hopf分支出现的条件

定理1：模型(2)的满足初始条件(3)的解是非负的，并且最终有界。

证明：解的非负性显然，下面证明其有界性。令 $k\left(x\right)=x{\text{e}}^{-q_{1}x}$，则 $\max \left\{f\left(x\right):x\ge 0\right\}=\frac{1}{q_{1}e}$。因此有

$\frac{\text{d}x}{\text{d}t}\le p_{1}x\left(t-\tau \right){\text{e}}^{-q_{1}x(t-\tau )}-{\delta }_{1}x\left(t\right)\le {\delta }_{1}x\left(t\right)+\frac{p_1}{q_{1}e}$

那么

$x\left(t\right)\le \left(x\left(0\right)-\frac{p_1}{{\delta }_1{q}_{1}e}\right){\text{e}}^{-{\delta }_{1}t}+\frac{p_1}{q_1{\delta }_{1}e}$

即有

$\underset{t\to \infty }{\lim }\sup x\left(t\right)\le \frac{p_1}{q_1{\delta }_{1}e}$

同理可以得到

$\frac{\text{d}y}{\text{d}t}\le p_{2}y\left(t-\tau \right){\text{e}}^{-q_{2}y(t-\tau )}-{\delta }_{2}y(\; t\; )$

即有

$\underset{t\to \infty }{\lim }\sup y\left(t\right)\le \frac{p_2}{q_2{\delta }_{\text{2}}e}$。

下面求模型的平衡点。由

$\{\begin{array}{l}x\left(p_1{\text{e}}^{-{\delta }_1\tau }{\text{e}}^{-q_1(x+y)}-{\delta }_1\right)=0\\ y\left(p_2{\text{e}}^{-{\delta }_2\tau }\frac{y}{x+y}\cdot {\text{e}}^{-q_2(x+y)}-{\delta }_2\right)=0\end{array}$ (4)

可知系统有四个平衡点 ${\tilde{E}}_0\left(0,0\right)$ 、 ${\tilde{E}}_1\left(\tilde{x},0\right)$ 、 ${\tilde{E}}_2\left(0,\tilde{y}\right)$ 以及 ${\tilde{E}}^{*}\left({\tilde{x}}^{*},{\tilde{y}}^{*}\right)$。易知

$\tilde{x}=-\frac{{\delta }_1\tau }{q_1}+\frac{1}{q_1}\ln \frac{p_1}{{\delta }_1}$， $\tilde{y}=-\frac{{\delta }_{\text{2}}\tau }{q_{\text{2}}}+\frac{1}{q_{\text{2}}}\ln \frac{p_{\text{2}}}{{\delta }_{\text{2}}}$

${\tilde{x}}^{*}=\left(-\frac{{\delta }_1\tau }{q_1}+\frac{1}{q_1}\ln \frac{p_1}{{\delta }_1}\right)\left[1-\frac{{\delta }_2}{p_2}{\left(\frac{p_1}{{\delta }_1}\right)}^{\frac{q_2}{q_1}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_1}}\right]$

${\tilde{y}}^{*}=\frac{{\delta }_2}{p_2}{\left(\frac{p_1}{{\delta }_1}\right)}^{\frac{q_2}{q_1}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_1}}\left(-\frac{{\delta }_1\tau }{q_1}+\frac{1}{q_1}\ln \frac{p_1}{{\delta }_1}\right)$

显然， ${\tilde{x}}^{*}+{\tilde{y}}^{*}=-\frac{{\delta }_1\tau }{q_1}+\frac{1}{q_1}\ln \frac{p_1}{{\delta }_1}$。当 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau }$ 时 ${\tilde{E}}_1$ 存在，当 $\frac{p_2}{{\delta }_2}>{\text{e}}^{{\delta }_2\tau }$ 时 ${\tilde{E}}_{\text{2}}$ 存在。关于正平衡点存在的条件，注意到必有 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau }$，且 $\left(-\frac{{\delta }_1\tau }{q_1}+\frac{1}{q_1}\ln \frac{p_1}{{\delta }_1}\right)\left[1-\frac{{\delta }_2}{p_2}{\left(\frac{p_1}{{\delta }_1}\right)}^{\frac{q_2}{q_1}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_1}}\right]>0$，即 $\frac{p_2}{{\delta }_2}>{\left(\frac{p_1}{{\delta }_1}\right)}^{\frac{q_2}{q_1}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_1}}$。等价地有 $\tilde{y}>\tilde{x}$。因此当 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau }$， $\tilde{y}>\tilde{x}$ 时，系统存在正平衡点。

下面分析各平衡点的稳定性。

定理2：对于模型(2)，我们有以下结论：

(1) 当 $\frac{p_1}{{\delta }_1}\le {\text{e}}^{{\delta }_1\tau }$， $\frac{p_{\text{2}}}{{\delta }_{\text{2}}}\le {\text{e}}^{{\delta }_{\text{2}}\tau }$ 时， ${\tilde{E}}_0$ 是全局渐近稳定的。

(2) 当 ${\text{e}}^{{\delta }_{\text{2}}\tau }<\frac{p_{\text{2}}}{{\delta }_{\text{2}}}\le {\text{e}}^{{\delta }_{\text{2}}\tau +2}$， $\tilde{y}>\tilde{x}$ 时， ${\tilde{E}}_{\text{2}}$ 是局部渐近稳定的。

(3) 当 $e^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\mathrm{e}}^{{\delta }_{\text{1}}\tau +2}$ 时， ${\tilde{E}}_{\text{1}}$ 是局部渐近稳定的。

证明：(1) 当 $\frac{p_1}{{\delta }_1}\le {\text{e}}^{{\delta }_1\tau }$， $\frac{p_{\text{2}}}{{\delta }_{\text{2}}}\le {\text{e}}^{{\delta }_{\text{2}}\tau }$ 时，只存在零平衡点 ${\tilde{E}}_0$。

$\frac{\text{d}}{\text{d}t}\left[x\left(t\right)+p_1{\text{e}}^{-{\delta }_1\tau }{\displaystyle {\int }_{t-\tau }^{t}x\left(s\right)}{\text{e}}^{-q_1(x(s)+y(s))}\text{d}s\right]=p_1{\text{e}}^{-{\delta }_1\tau }x\left(t\right){\text{e}}^{-q_1[x(t)+y(t)]}-{\delta }_{1}x(\; t\; )$

$\frac{\text{d}}{\text{d}t}\left[y\left(t\right)+p_2{\text{e}}^{-{\delta }_2\tau }{\displaystyle {\int }_{t-\tau }^t\frac{y^2\left(s\right)}{x\left(s\right)+y\left(s\right)}}{\text{e}}^{-q_2(x(s)+y(s))}\text{d}s\right]=p_2{\text{e}}^{-{\delta }_2\tau }\frac{y^2\left(t\right)}{x\left(t\right)+y\left(t\right)}{\text{e}}^{-q_2[x(t)+y(t)]}-{\delta }_{2}y(\; t\; )$

取 $V=x\left(t\right)+p_1{\text{e}}^{-{\delta }_1\tau }{\displaystyle {\int }_{t-\tau }^{t}x\left(s\right)}{\text{e}}^{-q_1(x(s)+y(s))}\text{d}s+y\left(t\right)+p_2{\text{e}}^{-{\delta }_{\text{2}}\tau }{\displaystyle {\int }_{t-\tau }^t\frac{y^2\left(s\right)}{x\left(s\right)+y\left(s\right)}}{\text{e}}^{-q_2(x(s)+y(s))}\text{d}s$

那么对 $t\ge 0$， $V>0$，有：

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\left[p_1{\text{e}}^{-{\delta }_1\tau }{\text{e}}^{-q_1[x(t)+y(t)]}-{\delta }_1\right]x\left(t\right)+\left[p_2{\text{e}}^{-{\delta }_{\text{2}}\tau }\frac{y\left(t\right)}{x\left(t\right)+y\left(t\right)}{\text{e}}^{-q_2[x(t)+y(t)]}-{\delta }_2\right]y\left(t\right)\\ \le \left(p_1{\text{e}}^{-{\delta }_1\tau }-{\delta }_1\right)x\left(t\right)+\left(p_2{\text{e}}^{-{\delta }_{\text{2}}\tau }-{\delta }_2\right)y\left(t\right)\le 0\end{array}$

且 $\frac{\text{d}V}{\text{d}t}=0$ 当且仅当 $x=y=0$。因此，零解 ${\tilde{E}}_0$ 全局渐近稳定。

(2) 将系统(2)在平衡点 $\tilde{E}=\left(x\left(t\right),y\left(t\right)\right)$ 进行线性化，得到如下形式：

$T^{\prime }\left(t\right)=A_{1}T\left(t\right)+BT\left(t-\tau \right)$ (5)

其中 $T\left(t\right)={\left(x\left(t\right),y\left(t\right)\right)}^{{\rm T}}$

$A_{\text{1}}=\left(\begin{array}{cc}-{\delta }_1& 0\\ 0& -{\delta }_2\end{array}\right)$ (6)

$B=\left(\begin{array}{cc}{p}_1{\text{e}}^{-{\delta }_1\tau }{\text{e}}^{-q_1(x+y)}-p_1{q}_1{\text{e}}^{-{\delta }_1\tau }x{\text{e}}^{-q_1(x+y)}& -p_1{q}_1{\text{e}}^{-{\delta }_1\tau }x{\text{e}}^{-q_1(x+y)}\\ -p_2{\text{e}}^{-{\delta }_2\tau }\frac{y^2}{{\left(x+y\right)}^2}{\text{e}}^{-q_2(x+y)}-p_2{q}_2{\text{e}}^{-{\delta }_2\tau }\frac{y^2}{x+y}{\text{e}}^{-q_2(x+y)}& p_2{\text{e}}^{-{\delta }_2\tau }\frac{2xy+y^2}{{\left(x+y\right)}^2}{\text{e}}^{-q_2(x+y)}-p_2{q}_2{\text{e}}^{-{\delta }_2\tau }\frac{y^2}{x+y}{\text{e}}^{-q_2(x+y)}\end{array}\right)$ (7)

我们可以通过如下公式来求解特征方程:

$\det \left(\lambda I-A-B{\text{e}}^{-\lambda \tau }\right)=0$ (8)

在 ${\tilde{E}}_{\text{2}}$ 处，其中

$A_{\text{1}}=\left(\begin{array}{cc}-{\delta }_1& 0\\ 0& -{\delta }_2\end{array}\right)B=\left(\begin{array}{cc}{p}_1{\left(\frac{{\delta }_2}{p_2}\right)}^{\frac{q_1}{q_2}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_2}}& 0\\ -{\delta }_2+{\delta }_2\left({\delta }_2\tau -\ln \frac{p_2}{{\delta }_2}\right)& {\delta }_2+{\delta }_2\left({\delta }_2\tau -\ln \frac{p_2}{{\delta }_2}\right)\end{array}\right)$

特征方程为

$\left[\lambda +{\delta }_1-p_1{\left(\frac{{\delta }_2}{p_2}\right)}^{\frac{q_1}{q_2}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_2}}{\text{e}}^{-\lambda \tau }\right]\left[\lambda +{\delta }_2-{\delta }_2\left(1+{\delta }_2\tau -\ln \frac{p_2}{{\delta }_2}\right){\text{e}}^{-\lambda \tau }\right]=0$ (9)

式子的前半部分， $a=-{\delta }_1$， $b=p_1{\left(\frac{{\delta }_2}{p_2}\right)}^{\frac{q_1}{q_2}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_2}}$，根据定理，当 $a+b<0$ 时， $\frac{p_2}{{\delta }_2}>{\left(\frac{p_1}{{\delta }_1}\right)}^{\frac{q_2}{q_1}}{\text{e}}^{\frac{(q_1{\delta }_2-q_2{\delta }_1)\tau }{q_1}}$，即 $\tilde{y}>\tilde{x}$ 时，满足前半部分为负实部的根。式子的后半部分， $a=-{\delta }_{\text{2}}$， $b={\delta }_2\left(1+{\delta }_2\tau -\ln \frac{p_2}{{\delta }_2}\right)$，根据文献 [14] 中的定理，当 $b\ge a$ 时，即 $e^{{\delta }_{\text{2}}\tau }<\frac{p_{\text{2}}}{{\delta }_{\text{2}}}\le {\text{e}}^{{\delta }_{\text{2}}\tau +2}$ 时，满足后半部分为负实部的根。

因此， $e^{{\delta }_{\text{2}}\tau }<\frac{p_{\text{2}}}{{\delta }_{\text{2}}}\le {\text{e}}^{{\delta }_{\text{2}}\tau +2}$， $\tilde{y}>\tilde{x}$ 时， ${\tilde{E}}_2$ 是局部渐近稳定的。

(3) 在 ${\tilde{E}}_1$ 处，其中

$A_1=\left(\begin{array}{cc}-{\delta }_1& 0\\ 0& -{\delta }_2\end{array}\right)B=\left(\begin{array}{cc}{\delta }_1-{\delta }_1\left({\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)& {\delta }_1\left({\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)\\ 0& 0\end{array}\right)$

特征方程为

$\left[\lambda +{\delta }_1-{\delta }_1\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right){\text{e}}^{-\lambda \tau }\right]\left(\lambda +{\delta }_2\right)=0$ (10)

由于式子的后半部分必为负实部的根，因此只需要看式子的前半部分即可。其中 $a=-{\delta }_1$，

$b={\delta }_1\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)$。同样根据定理，当 $b\ge a$ 时，即 ${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$ 时， ${\tilde{E}}_{\text{1}}$ 是局部渐近稳定的。

下面，我们将确定平衡点 ${\tilde{E}}_1$ 何时变得不稳定并发生Hopf分支。特征方程(10)可以写成

${\lambda }^{\text{2}}+A\lambda +B+\left(A_1\lambda +B_1\right){\text{e}}^{-\lambda \tau }=0$ (11)

其中， $A={\delta }_1+{\delta }_2$， $B={\delta }_1{\delta }_2$， $A_1=-{\delta }_1\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)$， $B_1=-{\delta }_1{\delta }_2\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)$

假设存在 $\tau ={\tau }^{*}$，使得(11)有一对纯虚根，用 $\lambda =\pm i\omega$， $\left(\omega >0\right)$ 表示。将 $\lambda =i\omega$ 带入，得到

$-{\omega }^2+iA\omega +B-\left(i{A}_1\omega +B_1\right)\left(\cos \omega \tau -i\sin \omega \tau \right)=0$

分离实部与虚部，有

$\{\begin{array}{l}{B}_1\cos \omega \tau +A_1\omega \sin \omega \tau ={\omega }^2-B\\ A_1\omega \cos \omega \tau -B_1\sin \omega \tau =-A\omega \end{array}$

两式平方相加得到

${\omega }^4+\left(2A^2-2B-A_1^2\right){\omega }^2+B^2-B_1^2=0$ (12)

令 $z={\omega }^2$，有

$z^2+\left(2A^2-2B-A_1^2\right)z+B^2-B_1^2=0$ (13)

由于(10)至少有一个负实部的特征根，即最多有一个正实部的特征根。因此(13)最多有一个正根。当

$B^2-B_1^2<0$，即 ${\delta }_1^2{\delta }_2^2\left[-2\left({\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)-{\left({\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^2\right]<0$ 时出现纯虚根的情况，即 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$ 时。

对于 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$ 时，有

${\tau }_k=\frac{\text{1}}{\omega }\left(\arccos \frac{B_1{\omega }^2-B{B}_1-A{A}_1{\omega }^2}{A_1^2{\omega }^2+B_1^2}+2k\pi \right)$

其中， $k=\text{0,1,2,}\cdots$。

令 ${\lambda }_k\left(\tau \right)={\alpha }_k\left(\tau \right)+i{\omega }_k\left(\tau \right)$ 是方程(11) $\tau ={\tau }_{\text{k}}$ 时附近的特征值，那么有 $\alpha \left({\tau }_k\right)=0$ 和 $\omega \left({\tau }_k\right)={\omega }_0$。因此我们有下面的结论。

定理3：如果 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$， ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$，存在 ${\tau }_k>0$，当 $\tau ={\tau }_k$ 时，模型(2)在 ${\tilde{E}}_1$ 处发生Hopf分支。

证明： ${\lambda }_k\left(\tau \right)={\alpha }_k\left(\tau \right)+i{\omega }_k\left(\tau \right)$ 是方程(11)在 $\tau ={\tau }_k$ 附近的特征值，将方程(11)对 $\tau$ 求导，我们得到

$\left(2\lambda +A\right)\frac{\text{d}\lambda }{\text{d}\tau }+A_1{\text{e}}^{-\lambda \tau }\frac{\text{d}\lambda }{\text{d}\tau }-\left(A_1\lambda +B_1\right){\text{e}}^{-\lambda \tau }\left(\lambda +\tau \frac{\text{d}\lambda }{\text{d}\tau }\right)=0$

即有

$\frac{\text{d}\lambda }{\text{d}\tau }=\frac{\lambda \left(A_1\lambda +B_1\right){\text{e}}^{-\lambda \tau }}{2\lambda +A+\left[A_1-\tau \left(A_1\lambda +B_1\right)\right]{\text{e}}^{-\lambda \tau }}$

由于 $\alpha \left({\tau }_k\right)=0$ 和 $\omega \left({\tau }_k\right)={\omega }_0$，有

$\begin{array}{c}{\left[\frac{\text{d}\mathrm{Re}\lambda \left({\tau }_k\right)}{\text{d}\tau }\right]}^{-1}=\mathrm{Re}{\left[\frac{\left(2\lambda +A\right){\text{e}}^{\lambda \tau }+A_1}{\lambda \left(A_1\lambda +B_1\right)}\right]}_{\tau ={\tau }_k}\\ =\frac{1}{Z}\left[-2\omega \left(A_1\omega \cos \omega \tau -B_1\sin \omega \tau \right)-A\left(B_1\cos \omega \tau +A_1\omega \sin \omega \tau \right)-A_1{B}_1\right]\\ =\frac{1}{Z}\left(A{\omega }^2+AB-A_1{B}_1\right)\end{array}$

其中， $Z={\omega }^2\left(A_1^2{\omega }^2+B_1^2\right)$。

当 $AB-A_1{B}_1>0$，即 ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$ 时，有 $\frac{\text{d}\mathrm{Re}\lambda \left({\tau }_k\right)}{\text{d}\tau }>0$。

3. Hopf分支的方向与稳定性

我们已经证明了Hopf分支的存在。接下来，我们由文献 [15] 和 [16] 中的Hopf分支定理，利用中心流形定理和规范型理论给出确定Hopf分支的方向和周期解的稳定性的计算公式。

令 $m\left(t\right)=x\left(\tau t-\tilde{x}\right)$， $n\left(t\right)=y\left(\tau t\right)$，则有：

$\begin{array}{c}{m}^{\prime }\left(t\right)=\tau [-{\delta }_{1}m\left(t\right)+\left({\delta }_1-q_1{\delta }_1\tilde{x}\right)m\left(t-1\right)-q_1{\delta }_1\tilde{x}n\left(t-1\right)+\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)m^2\left(t-1\right)\\ +\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}{n}^2\left(t-1\right)+\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)m\left(t-1\right)n\left(t-1\right)+o\left(\rho \right)]\end{array}$

$n^{\prime }\left(t\right)=\tau \left[-{\delta }_{2}n\left(t\right)+p_2{\text{e}}^{-{\delta }_2\tau }\frac{n^2\left(t-1\right)}{m\left(t-1\right)+n\left(t-1\right)+\tilde{x}}-\frac{p_2{e}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}{n}^2\left(t-1\right)+o\left(\rho \right)\right]$

其中， $o\left(\rho \right)=\sqrt{m^2\left(t-1\right)+n^2\left(t-1\right)}$。

系统(2)可以写成：

$\left(\begin{array}{l}{m}^{\prime }\left(t\right)\\ n^{\prime }\left(t\right)\end{array}\right)=\tau A_1\left(\begin{array}{l}m\left(t\right)\\ n\left(t\right)\end{array}\right)+\tau B\left(\begin{array}{l}m\left(t-1\right)\\ n\left(t-1\right)\end{array}\right)+F\left(m_t,n_t,\tau \right)$ (14)

其中，

$A_{\text{1}}=\left(\begin{array}{cc}-{\delta }_1& 0\\ 0& -{\delta }_2\end{array}\right)$， $B=\left(\begin{array}{cc}{\delta }_1-q_1{\delta }_1\tilde{x}& -q_1{\delta }_1\tilde{x}\\ 0& 0\end{array}\right)$， $F=\left(\begin{array}{c}\tau K_1\\ \tau K_2\end{array}\right)$，

$K_1=\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)m^2\left(t-1\right)+\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}{n}^2\left(t-1\right)+\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)m\left(t-1\right)n\left(t-1\right)+o(\; \rho \; )$

$K_2=p_2{\text{e}}^{-{\delta }_2\tau }\frac{n^2\left(t-1\right)}{m\left(t-1\right)+n\left(t-1\right)+\tilde{x}}-\frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}{n}^2\left(t-1\right)+o(\; \rho \; )$

为了简便，设 $\tau =\hat{\tau }$ 时系统在 ${\tilde{E}}_1$ 处产生Hopf分支， $\pm i\hat{\omega }$ 是其对应的纯虚根。令 $\tau =\hat{\tau }+\mu$，则 $\mu =0$ 是系统的Hopf分支值。

记 $C^k\left[-1,0\right]=\left\{\varphi |\varphi =\left({\varphi }_1,{\varphi }_2\right):\left[-1,0\right]\to \text{}{R}^2\right\}$，其中 ${\varphi }_{\text{1}}$， ${\varphi }_{\text{2}}$ 在 $\left[-1,0\right]$ 上 $k$ 次连续可微。

记 $C\left[-1,0\right]：=C^{\text{0}}\left[-1,0\right]$，对 $\varphi ={\left({\varphi }_1,{\varphi }_2\right)}^{{\rm T}}\in C\left[-1,0\right]$。定义 $L_{\mu }:C\left[-1,0\right]\to \text{}{R}^2$ 为：

$L_{\mu }\varphi =\left(\hat{\tau }+\mu \right)A_1\varphi \left(0\right)+\left(\hat{\tau }+\mu \right)B\varphi (\; -\; 1\; )$

由Riesz表示定理，存在有界变差二阶函数矩阵 $\eta \left(\theta ,\mu \right):\left[-1,0\right]\to \text{}{R}^{2\times 2}$， $\theta \in \left[-1,0\right]$。使得：

$L_{\mu }\varphi ={\displaystyle {\int }_{-1}^0\text{d}}\eta \left(\theta ,\mu \right)\varphi \left(\theta \right)\text{,}\varphi \in C$

选取 $\eta \left(\theta ,\mu \right)=\left(\hat{\tau }+\mu \right)A_1\delta \left(\theta \right)-\left(\hat{\tau }+\mu \right)B\delta \left(\theta +1\right)$，其中

$\delta \left(\theta \right)=\{\begin{array}{l}1,\theta =0\\ 0,\theta \ne 0\end{array}$

对于 $\varphi \in C^1\left(\left[-1,0\right],R^2\right)$，定义

$A\left(\mu \right)\varphi =\{\begin{array}{l}\frac{\text{d}\varphi \left(\theta \right)}{\text{d}\theta },\theta \in \left[-1,0\right)\\ {\displaystyle {\int }_{-1}^0\text{d}}\eta \left(t,\mu \right)\varphi \left(t\right),\theta =0\end{array}$

$R\left(\mu \right)\varphi =\{\begin{array}{l}0,\theta \in \left[-1,0\right)\\ F\left(\varphi ,\hat{\tau }+\mu \right),\theta =0\end{array}$

则(14)可以等价的表示成：

${v^{\prime }}_t=A\left(\mu \right)v_t+R\left(\mu \right)v_t$， (15)

其中 $v_t={\left(x_t,y_t\right)}^{{\rm T}}$， $v_t\left(\theta \right)=v\left(t+\theta \right)$， $\theta \in \left[-1,0\right)$。

对于 $\phi \in C^1\left(\left[0,1\right],{\left(C^2\right)}^{*}\right)$，其中 ${\left(C^2\right)}^{*}$ 是2维复行向量空间，定义

$A^{*}\phi \left(s\right)=\{\begin{array}{l}-\text{d}\phi \left(s\right)\text{d}s,s\in \left(0,1\right]\\ {\displaystyle {\int }_{-1}^0\text{d}\eta \left(t,0\right)\phi \left(-t\right)},s=0\end{array}$

进一步，对于 $\varphi \in C^1\left(\left[-1,0\right],R^2\right)$， $\phi \in C^1\left(\left[0,1\right],{\left(C^2\right)}^{*}\right)$，定义双线性形式如下：

$\langle \phi ,\varphi \rangle =\bar{\phi }\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-1}^0{\displaystyle {\int }_{\xi =0}^{\theta }\bar{\phi }\left(\xi -\theta \right)\text{d}\eta }}\left(\theta ,0\right)\varphi \left(\xi \right)\text{d}\xi$

那么 $A^{*}$ 与 $A\left(0\right)$ 是伴随算子，使得 $\langle \phi ,A\varphi \rangle =\langle A^{*}\phi ,\varphi \rangle$。令 $q\left(\theta \right)$ 和 $q^{*}\left(s\right)$ 是 $A$ 和 $A^{*}$ 对应特征值 $i\hat{\omega }\hat{\tau }$ 和 $-i\hat{\omega }\hat{\tau }$ 的特征向量，经过计算，我们得到下面的结果。

引理1：向量 $q\left(\theta \right)={\left(1,{\alpha }_1\right)}^{{\rm T}}{\text{e}}^{i\hat{\omega }\hat{\tau }\theta }$ 为算子 $A$ 关于特征值 $i\hat{\omega }\hat{\tau }$ 的特征向量， $q^{*}\left(s\right)=\bar{D}{\left(1,{\alpha }_2\right)}^{{\rm T}}{\text{e}}^{i\hat{\omega }\hat{\tau }s}$ 为算子

$A^{*}$ 关于特征值 $-i\hat{\omega }\hat{\tau }$ 的特征向量，并且 $\langle q^{*},q\rangle =1$， $\langle q^{*},\bar{q}\rangle =0$，其中， ${\alpha }_1=\frac{-{\delta }_1-i\hat{\omega }+\left({\delta }_1-q_1{\delta }_1\tilde{x}\right){\text{e}}^{-i\hat{\omega }\hat{\tau }}}{q_1{\delta }_1\tilde{x}{\text{e}}^{-i\hat{\omega }\hat{\tau }}}$， ${\alpha }_2=\frac{q_1{\delta }_1\tilde{x}{\text{e}}^{i\hat{\omega }\hat{\tau }}}{i\hat{\omega }-{\delta }_2}$，

$D=\frac{1}{1+{\alpha }_1{\bar{\alpha }}_2+\hat{\tau }{\text{e}}^{-i\hat{\omega }\hat{\tau }}\left({\delta }_1-q_1{\delta }_1\tilde{x}-q_1{\delta }_1\tilde{x}{\alpha }_1\right)}$。

证明：设算子 $A$ 关于特征值 $i\hat{\omega }\hat{\tau }$ 的特征向量为 $q\left(\theta \right)={\left(1,{\alpha }_1\right)}^{{\rm T}}{\text{e}}^{i\hat{\omega }\hat{\tau }\theta }$，故由 $Aq\left(0\right)=i\hat{\omega }\hat{\tau }q\left(0\right)$， ${\alpha }_1$ 满足

$\hat{\tau }{A}_1\left(\begin{array}{l}1\\ {\alpha }_1\end{array}\right)+\hat{\tau }B\left(\begin{array}{l}1\\ {\alpha }_1\end{array}\right){\text{e}}^{-i\hat{\omega }\hat{\tau }}=i\hat{\omega }\hat{\tau }\left(\begin{array}{l}1\\ {\alpha }_1\end{array}\right)$

通过计算得到 ${\alpha }_1=\frac{-{\delta }_1-i\hat{\omega }+\left({\delta }_1-q_1{\delta }_1\tilde{x}\right){\text{e}}^{-i\hat{\omega }\hat{\tau }}}{q_1{\delta }_1\tilde{x}{\text{e}}^{-i\hat{\omega }\hat{\tau }}}$。同理可设算子 $A^{*}$ 关于特征值 $-i\hat{\omega }\hat{\tau }$ 的特征向量为 $q^{*}\left(s\right)=\bar{D}{\left(1,{\alpha }_2\right)}^{{\rm T}}{\text{e}}^{i\hat{\omega }\hat{\tau }s}$，由 $A^{*}{q}^{*}\left(0\right)=-i\hat{\omega }\hat{\tau }{q}^{*}\left(0\right)$，得到 ${\alpha }_2=\frac{q_1{\delta }_1\tilde{x}{\text{e}}^{i\hat{\omega }\hat{\tau }}}{i\hat{\omega }-{\delta }_2}$。

因为

$\begin{array}{c}\langle q^{*},q\rangle ={\bar{q}}^{*}\left(0\right)q\left(0\right)-{\displaystyle {\int }_{-1}^0{\displaystyle {\int }_{\xi =0}^{\theta }{\bar{q}}^{*}\left(\xi -\theta \right)\text{d}\eta }}\left(\theta ,0\right)q\left(\xi \right)\text{d}\xi \\ =D\left(1+{\alpha }_1{\bar{\alpha }}_2\right)-{\displaystyle {\int }_{-1}^0{\displaystyle {\int }_{\xi =0}^{\theta }D\left(1,{\bar{\alpha }}_2\right){\text{e}}^{-i\hat{\omega }\hat{\tau }(\xi -\theta )}}}\text{d}\eta \left(\theta ,0\right){\left(1,{\bar{\alpha }}_1\right)}^{{\rm T}}{\text{e}}^{i\hat{\omega }\hat{\tau }\xi }\text{d}\xi \\ =D\left(1+{\alpha }_1{\bar{\alpha }}_2\right)-D\left(1,{\bar{\alpha }}_2\right){\displaystyle {\int }_{-1}^0\text{d}\eta \left(\theta ,0\right)\theta }{\text{e}}^{i\hat{\omega }\hat{\tau }\theta }{\left(1,{\bar{\alpha }}_1\right)}^{{\rm T}}\end{array}$

$\begin{array}{l}=D\left(1+{\alpha }_1{\bar{\alpha }}_2\right)-D\left(1,{\bar{\alpha }}_2\right)\left[\hat{\tau }B\left(-1\right){\text{e}}^{-i\hat{\omega }\hat{\tau }}{\left(1,{\bar{\alpha }}_1\right)}^{{\rm T}}\right]\\ =D\left(1+{\alpha }_1{\bar{\alpha }}_2\right)+D\hat{\tau }{\text{e}}^{-i\hat{\omega }\hat{\tau }}\left({\delta }_1-q_1{\delta }_1\tilde{x}-q_1{\delta }_1\tilde{x}{\alpha }_1\right)\\ =D\left[1+{\alpha }_1{\bar{\alpha }}_2+\hat{\tau }{\text{e}}^{-i\hat{\omega }\hat{\tau }}\left({\delta }_1-q_1{\delta }_1\tilde{x}-q_1{\delta }_1\tilde{x}{\alpha }_1\right)\right]\end{array}$

要使得 $\langle q^{*},q\rangle =1$，只需要取

$D=\frac{1}{1+{\alpha }_1{\bar{\alpha }}_2+\hat{\tau }{\text{e}}^{-i\hat{\omega }\hat{\tau }}\left({\delta }_1-q_1{\delta }_1\tilde{x}-q_1{\delta }_1\tilde{x}{\alpha }_1\right)}$

即可。文献已证得 $\langle q^{*},\bar{q}\rangle =0$。

下面运用文献 [15] 中的方法，计算在 $\mu =0$ 时系统在中心流形 $C_0$ 的坐标。令 $v_t={\left(m_t,n_t\right)}^{{\rm T}}$ 为(14)在 $\tau =\hat{\tau }$ 时的解，定义

$z\left(t\right)=\langle q^{*},v_t\rangle ,W\left(t,\theta \right)=v_t\left(\theta \right)-2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\},$ (16)

在中心流形 $C_0$ 上，有 $W\left(t,\theta \right)=W\left(z\left(t\right),\bar{z}\left(t\right),\theta \right)$，其中

$W\left(z\left(t\right),\bar{z}\left(t\right),\theta \right)=W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\bar{z}+W_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+W_{30}\left(\theta \right)\frac{z^3}{6}+\cdots$

$z$ 和 $\bar{z}$ 是中心流形 $C_0$ 在 $q^{*}$ 和 ${\bar{q}}^{*}$ 方向上的局部坐标，分析在中心流形 $C_0$ 上抽象方程(15)的解 $v_t$ 可得

$\begin{array}{c}{z}^{\prime }\left(t\right)=i\hat{\omega }\hat{\tau }z\left(t\right)+\langle q^{*}\left(s\right),F\left(W\left(t,\cdot \right)+2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\}\right)\rangle \\ =i\hat{\omega }\hat{\tau }z\left(t\right)+{\bar{q}}^{*}\left(0\right)F\left(W\left(z,\bar{z},0\right)+2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\}\right)\\ :=i\hat{\omega }\hat{\tau }z\left(t\right)+{\bar{q}}^{*}\left(0\right)F_0\left(z,\bar{z}\right)\end{array}$

将上式重新表示为

$z^{\prime }\left(t\right)=i\hat{\omega }\hat{\tau }z\left(t\right)+g\left(z,\bar{z}\right)$，(17)

其中

$g\left(z,\bar{z}\right)={\bar{q}}^{*}\left(0\right)F_0\left(z,\bar{z}\right)=g_{20}\frac{z^2}{2}+g_{11}z\bar{z}+g_{02}\frac{{\bar{z}}^2}{2}+g_{21}\frac{z^2\bar{z}}{2}+\cdots$， (18)

从(15)和(17)中，我们得到

$\begin{array}{c}{W}^{\prime }={v^{\prime }}_t-z^{\prime }q-{\bar{z}}^{\prime }q=\{\begin{array}{l}AW-2\mathrm{Re}\left\{{\bar{q}}^{*}\left(0\right)F_{0}q\left(\theta \right)\right\},\theta \in \left[-1,0\right)\\ AW-2\mathrm{Re}\left\{{\bar{q}}^{*}\left(0\right)F_{0}q\left(\theta \right)+F_0\right\},\theta =0\end{array}\\ :=AW+H\left(z,\bar{z},\theta \right)\end{array}$

其中

$H\left(z,\bar{z},\theta \right)=H_{20}\frac{z^2}{2}+H_{11}z\bar{z}+H_{02}\frac{{\bar{z}}^2}{2}+\cdots$， (19)

展开上述级数并对比相应系数，我们得到

$\left(A-2i\hat{\omega }\hat{\tau }\right)W_{20}\left(\theta \right)=-H_{20}\left(\theta \right)$，

$A{W}_{11}\left(\theta \right)=-H_{11}\left(\theta \right)$，

$\left(A+2i\hat{\omega }\hat{\tau }\right)W_{02}\left(\theta \right)=-H_{02}\left(\theta \right)$，... (20)

由(16)式有

$\begin{array}{c}{v}_t\left(\theta \right)=W\left(t,\theta \right)+2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\}\\ =z\left(t\right)q\left(\theta \right)+\bar{z}\left(t\right)\bar{q}\left(\theta \right)+W_{20}\frac{z^2}{2}+W_{11}z\bar{z}+W_{02}\frac{{\bar{z}}^2}{2}+\cdots \end{array}$

从而有

${\bar{q}}^{*}\left(0\right)=D\left(1,{\bar{\alpha }}_2\right)$，

$m\left(t\right)=z+\bar{z}+W_{20}^{(1)}\left(0\right)\frac{z^2}{2}+W_{11}^{(1)}\left(0\right)z\bar{z}+W_{02}^{(1)}\left(0\right)\frac{{\bar{z}}^2}{2}+\cdots$，

$n\left(t\right)=z{\alpha }_1+\bar{z}{\bar{\alpha }}_1+W_{20}^{(2)}\left(0\right)\frac{z^2}{2}+W_{11}^{(2)}\left(0\right)z\bar{z}+W_{02}^{(2)}\left(0\right)\frac{{\bar{z}}^2}{2}+\cdots$，

$m\left(t-1\right)=z{\text{e}}^{-i\hat{\omega }\hat{\tau }}+\bar{z}{\text{e}}^{i\hat{\omega }\hat{\tau }}+W_{20}^{(1)}\left(-1\right)\frac{z^2}{2}+W_{11}^{(1)}\left(-1\right)z\bar{z}+W_{02}^{(1)}\left(-1\right)\frac{{\bar{z}}^2}{2}+\cdots$，

$n\left(t-1\right)=z{\alpha }_1{\text{e}}^{-i\hat{\omega }\hat{\tau }}+\bar{z}{\bar{\alpha }}_1{\text{e}}^{i\hat{\omega }\hat{\tau }}+W_{20}^{(2)}\left(-1\right)\frac{z^2}{2}+W_{11}^{(2)}\left(-1\right)z\bar{z}+W_{02}^{(2)}\left(-1\right)\frac{{\bar{z}}^2}{2}+\cdots$，

$F_0=\left(\begin{array}{l}\hat{\tau }\left[\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)m^2\left(t-1\right)+\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}{n}^2\left(t-1\right)+\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)m\left(t-1\right)n\left(t-1\right)+o\left(\rho \right)\right]\\ \hat{\tau }\left[p_2{\text{e}}^{-{\delta }_2\tau }\frac{n^2\left(t-1\right)}{m\left(t-1\right)+n\left(t-1\right)+\tilde{x}}-\frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}{n}^2\left(t-1\right)+o\left(\rho \right)\right]\end{array}\right)$

与(18)对比系数得到：

$\frac{g_{20}}{2}=\left[\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)+{\alpha }_1\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)+\frac{{\alpha }_1{}^2{\delta }_1\tilde{x}{q}_1{}^2}{2}+{\bar{\alpha }}_2{\alpha }_1{}^2\frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}\right]{\text{e}}^{-2i\hat{\omega }\hat{\tau }}D\hat{\tau }$

$g_{11}=\left[\left({\delta }_1\tilde{x}{q}_1{}^2-2{\delta }_1{q}_1\right)+{\alpha }_1{\bar{\alpha }}_1{\delta }_1\tilde{x}{q}_1{}^2+\left({\alpha }_1+{\bar{\alpha }}_1\right)\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)+2{\alpha }_1{\bar{\alpha }}_1{\bar{\alpha }}_2\frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}\right]D\hat{\tau }$

$\frac{g_{02}}{2}=\left[\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)+{\bar{\alpha }}_1\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)+\frac{{\bar{\alpha }}_1{}^2{\delta }_1\tilde{x}{q}_1{}^2}{2}+{\bar{\alpha }}_2{\bar{\alpha }}_1{}^2\frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}\right]{\text{e}}^{2i\hat{\omega }\hat{\tau }}D\hat{\tau }$

$\begin{array}{c}\frac{g_{21}}{2}=[\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)\left(2{\text{e}}^{-i\hat{\omega }\hat{\tau }}{W}_{11}^{(1)}\left(-1\right)+{\text{e}}^{i\hat{\omega }\hat{\tau }}{W}_{20}^{(1)}\left(-1\right)\right)+\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}\left(2{\alpha }_1{\text{e}}^{-i\hat{\omega }\hat{\tau }}{W}_{11}^{(2)}\left(-1\right)+{\bar{\alpha }}_1{\text{e}}^{i\hat{\omega }\hat{\tau }}{W}_{20}^{(2)}\left(-1\right)\right)\\ +\left(\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1\right)\left({\text{e}}^{-i\hat{\omega }\hat{\tau }}{W}_{11}^{(2)}\left(-1\right)+\frac{1}{2}{\text{e}}^{i\hat{\omega }\hat{\tau }}{W}_{20}^{(2)}\left(-1\right)+\frac{1}{2}{\bar{\alpha }}_1{\text{e}}^{i\hat{\omega }\hat{\tau }}{W}_{20}^{(1)}\left(-1\right)+{\alpha }_1{\text{e}}^{-i\hat{\omega }\hat{\tau }}{W}_{11}^{(1)}\left(-1\right)\right)\\ +{\bar{\alpha }}_2\frac{p_2{e}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}\left(2{\alpha }_1{\text{e}}^{-i\hat{\omega }\hat{\tau }}{W}_{11}^{(2)}\left(-1\right)+{\bar{\alpha }}_1{\text{e}}^{i\hat{\omega }\hat{\tau }}{W}_{20}^{(2)}\left(-1\right)\right)]D\hat{\tau }\end{array}$

接下来要计算 $g_{21}$ 的值，还需要计算 $W_{20}\left(\theta \right)$ 和 $W_{11}\left(\theta \right)$ 的值。对于 $\theta \in \left[-1,0\right)$，有

$\begin{array}{c}H\left(z,\bar{z},\theta \right)=-2\mathrm{Re}\left\{{\bar{z}}^{*}\left(0\right)F_{0}q\left(\theta \right)\right\}=-gq\left(\theta \right)-\bar{g}\bar{q}\left(\theta \right)\\ =-\left(g_{20}\frac{z^2}{2}+g_{11}z\bar{z}+g_{02}\frac{{\bar{z}}^2}{2}+\cdots \right)q\left(\theta \right)-\left({\bar{g}}_{20}\frac{{\bar{z}}^2}{2}+g_{11}z\bar{z}+{\bar{g}}_{02}\frac{z^2}{2}+\cdots \right)\bar{q}(\; \theta \; )\end{array}$

与(19)对比相应系数，得到

$H_{20}\left(\theta \right)=-g_{20}q\left(\theta \right)-{\bar{g}}_{02}\bar{q}\left(\theta \right),H_{11}\left(\theta \right)=-g_{11}q\left(\theta \right)-{\bar{g}}_{11}\bar{q}(\; \theta \; )$

由(20)与 $A$ 的定义有

${W^{\prime }}_{20}=2i\hat{\omega }\hat{\tau }{W}_{20}\left(\theta \right)+g_{20}q\left(0\right){\text{e}}^{i\hat{\omega }\hat{\tau }\theta }+{\bar{g}}_{02}\bar{q}\left(0\right){\text{e}}^{-i\hat{\omega }\hat{\tau }\theta }$

由常数变易法，解得

$\begin{array}{c}{W}_{20}\left(\theta \right)={\text{e}}^{{\displaystyle \int 2i\hat{\omega }\hat{\tau }\theta }}\left({\displaystyle \int \left(g_{20}q\left(\theta \right)+{\bar{g}}_{02}\bar{q}\left(\theta \right)\right){\text{e}}^{-{\displaystyle \int 2i\hat{\omega }\hat{\tau }\theta }}}\text{d}\theta +E_1\right)\\ =\frac{i{g}_{20}q\left(0\right){\text{e}}^{i\hat{\omega }\hat{\tau }\theta }}{\hat{\omega }\hat{\tau }}+\frac{i{\bar{g}}_{02}\bar{q}\left(0\right){\text{e}}^{-i\hat{\omega }\hat{\tau }\theta }}{3\hat{\omega }\hat{\tau }}+E_1{\text{e}}^{2i\hat{\omega }\hat{\tau }\theta }\end{array}$ (21)

同理有

$W_{11}\left(\theta \right)=-\frac{i{g}_{11}q\left(0\right){\text{e}}^{i\hat{\omega }\hat{\tau }\theta }}{\hat{\omega }\hat{\tau }}+\frac{i{\bar{g}}_{11}\bar{q}\left(0\right){\text{e}}^{-i\hat{\omega }\hat{\tau }\theta }}{\hat{\omega }\hat{\tau }}+E_2$ (22)

其中 $E_1$， $E_2$ 是二维向量，它们的值可以通过 $\theta =0$ 时 $H$ 的值来确定。当 $\theta =0$ 时，有

$H\left(z,\bar{z},0\right)=-2\mathrm{Re}\left\{{\bar{z}}^{*}\left(0\right)F_{0}q\left(0\right)\right\}+F_0$

因此

$H_{20}\left(0\right)=-g_{20}q\left(0\right)-{\bar{g}}_{02}\bar{q}\left(0\right)+2\hat{\tau }\left(\begin{array}{c}\left[\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}-{\delta }_1{q}_1+\frac{{\delta }_1\tilde{x}{q}_1{}^2}{2}{\alpha }_1^2+{\alpha }_1\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)\right]{\text{e}}^{-2i\hat{\omega }\hat{\tau }}\\ \frac{p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}{\alpha }_1^2{\text{e}}^{-2i\hat{\omega }\hat{\tau }}\end{array}\right)$

$H_{11}\left(0\right)=-g_{11}q\left(0\right)-{\bar{g}}_{11}\bar{q}\left(0\right)+\hat{\tau }\left(\begin{array}{c}{\delta }_1\tilde{x}{q}_1{}^2-2{\delta }_1{q}_1+{\delta }_1\tilde{x}{q}_1{}^2{\alpha }_1{\bar{\alpha }}_1+\left({\delta }_1\tilde{x}{q}_1{}^2-{\delta }_1{q}_1\right)\left({\alpha }_1+{\bar{\alpha }}_1\right)\\ \frac{2p_2{\text{e}}^{-{\delta }_2\tau }\left({\text{e}}^{-q_2\tilde{x}}-1\right)}{\tilde{x}}{\alpha }_1{\bar{\alpha }}_1\end{array}\right)$

从(20)和 $A$ 的定义可得

$\hat{\tau }{A}_1{W}_{20}\left(0\right)+\hat{\tau }B{W}_{20}\left(-1\right)=2i\hat{\omega }\hat{\tau }{W}_{20}\left(0\right)-H_{20}(\; 0\; )$

$\hat{\tau }{A}_1{W}_{\text{11}}\left(0\right)+\hat{\tau }B{W}_{\text{11}}\left(-1\right)=-H_{\text{11}}(\; 0\; )$

将(21)和(22)带入上面两个方程，可以分别得到

$\begin{array}{c}{E}_1=-\frac{1}{\hat{\tau }}{\left(A_1+B{\text{e}}^{-2i\hat{\omega }\hat{\tau }}-2i\hat{\omega }I\right)}^{-1}[H_{20}\left(0\right)+\frac{i}{\hat{\omega }}{g}_{20}{A}_{1}q\left(0\right)+\frac{i}{3\hat{\omega }}{\bar{g}}_{02}{A}_1\bar{q}\left(0\right)+\frac{i}{\hat{\omega }}{g}_{20}{\text{e}}^{-i\hat{\omega }\hat{\tau }}Bq\left(0\right)\\ +\frac{i}{3\hat{\omega }}{\bar{g}}_{02}{\text{e}}^{i\hat{\omega }\hat{\tau }}B\bar{q}\left(0\right)+2g_{20}q\left(0\right)+\frac{2}{3}{\bar{g}}_{02}\bar{q}\left(0\right)]\end{array}$

$E_2=-\frac{1}{\hat{\tau }}{\left(A_1+B\right)}^{-1}\left[H_{11}\left(0\right)+\frac{2}{\hat{\omega }}{A}_1\mathrm{Im}\left(g_{11}q\left(0\right)\right)+\frac{2}{\hat{\omega }}B\mathrm{Im}\left(g_{11}q\left(0\right)\right){\text{e}}^{-i\hat{\omega }\hat{\tau }}\right]$

其中， $I$ 是三阶单位矩阵，从而可以求出 $g_{21}$ 的值，下面式子的值也可以通过计算得到：

$C_1\left(0\right)=\frac{i}{2\hat{\omega }\hat{\tau }}\left(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{{\left|g_{02}\right|}^2}{3}\right)+\frac{g_{21}}{2}$

${\mu }_2=-\frac{\mathrm{Re}\left(C_1\left(0\right)\right)}{\mathrm{Re}{\lambda }^{\prime }(\; \tau \; ^\; )}$

${\beta }_2=2\mathrm{Re}\left(C_1(\; 0\; )\right)$

定理4：分支方向由 ${\mu }_2$ 决定，分支周期解的稳定性由 ${\beta }_2$ 决定，具体有：

(1) 如果 ${\mu }_2>0$ ( ${\mu }_2<0$ )，则系统的Hopf分支是上临界(下临界)的。

(2) 如果 ${\beta }_2<0$ ( ${\beta }_2>0$ )，则分支周期解是稳定(不稳定)的。

4. 数值模拟与总结

为了更好的展示模型的动力学稳定性，通过选取不同的参数应用MATLAB对平衡点 ${\tilde{E}}_1$ 的稳定性进行数值模拟。当 ${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$ 时，分别选取不同的参数，它们的图像如图1所示。

图1所选的参数值为 $p_1=\text{0}\text{.8}$， $p_2=\text{0}\text{.8}$， ${\delta }_{\text{1}}=\text{0}\text{.03}$， ${\delta }_{\text{2}}=\text{0}\text{.03}$， $q_1=1\times {10}^{-8}$， $q_2=1\times {10}^{-\text{8}}$， $\tau =1\text{00}$。分别取初值 $x_0=\text{5}\times {10}^6$， $y_0=1\times {10}^{\text{6}}$ ； $x_0=\text{1}\times {10}^{\text{5}}$， $y_0=1\times {10}^{\text{6}}$ 得到图中的解曲线。图1表明：当

${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$ 时， ${\tilde{E}}_1$ 是局部渐近稳定的。

当 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$， ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$ 时，分别选取不同的参数，它们的图像如图2所示：

Figure 1. The solution curves of ${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$

图1. ${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$ 时的解曲线

Figure 2. The solution curves of $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$, ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$

图2. $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$， ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$ 时的解曲线

图2所选的参数值为 $p_1=\text{0}\text{.8}$， $p_2=\text{0}\text{.8}$， ${\delta }_{\text{1}}=\text{0}\text{.01}$， ${\delta }_{\text{2}}=\text{0}\text{.05}$， $q_1=1\times {10}^{-8}$， $q_2=1\times {10}^{-\text{8}}$， $\tau =1\text{00}$。

分别取初值 $x_0=5\times {10}^6$， $y_0=1\times {10}^{\text{6}}$ ； $x_0=\text{1}\times {10}^{\text{5}}$， $y_0=1\times {10}^{\text{7}}$。图2表明：当 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$， ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$ 时， ${\tilde{E}}_1$ 处出现了周期解。

5. 结论

本文建立了时滞微分方程模型研究Wolbachia氏菌在蚊群中的传播，首先考虑了系统的非负性和有

界性，讨论了平衡点的存在性和在各种不同的参数条件下的稳定性态。当 ${\text{e}}^{{\delta }_{\text{1}}\tau }<\frac{p_{\text{1}}}{{\delta }_{\text{1}}}\le {\text{e}}^{{\delta }_{\text{1}}\tau +2}$ 时， ${\tilde{E}}_1$ 是局部

渐近稳定的。如果初值的选取能够使得解趋于 ${\tilde{E}}_1$，则感染蚊群会持续存在，未感染蚊群会灭绝，Wolbachia

能够成功入侵到整个蚊群当中。当 $\frac{p_1}{{\delta }_1}>{\text{e}}^{{\delta }_1\tau +2}$， ${\left(1+{\delta }_1\tau -\ln \frac{p_1}{{\delta }_1}\right)}^{\text{2}}<1+\frac{{\delta }_2}{{\delta }_1}$ 时， ${\tilde{E}}_1$ 处出现了周期解，感染蚊

群和未感染蚊群的数量最终会趋于动态平衡。在这种情况下，未感染蚊群会灭绝，感染蚊群会持续存在，Wolbachia能够成功入侵到整个蚊群当中。

基金项目

国家自然科学基金(11771104)资助项目。