成熟时滞与消化时滞共同作用对一类捕食模型动力学行为的影响

The joint impacts of maturation delay and digestion delay on the dynamics of a predator⁃prey model

摘要

研究了一类具有成熟时滞和消化时滞的捕食者⁃食饵模型的动力学。首先，运用几何方法对具有参数依赖时滞的特征方程进行了分析，得到了时滞参数平面上的Hopf分支曲线和不稳定区域。其次，运用中心流形和规范型理论推导出了判断分支方向和周期轨道稳定性的准则。最后，结合数值模拟发现，随着消化时滞的增大，模型将失稳；而随着成熟时滞增加，模型经过有限次稳定性切换后趋于稳定状态。

Abstract

The dynamics of a predator⁃prey model with maturation and digestion delays were studied. The characteristic equation with delay-dependent parameters was first analyzed using the geometric method. Accordingly， the Hopf bifurcation curve and unstable region on the two-delay plane were obtained. The center manifold and normal form theory were then applied to deduce the criteria for judging the bifurcation direction and the stability of the bifurcation periodic orbit. Combined with numerical examples， it was found that， as the digestion delay increases， the model would lose its stability. In contrast， as the maturation lag increases， the model would tend to be stable after finite stability switches.

Graphical abstract

关键词

时滞 / Hopf分支 / 稳定性开关 / 规范型

Key words

delays / Hopf bifurcation / stability switches / normal form

0 引言

捕食者⁃食饵模型可用来描述生物群落中重要的捕食关系。捕食者与猎物在自然界相互作用中伴随一些复杂的生物过程，其中恐惧效应指的是猎物面临被捕食的风险而产生的反应，其比直接捕食更能影响猎物的行为和心理，可通过改变繁殖、觅食和栖息等因素，最终影响捕食者和猎物的种群动态^［1］。WANG等^［2］首次基于经典的捕食者⁃食饵模型考虑了恐惧效应，并构造如下模型：

