一类具有藻毒素动态影响的水域生态模型的动力学问题研究

Study on the Dynamics of an Aquatic Ecological Model with Dynamic Effects of Algo Toxins

摘要: 在考虑水域中的藻毒素是变量的情况下，构建一类具有藻毒素动态影响的水域生态模型，探索藻毒素对产毒藻类与鱼类生长共存模式的影响机制。数学理论工作主要调查了模型平衡点的存在性、稳定性和Hopf分岔动力学行为，推导出模型诱发Hopf分岔的临界条件。数值模拟工作既验证了理论推导结果的有效性，又阐明模型分岔动力学行为的演化特性，并揭示了过高或过低的藻毒素接触率参数将会导致产毒藻类种群和鱼类种群从常稳态持久生存模式变为周期振荡共存模式。总而言之，希望这些研究结果既有利于促进藻鱼生态系统研究的快速发展，又可为生物操纵控藻技术的应用提供一定的理论基础。

Abstract: Considering that algal toxins in water were variable, an aquatic ecological model with dynamic effects of algal toxins was constructed, which could explore the impact mechanism of algal toxins on the growth coexistence mode of toxin producing algae and fish. The mathematical theory work mainly investigated the existence, stability, and Hopf bifurcation dynamic behavior of model equilibrium points, and derived the critical conditions for the model to induce Hopf bifurcation. The numerical simulation work not only verified the effectiveness of the theoretical derivation results, but also elucidated the evolutionary characteristics of the bifurcation dynamics behavior of the model, and revealed that excessively high or low algal toxin exposure rates would cause the toxic algal population and fish population to shift from a stable coexistence mode to a cyclic oscillation mode. In conclusion, it was hoped that these research results would not only promote the rapid development of algae fish ecosystem research, but also provide a theoretical basis for the application of biological manipulation algae control technology.

1. 引言

在河流、湖泊和海洋等水域中，生活着大量缺乏有效移动能力，随水流而动的漂浮生物。它们种类繁多，形态各异，隶属于各个门类，统称为浮游生物。几乎所有的水生生物都以浮游生物为基础，浮游生物是在所有水生环境(即湖泊、河流、河口和海洋)表面附近自由漂浮的最丰富的生命形式 [1] [2] 。

浮游生物中的浮游植物不仅是所有水生食物链的基础，它们在吸收周围环境中的二氧化碳后，为人类和其他活体动物产生大量氧气，从而对其他生物拥有非常多有利的影响 [3] 。但此外，浮游植物也会造成不好的影响。浮游植物种群最常见的特征是，在一定时间周期内，细胞快速增殖，生物量迅速增加，随后一段时间内种群数量几乎同样迅速减少。在这种被称为“开花”的浮游植物种群密度的快速变化现象下，其中一些水华由于高生物量的积累或毒性的存在，更恰当地来说应当被称为“有害藻类水华”(HABs, [4] )，对海洋生态系统或人类健康有害，并可能造成巨大的社会经济损害。

为了人为地控制并缓解“有害藻类水华”对生态环境的影响，建立有效的生态系统动力学模型并加以分析，是我们能够做到的有效手段。在这之中，猎物和捕食物种之间的相互作用是生态系统数学模型中的一个重点研究问题。生态系统通常会受到人类开发活动的干扰，如渔业、野生动物和林业管理等领域中常见的生物资源开发和收获 [5] [6] [7] [8] 。于是，为了探索各种干扰对生态系统的影响，大多数研究人员都专注于生物系统的分析和建模以及收获 [9] 。此外，与其他生物资源一样，捕食–被捕食系统是生态系统中的重要生物资源。因此，利用数学模型研究捕食–被捕食系统是本世纪的一个重要研究领域。

减少浮游植物释放有毒物质导致的浮游动物的放牧压力是过去几十年中最有趣的研究课题之一。Chattopadhyay等人在论文 [10] 中提出了一类具有毒浮游植物和浮游动物相互作用的数学模型，探讨了有毒浮游植物在有害藻华环境中的作用机制。基于Chattopadhyay等人的想法，Tapan Saha在论文 [11] 中研究了类似的模型，并假设捕食功能反应函数和毒素影响函数由相同类型的函数描述，即Holling II型函数。在大多数文章中，作者都假设捕食功能反应函数和毒素影响函数是不同的类型。但由于浮游植物是水生环境中浮游动物最有利的食物来源，Holling II型功能形式是描述捕食规律的合理假设。所以无论是有助于浮游动物物种的生长，还是由于有毒物质的存在而抑制放牧速率，Tapan Saha的假设都是很合理的。这之后，Rana等人在论文 [12] 中研究了纳米颗粒对浮游植物–浮游动物模型动力学性态的影响，通过假设纳米颗粒会导致猎物物种(浮游植物)的生长速度降低，将纳米颗粒引入模型。M. Hossain等人 [13] 研究了纳米颗粒对一类(浮游植物–浮游动物–鱼类)三营养食物链模型动力学的影响。他们将纳米颗粒的影响纳入三物种食物链模型中，并假设猎物种群的增长由于与纳米颗粒的相互作用而减少，揭示模型动力学性态的演化规律。此外，这些研究 [14] - [20] 也关注并探讨了毒素对生态系统稳定性和生物多样性的影响。

以上各位研究人员的探索过程对含藻毒素环境下的生物生态系统的研究都起到了很大的推动作用，但遗憾的是，研究内容仅限于常量毒素条件下以及固定输入速率的环境毒素影响下的生态模型，因此我们对其中模型的优缺点进行分析，然后选择并自主构建一类适合藻毒素自然变化规律的水域生态模型，对其动力学性态进行定性分析与数值模拟研究。

2. 生态模型的构建

基于各种参考文献中所建的模型，我们引入变量 $p\left(t\right)$ 表示水环境中藻毒素浓度，提出一类具有藻毒素动态影响的水域生态模型，其表示如下：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{K}\right)-\frac{axy}{1+b_{1}x},\\ \frac{\text{d}y}{\text{d}t}=\frac{eaxy}{1+b_{1}x}-dy-\frac{\beta py}{1+b_{2}x},\\ \frac{\text{d}p}{\text{d}t}=A-\beta py-mp\text{.}\end{array}$ (1)

其中 $x\left(t\right)$ 和 $y\left(t\right)$ 分别代表产毒藻类和鱼类的种群密度， $p\left(t\right)$ 表示水环境中的藻毒素浓度。其他参数所代表的意义见表1。

Table 1. The meanings of various parameters

表1. 关于各参数的含义

我们可以看到，在该生态数学模型中，代表环境藻毒素浓度 $p\left(t\right)$ 的增长速率为固定的输入速率A，这显然不符合我们所要研究的由藻类产生毒素影响下的环境。于是我们将之修正为与产毒藻类种群密度相关的 $\alpha x$ ，其中 $\alpha$ 为藻毒素的生产速率。于是我们得到以下模型：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{K}\right)-\frac{axy}{1+b_{1}x},\\ \frac{\text{d}y}{\text{d}t}=\frac{eaxy}{1+b_{1}x}-dy-\frac{\beta py}{1+b_{2}x},\\ \frac{\text{d}p}{\text{d}t}=\alpha x-\beta py-mp\text{.}\end{array}$ (2)

针对所建动态模型(2)，我们主要研究模型平衡点的存在性和稳定性，探索Hopf分岔的存在性、方向性和稳定性，并给出相应的阈值条件。同时，对模型特定动力学行为进行数值模拟试验，验证理论结果的可行性，揭示关键参数对模型动力学性态演化的影响机制。

3. 平衡点的存在性与稳定性

平衡点是模型(2)的特殊解，它将表现出模型(2)的丰富动力学性质，因此，本节将讨论模型(2)中所有可能平衡点的存在性和稳定性。