\left\{\begin{array}{l}\frac{\mathrm{d}N(t)}{\mathrm{d}t}=\frac{rN}{\mathrm{1}+kP}-d_{\mathrm{1}}N(t)-\\ d_{\mathrm{2}}{N}^{\mathrm{2}}(t)-\frac{\alpha NP}{\mathrm{1}+hN},\\ \frac{\mathrm{d}P(t)}{\mathrm{d}t}=\beta \frac{\alpha NP}{\mathrm{1}+hN}-mP(t),\end{array}\right.(1)

式中：N（t）和P（t）分别表示食饵和捕食者在t时刻的生物量；r\mathrm{、}{d}_{\mathrm{1}}和d_{\mathrm{2}}分别为食饵的出生率、自然死亡率和种内竞争诱导的死亡率；k刻画恐惧效应的强度；α为捕食者的攻击率；h为半饱和常数；β为食饵的转化率；m为捕食者的死亡率。研究表明，恐惧效应可增加模型（1）的稳定性。

许多生物过程都不是瞬时完成的，都存在时间延迟^［3］。例如，猎物从出生到能够繁衍后代需要经历成熟期。再如，捕食者需要时间去消化食物，然后才能进行下一步的活动和反应。KUANG在其专著^［4］中指出：可以用时滞微分方程来表征成熟期和消化时间。MA等^［5］通过沿特征线求解具有年龄结构模型推导出下面可以描述食饵成熟期的模型：

\left\{\begin{array}{l}\frac{\mathrm{d}N(t)}{\mathrm{d}t}={\mathrm{e}}^{-d_J\tau }\frac{rN(t-\tau )}{\mathrm{1}+kP(t-\tau )}-d_{\mathrm{1}}N(t)-\\ d_{\mathrm{2}}{N}^{\mathrm{2}}(t)-\frac{\alpha N(t)P(t)}{\mathrm{1}+hN(t)},\\ \frac{\mathrm{d}P(t)}{\mathrm{d}t}=\beta \frac{\alpha N(t)P(t)}{\mathrm{1}+hN(t)}-mP(t),\end{array}\right.(2)

式中：τ和d_J分别表示食饵的成熟期和幼年食饵的死亡率。

在模型（2）的基础上，进一步考虑捕食者的消化时间，得到下面双时滞模型：

\left\{\begin{array}{l}\frac{\mathrm{d}N(t)}{\mathrm{d}t}={\mathrm{e}}^{-d_J\tau }\frac{rN(t-\tau )}{\mathrm{1}+kP(t-\tau )}-d_{\mathrm{1}}N(t)-\\ d_{\mathrm{2}}{N}^{\mathrm{2}}(t)-\frac{\alpha N(t)P(t)}{\mathrm{1}+hN(t)},\\ \frac{\mathrm{d}P(t)}{\mathrm{d}t}=\beta \frac{\alpha N(t-\sigma )P(t-\sigma )}{\mathrm{1}+hN(t-\sigma )}-mP(t),\end{array}\right.(3)

式中：σ为消化时滞。

模型（3）涉及两个多时滞且参数依赖时滞，研究其动力学行为具有一定的难度。本文将运用可使两时滞同时变化的几何方法，获得时滞平面上的Hopf分支曲线，进而计算Hopf分支的规范型来判定模型失稳后周期轨道的稳定性。研究结果可揭示成熟期和消化时滞共同作用对捕食⁃被捕食者动力学的影响。

1 模型（3）的线性分析

下面将采用几何方法分析模型（3）在正平衡点处的特征方程，从而确定时滞参数平面上线性稳定区域与Hopf分支产生的条件。

从可持续发展角度，重点关注捕食者和猎物的共存，因此主要考虑模型（3）的正平衡点。

若r{\mathrm{e}}^{-d_J\tau }>d_{\mathrm{1}}和

(r{\mathrm{e}}^{-d_J\tau }-d_{\mathrm{1}})(\alpha \beta -mh)>d_{\mathrm{2}}m，

则模型（3）存在唯一正稳态解E^* =（N_*，P_* （τ）），其中

N_* =\frac{m}{\alpha \beta -mh}，

P_* （τ）满足下面方程：

\begin{array}{l}\alpha k{P}_{*}^{\mathrm{2}}(\tau )+(k(\mathrm{1}+h{N}_{*})(d_{\mathrm{1}}+\\ d_{\mathrm{2}}{N}_{*})+\alpha )P_{*}(\tau )+\\ (\mathrm{1}+h{N}_{*})(d_{\mathrm{2}}{N}_{*}-r{\mathrm{e}}^{-d_J\tau }+d_{\mathrm{1}})=\mathrm{0}\mathrm{。}\end{array}

模型（3）在E^* 处的线性化系统为：

\left\{\begin{array}{l}\frac{\mathrm{d}N(t)}{\mathrm{d}t}=a_{\mathrm{11}}(\tau )N(t)+a_{\mathrm{12}}P(t)+\\ b_{\mathrm{11}}(\tau )N(t-\tau )+b_{\mathrm{12}}(\tau )P(t-\tau ),\\ \frac{\mathrm{d}P(t)}{\mathrm{d}t}=a_{\mathrm{22}}P(t)+c_{\mathrm{21}}(\tau )N(t-\sigma )+\\ c_{\mathrm{22}}P(t-\sigma ),\end{array}\right.(4)

式中：

a_{\mathrm{11}}(\tau )=-(d_{\mathrm{1}}+\mathrm{2}{d}_{\mathrm{2}}{N}_{*}+\frac{\alpha P_{*}(\tau )}{(\mathrm{1}+h{N}_{*}{)}^{\mathrm{2}}}),

a_{\mathrm{12}}=-\frac{\alpha N_{*}}{\mathrm{1}+h{N}_{*}},a_{\mathrm{22}}=-m,

b_{\mathrm{11}}(\tau )=\frac{r{\mathrm{e}}^{-d_J\tau }}{\mathrm{1}+k{P}_{*}(\tau )},

b_{\mathrm{12}}(\tau )=-\frac{rk{N}_{*}{\mathrm{e}}^{-d_J\tau }}{(\mathrm{1}+k{P}_{*}{(\tau ))}^{\mathrm{2}}},

c_{\mathrm{21}}(\tau )=\frac{\alpha \beta P_{*}(\tau )}{(\mathrm{1}+h{N}_{*}{)}^{\mathrm{2}}},c_{\mathrm{22}}=\frac{\alpha \beta N_{*}}{\mathrm{1}+h{N}_{*}}\mathrm{。}

将N（t）=ξ₁e ^λt，P（t）=ξ₂e ^λt 代入系统（4）得

\left\{\begin{array}{l}\lambda {\xi }_{\mathrm{1}}=\left(a_{\mathrm{11}}(\tau )+b_{\mathrm{11}}(\tau ){\mathrm{e}}^{-\lambda \tau }\right){\xi }_{\mathrm{1}}+\\ \left(a_{\mathrm{12}}+b_{\mathrm{12}}(\tau ){\mathrm{e}}^{-\lambda \tau }\right){\xi }_{\mathrm{2}},\\ \lambda {\xi }_{\mathrm{2}}=c_{\mathrm{21}}(\tau ){\mathrm{e}}^{-\lambda \sigma }{\xi }_{\mathrm{1}}+\left(a_{\mathrm{22}}+c_{\mathrm{22}}{\mathrm{e}}^{-\lambda \sigma }\right){\xi }_{\mathrm{2}}\mathrm{。}\end{array}\right.(5)

若

\left|\begin{array}{ll}\lambda -a_{\mathrm{11}}(\tau )-b_{\mathrm{11}}(\tau ){\mathrm{e}}^{-\lambda \tau }& -a_{\mathrm{12}}-b_{\mathrm{12}}(\tau ){\mathrm{e}}^{-\lambda \tau }\\ -c_{\mathrm{21}}(\tau ){\mathrm{e}}^{-\lambda \sigma }& \lambda -a_{\mathrm{22}}-c_{\mathrm{22}}{\mathrm{e}}^{-\lambda \sigma }\end{array}\right|=\mathrm{0},(6)

则方程（5）存在非平凡解（ξ₁，ξ₂）。

通过整理方程（6），得到系统（4）的特征方程：

\begin{array}{l}P(\lambda ;\tau ,\sigma )=p_{\mathrm{0}}(\lambda ,\tau )+p_{\mathrm{1}}(\lambda ,\tau ){\mathrm{e}}^{-\mathrm{\lambda }\mathrm{\tau }}+\\ p_{\mathrm{2}}(\lambda ,\tau ){\mathrm{e}}^{-\mathrm{\lambda }\mathrm{\sigma }}+p_{\mathrm{3}}(\lambda ,\tau ){\mathrm{e}}^{-\mathrm{\lambda }(\mathrm{\tau }+\mathrm{\sigma })}=\mathrm{0},\end{array}(7)

式中：

p_{\mathrm{0}}(\lambda ,\tau )={\lambda }^{\mathrm{2}}-\left(a_{\mathrm{11}}(\tau )+a_{\mathrm{22}}\right)\lambda +a_{\mathrm{11}}(\tau )a_{\mathrm{22}},

p_{\mathrm{1}}(\lambda ,\tau )=-b_{\mathrm{11}}(\tau )\lambda +a_{\mathrm{22}}{b}_{\mathrm{11}}(\tau ),

p_{\mathrm{2}}(\lambda ,\tau )=-c_{\mathrm{22}}\lambda +a_{\mathrm{11}}(\tau )c_{\mathrm{22}}-a_{\mathrm{12}}{c}_{\mathrm{21}}(\tau ),

p_{\mathrm{3}}(\lambda ,\tau )=b_{\mathrm{11}}(\tau )c_{\mathrm{22}}-b_{\mathrm{12}}(\tau )c_{\mathrm{21}}(\tau )\mathrm{。}

将采用文献［6］中的方法处理特征方程（7），进而确定（τ，σ）平面上的穿越曲线以及稳定与不稳定区域。把\lambda =\mathrm{i}\omega \left(\omega >\mathrm{0}\right)代入方程（7）有

\begin{array}{l}\left(p_{\mathrm{2}}\left(\mathrm{i}\right.\omega ,\tau )+p_{\mathrm{3}}\left(\mathrm{i}\right.\omega ,\tau ){\mathrm{e}}^{-\mathrm{i}\omega \tau }\right){\mathrm{e}}^{-\mathrm{i}\mathrm{\omega }\mathrm{\sigma }}+\\ p_{\mathrm{0}}(\mathrm{i}\omega ,\tau )+p_{\mathrm{1}}(\mathrm{i}\omega ,\tau ){\mathrm{e}}^{-\mathrm{i}\mathrm{\omega }\mathrm{\tau }}=\mathrm{0}\mathrm{。}\end{array}(8)

对方程（8）两端取模并注意到

|{\mathrm{e}}^{-\mathrm{i}\mathrm{\omega }\mathrm{\sigma }}|=|\mathrm{c}\mathrm{o}\mathrm{s}(\omega \sigma )+\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}(\omega \sigma )|=\mathrm{1},

可得

\begin{array}{l}|p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}+|p_{\mathrm{1}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-\\ |p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-|p_{\mathrm{3}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}=\\ \mathrm{2}\left(E(\omega ,\tau )\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )-F(\omega ,\tau )\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )\right),\end{array}(9)

式中：

\begin{array}{l}E(\omega ,\tau )=\mathrm{R}\mathrm{e}(p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{3}}(\mathrm{i}\omega ,\tau ))-\\ \mathrm{R}\mathrm{e}(p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{1}}(\mathrm{i}\omega ,\tau )),\end{array}

\begin{array}{l}F(\omega ,\tau )=\mathrm{I}\mathrm{m}(p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{3}}(\mathrm{i}\omega ,\tau ))-\\ \mathrm{I}\mathrm{m}(p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{1}}(\mathrm{i}\omega ,\tau ))\mathrm{。}\end{array}

令

\begin{array}{l}\phi (\omega ,\tau )=\mathrm{a}\mathrm{r}\mathrm{g}(p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{3}}(\mathrm{i}\omega ,\tau )-\\ p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){\overline{p}}_{\mathrm{1}}(\mathrm{i}\omega ,\tau )),\end{array}

那么方程（9）可写成

\begin{array}{l}|p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}+|p_{\mathrm{1}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-\\ |p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-|p_{\mathrm{3}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}=\\ \mathrm{2}\sqrt[ ]{E^{\mathrm{2}}(\omega ,\tau )+F^{\mathrm{2}}(\omega ,\tau )}\mathrm{c}\mathrm{o}\mathrm{s}(\phi (\omega ,\tau )+\omega \tau )\mathrm{。}\end{array}（10）

由式（10）可知，特征方程（7）存在纯虚根的必要条件是：（τ，ω）满足下面不等式

\begin{array}{l}\left(|p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}+|p_{\mathrm{1}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-\right.\\ |p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-{\left.|p_{\mathrm{3}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}\right)}^{\mathrm{2}}\le \\ \mathrm{4}{\left(E^{\mathrm{2}}(\omega ,\tau )+F^{\mathrm{2}}(\omega ,\tau )\right)}^{\mathrm{2}},\end{array}

称\left(\tau ,\omega \right)满足上述不等式的区域D为可行区域。

此外，记

\begin{array}{l}\phi (\omega ,\tau )=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}((|p_{\mathrm{0}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}+\\ |p_{\mathrm{1}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-|p_{\mathrm{2}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}}-|p_{\mathrm{3}}(\mathrm{i}\omega ,\tau ){|}^{\mathrm{2}})/\\ \mathrm{2}\sqrt[ ]{E^{\mathrm{2}}(\omega ,\tau )+F^{\mathrm{2}}(\omega ,\tau )})\mathrm{。}\end{array}

由此可定义函数f_{m_{\mathrm{1}}}^{\pm }:D\to \mathbf{R}为

\begin{array}{l}{f}_{m_{\mathrm{1}}}^{\pm }\left(\tau ,\omega \right)=\phi (\omega ,\tau )\mp \phi (\omega ,\tau )-\\ \mathrm{2}{m}_{\mathrm{1}}\mathrm{\pi }+\omega \tau ,m_{\mathrm{1}}\in \mathbf{Z}\mathrm{。}\end{array}

若\left({\tau }_{m_{\mathrm{1}}}^{\pm },\omega \right)满足f_{m_{\mathrm{1}}}^{\pm }\left(\tau ,\omega \right)=\mathrm{0},代入式（8）可得到σ的表达式

\begin{array}{l}{\sigma }_{m_{\mathrm{2}}}^{\pm }=\frac{\mathrm{1}}{\omega }\mathrm{a}\mathrm{r}\mathrm{g}(-\frac{p_{\mathrm{2}}\left(\mathrm{i}\omega ,{\tau }_{n_{\mathrm{1}}}^{\pm }\right)+p_{\mathrm{3}}\left(\mathrm{i}\omega ,{\tau }_{n_{\mathrm{1}}}^{\pm }\right){\mathrm{e}}^{-\mathrm{i}\omega {\tau }_{m_{\mathrm{1}}}^{\pm }}}{p_{\mathrm{0}}\left(\mathrm{i}\omega ,{\tau }_{m_{\mathrm{1}}}^{\pm }\right)+p_{\mathrm{1}}\left(\mathrm{i}\omega ,{\tau }_{m_{\mathrm{1}}}^{\pm }\right){\mathrm{e}}^{-\mathrm{i}\omega {\tau }_{m_{\mathrm{1}}}^{\pm }}})+\\ \mathrm{2}{m}_{\mathrm{2}}\mathrm{\pi },m_{\mathrm{2}}\in \mathbf{Z}\mathrm{。}\end{array}

于是得到时滞平面上的穿越曲线：

S=\left({\tau }_{m_{\mathrm{1}}}^{\pm },{\sigma }_{m_{\mathrm{2}}}^{\pm }\right)\in {\mathbf{R}}_{+}^{\mathrm{2}},m_{\mathrm{1}},m_{\mathrm{2}}\in \mathbf{Z}\mathrm{。}(11)

若({\tau }_{*},{\sigma }_{*})\in S且相应横截条件成立时，则存在{\omega }_{*}满足P(\mathrm{i}{\omega }_{*};{\tau }_{*},{\sigma }_{*})=\mathrm{0}以及在({\tau }_{*},{\sigma }_{*})附近存在{\lambda }^{\pm }(\tau ,\sigma )=\alpha (\tau ,\sigma )\pm \mathrm{i}\omega (\tau ,\sigma )满足\alpha ({\tau }_{*},{\sigma }_{*})=\mathrm{0},\omega ({\tau }_{*},{\sigma }_{*})={\omega }_{*}。

下面将给出时滞(\tau ,\sigma )偏离({\tau }_{*},{\sigma }_{*})时的特征值穿越方向。由文献［7］的定义知：沿穿越曲线ω增加的方向为正向；沿穿越曲线正向走，穿越曲线的左侧（右侧）区域为左侧（右侧）区域。记

R_{\mathrm{1}}(\tau ,\sigma )=\mathrm{R}\mathrm{e}(\frac{\partial P\left(\lambda ;\tau ,\sigma \right)}{\partial \tau }),

I_{\mathrm{1}}(\tau ,\sigma )=\mathrm{I}\mathrm{m}(\frac{\partial P\left(\lambda ;\tau ,\sigma \right)}{\partial \tau }),

R_{\mathrm{2}}(\tau ,\sigma )=\mathrm{R}\mathrm{e}(\frac{\partial P\left(\lambda ;\tau ,\sigma \right)}{\partial \sigma }),

I_{\mathrm{2}}(\tau ,\sigma )=\mathrm{I}\mathrm{m}(\frac{\partial P\left(\lambda ;\tau ,\sigma \right)}{\partial \sigma })\mathrm{。}

由隐函数定理和文献［7］可知下面定理成立：

定理1　当({\tau }_{*},{\sigma }_{*})\in S时，特征方程（7）存在纯虚根。

此外，若

则当（τ，σ）穿过({\tau }_{*},{\sigma }_{*})到穿越曲线右侧（左侧）时，特征方程（7）有两对具有正实部的特征值。

2 Hopf分支的规范型

下面运用中心流形定理和规范型理论推导模型（3）在Hopf分支点附近的动力学行为。

由定理（1）知，固定τ，存在Hopf分支值σ_* 使得特征方程（7）存在纯虚根\mathrm{i}{\omega }_{*}且其余特征根均具有负实部。采用文献［8］的方法，计算模型（3）在σ_* 处Hopf分支的规范型，并确定分支的方向和周期轨的稳定性。

记\mathit{U}(t)=(U_{\mathrm{1}}(t),U_{\mathrm{2}}{(t))}^{\mathrm{T}}，U_{\mathrm{1}}(t)=N(t\sigma )-N_{*},U_{\mathrm{2}}(t)=P(t\sigma )-P_{*}(\tau ),则模型（3）转化为

\frac{\mathrm{d}\mathit{U}(t)}{\mathrm{d}t}=L\left(\sigma \right)\left({\mathit{U}}_t\right)+g\left({\mathit{U}}_t,\sigma \right),(12)

式中：

L(\sigma )(\mathit{\psi })=\sigma {\mathit{L}}_{\mathrm{1}}\mathit{\psi }(\mathrm{0})+\sigma {\mathit{L}}_{\mathrm{2}}\mathit{\psi }(\frac{-\tau }{\sigma })+\sigma {\mathit{L}}_{\mathrm{3}}\mathit{\psi }(-\mathrm{1}),

g(\mathit{\psi },\sigma )=\sigma (g_{\mathrm{1}},g_{\mathrm{2}}{)}^{\mathrm{T}},\mathit{\psi }=({\psi }_{\mathrm{1}},{\psi }_{\mathrm{2}}{)}^{\mathrm{T}},(13)

式中：

\begin{array}{l}{g}_{\mathrm{1}}=\sum _{i+j+k+l\ge \mathrm{2}}\frac{\mathrm{1}}{i!j!k!l!}{g}_{\mathrm{1}}^{ijkl}{\psi }_{\mathrm{1}}^i\left(\mathrm{0}\right){\psi }_{\mathrm{2}}^j\left(\mathrm{0}\right)·\\ {\psi }_{\mathrm{1}}^k(\frac{-\tau }{\sigma }){\psi }_{\mathrm{2}}^l(\frac{-\tau }{\sigma }),\end{array}

g_{\mathrm{2}}=\sum _{m+n\ge \mathrm{2}}\frac{\mathrm{1}}{m!n!}{g}_{\mathrm{2}}^{mn}{\psi }_{\mathrm{1}}^m(-\mathrm{1})·{\psi }_{\mathrm{2}}^n(-\mathrm{1}),