3.1. 平衡点的存在性

显然，模型(2)最多有四个可能的平衡点：平凡的灭绝平衡点 $E_0\left(0,0,0\right)$ ，一个无鱼类平衡点 $E_1\left(K,0,\frac{\alpha K}{m}\right)$ ，两个共存平衡点 $E_1^{*}\left(x_1,y_1,p_1\right)$ 和 $E_2^{*}\left(x_2,y_2,p_2\right)$ 。

其中平衡点E₀和E₁无条件存在，而对于共存平衡点 $E^{*}$ ，x₁和x₂是以下方程的正实数根：

$\frac{\beta r{b}_1}{eaK}\left(ea-d{b}_1\right)x^3+{\rho }_1{x}^2+{\rho }_{2}x+{\rho }_3=0,$ (3)

其中，

${\rho }_1=\alpha \left(d{b}_1{b}_2+\beta b_1-ea{b}_2\right)+\frac{\beta r}{ea}\left(\frac{ea-2d{b}_1}{K}-ea{b}_1+d{b}_1^2\right),$ (4)

${\rho }_2=\alpha \left(d{b}_2+\beta \right)-m\left(ea-d{b}_1\right)+\frac{\beta r}{ea}\left(2d{b}_1-\frac{d}{K}-ea\right),$ (5)

${\rho }_3=\frac{\beta rd}{ea}+md\text{.}$ (6)

相对应的y₁，y₂，p₁，p₂为，

$\{\begin{array}{l}{y}_1=\frac{r}{ea}\left(1+b_1{x}_1\right)\left(1-\frac{x_1}{K}\right),\\ p_1=\frac{\left(ea-d{b}_1\right)x_1-d}{\left(d{b}_1{b}_2+\beta b_1-ea{b}_2\right)x_1+d{b}_2\beta }\text{.}\end{array}$ (7)

$\{\begin{array}{l}{y}_2=\frac{r}{ea}\left(1+b_1{x}_2\right)\left(1-\frac{x_2}{K}\right),\\ p_2=\frac{\left(ea-d{b}_1\right)x_2-d}{\left(d{b}_1{b}_2+\beta b_1-ea{b}_2\right)x_2+d{b}_2\beta }\text{.}\end{array}$ (8)

共存平衡点的存在显然是有条件的，基于对生物学意义的考虑，我们必须满足条件 $0<x_i\le K$ ， $y_i>0$ ， $p_i>0$ ， $i=1,2$ 。则根据(7) (8)可得， $ea-d{b}_1>0$ 。共存平衡点的具体存在条件见表2。共存平衡点对模型的动态行为具有决定性的影响，例如，如果不存在共存平衡点，鱼类种群最终会灭绝，模型(2)将不稳定。

Table 2. The existence conditions for equilibrium points

表2. 平衡点的存在条件

表2为模型平衡点存在的详细条件，表格中对应的各参数详细如下：

$B=\frac{\beta r{b}_1}{eaK}\left(ea-d{b}_1\right),$ (9)

$\Delta =27B^2{\rho }_3^2-18B{\rho }_1{\rho }_2{\rho }_3+4B{\rho }_2^3+4{\rho }_1^3{\rho }_3-{\rho }_1^2{\rho }_2^2.$ (10)

3.2. 平衡点的稳定性

通过雅可比矩阵特征值实部的符号可以得到平衡点的稳定性，首先得到模型(2)的雅可比矩阵表达式为：

$J_{E\left(x,y,p\right)}=\left(\begin{array}{ccc}r-\frac{2r}{K}x-\frac{eay}{{\left(1+b_{1}x\right)}^2}& -\frac{eax}{1+b_{1}x}& 0\\ \frac{eay}{{\left(1+b_{1}x\right)}^2}& \frac{eax}{1+b_{1}x}-d-\frac{\beta p}{1+b_{2}p}& -\frac{\beta y}{{\left(1+b_{2}p\right)}^2}\\ \alpha & -\beta p& -\beta y-m\end{array}\right)$ ,

由此我们得到了一些关于平衡点稳定性的定理。

定理1.平衡点 $E_0\left(0,0,0\right)$ 始终存在，且始终不稳定。

证明：将平衡点 $E_0\left(0,0,0\right)$ 带入雅可比矩阵，我们可以得到：

$J_{E_0\left(0,0,0\right)}=\left(\begin{array}{ccc}r& 0& 0\\ 0& -d& 0\\ \alpha & 0& -m\end{array}\right)$ ,

显然该矩阵为下三角矩阵，我们易得对角线上的三个特征根，分别为r，-d，-m。由于其中一个特征根 $r>0$ ，则模型(2)在平衡点 $E_0\left(0,0,0\right)$ 附近始终不稳定。

定理2.平衡点 $E_1\left(K,0,\frac{\alpha K}{m}\right)$ 始终存在，且当 $d>\frac{eaK}{1+b_{1}K}-\frac{\alpha \beta K}{m+\alpha b_{2}K}$ 时，平衡点E₁是局部渐近稳定的，当 $d<\frac{eaK}{1+b_{1}K}-\frac{\alpha \beta K}{m+\alpha b_{2}K}$ 时，平衡点E₁是不稳定的。

证明：将平衡点 $E_1\left(K,0,\frac{\alpha K}{m}\right)$ 带入雅可比矩阵，我们可以得到：

$J_{E_1\left(K,0,\frac{\alpha K}{m}\right)}=\left(\begin{array}{ccc}-r& -\frac{eaK}{1+b_{1}K}& 0\\ 0& \frac{eaK}{1+b_{1}K}-d-\frac{\alpha \beta K}{m+\alpha b_{2}K}& 0\\ \alpha & \frac{-\alpha \beta K}{m}& -m\end{array}\right)$ ,

我们可以容易解得该矩阵的三个特征根，分别为 ${\lambda }_1=-r,{\lambda }_2=\frac{eaK}{1+b_{1}K}-m-\frac{\alpha \beta K}{m+\alpha b_{2}K},{\lambda }_3=-m$ 。由于 ${\lambda }_1<0,{\lambda }_3<0$ ，为了使平衡点 $E_1\left(K,0,\frac{\alpha K}{m}\right)$ 是局部渐近稳定的，则应使 ${\lambda }_2<0$ ，即 $d>\frac{eaK}{1+b_{1}K}-\frac{\alpha \beta K}{m+\alpha b_{2}K}$ 。

定理3.当共存平衡点 $E_1^{*}\left(x_1,y_1,p_1\right)$ 和 $E_2^{*}\left(x_2,y_2,p_2\right)$ 存在时，且当 ${\sigma }_1>0,{\sigma }_3>0,{\sigma }_1{\sigma }_2>{\sigma }_3$ 时，平衡点具有局部渐近稳定性。

证明：将平衡点 $E^{*}\left(x^{*},y^{*},p^{*}\right)$ 带入雅可比矩阵，我们可以得到：

$J_{E^{\text{*}}\left(x^{\text{*}},y^{\text{*}},p^{\text{*}}\right)}=\left(\begin{array}{ccc}r-\frac{2r}{K}{x}^{\text{*}}-\frac{ea{y}^{\text{*}}}{{\left(1+b_1{x}^{\text{*}}\right)}^2}& -\frac{ea{x}^{\text{*}}}{1+b_1{x}^{\text{*}}}& 0\\ \frac{ea{y}^{\text{*}}}{{\left(1+b_1{x}^{\text{*}}\right)}^2}& \frac{ea{x}^{\text{*}}}{1+b_1{x}^{\text{*}}}-d-\frac{\beta p^{\text{*}}}{1+b_2{p}^{\text{*}}}& -\frac{\beta y^{\text{*}}}{{\left(1+b_2{p}^{\text{*}}\right)}^2}\\ \alpha & -\beta p^{\text{*}}& -\beta y^{\text{*}}-m\end{array}\right)$ ,

再将平衡点带入模型(2)中，我们能够得到 $\frac{ea{x}^{*}}{1+b_1{x}^{*}}-d-\frac{\beta p^{*}}{1+b_2{p}^{*}}=0$ ，因此我们将雅可比矩阵化简为：

$J_{E^{*}\left(x^{*},y^{*},p^{*}\right)}=\left(\begin{array}{ccc}r-\frac{2r}{K}{x}^{*}-\frac{ea{y}^{*}}{{\left(1+b_1{x}^{*}\right)}^2}& -\frac{ea{x}^{*}}{1+b_1{x}^{*}}& 0\\ \frac{ea{y}^{*}}{{\left(1+b_1{x}^{*}\right)}^2}& 0& -\frac{\beta y^{*}}{{\left(1+b_2{p}^{*}\right)}^2}\\ \alpha & -\beta p^{*}& -\beta y^{*}-m\end{array}\right)$ ,

下记为：

$J_{E^{\text{*}}\left(x^{\text{*}},y^{\text{*}},p^{\text{*}}\right)}=\left(\begin{array}{ccc}{M}_1& M_2& 0\\ M_3& 0& M_4\\ M_5& M_6& M_7\end{array}\right)$ ,

求解该矩阵的特征根则，

$\left(M_1-\lambda \right)\lambda \left(M_7-\lambda \right)+\left(M_1-\lambda \right)M_4{M}_6+M_2{M}_3\left(M_7-\lambda \right)-M_2{M}_4{M}_5=0$

化简得，

$\begin{array}{l}{\lambda }^3-\left(M_1+M_7\right){\lambda }^2+\left(M_1{M}_7-M_4{M}_6-M_2{M}_3\right)\lambda \\ +\left(M_1{M}_4{M}_6+M_2{M}_3{M}_7-M_2{M}_4{M}_5\right)=0\end{array}$

令

${\sigma }_1=-\left(M_1+M_7\right),{\sigma }_2=M_1{M}_7-M_4{M}_6-M_2{M}_3,$

${\sigma }_3=M_1{M}_4{M}_6+M_2{M}_3{M}_7-M_2{M}_4{M}_5,$

因此，根据劳斯判据，为使其平衡点具有局部渐近稳定性，则应使 ${\sigma }_1>0,{\sigma }_3>0,{\sigma }_1{\sigma }_2>{\sigma }_3$ ，此时该平衡点具有局部渐近稳定性。

4. 局部分岔动力学分析

本节将详细讨论模型(2)的局部分岔动力学行为。由于是在藻毒素环境影响下的生态模型，显然作为鱼类与藻毒素的接触率参数 $\beta$ 非常关键，于是我们将主要讨论基于参数 $\beta$ 发生的Hopf分岔对模型动力学行为带来的影响，同时也对消耗率参数m做出了讨论。

4.1. Hopf分岔的存在性

定理4.当接触率 $\beta$ 超过临界值 ${\beta }^{*}$ 时，模型(2)在内部平衡点 $E^{*}\left(x^{*},y^{*},z^{*}\right)$ 周围经历Hopf分岔，且满足以下条件：

(1) $H\left({\beta }^{\text{*}}\right)\equiv {\sigma }_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{\text{*}}\right)-{\sigma }_3\left({\beta }^{\text{*}}\right)=0$ (11)

(2) ${{\sigma }^{\prime }}_3\left({\beta }^{\text{*}}\right)-{{\sigma }^{\prime }}_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{\text{*}}\right)-{\sigma }_1\left({\beta }^{\text{*}}\right){{\sigma }^{\prime }}_2\left({\beta }^{\text{*}}\right)\ne 0,{\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Re}\left(\lambda \left(\beta \right)\right)\right)|}_{\beta ={\beta }^{\text{*}}}\ne 0$ (12)