{\mathit{L}}_{\mathrm{1}}=\left(\begin{array}{ll}{a}_{\mathrm{11}}& a_{\mathrm{12}}\\ \mathrm{0}& a_{\mathrm{22}}\end{array}\right),{\mathit{L}}_{\mathrm{2}}=\left(\begin{array}{ll}{b}_{\mathrm{11}}& b_{\mathrm{11}}\\ \mathrm{0}& \mathrm{0}\end{array}\right),{\mathit{L}}_{\mathrm{3}}=\left(\begin{array}{ll}\mathrm{0}& \mathrm{0}\\ c_{\mathrm{21}}& c_{\mathrm{22}}\end{array}\right)\mathrm{。}

简单计算可得g_{\mathrm{1}}^{ijkl}(\mathrm{2}\le i+j+k+l\le \mathrm{3})和g_{\mathrm{2}}^{mn}(\mathrm{2}\le m+n\le \mathrm{3})如下：

g_{\mathrm{1}}^{\mathrm{2000}}=-\mathrm{2}\frac{h(d_{\mathrm{2}}{N}_{*}^{\mathrm{3}}{h}^{\mathrm{2}}+\mathrm{3}{d}_{\mathrm{2}}{N}_{*}^{\mathrm{2}}h+\mathrm{3}{d}_{\mathrm{2}}{N}_{*}-P_{*}\alpha )+d_{\mathrm{2}}}{(\mathrm{1}+N_{*}{h)}^{\mathrm{3}}},

g_{\mathrm{1}}^{\mathrm{1100}}=-\frac{\alpha }{(\mathrm{1}+N_{*}{h)}^{\mathrm{2}}},g_{\mathrm{1}}^{\mathrm{0011}}=-\frac{r{\mathrm{e}}^{-d_J\tau }k}{(\mathrm{1}+P_{*}{k)}^{\mathrm{2}}},

g_{\mathrm{1}}^{\mathrm{0002}}=\mathrm{2}\frac{{\mathrm{e}}^{-d_J\tau }r{N}_{*}{k}^{\mathrm{2}}}{(\mathrm{1}+P_{*}{k)}^{\mathrm{3}}},g_{\mathrm{1}}^{\mathrm{3000}}=-\mathrm{6}\frac{\alpha P_{*}{h}^{\mathrm{2}}}{(\mathrm{1}+N_{*}{h)}^{\mathrm{4}}},

g_{\mathrm{1}}^{\mathrm{2100}}=\mathrm{2}\frac{\alpha h}{(\mathrm{1}+N_{*}{h)}^{\mathrm{3}}},g_{\mathrm{1}}^{\mathrm{0012}}=\mathrm{2}\frac{{\mathrm{e}}^{-d_J\tau }r{k}^{\mathrm{2}}}{(\mathrm{1}+P_{*}{k)}^{\mathrm{3}}},

g_{\mathrm{1}}^{\mathrm{0003}}=-\mathrm{6}\frac{{\mathrm{e}}^{-d_J\tau }r{N}_{*}{k}^{\mathrm{3}}}{(\mathrm{1}+P_{*}{k)}^{\mathrm{4}}},

g_{\mathrm{2}}^{\mathrm{20}}=-\mathrm{2}\frac{\beta \alpha P_{*}h}{(\mathrm{1}+N_{*}{h)}^{\mathrm{3}}},g_{\mathrm{2}}^{\mathrm{11}}=\frac{\beta \alpha }{(\mathrm{1}+N_{*}{h)}^{\mathrm{2}}},

g_{\mathrm{2}}^{\mathrm{30}}=\mathrm{6}\frac{\beta \alpha P_{*}{h}^{\mathrm{2}}}{(\mathrm{1}+N_{*}{h)}^{\mathrm{4}}},g_{\mathrm{2}}^{\mathrm{21}}=-\mathrm{2}\frac{\beta \alpha h}{(\mathrm{1}+N_{*}{h)}^{\mathrm{3}}}\mathrm{。}

令\sigma ={\sigma }_{*}+\epsilon，那么系统（12）可写成

\frac{\mathrm{d}\mathit{U}(t)}{\mathrm{d}t}=L\left({\sigma }_{*}\right)\left({\mathit{U}}_t\right)+G\left({\mathit{U}}_t,\epsilon \right),(14)

式中：

\begin{array}{l}L({\sigma }_{*})(\mathit{\psi })={\sigma }_{*}{L}_{\mathrm{1}}\mathit{\psi }(\mathrm{0})+\\ {\sigma }_{*}{L}_{\mathrm{2}}\mathit{\psi }(\frac{-\tau }{{\sigma }_{*}})+{\sigma }_{*}{L}_{\mathrm{3}}\mathit{\psi }(-\mathrm{1}),\end{array}

\begin{array}{l}G(\mathit{\psi },\epsilon )=(L({\sigma }_{*}+\epsilon )-L({\sigma }_{*}))(\mathit{\psi })+\\ g(\mathit{\psi },{\sigma }_{*}+\epsilon ),\end{array}

易知\epsilon =\mathrm{0}是系统（14）的Hopf分支点，\pm \mathrm{i}{\omega }_{*}{\sigma }_{*}是下面线性系统的特征值：

\frac{\mathrm{d}\mathit{U}(t)}{\mathrm{d}t}=L\left({\sigma }_{*}\right)\left({\mathit{U}}_t\right)\mathrm{。}

由Riesz表示定理知，存在有界变差函数V(\theta )(-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}\le \theta \le \mathrm{0}),使得对\psi \in C([-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0}],{\mathbf{R}}^{\mathrm{2}})有

L\left({\sigma }_{*}\right)(\mathit{\psi })={\int }_{-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}}^{\mathrm{0}}\mathrm{d}V(\theta )\mathit{\psi }(\theta )\mathrm{。}(15)

事实上，式（15）中的V(\theta )可取为

\begin{array}{l}V(\theta )={\sigma }_{*}(L_{\mathrm{1}}\delta (\theta )+L_{\mathrm{2}}\delta (\theta +\frac{\tau }{{\sigma }_{*}})+\\ L_{\mathrm{2}}\delta \left(\theta +\mathrm{1}\right)),\end{array}

对\psi (\theta )\in C([-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0}],{\mathbf{R}}^{\mathrm{2}})，定义

A\left(\mathit{\psi }(\theta )\right)=\left\{\begin{array}{l}\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\theta },\theta \in [-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0}),\\ {\int }_{-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}}^{\mathrm{0}}\mathrm{d}V(\theta )\mathit{\psi }(\theta ),\theta =\mathrm{0},\end{array}\right.

R(\epsilon )(\mathit{\psi }(\theta ))=\left\{\begin{array}{l}\mathrm{0},\theta \in [-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0}),\\ G\left(\mathit{\psi }(\theta ),\epsilon \right),\theta =\mathrm{0}\mathrm{。}\end{array}\right.

那么系统（14）可写成

\frac{\mathrm{d}{\mathit{U}}_t}{\mathrm{d}t}=\mathit{A}{\mathit{U}}_t+R(\epsilon )\left({\mathit{U}}_t\right),(16)

对\mathit{\rho }(s)\in C(\left[\mathrm{0}\right.,\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}],({\mathbf{R}}^{\mathrm{2}}{)}^{*})，定义

{\mathit{A}}^{*}\left(\mathit{\rho }(s)\right)=\left\{\begin{array}{l}-\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}s},s\in (\mathrm{0},\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}],\\ {\int }_{-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}}^{\mathrm{0}}\mathit{\rho }\left(-s\right)\mathrm{d}V(s),s=\mathrm{0}\mathrm{。}\end{array}\right.

此外定义下面的双线性形式：

\begin{array}{l}〈\mathit{\rho }(s),\mathit{\psi }(\theta )〉=\overline{\mathit{\rho }}\left(\mathrm{0}\right)·\mathit{\psi }(\mathrm{0})-\\ {\int }_{-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\}}^{\mathrm{0}}{\int }_{\mathrm{0}}^{\theta }\overline{\mathit{\rho }}\left(\xi -\theta \right)\mathrm{d}V(\theta )\mathit{\psi }\left(\xi \right)\mathrm{d}\xi ,\end{array}

则\mathit{A}与{\mathit{A}}^{*}满足下面的伴随关系：

〈\mathit{\rho }(s),\mathit{A}\mathit{\psi }(\theta )〉=〈{\mathit{A}}^{*}\mathit{\rho }(s),\mathit{\psi }(\theta )〉\mathrm{。}

由前可知，\pm \mathrm{i}{\omega }_{*}{\sigma }_{*}是\mathit{A}和{\mathit{A}}^{*}的特征值。取

\psi (\theta )=\left(\begin{array}{l}\mathrm{1}\\ {\psi }_{\mathrm{2}}\end{array}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{*}{\sigma }_{*}\theta },

式中：

{\psi }_{\mathrm{2}}=\frac{c_{\mathrm{21}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}}{\mathrm{i}{\omega }_{*}-a_{\mathrm{22}}-c_{\mathrm{22}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}},

\mathit{\rho }(s)=D\left(\mathrm{1},{\rho }_{\mathrm{2}}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{*}{\sigma }_{*}s},

{\rho }_{\mathrm{2}}=-\frac{a_{\mathrm{11}}+\mathrm{i}{\omega }_{*}+b_{\mathrm{11}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}\tau }}{c_{\mathrm{21}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}},

\overline{D}=(\mathrm{1}+{\overline{\rho }}_{\mathrm{2}}{\psi }_{\mathrm{2}}+\left(b_{\mathrm{11}}+b_{\mathrm{12}}{\psi }_{\mathrm{2}}\right)\tau {\mathrm{e}}^{-\mathrm{i}{\omega }_{*}\tau }+

\left({\overline{\rho }}_{\mathrm{2}}{c}_{\mathrm{21}}+{\overline{\rho }}_{\mathrm{2}}{c}_{\mathrm{22}}{\psi }_{\mathrm{2}}\right){\sigma }_{*}{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}{)}^{-\mathrm{1}},

则

\mathit{A}(\psi )=\mathrm{i}{\omega }_{*}{\sigma }_{*}\mathit{\psi },{\mathit{A}}^{*}(\mathit{\rho })=-\mathrm{i}{\omega }_{*}{\sigma }_{*}\mathit{\rho },

〈\mathit{\rho }(s),\mathit{\psi }(\theta )〉=\mathrm{1},〈\mathit{\rho }(s),\overline{\mathit{\psi }}(\theta )〉=\mathrm{0}\mathrm{。}

令{\mathit{U}}_t是系统（16）当ε=0时的解。定义

z(t)=〈\mathit{\rho },{\mathit{U}}_t〉,

\mathit{W}(t,\theta )={\mathit{U}}_t-\mathrm{2}\mathrm{R}\mathrm{e}(z(t)\mathit{\psi }(\theta ))\mathrm{。}(17)

在中心流形C₀上有

W(t,\theta )=w_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+w_{\mathrm{11}}(\theta )z\overline{z}+w_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,

\begin{array}{l}{w}_{\mathrm{20}}(\theta )=\left(\begin{array}{l}{w}_{\mathrm{1}}^{\mathrm{20}}(\theta )\\ w_{\mathrm{2}}^{\mathrm{20}}(\theta )\end{array}\right),w_{\mathrm{11}}(\theta )=\left(\begin{array}{l}{w}_{\mathrm{1}}^{\mathrm{11}}(\theta )\\ w_{\mathrm{2}}^{\mathrm{11}}(\theta )\end{array}\right),\\ w_{\mathrm{02}}(\theta )=\left(\begin{array}{l}{w}_{\mathrm{1}}^{\mathrm{02}}(\theta )\\ w_{\mathrm{2}}^{\mathrm{02}}(\theta )\end{array}\right)\mathrm{。}\end{array}

进而可得

\begin{array}{l}{U}_{\mathrm{1}\mathrm{t}}(\mathrm{0})=z(t)+\overline{z}(t)+w_{\mathrm{1}}^{\mathrm{20}}(\mathrm{0})\frac{z^{\mathrm{2}}}{\mathrm{2}}+\\ w_{\mathrm{1}}^{\mathrm{11}}(\mathrm{0})z\overline{z}+w_{\mathrm{1}}^{\mathrm{02}}(\mathrm{0})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}

\begin{array}{l}{U}_{\mathrm{2}\mathrm{t}}(\mathrm{0})={\psi }_{\mathrm{2}}z(t)+{\overline{\psi }}_{\mathrm{2}}\overline{z}(t)+w_{\mathrm{2}}^{\mathrm{20}}(\mathrm{0})\frac{z^{\mathrm{2}}}{\mathrm{2}}+\\ w_{\mathrm{2}}^{\mathrm{11}}(\mathrm{0})z\overline{z}+w_{\mathrm{2}}^{\mathrm{02}}(\mathrm{0})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}

\begin{array}{l}{U}_{\mathrm{1}\mathrm{t}}(-\mathrm{1})=z(t){\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}+\overline{z}(t){\mathrm{e}}^{\mathrm{i}{\omega }_{*}{\sigma }_{*}}+w_{\mathrm{1}}^{\mathrm{20}}(-\mathrm{1})\frac{z^{\mathrm{2}}}{\mathrm{2}}+\\ w_{\mathrm{1}}^{\mathrm{11}}(-\mathrm{1})z\overline{z}+w_{\mathrm{1}}^{\mathrm{02}}(-\mathrm{1})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}

\begin{array}{l}{U}_{\mathrm{2}\mathrm{t}}(-\mathrm{1})={\psi }_{\mathrm{2}}z(t){\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}+{\overline{\psi }}_{\mathrm{2}}\overline{z}(t){\mathrm{e}}^{\mathrm{i}{\omega }_{*}{\sigma }_{*}}+\\ w_{\mathrm{2}}^{\mathrm{20}}(-\mathrm{1})\frac{z^{\mathrm{2}}}{\mathrm{2}}+w_{\mathrm{2}}^{\mathrm{11}}(-\mathrm{1})z\overline{z}+w_{\mathrm{2}}^{\mathrm{02}}(-\mathrm{1})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}

\begin{array}{l}{U}_{\mathrm{1}\mathrm{t}}(-\frac{\tau }{{\sigma }_{*}})=z(t){\mathrm{e}}^{-\mathrm{i}{\omega }_{*}\tau }+\overline{z}(t){\mathrm{e}}^{\mathrm{i}{\omega }_{*}\tau }+w_{\mathrm{1}}^{\mathrm{20}}(-\frac{\tau }{{\sigma }_{*}})\frac{z^{\mathrm{2}}}{\mathrm{2}}+\\ w_{\mathrm{1}}^{\mathrm{11}}(-\frac{\tau }{{\sigma }_{*}})z\overline{z}+w_{\mathrm{1}}^{\mathrm{02}}(-\frac{\tau }{{\sigma }_{*}})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}