其中， $\lambda$ 为内部平衡点对应的雅可比矩阵的特征根。

证明：我们首先验证在条件(11)的情况下，雅可比矩阵在内部平衡点处特征根的情况。

当 $\beta ={\beta }^{*}$ 时，由 $H\left({\beta }^{*}\right)\equiv {\sigma }_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{\text{*}}\right)-{\sigma }_3\left({\beta }^{\text{*}}\right)=0$ ，则求特征根的方程可化为：

${\lambda }^3+{\sigma }_1{\lambda }^2+{\sigma }_2\lambda +{\sigma }_3={\lambda }^3+{\sigma }_1{\lambda }^2+{\sigma }_2\lambda +{\sigma }_1{\sigma }_2=\left({\lambda }^2+{\sigma }_2\right)\left(\lambda +{\sigma }_1\right)$

因此，可以得到 $\left({\lambda }^2+{\sigma }_2\right)\left(\lambda +{\sigma }_1\right)=0$ ，

解得三个特征根分别为 ${\lambda }_1=i\sqrt{{\sigma }_2},{\lambda }_2=-i\sqrt{{\sigma }_2},{\lambda }_3=-{\sigma }_1$ ，

显然此时拥有一对纯虚根和一个非零负实根。

接下来我们验证其横截性条件(12)。

对于一般的 $\beta$ ，特征方程的根有如下形式：

${\lambda }_1\left(\beta \right)={\varphi }_1\left(\beta \right)+i{\varphi }_2\left(\beta \right),$ (13)

${\lambda }_2\left(\beta \right)={\varphi }_1\left(\beta \right)-i{\varphi }_2\left(\beta \right),$ (14)

${\lambda }_3\left(\beta \right)=-{\sigma }_1\text{.}$ (15)

将(13) (14)分别带入特征根方程并求导得，

$\begin{array}{l}A\left(\beta \right){{\varphi }^{\prime }}_1\left(\beta \right)-B\left(\beta \right){{\varphi }^{\prime }}_2\left(\beta \right)+C\left(\beta \right)=0,\\ B\left(\beta \right){{\varphi }^{\prime }}_1\left(\beta \right)+A\left(\beta \right){{\varphi }^{\prime }}_2\left(\beta \right)+D\left(\beta \right)=0.\end{array}$

这里 $A\left(\beta \right),B\left(\beta \right),C\left(\beta \right),D\left(\beta \right)$ 分别为，

$\begin{array}{l}A\left(\beta \right)=3{\varphi }_1^2\left(\beta \right)+2{\sigma }_1\left(\beta \right){\varphi }_1\left(\beta \right)+{\sigma }_2\left(\beta \right)-3{\varphi }_2^2\left(\beta \right),\\ B\left(\beta \right)=6{\varphi }_1\left(\beta \right){\varphi }_2\left(\beta \right)+2{\sigma }_1\left(\beta \right){\varphi }_2\left(\beta \right),\\ C\left(\beta \right)={\varphi }_1^2\left(\beta \right){{\sigma }^{\prime }}_1\left(\beta \right)+{{\sigma }^{\prime }}_2\left(\beta \right){\varphi }_1\left(\beta \right)+{{\sigma }^{\prime }}_3\left(\beta \right)-{{\sigma }^{\prime }}_1\left(\beta \right){\varphi }_2^2\left(\beta \right),\\ D\left(\beta \right)=2{\varphi }_1\left(\beta \right){\varphi }_2\left(\beta \right){{\sigma }^{\prime }}_1\left(\beta \right)+{{\sigma }^{\prime }}_2\left(\beta \right){\varphi }_2\left(\beta \right)\text{.}\end{array}$

由 ${\varphi }_1\left({\beta }^{*}\right)=0,{\varphi }_2\left({\beta }^{*}\right)=\sqrt{{\sigma }_2\left({\beta }^{*}\right)}$ ，则