\begin{array}{l}{U}_{\mathrm{2}\mathrm{t}}(-\frac{\tau }{{\sigma }_{*}})={\psi }_{\mathrm{2}}z(t){\mathrm{e}}^{-\mathrm{i}{\omega }_{*}\tau }+{\overline{\psi }}_{\mathrm{2}}\overline{z}(t){\mathrm{e}}^{\mathrm{i}{\omega }_{*}\tau }+\\ w_{\mathrm{2}}^{\mathrm{20}}(-\frac{\tau }{{\sigma }_{*}})\frac{z^{\mathrm{2}}}{\mathrm{2}}+w_{\mathrm{2}}^{\mathrm{11}}(-\frac{\tau }{{\sigma }_{*}})z\overline{z}+\\ w_{\mathrm{2}}^{\mathrm{02}}(-\frac{\tau }{{\sigma }_{*}})\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots \mathrm{。}\end{array}

由式（17）可知，当\epsilon =\mathrm{0}时，

\frac{\mathrm{d}z(t)}{\mathrm{d}t}=\mathrm{i}{\omega }_{*}{\sigma }_{*}z+g\left(z,\overline{z}\right),(18)

式中：

\begin{array}{l}g(z,\overline{z})=\overline{\rho }(\mathrm{0})G(W(t,\theta )+\\ \mathrm{2}\mathrm{R}\mathrm{e}(z(t)\psi (\theta )),\mathrm{0})=\end{array}

g_{\mathrm{20}}\frac{z^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{11}}z\overline{z}+g_{\mathrm{02}}\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{21}}\frac{z^{\mathrm{2}}\overline{z}}{\mathrm{2}}+\dots \mathrm{。}

通过计算可得，

\begin{array}{l}G(\mathit{W}(t,\theta )+\mathrm{2}\mathrm{R}\mathrm{e}(z(t)\mathit{\psi }(\theta )),\mathrm{0})=\\ {\sigma }_{*}({\mathit{G}}^{zz}\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\mathit{G}}^{z\overline{z}}z\overline{z}+{\mathit{G}}^{\overline{z}\overline{z}}\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+{\mathit{G}}^{zz\overline{z}}\frac{z^{\mathrm{2}}\overline{z}}{\mathrm{2}}),\end{array}

式中：

{\mathit{G}}^{zz}=\left(\begin{array}{l}{G}_{\mathrm{1}}^{zz}\\ G_{\mathrm{2}}^{zz}\end{array}\right),{\mathit{G}}^{z\overline{z}}=\left(\begin{array}{l}{G}_{\mathrm{1}}^{z\overline{z}}\\ G_{\mathrm{2}}^{z\overline{z}}\end{array}\right),{\mathit{G}}^{\overline{z}\overline{z}}=\left(\begin{array}{l}{G}_{\mathrm{1}}^{\overline{z}\overline{z}}\\ G_{\mathrm{2}}^{\overline{z}\overline{z}}\end{array}\right),{\mathit{G}}^{zz\overline{z}}=\left(\begin{array}{l}{G}_{\mathrm{1}}^{zz\overline{z}}\\ G_{\mathrm{2}}^{zz\overline{z}}\end{array}\right),

\begin{array}{l}{G}_{\mathrm{1}}^{zz}={\psi }_{\mathrm{2}}(g_{\mathrm{1}}^{\mathrm{0002}}{\psi }_{\mathrm{2}}+\mathrm{2}{g}_{\mathrm{1}}^{\mathrm{0011}}){\mathrm{e}}^{-\mathrm{2}\mathrm{i}{\omega }_{*}\tau }+\\ \mathrm{2}{\psi }_{\mathrm{2}}{g}_{\mathrm{1}}^{\mathrm{1100}}+g_{\mathrm{1}}^{\mathrm{2000}},\end{array}

G_{\mathrm{2}}^{zz}={\mathrm{e}}^{-\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}}\left(\mathrm{2}{g}_{\mathrm{2}}^{\mathrm{11}}{\psi }_{\mathrm{2}}+g_{\mathrm{2}}^{\mathrm{20}}\right),

\begin{array}{l}{G}_{\mathrm{1}}^{z\overline{z}}=(g_{\mathrm{1}}^{\mathrm{0011}}+g_{\mathrm{1}}^{\mathrm{1100}}){\overline{\psi }}_{\mathrm{2}}+|{\psi }_{\mathrm{2}}{|}^{\mathrm{2}}{g}_{\mathrm{1}}^{\mathrm{0002}}+\\ g_{\mathrm{1}}^{\mathrm{0011}}{\psi }_{\mathrm{2}}+g_{\mathrm{1}}^{\mathrm{1100}}{\psi }_{\mathrm{2}}+g_{\mathrm{1}}^{\mathrm{2000}},\end{array}

G_{\mathrm{2}}^{z\overline{z}}=\mathrm{2}\mathrm{R}\mathrm{e}({\psi }_{\mathrm{2}})g_{\mathrm{2}}^{\mathrm{11}}+g_{\mathrm{2}}^{\mathrm{20}},

G_{\mathrm{1}}^{\overline{z}\overline{z}}=\overline{G_{\mathrm{1}}^{zz}},G_{\mathrm{2}}^{\overline{z}\overline{z}}=\overline{G_{\mathrm{2}}^{zz}},

\begin{array}{l}{G}_{\mathrm{1}}^{zz\overline{z}}=[g_{\mathrm{1}}^{\mathrm{0003}}{\overline{\psi }}_{\mathrm{2}}{\psi }_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{g}_{\mathrm{1}}^{\mathrm{0012}}|{\psi }_{\mathrm{2}}{|}^{\mathrm{2}}+g_{\mathrm{1}}^{\mathrm{0012}}{\psi }_{\mathrm{2}}^{\mathrm{2}}+\\ (\mathrm{2}{w}_{\mathrm{1}}^{\mathrm{11}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0011}}+\\ \mathrm{2}{w}_{\mathrm{2}}^{\mathrm{11}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0002}}){\psi }_{\mathrm{2}}+\\ \mathrm{2}{w}_{\mathrm{2}}^{\mathrm{11}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0011}}]{\mathrm{e}}^{-\mathrm{i}{\omega }_{*}\tau }+\\ [(w_{\mathrm{1}}^{\mathrm{20}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0011}}+\\ w_{\mathrm{2}}^{\mathrm{20}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0002}}){\overline{\psi }}_{\mathrm{2}}+\\ w_{\mathrm{2}}^{\mathrm{20}}(-\frac{\tau }{{\sigma }_{*}})g_{\mathrm{1}}^{\mathrm{0011}}]{\mathrm{e}}^{\mathrm{i}{\omega }_{*}\tau }+\\ \left(w_{\mathrm{1}}^{\mathrm{20}}\left(\mathrm{0}\right)g_{\mathrm{1}}^{\mathrm{1100}}+g_{\mathrm{1}}^{\mathrm{2100}}\right){\overline{\psi }}_{\mathrm{2}}+\\ \left(\mathrm{2}{w}_{\mathrm{1}}^{\mathrm{11}}\left(\mathrm{0}\right)g_{\mathrm{1}}^{\mathrm{1100}}+\mathrm{2}{g}_{\mathrm{1}}^{\mathrm{2100}}\right){\psi }_{\mathrm{2}}+\\ \left(w_{\mathrm{2}}^{\mathrm{20}}(\mathrm{0})+\mathrm{2}{w}_{\mathrm{2}}^{\mathrm{11}}\left(\mathrm{0}\right)\right)g_{\mathrm{1}}^{\mathrm{1100}}+g_{\mathrm{1}}^{\mathrm{3000}}+\\ \left(\mathrm{2}{w}_{\mathrm{1}}^{\mathrm{11}}(\mathrm{0})+w_{\mathrm{1}}^{\mathrm{20}}\left(\mathrm{0}\right)\right)g_{\mathrm{1}}^{\mathrm{2000}},\end{array}

\begin{array}{l}{G}_{\mathrm{2}}^{zz\overline{z}}=(g_{\mathrm{2}}^{\mathrm{21}}{\overline{\psi }}_{\mathrm{2}}+(\mathrm{2}{w}_{\mathrm{1}}^{\mathrm{11}}(-\mathrm{1}){\psi }_{\mathrm{2}}+\\ \mathrm{2}{w}_{\mathrm{2}}^{\mathrm{11}}(-\mathrm{1}))g_{\mathrm{2}}^{\mathrm{11}}+\mathrm{2}{g}_{\mathrm{2}}^{\mathrm{20}}{w}_{\mathrm{1}}^{\mathrm{11}}(-\mathrm{1})+\\ \mathrm{2}{g}_{\mathrm{2}}^{\mathrm{21}}{\psi }_{\mathrm{2}}+g_{\mathrm{2}}^{\mathrm{30}}){\mathrm{e}}^{-\mathrm{i}{\omega }_{*}{\sigma }_{*}}+\\ (g_{\mathrm{2}}^{\mathrm{11}}{w}_{\mathrm{1}}^{\mathrm{20}}(-\mathrm{1}){\overline{\psi }}_{\mathrm{2}}+\\ g_{\mathrm{2}}^{\mathrm{20}}{w}_{\mathrm{1}}^{\mathrm{20}}(-\mathrm{1})+\\ g_{\mathrm{2}}^{\mathrm{11}}{w}_{\mathrm{2}}^{\mathrm{20}}(-\mathrm{1})){\mathrm{e}}^{\mathrm{i}{\omega }_{*}{\sigma }_{*}}\mathrm{。}\end{array}

因此可得，

\begin{array}{l}{g}_{\mathrm{20}}={\sigma }_{*}\overline{D}(G_{\mathrm{1}}^{zz}+{\overline{\rho }}_{\mathrm{2}}{G}_{\mathrm{2}}^{zz}),\\ g_{\mathrm{11}}={\sigma }_{*}\overline{D}(G_{\mathrm{1}}^{z\overline{z}}+{\overline{\rho }}_{\mathrm{2}}{G}_{\mathrm{2}}^{z\overline{z}}),\\ g_{\mathrm{02}}={\sigma }_{*}\overline{D}(G_{\mathrm{1}}^{\overline{z}\overline{z}}+{\overline{\rho }}_{\mathrm{2}}{G}_{\mathrm{2}}^{\overline{z}\overline{z}}),\\ g_{\mathrm{21}}={\sigma }_{*}\overline{D}(G_{\mathrm{1}}^{zz\overline{z}}+{\overline{\rho }}_{\mathrm{2}}{G}_{\mathrm{2}}^{zz\overline{z}})\mathrm{。}\end{array}(19)

现计算g_{\mathrm{21}}中的{\mathit{w}}_{\mathrm{20}}(\theta )和{\mathit{w}}_{\mathrm{11}}(\theta )。由式（16）和式（18）可得，

\begin{array}{l}\frac{\mathrm{d}\mathit{W}(t,\theta )}{\mathrm{d}t}=\frac{\mathrm{d}{\mathit{U}}_t}{\mathrm{d}t}-\frac{\mathrm{d}z}{\mathrm{d}t}\mathit{\psi }(\theta )-\frac{\mathrm{d}\overline{z}}{\mathrm{d}t}\overline{\mathit{\psi }}(\theta )=\\ \left\{\begin{array}{l}\mathit{A}\mathit{W}-\mathrm{2}\mathrm{R}\mathrm{e}\left(g\left(z,\overline{z}\right)\mathit{\psi }(\theta )\right),\\ \theta \in [-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0});\\ \mathit{A}\mathit{W}-\mathrm{2}\mathrm{R}\mathrm{e}\left(g\left(z,\overline{z}\right)\mathit{\psi }\left(\mathrm{0}\right)\right)+\\ G\left(\mathit{W}\left(t,\theta \right)+\mathrm{2}\mathrm{R}\mathrm{e}\left(z(t)\mathit{\psi }(\theta )\right),\mathrm{0}\right),\\ \theta =\mathrm{0},\end{array}\right.\end{array}

并记

\begin{array}{l}\frac{\mathrm{d}\mathit{W}(t,\theta )}{\mathrm{d}t}=\mathit{A}\mathit{W}+\mathit{H}\left(z,\overline{z},\theta \right)=\\ \mathit{A}\mathit{W}+{\mathit{H}}_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\mathit{H}}_{\mathrm{11}}(\theta )z\overline{z}+\\ {\mathit{H}}_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\dots ,\end{array}(20)

式中：

{\mathit{H}}_{\mathrm{20}}(\theta )=\left\{\begin{array}{l}-\left(g_{\mathrm{20}}\mathit{\psi }(\theta )+{\overline{g}}_{\mathrm{20}}\overline{\mathit{\psi }}(\theta )\right),\\ \theta \in [-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0});\\ -\left(g_{\mathrm{20}}\mathit{\psi }(\theta )+{\overline{g}}_{\mathrm{20}}\overline{\mathit{\psi }}(\theta )\right)+{\sigma }_{*}{G}_{zz},\\ \theta =\mathrm{0},\end{array}\right.

{\mathit{H}}_{\mathrm{11}}(\theta )=\left\{\begin{array}{l}-\left(g_{\mathrm{11}}\mathit{\psi }(\theta )+{\overline{g}}_{\mathrm{11}}\overline{\mathit{\psi }}(\theta )\right),\\ \theta \in [-\mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\tau }{{\sigma }_{*}},\mathrm{1}\},\mathrm{0});\\ -\left(g_{\mathrm{11}}\mathit{\psi }(\theta )+{\overline{g}}_{\mathrm{11}}\overline{\mathit{\psi }}(\theta )\right)+{\sigma }_{*}{G}_{z\overline{z}},\\ \theta =\mathrm{0}\mathrm{。}\end{array}\right.

另一方面，

\begin{array}{l}\frac{\mathrm{d}\mathit{W}\left(t,\theta \right)}{\mathrm{d}t}={\mathit{w}}_{\mathrm{20}}(\theta )z\frac{\mathrm{d}z}{\mathrm{d}t}+{\mathit{w}}_{\mathrm{11}}(\theta )\frac{\mathrm{d}z}{\mathrm{d}t}\overline{z}+\\ {\mathit{w}}_{\mathrm{11}}(\theta )\frac{\mathrm{d}z}{\mathrm{d}t}\overline{z}+{\mathit{w}}_{\mathrm{02}}(\theta )\frac{\mathrm{d}\overline{z}}{\mathrm{d}t}\overline{z}+\dots \mathrm{。}\end{array}(21)

比较式（20）和式（21）的系数有，

\begin{array}{l}\left(\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}{\mathit{I}}_{\mathrm{2}}-\mathit{A}\right){\mathit{w}}_{\mathrm{20}}(\theta )={\mathit{H}}_{\mathrm{20}}(\theta ),\\ -\mathit{A}{\mathit{w}}_{\mathrm{11}}(\theta )={\mathit{H}}_{\mathrm{11}}(\theta ),\dots \end{array}(22)

从式（22）可解得

\begin{array}{l}{\mathit{w}}_{\mathrm{20}}(\theta )=(\frac{\mathrm{i}{g}_{\mathrm{20}}}{{\omega }_{*}{\sigma }_{*}}\mathit{\psi }(\theta )+\frac{\mathrm{i}{\overline{g}}_{\mathrm{02}}}{\mathrm{3}{\omega }_{*}{\sigma }_{*}}\overline{\mathit{\psi }}(\theta ))+\\ {\mathit{E}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}\theta },\end{array}