$\begin{array}{l}A\left({\beta }^{\text{*}}\right)=-2{\sigma }_2\left({\beta }^{\text{*}}\right),\\ B\left({\beta }^{\text{*}}\right)=2{\sigma }_1\left({\beta }^{\text{*}}\right)\sqrt{{\sigma }_2\left({\beta }^{\text{*}}\right)},\\ C\left({\beta }^{\text{*}}\right)={{\sigma }^{\prime }}_3\left({\beta }^{\text{*}}\right)-{{\sigma }^{\prime }}_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{\text{*}}\right),\\ D\left({\beta }^{\text{*}}\right)={{\sigma }^{\prime }}_2\left({\beta }^{\text{*}}\right)\sqrt{{\sigma }_2\left({\beta }^{\text{*}}\right)}\text{.}\end{array}$

将之带入横截条件(12)，可得

$\begin{array}{l}{\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Re}\left(\lambda \left(\beta \right)\right)\right)|}_{\beta ={\beta }^{*}}=-\frac{A\left({\beta }^{*}\right)C\left({\beta }^{*}\right)+B\left({\beta }^{*}\right)D\left({\beta }^{*}\right)}{A{\left({\beta }^{*}\right)}^2+B{\left({\beta }^{*}\right)}^2}\\ =-\frac{-2{\sigma }_2\left({\beta }^{*}\right)\left({{\sigma }^{\prime }}_3\left({\beta }^{\text{*}}\right)-{{\sigma }^{\prime }}_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{*}\right)\right)+2{\sigma }_1\left({\beta }^{*}\right)\sqrt{{\sigma }_2\left({\beta }^{*}\right)}{{\sigma }^{\prime }}_2\left({\beta }^{\text{*}}\right)\sqrt{{\sigma }_2\left({\beta }^{*}\right)}}{4{\sigma }_2^2\left({\beta }^{*}\right)+4{\sigma }_1^2\left({\beta }^{*}\right){\sigma }_2\left({\beta }^{*}\right)}\\ =\frac{{{\sigma }^{\prime }}_3\left({\beta }^{\text{*}}\right)-{{\sigma }^{\prime }}_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{*}\right)-{\sigma }_1\left({\beta }^{*}\right){{\sigma }^{\prime }}_2\left({\beta }^{\text{*}}\right)}{2{\sigma }_2\left({\beta }^{*}\right)+2{\sigma }_1^2\left({\beta }^{*}\right)}\ne 0\end{array}$

因此当 ${{\sigma }^{\prime }}_3\left({\beta }^{\text{*}}\right)-{{\sigma }^{\prime }}_1\left({\beta }^{\text{*}}\right){\sigma }_2\left({\beta }^{\text{*}}\right)-{\sigma }_1\left({\beta }^{\text{*}}\right){{\sigma }^{\prime }}_2\left({\beta }^{\text{*}}\right)\ne 0$ 时，横截性条件(12)成立。所以，模型(2)将在 $\beta ={\beta }^{*}$ 处发生Hopf分岔。

4.2. Hopf分岔的方向性和稳定性

本部分将利用规范型理论，推导出模型(2)Hopf分岔的方向性、稳定性和分岔周期解的周期。根据规范型理论，我们首先利用内部平衡点 $E^{*}$ 处的雅可比矩阵，计算出当 $\beta ={\beta }^{*}$ 时，分别对应 ${\lambda }_1,{\lambda }_3$ 的特征向量，其中 ${\lambda }_1=i\omega ,{\lambda }_3=-{\sigma }_1,\omega =\sqrt{{\sigma }_2}$ 。

我们得到对应的特征向量如下：

$v_1=\left(\begin{array}{c}{b}_{11}-i{b}_{12}\\ b_{21}-i{b}_{22}\\ b_{31}-i{b}_{32}\end{array}\right),v_3=\left(\begin{array}{c}{b}_{13}\\ b_{23}\\ b_{33}\end{array}\right).$

其中，

$b_{11}=1,b_{12}=0,b_{13}=1,b_{21}=-\frac{M_1}{M_2},b_{22}=-\frac{\omega }{M_2},b_{23}=-\frac{{\sigma }_1+M_1}{M_2},$

$b_{31}=-\frac{{\omega }^2+M_2{M}_3}{M_4},b_{32}=\frac{\omega M_1}{M_2{M}_4},b_{33}=\frac{{\sigma }_1^2+{\sigma }_1{M}_1-M_2{M}_3}{M_2{M}_4},$

接下来我们对原模型做出以下变换，令

$\begin{array}{l}x=x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3,\\ y=y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3,\\ p=p^{\text{*}}+b_{31}{s}_1+b_{32}{s}_2+b_{33}{s}_3,\end{array}$

则通过以上变换，我们得到新模型：

$\begin{array}{l}\frac{\text{d}{s}_1}{\text{d}t}=\frac{E_1\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_2{b}_{32}-E_3{b}_{22}}{D}\equiv Q^1,\\ \frac{\text{d}{s}_2}{\text{d}t}=\frac{E_1\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_2\left(b_{31}-b_{33}\right)+E_3\left(b_{21}-b_{23}\right)}{D}\equiv Q^2,\\ \frac{\text{d}{s}_3}{\text{d}t}=\frac{E_1\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)-E_2{b}_{32}+E_3{b}_{22}}{D}\equiv Q^3\text{.}\end{array}$ (16)

其中，

$D=b_{21}{b}_{32}-b_{22}{b}_{31}+b_{22}{b}_{33}-b_{23}{b}_{32},$

$\begin{array}{c}{E}_1=r\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)\left(1-\frac{x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3}{K}\right)\\ -\frac{ea\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)\left(y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3\right)}{1+b_1\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)},\end{array}$

$\begin{array}{c}{E}_2=\frac{ea\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)\left(y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3\right)}{1+b_1\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)}\\ -d\left(y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3\right)-\frac{\beta \left(p^{\text{*}}+b_{31}{s}_1+b_{32}{s}_2+b_{33}{s}_3\right)\left(y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3\right)}{1+b_2\left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)},\end{array}$

$\begin{array}{c}{E}_3=\alpha \left(x^{\text{*}}+b_{11}{s}_1+b_{12}{s}_2+b_{13}{s}_3\right)-\beta \left(p^{\text{*}}+b_{31}{s}_1+b_{32}{s}_2+b_{33}{s}_3\right)\left(y^{\text{*}}+b_{21}{s}_1+b_{22}{s}_2+b_{23}{s}_3\right)\\ -m\left(p^{\text{*}}+b_{31}{s}_1+b_{32}{s}_2+b_{33}{s}_3\right)\text{.}\end{array}$

我们将新模型(16)的雅可比矩阵写作：

$J=\left(\begin{array}{ccc}\frac{\partial Q^1}{\partial s_1}& \frac{\partial Q^1}{\partial s_2}& \frac{\partial Q^1}{\partial s_3}\\ \frac{\partial Q^2}{\partial s_1}& \frac{\partial Q^2}{\partial s_2}& \frac{\partial Q^2}{\partial s_3}\\ \frac{\partial Q^3}{\partial s_1}& \frac{\partial Q^3}{\partial s_2}& \frac{\partial Q^3}{\partial s_3}\end{array}\right),$

显然， $\left(0,0,0\right)$ 为模型(16)的平衡点。

由此我们可以将新模型的雅可比矩阵化简为：

$J=\left(\begin{array}{ccc}0& \omega & 0\\ \omega & 0& 0\\ 0& 0& D_1\end{array}\right),$

通过上述矩阵，我们计算表达式 $g_{11},g_{02},g_{20},G_{101},G_{110},G_{21},h_{11},h_{20},w,w_{20},w_{11},G_{21}$ ，如下：

$\begin{array}{c}{g}_{11}=\frac{1}{4}\left[\left(\frac{{\partial }^2{Q}^1}{\partial s_1^2}+\frac{{\partial }^2{Q}^2}{\partial s_2^2}\right)+i\left(\frac{{\partial }^2{Q}^2}{\partial s_1^2}+\frac{{\partial }^2{Q}^1}{\partial s_2^2}\right)\right]\\ =\frac{1}{4D}[E_{1xx}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xx}{b}_{32}-E_{3xx}{b}_{22}+E_{1yy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2yy}\left(b_{31}-b_{33}\right)\\ +E_{3yy}\left(b_{21}-b_{23}\right)+i(E_{1xx}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xx}\left(b_{31}-b_{33}\right)+E_{3xx}\left(b_{21}-b_{23}\right)\\ +E_{1yy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2yy}{b}_{32}-E_{3yy}{b}_{22})],\end{array}$

$\begin{array}{c}{g}_{02}=\frac{1}{4}\left[\left(\frac{{\partial }^2{Q}^1}{\partial s_1^2}-\frac{{\partial }^2{Q}^1}{\partial s_2^2}-2\frac{{\partial }^2{Q}^2}{\partial s_1\partial s_2}\right)+i\left(\frac{{\partial }^2{Q}^2}{\partial s_1^2}-\frac{{\partial }^2{Q}^2}{\partial s_2^2}+2\frac{{\partial }^2{Q}^1}{\partial s_1\partial s_2}\right)\right]\\ =\frac{1}{4D}[E_{1xx}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xx}{b}_{32}-E_{3xx}{b}_{22}-E_{1yy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-E_{2yy}{b}_{32}+E_{3yy}{b}_{22}\\ -2E_{1xy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)+2E_{2xy}\left(b_{31}-b_{33}\right)-2E_{3xy}\left(b_{21}-b_{23}\right)+i(E_{1xx}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)\\ -E_{2xx}\left(b_{31}-b_{33}\right)+E_{3xx}\left(b_{21}-b_{23}\right)-E_{1yy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)+E_{2yy}\left(b_{31}-b_{33}\right)\\ -E_{3yy}\left(b_{21}-b_{23}\right)+2E_{1xy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+2E_{2xy}{b}_{32}-2E_{3xy}{b}_{22})],\end{array}$

$\begin{array}{c}{g}_{20}=\frac{1}{4}\left[\left(\frac{{\partial }^2{Q}^1}{\partial s_1^2}-\frac{{\partial }^2{Q}^1}{\partial s_2^2}+2\frac{{\partial }^2{Q}^2}{\partial s_1\partial s_2}\right)+i\left(\frac{{\partial }^2{Q}^2}{\partial s_1^2}-\frac{{\partial }^2{Q}^2}{\partial s_2^2}-2\frac{{\partial }^2{Q}^1}{\partial s_1\partial s_2}\right)\right]\\ =\frac{1}{4D}[E_{1xx}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xx}{b}_{32}-E_{3xx}{b}_{22}-E_{1yy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-E_{2yy}{b}_{32}+E_{3yy}{b}_{22}\\ +2E_{1xy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-2E_{2xy}\left(b_{31}-b_{33}\right)+2E_{3xy}\left(b_{21}-b_{23}\right)+i(E_{1xx}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)\\ -E_{2xx}\left(b_{31}-b_{33}\right)+E_{3xx}\left(b_{21}-b_{23}\right)-E_{1yy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)+E_{2yy}\left(b_{31}-b_{33}\right)\\ -E_{3yy}\left(b_{21}-b_{23}\right)-2E_{1xy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-2E_{2xy}{b}_{32}+2E_{3xy}{b}_{22})],\end{array}$

$\begin{array}{c}{G}_{21}=\frac{1}{8}\left[\left(\frac{{\partial }^3{Q}^1}{\partial s_1^3}+\frac{{\partial }^3{Q}^1}{\partial s_1\partial s_2^2}+\frac{{\partial }^3{Q}^2}{\partial s_1^2\partial s_2}+\frac{{\partial }^3{Q}^2}{\partial s_2^3}\right)+i\left(\frac{{\partial }^3{Q}^2}{\partial s_1^3}+\frac{{\partial }^3{Q}^2}{\partial s_1\partial s_2^2}-\frac{{\partial }^3{Q}^1}{\partial s_1^2\partial s_2}-\frac{{\partial }^3{Q}^1}{\partial s_2^3}\right)\right]\\ =\frac{1}{8D}[E_{1xxx}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xxx}{b}_{32}-E_{3xxx}{b}_{22}+E_{1xyy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xyy}{b}_{32}-E_{3xyy}{b}_{22}\\ +E_{1xxy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xxy}\left(b_{31}-b_{33}\right)+E_{3xxy}\left(b_{21}-b_{23}\right)+E_{1yyy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)\\ -E_{2yyy}\left(b_{31}-b_{33}\right)+E_{3yyy}\left(b_{21}-b_{23}\right)]+\frac{i}{8D}[E_{1xxx}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xxx}\left(b_{31}-b_{33}\right)\\ +E_{3xxx}\left(b_{21}-b_{23}\right)+E_{1xyy}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xyy}\left(b_{31}-b_{33}\right)+E_{3xyy}\left(b_{21}-b_{23}\right)\\ -E_{1xxy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-E_{2xxy}{b}_{32}+E_{3xxy}{b}_{22}-E_{1yyy}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-E_{2yyy}{b}_{32}+E_{3yyy}{b}_{22}].\end{array}$

同时，我们可以得到，

$\omega =-\frac{\partial Q^1}{\partial s_2}=-\frac{1}{D}\left[E_{1y}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2y}{b}_{32}-E_{3y}{b}_{22}\right],$

$\begin{array}{c}{h}_{11}=\frac{1}{4}\left(\frac{{\partial }^2{Q}^3}{\partial s_1^2}+\frac{{\partial }^2{Q}^3}{\partial s_2^2}\right)\\ =\frac{1}{4D}[E_{1xx}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)-E_{2xx}{b}_{32}+E_{3xx}{b}_{22}\\ +E_{1yy}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)-E_{2yy}{b}_{32}+E_{3yy}{b}_{22}],\end{array}$

$\begin{array}{c}{h}_{20}=\frac{1}{4}\left(\frac{{\partial }^2{Q}^3}{\partial s_1^2}-\frac{{\partial }^2{Q}^3}{\partial s_2^2}-2i\frac{{\partial }^2{Q}^3}{\partial s_1\partial s_2}\right)\\ =\frac{1}{4D}[E_{1xx}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)-E_{2xx}{b}_{32}+E_{3xx}{b}_{22}-E_{1yy}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)\\ +E_{2yy}{b}_{32}-E_{3yy}{b}_{22}-E_{1xy}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)+E_{2xy}{b}_{32}-E_{3xy}{b}_{22}].\end{array}$

又有，

$D_1=\frac{\partial Q^3}{\partial s_3}=\frac{1}{D}\left[E_{1z}\left(b_{21}{b}_{32}-b_{22}{b}_{31}\right)-E_{2z}{b}_{32}+E_{3z}{b}_{22}\right]$

由上可以求得，

$w_{11}=-\frac{h_{11}}{D_1},w_{20}=-\frac{h_{20}}{D-2i\omega }.$

此外，

$\begin{array}{c}{G}_{110}=\frac{1}{2}\left[\left(\frac{{\partial }^2{Q}^1}{\partial s_1\partial s_3}+\frac{{\partial }^2{Q}^2}{\partial s_2\partial s_3}\right)+i\left(\frac{{\partial }^2{Q}^2}{\partial s_1\partial s_3}-\frac{{\partial }^2{Q}^1}{\partial s_2\partial s_3}\right)\right]\\ =\frac{1}{2D}[E_{1xz}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xz}{b}_{32}-E_{3xz}{b}_{22}+E_{1yz}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2yz}\left(b_{31}-b_{33}\right)\\ +E_{3yz}\left(b_{21}-b_{23}\right)]+\frac{i}{2D}[E_{1xz}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xz}\left(b_{31}-b_{33}\right)+E_{3xz}\left(b_{21}-b_{23}\right)\\ -E_{1yz}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)-E_{2yz}{b}_{32}+E_{3yz}{b}_{22}],\end{array}$