{\mathit{w}}_{\mathrm{11}}(\theta )=(\frac{-\mathrm{i}{g}_{\mathrm{11}}}{{\omega }_{*}{\sigma }_{*}}\mathit{\psi }(\theta )+\frac{\mathrm{i}{\overline{g}}_{\mathrm{11}}}{{\omega }_{*}{\sigma }_{*}}\overline{\mathit{\psi }}(\theta ))+{\mathit{E}}_{\mathrm{2}}\mathrm{。}(23)

令\theta =\mathrm{0}，可得

\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}{E}_{\mathrm{1}}-A\left(E_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}\theta }\right)={\sigma }_{*}{G}_{zz},

-A(E_{\mathrm{2}})={\sigma }_{*}{G}_{z\overline{z}}\mathrm{。}

进而求得式（23）中的E₁和E₂：

\begin{array}{l}{E}_{\mathrm{1}}=(\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}{I}_{\mathrm{2}}-L_{\mathrm{0}}-L_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}{\omega }_{*}\tau }-\\ L_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}{\omega }_{*}{\sigma }_{*}}{)}^{-\mathrm{1}}·{\sigma }_{*}{G}_{zz},\end{array}

E_{\mathrm{2}}={\left(-L_{\mathrm{0}}-L_{\mathrm{1}}-L_{\mathrm{2}}\right)}^{-\mathrm{1}}{\sigma }_{*}{G}_{z\overline{z}}\mathrm{。}

将式（23）代入式（19）可得到g_{\mathrm{21}}的显式表达式。进而可得，

c_{\mathrm{1}}(\mathrm{0})=\frac{\mathrm{i}}{\mathrm{2}{\omega }_{*}{\sigma }_{*}}(g_{\mathrm{11}}{g}_{\mathrm{20}}-\mathrm{2}|g_{\mathrm{11}}{|}^{\mathrm{2}}-\frac{|g_{\mathrm{02}}{|}^{\mathrm{2}}}{\mathrm{3}})+\frac{g_{\mathrm{21}}}{\mathrm{2}},

{\mu }_{\mathrm{2}}=-\frac{\mathrm{R}\mathrm{e}(c_{\mathrm{1}}\left(\mathrm{0}\right))}{\mathrm{R}\mathrm{e}({\lambda }^{\text{'}}\left({\sigma }_{*}\right))},{\beta }_{\mathrm{2}}=\mathrm{2}\mathrm{R}\mathrm{e}(c_{\mathrm{1}}\left(\mathrm{0}\right)),

{\tau }_{\mathrm{2}}=-\frac{\mathrm{I}\mathrm{m}(c_{\mathrm{1}}\left(\mathrm{0}\right))+{\mu }_{\mathrm{2}}\mathrm{I}\mathrm{m}({\lambda }^{\text{'}}\left({\sigma }_{*}\right))}{{\omega }_{*}{\sigma }_{*}}\mathrm{。}

当时，模型（3）将在时出现周期轨；当时，分支出的周期轨是轨道（不）稳定的；当时，周期轨的周期随着σ的增加而增加（减少）。

3 数值例子

通过一个数值例子佐证上述理论结果的正确性，并探究成熟时滞与消化时滞对模型（3）动力学的影响。参数取自文献［5］：r=\mathrm{0.4}，k=\mathrm{20}，d_J=\mathrm{0.1}，d_{\mathrm{1}}=\mathrm{0.01}，d_{\mathrm{2}}=\mathrm{0.05}，\alpha =\mathrm{0.5}，h=\mathrm{0.6}，\beta =\mathrm{0.4}，m=\mathrm{0.05}。可确定的可行域和S_{\mathrm{0}}^{\pm }(\tau ,\omega )曲线如图1所示，（τ，σ）平面上的分支曲线如图2所示。图2中箭头所指向的区域为不稳定区域。

由图2可知，成熟时滞τ可引起多次稳定性切换，而较大的消化时滞σ具有失稳效应。例如，固定τ=15，可计算出\omega =\mathrm{0}.\mathrm{029}\mathrm{8}，{\sigma }_{*}=\mathrm{2}.\mathrm{913}\mathrm{5}，通过计算规范型可得c_{\mathrm{1}}(\mathrm{0})=-\mathrm{0.546}\mathrm{4}-\mathrm{0.286}\mathrm{1}\mathrm{i}，{\mu }_{\mathrm{2}}=\mathrm{2}\mathrm{288.778}\mathrm{2}，{\beta }_{\mathrm{2}}=-\mathrm{1.092}\mathrm{8}，{\tau }_{\mathrm{2}}=\mathrm{21.920}\mathrm{6}。进而可得，当σ>σ_* =2.913 5时，E^* 失稳且分支出的周期轨道是渐近稳定的。数值仿真结果（图3）验证了该结论。

4 结束语

探究了成熟时滞与消化时滞共同作用对一类捕食模型动力学的影响。首先，运用能够让双时滞同步变化的几何方法处理具有参数依赖型时滞的特征方程，由此得到时滞平面内的Hopf分支曲线。其次，运用中心流形定理和规范型理论，得到模型在Hopf分支点处的等价系统，确定了分支方向与周期轨道的稳定性。最后，结合数值算例发现，随着消化时滞参数的增大，模型将失去稳定性；随着成熟时滞增大，模型经过有限次稳定性开关后最终趋于稳定。理论结果揭示了生物过程非瞬时性可显著影响捕食和猎物的动力学性质。

参考文献

[1] 肖江龙，宋永利，夏永辉. 一个扩散模型中恐惧效应与集群行为协同诱导的时空动力学研究［J］. 数学物理学报， 2024， 44 （6）： 1577⁃1594.

[2] XIAO Jianglong， SONG Yongli， XIA Yonghui. Spatiotemporal dynamics induced by the interaction between fear and schooling behavior in a diffusive model［J］. Acta Mathematica Scientia， 2024， 44 （6）： 1577⁃1594.

[3] WANG Xiaoying， ZANETTE Liana， ZOU Xingfu. Modelling the fear effect in predator⁃prey interactions［J］. Journal of Mathematical Biology， 2016， 73（5）： 1179⁃1204.

[4] 宋鸽， 甘静雯. 一类具有避难所和时滞的非自治阶段结构捕食系统的动力学分析［J］. 信阳师范学院学报（自然科学版）， 2024， 37（2）： 203⁃209.

[5] SONG Ge， GAN Jingwen. Dynamical analyses of a nonautonomous prey-predator stage structure system with refuges and time delay［J］. Journal of Xinyang Normal University （Natural Science Edition）， 2024， 37（2）： 203⁃209.

[6] KUANG Yang. Delay differential equations with applications in population dynamics ［M］. New York： Academic Press， 1993.

[7] MA Xiaoke， SU Ying， ZOU Xingfu. Joint impact of maturation delay and fear effect on the population dynamics of a predator⁃prey system［J］. SIAM Journal on Applied Mathematics， 2024， 84（4）： 1557⁃1579.

[8] LI Shuai， SONG Xinyu， HUANG Chengdai. Further study on the crossing curves in two⁃delay differential equations with delay⁃dependent coefficients［J］. Applied Mathematics Letters， 2024， 158： 109264.

[9] AN Qi， BERETTA E， KUANG Yang， et al. Geometric stability switch criteria in delay differential equations with two delays and delay dependent parameters［J］. Journal of Differential Equations， 2019， 266（11）： 7073⁃7100.

[10] HASSARD B D， KAZARINOFF N D， WAN Yieh‑Hei. Theory and applications of Hopf bifurcation［M］. New York： Cambridge University Press， 1981.

基金资助

国家自然科学基金项目(12401651)

信阳师范大学“南湖学者奖励计划”青年项目