$\begin{array}{c}{G}_{101}=\frac{1}{2}\left[\left(\frac{{\partial }^2{Q}^1}{\partial s_1\partial s_3}-\frac{{\partial }^2{Q}^2}{\partial s_2\partial s_3}\right)+i\left(\frac{{\partial }^2{Q}^2}{\partial s_1\partial s_3}+\frac{{\partial }^2{Q}^1}{\partial s_2\partial s_3}\right)\right]\\ =\frac{1}{2D}[E_{1xz}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2xz}{b}_{32}-E_{3xz}{b}_{22}-E_{1yz}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)+E_{2yz}\left(b_{31}-b_{33}\right)\\ -E_{3yz}\left(b_{21}-b_{23}\right)]+\frac{i}{2D}[E_{1xz}\left(-b_{21}{b}_{33}+b_{23}{b}_{31}\right)-E_{2xz}\left(b_{31}-b_{33}\right)+E_{3xz}\left(b_{21}-b_{23}\right)\\ +E_{1yz}\left(b_{22}{b}_{33}-b_{23}{b}_{32}\right)+E_{2yz}{b}_{32}-E_{3yz}{b}_{22}],\end{array}$

$g_{21}=G_{21}+2G_{110}{w}_{11}+G_{101}{w}_{20}.$

其中，

$E_{1xx}=-\frac{2r}{K}-\frac{2ea{b}_{21}}{b_1{x}^{*}+1}+\frac{2ea{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{2ea{x}^{*}{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{2a{x}^{*}{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3},$

$\begin{array}{c}{E}_{2xx}=\frac{2ea{b}_{21}}{b_1{x}^{*}+1}-\frac{2ea{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{2ea{x}^{*}{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{2ea{x}^{*}{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}-\frac{2\beta b_{31}{b}_{21}}{b_2{p}^{*}+1}\\ +\frac{2\beta b_{31}^2{y}^{*}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{2\beta p^{*}{b}_{21}{b}_2{b}_{31}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{y}^{*}{b}_2^2{b}_{31}^2}{{\left(b_2{p}^{*}+1\right)}^3},\end{array}$

$E_{3xx}=-2\beta b_{31}{b}_{21},$

$E_{1yy}=0,$

$E_{2yy}=-\frac{2\beta b_{32}{b}_{22}}{b_2{p}^{*}+1}+\frac{2\beta b_{32}^2{y}^{*}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{2\beta p^{*}{b}_{22}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{y}^{*}{b}_2^2{b}_{32}^2}{{\left(b_2{p}^{*}+1\right)}^3},$

$E_{3yy}=-2\beta b_{32}{b}_{22},$

$E_{1xy}=-\frac{ea{b}_{22}}{b_1{x}^{*}+1}+\frac{ea{x}^{*}{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2},$

$\begin{array}{c}{E}_{2xy}=\frac{ea{b}_{22}}{b_1{x}^{*}+1}-\frac{ea{x}^{*}{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{\beta b_{31}{b}_{22}}{b_2{p}^{*}+1}+\frac{2\beta b_{31}{y}^{*}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{\beta b_{32}{b}_{21}}{b_2{p}^{*}+1}\\ +\frac{\beta p^{*}{b}_{21}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{\beta p^{*}{b}_{22}{b}_2{b}_{31}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{y}^{*}{b}_2^2{b}_{31}{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^3},\end{array}$

$E_{3xy}=-\beta b_{32}{b}_{21}-\beta b_{31}{b}_{22},$

$E_{1xxx}=\frac{6ea{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{6ea{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}-\frac{6ea{x}^{*}{b}_{21}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}+\frac{6ea{x}^{*}{y}^{*}{b}_1^3}{{\left(b_1{x}^{*}+1\right)}^4},$

$\begin{array}{c}{E}_{2xxx}=-\frac{6ea{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{6ea{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}+\frac{6ea{x}^{*}{b}_{21}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}-\frac{6ea{x}^{*}{y}^{*}{b}_1^3}{{\left(b_1{x}^{*}+1\right)}^4}+\frac{6\beta b_{31}^2{b}_{21}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}\\ -\frac{6\beta b_{31}^3{y}^{*}{b}_2^2}{{\left(b_2{p}^{*}+1\right)}^3}-\frac{6\beta p^{*}{b}_{21}{b}_2^2{b}_{31}^2}{{\left(b_2{p}^{*}+1\right)}^3}+\frac{6\beta p^{*}{y}^{*}{b}_2^3{b}_{31}^3}{{\left(b_2{p}^{*}+1\right)}^4},\end{array}$

$E_{3xxx}=0,$

$E_{1yyy}=0,$

$E_{2yyy}=\frac{6\beta b_{32}^2{b}_{22}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{6\beta b_{32}^3{y}^{*}{b}_2^2}{{\left(b_2{p}^{*}+1\right)}^3}-\frac{6\beta p^{*}{b}_{22}{b}_2^2{b}_{32}^2}{{\left(b_2{p}^{*}+1\right)}^3}+\frac{6\beta p^{*}{y}^{*}{b}_2^3{b}_{32}^3}{{\left(b_2{p}^{*}+1\right)}^4},$

$E_{3yyy}=0,$

$E_{1xxy}=\frac{2ea{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{2ea{x}^{*}{b}_{22}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3},$

$\begin{array}{c}{E}_{2xxy}=-\frac{2ea{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{2ea{x}^{*}{b}_{22}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}+\frac{4\beta b_{31}{b}_{21}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{2\beta b_{31}^2{b}_{22}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{6\beta b_{31}^2{y}^{*}{b}_2^2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^3}\\ -\frac{4\beta p^{*}{b}_{21}{b}_2^2{b}_{31}{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^3}-\frac{2\beta p^{*}{b}_{22}{b}_2^2{b}_{31}^2}{{\left(b_2{p}^{*}+1\right)}^3}+\frac{6\beta p^{*}{y}^{*}{b}_2^3{b}_{31}^2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^4},\end{array}$

$E_{3xxy}=0,$

$E_{1xyy}=0,$

$\begin{array}{c}{E}_{2xyy}=\frac{4\beta b_{31}{b}_{22}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{6\beta b_{31}{y}^{*}{b}_2^2{b}_{32}^2}{{\left(b_2{p}^{*}+1\right)}^3}+\frac{2\beta b_{32}^2{b}_{21}{b}_2}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{b}_{21}{b}_2^2{b}_{32}^2}{{\left(b_2{p}^{*}+1\right)}^3}\\ -\frac{4\beta p^{*}{b}_{22}{b}_2^2{b}_{31}{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^3}+\frac{6\beta p^{*}{y}^{*}{b}_2^3{b}_{31}{b}_{32}^2}{{\left(b_2{p}^{*}+1\right)}^4},\end{array}$

$E_{3xyy}=0,$

$E_{1y}=-\frac{ea{x}^{*}{b}_{22}}{b_1{x}^{*}+1},$

$E_{2y}=\frac{ea{x}^{*}{b}_{22}}{b_1{x}^{*}+1}-d{b}_{22}-\frac{\beta b_{32}{y}^{*}}{b_2{p}^{*}+1}-\frac{\beta p^{*}{b}_{22}}{b_2{p}^{*}+1}+\frac{\beta p^{*}{y}^{*}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2},$

$E_{3y}=-\beta p^{*}{b}_{22}-\beta b_{32}{y}^{*}-m{b}_{32},$

$E_{1z}=r\left(1-\frac{x^{*}}{K}\right)-\frac{r{x}^{*}}{K}-\frac{ea{y}^{*}}{b_1{x}^{*}+1}-\frac{ea{x}^{*}{b}_{23}}{b_1{x}^{*}+1}+\frac{ea{x}^{*}{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2},$

$E_{2z}=\frac{ea{y}^{*}}{b_1{x}^{*}+1}+\frac{ea{x}^{*}{b}_{23}}{b_1{x}^{*}+1}-\frac{ea{x}^{*}{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-d{b}_{23}-\frac{\beta b_{33}{y}^{*}}{b_2{p}^{*}+1}-\frac{\beta p^{*}{b}_{23}}{b_2{p}^{*}+1}+\frac{\beta p^{*}{y}^{*}{b}_2{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^2},$

$E_{3z}=-\beta p^{*}{b}_{23}-\beta b_{33}{y}^{*}-m{b}_{33}+\alpha ,$

$E_{1xz}=-\frac{2r}{K}-\frac{ea{b}_{23}}{b_1{x}^{*}+1}+\frac{2ea{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{ea{b}_{21}}{b_1{x}^{*}+1}+\frac{ea{x}^{*}{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{ea{x}^{*}{b}_{23}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{2ea{x}^{*}{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3},$

$\begin{array}{c}{E}_{2xz}=\frac{ea{b}_{23}}{b_1{x}^{*}+1}-\frac{2ea{y}^{*}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{ea{b}_{21}}{b_1{x}^{*}+1}-\frac{ea{x}^{*}{b}_{21}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{ea{x}^{*}{b}_{23}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}+\frac{2ea{x}^{*}{y}^{*}{b}_1^2}{{\left(b_1{x}^{*}+1\right)}^3}-\frac{\beta b_{31}{b}_{23}}{b_2{p}^{*}+1}\\ +\frac{2\beta b_{31}{y}^{*}{b}_2{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{\beta b_{33}{b}_{21}}{b_2{p}^{*}+1}+\frac{\beta p^{*}{b}_{21}{b}_2{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{\beta p^{*}{b}_{23}{b}_2{b}_{31}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{y}^{*}{b}_2^2{b}_{31}{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^3},\end{array}$

$E_{3xz}=-\beta b_{33}{b}_{21}-\beta b_{31}{b}_{23},$

$E_{1yz}=-\frac{ea{b}_{22}}{b_1{x}^{*}+1}+\frac{ea{x}^{*}{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2},$

$\begin{array}{c}{E}_{2yz}=\frac{ea{b}_{22}}{b_1{x}^{*}+1}-\frac{ea{x}^{*}{b}_{22}{b}_1}{{\left(b_1{x}^{*}+1\right)}^2}-\frac{\beta b_{32}{b}_{23}}{b_2{p}^{*}+1}+\frac{2\beta b_{32}{y}^{*}{b}_2{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{\beta b_{33}{b}_{22}}{b_2{p}^{*}+1}\\ +\frac{\beta p^{*}{b}_{22}{b}_2{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^2}+\frac{\beta p^{*}{b}_{23}{b}_2{b}_{32}}{{\left(b_2{p}^{*}+1\right)}^2}-\frac{2\beta p^{*}{y}^{*}{b}_2^2{b}_{32}{b}_{33}}{{\left(b_2{p}^{*}+1\right)}^3},\end{array}$

$E_{3yz}=-\beta b_{33}{b}_{22}-\beta b_{32}{b}_{23}.$

根据以上各式，我们可以得到以下数值表达式：

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

${\mu }_2=-\frac{\mathrm{Re}{C}_1\left(0\right)}{{\alpha }^{\prime }\left(0\right)},$ (18)

${\tau }_2=-\frac{\mathrm{Im}{C}_1\left(0\right)+{\mu }_2{\omega }^{\prime }\left(0\right)}{\omega },$ (19)

${\beta }_2=2\mathrm{Re}{C}_1\left(0\right)\text{.}$ (20)

其中， ${\alpha }^{\prime }\left(0\right)={\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Re}\left({\lambda }_1\left(\beta \right)\right)\right)|}_{\beta ={\beta }^{*}},{\omega }^{\prime }\left(0\right)={\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Im}\left({\lambda }_1\left(\beta \right)\right)\right)|}_{\beta ={\beta }^{*}}$ 。

于是我们得到以下结论：

1) ${\mu }_2$ 决定了Hopf分岔的方向。当 ${\mu }_2>0$ 时，Hopf分岔为超临界的，且当 $\beta >{\beta }^{*}$ 时，存在分岔周期解；当 ${\mu }_2<0$ 时，Hopf分岔为亚临界的，且当 $\beta <{\beta }^{*}$ 时，存在分岔周期解。

2) ${\beta }_2$ 决定了分岔周期解的稳定性。当 ${\beta }_2<0$ 时，分岔周期解是稳定的；当 ${\beta }_2>0$ 时，分岔周期解是不稳定的。

3) ${\tau }_2$ 决定了分岔周期解的周期，当 ${\tau }_2>0$ 时，分岔周期解的周期将增大；当 ${\tau }_2<0$ 时，分岔周期解的周期将减小。

基于数学理论分析，我们探索了模型(2)平衡点的存在性和稳定性，探究了Hopf分岔的存在性、方向性和稳定性，给出了一些关键参数制约下的阈值条件，这些理论结果为后续的数值模拟提供了一定的理论基础，也理论证明了模型(2)具有较为复杂的动力学性态，为藻鱼可以形成周期振荡生长共存模式提供了理论支撑。

5. 数值模拟

为了验证关于平衡点和Hopf分岔理论推导结果的有效性和可行性，我们使用MATLAB进行了一些数值计算。其中，我们在模型(2)中选择如下一组参数值：

$r=1,K=5,a=0.5,b_1=0.3,b_2=0.5,d=0.1,\alpha =1,e=0.5,\beta =\mathrm{0.1.}$ (21)

通过将以上参数取值代入模型(2)，我们得到了模型在生物学上有意义的唯一内平衡点 $E^{*}\left(0.2922437696,2.048199983,0.4146360351\right)$ 。

通过定理3，我们分别计算对应于该内平衡点的雅可比矩阵，矩特征根方程中的 ${\sigma }_1,{\sigma }_2,{\sigma }_3$ ，则有：

${\sigma }_1=0.6969582547>0,{\sigma }_3=0.1031576106>0,$

${\sigma }_1{\sigma }_2-{\sigma }_3=0.0270982788>0.$

显然，满足了该内平衡点的稳定性条件，因此模型(2)的平衡点 $E^{*}$ 是局部渐近稳定的，这说明产毒藻类种群和鱼类种群可以形成局部常稳态生长共存模式，详细结果如图1。

Figure 1. The phase diagram of the model (2) shows the presence of stable focal points with $\beta =0.1$

图1. 在 $\beta =0.1$ 时，模型(2)的相图展示存在稳定的焦点

如果我们逐渐增大或逐渐减小 $\beta$ 的值，模型(2)内平衡点都将失去稳定性，模型(2)将发生Hopf分岔动力学行为，出现极限f环，这说明模型(2)既可以发生亚临界Hopf分岔，也可以发生跨临界Hopf分岔。也就说当鱼类和藻毒素接触比较少时，产毒藻类种群和鱼类种群可以形成周期振荡生长共存模式，这个时候鱼类种群平均密度比较高，而产毒藻类种群密度比较低。当鱼类和藻毒素接触比较多时，产毒藻类种群和鱼类种群也可以形成周期振荡生长共存模式，但是鱼类种群平均密度偏低一些，而产毒藻类种群密度比较高，详细结果见图2。通过这个数值模拟结果可知，当鱼类和藻毒素接触比较多时，虽然不利于鱼类种群形成较高密度的生长共存模式，但是可以保存鱼类种群的持久生存。

(a) (b)

Figure 2.The phase diagram of the model (2) represents the limit cycle with (a) $\beta =0.03$ and (b) $\beta =0.5$

图2. 分别展示了在(a) $\beta =0.03$ 和(b) $\beta =0.5$ 时，模型(2)的相位图，表示了极限环

为了更加清晰地描述这一现象，下面我们给出了模型(2)控制参数为 $\beta$ 的分岔图3。通过定理4以及分岔图，我们可以发现模型在 ${\beta }_1^{*}\approx 0.05$ 以及 ${\beta }_2^{*}\approx 0.31$ 处均会发生Hopf分岔。

同时，我们也分别计算了其Hopf分岔的横截性条件：

${\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Re}\left(\lambda \left(\beta \right)\right)\right)|}_{\beta ={\beta }_1^{\text{*}}}=-0.09652214120<0$

${\frac{\text{d}}{\text{d}\beta }\left(\mathrm{Re}\left(\lambda \left(\beta \right)\right)\right)|}_{\beta ={\beta }_2^{*}}=0.1212853383>0$

显然均满足其Hopf分岔的横截性条件。

因此，当 ${\beta }_1^{*}<\beta <{\beta }_2^{*}$ 时，模型(2)在平衡态 $E^{*}$ 处渐近稳定；当 $\beta <{\beta }_1^{*}$ 或 $\beta >{\beta }_2^{*}$ 时，模型(2)在平衡态 $E^{*}$ 处失去稳定性，发生Hopf分岔动力学行为。由此我们得出，在藻毒素的环境影响下，不管鱼类种群接触藻毒素是多还是少，产毒藻类种群和鱼类种群都可以形成生长共存模式，在 $0.05<\beta <0.31$ ，产毒藻类种群和鱼类种群可以常稳态持久生存，在 $\beta <0.05$ 和在 $0.31<\beta$ ，产毒藻类种群和鱼类种群可以稳态周期振荡持久生存。因此，也就是说，不管鱼类种群接触藻毒素是多还是少，产毒藻类种群和鱼类种群都可以持久生存，仅仅是种群密度有所差异。

此外，通过推导结果(17) (18) (19) (20)，我们能够得到更多的结论，并验证我们模拟图像的准确性。我们分别计算对应于 ${\beta }_1^{*}$ 以及 ${\beta }_2^{*}$ 的以下数值：

1) 对于 ${\beta }_1^{*}=0.05$

$C_1\left(0\right)=-0.00790978273-0.2338730498\ast I,$

(a) (b) (c)

Figure 3. The bifurcation diagram of the model (2) with $\beta$ as the control parameter

图3. 模型(2)以 $\beta$ 为控制参数的分岔图

${\mu }_2=-0.03749667322<0,$

${\tau }_2=0.7142426584>0,$

${\beta }_2=-0.01581956546<0.$

2) 对于 ${\beta }_2^{*}=0.31$

$C_1\left(0\right)=-0.00324797872-0.09486681945\ast I,$

${\mu }_2=0.02677964843>0,$

${\tau }_2=0.3436759209>0,$

${\beta }_2=-0.00649595744<0.$

显然，对于 ${\beta }_1^{*}$ ，由于其对应的 ${\mu }_2<0$ ，所以在该处产生的Hopf分岔为亚临界的，且当 $\beta <{\beta }_1^{*}$ 时，存在分岔周期解。对于 ${\beta }_2^{*}$ ，由于其对应的 ${\mu }_2>0$ ，所以在该处产生的Hopf分岔为超临界的，且当 $\beta >{\beta }_2^{*}$ 时，存在分岔周期解。另外，对于 ${\beta }_1^{*}$ 以及 ${\beta }_2^{*}$ 其对应的 ${\beta }_2$ 参数均小于0，能够得到当 $\beta >{\beta }_2^{*}$ 或当 $\beta <{\beta }_1^{*}$ 时，其分岔周期解具有稳定性。同时，对于 ${\beta }_1^{*}$ 以及 ${\beta }_2^{*}$ 其对应的 ${\tau }_2$ 参数均小于0，能够得到当 $\beta >{\beta }_2^{*}$ 或当 $\beta <{\beta }_1^{*}$ 时，其分岔周期解的周期将逐渐增大。

为了探究人为控制藻毒素对该模型动力学和种群密度的影响，我们再次选取了另一参数，藻毒素的消耗率 ，并分别绘制了当 $\beta =0.5$ 和 $\beta =0.03$ 对于参数m的分岔图4、图5，其他参数与(21)保持一致。

从图4、图5中我们可以看到，当藻毒素接触率 $\beta$ 较大( $\beta =0.5$ )时，随着藻毒素消耗率m的增加，通过分岔点( $m^{*}\approx 0.10$ )后，模型(2)的动力学行为将从常稳态转变为周期振荡稳态，且随着参数取值的增加，鱼类种群密度具有缓慢增加趋势，产毒藻类种群和藻毒素都具有递减趋势。当关键参数取值超过关键值时，产毒藻类种群和鱼类种群可以形成稳态周期振荡生长共存模式，藻毒素的浓度平均减低，鱼类种群的密度平均升高。因此，当藻毒素与鱼类种群接触较多时，有效人为控制藻毒素，有利于产毒藻类种群和鱼类种群的稳态周期振荡持久生存。

(a) (b) (c)

Figure 4.The bifurcation diagram of the model (2) with m as the control parameter with $\beta =0.5$

图4. 当 $\beta =0.5$ 时，模型(2)以m为控制参数的分岔图

(a) (b) (c)

Figure 5.The bifurcation diagram of the model (2) with m as the control parameter with $\beta =0.03$

图5. 当 $\beta =0.03$ 时，模型(2)以m为控制参数的分岔图

当藻毒素接触率 $\beta$ 很小( $\beta =0.03$ )时，当 $m<m_1^{*}\approx 0.12$ 时，产毒藻类种群和鱼类种群形成周期振荡共存模式，而随着m的继续增加，产毒藻类种群和鱼类种群形成常稳态共存模式，在增加超过 $m_2^{*}\approx 0.30$ 后，产毒藻类种群和鱼类种群又形成周期振荡共存模式。同时可以发现，在常稳态共存模式下，增加人为控制藻毒素的消耗，对产毒藻类种群和鱼类种群的密度变化影响不大，也就是说在低藻毒素环境下，人为控制藻毒素对藻类种群和鱼类种群的种群密度变化影响不大，可以用自然消耗就可以了。

通过上述数值模拟结果分析可知，关键参数 $\beta$ 和m严重影响模型的动力学性态演化进程，导致模型动力学行为在稳态平衡点和稳态极限环之间转换，进而导致产毒藻类种群和鱼类种群的生长共存模式在常稳态和周期振荡稳态之间转换，但是相对来说，稳态周期振荡生长共存模式更有利于生态循环。

6. 总结

本论文把藻毒素当作一个变量，构建了一类具有藻毒素动态影响的水域生态模型，对其动力学行为进行理论分析与数值模拟研究。首先，理论推导出模型(2)平衡点的存在性和稳定性，给出模型(2)发生Hopf分岔的阈值条件，推理出Hopf分岔的方向性、稳定性和周期性。其次，对相关动力学行为进行数值模拟试验，阐明模型(2)动力学性态的演化进程，揭示关键参数取值对模型(2)动力学的影响机制。最后，从种群动力学的角度，阐明产毒藻类种群和鱼类种群可以形成常稳态生长共存模式和稳态周期振荡生长共存模式，并指出它们的内在驱动机制。总体来说，希望本研究结果有利于促进藻鱼生态系统的快速发展，可为藻毒素的有效控制与治理提供一定的理论支撑。