ABSTRACT

The periodic cyclical neutrophil oscillating model of hematopoietic is discussed. Mathematically, the model is a multi-delay DDEs with high nonlinearity. Varying time delay brings forth Hopf bifurcation in system dynamic behavior. By calculating coefficients of normal form, it is analyzed that Hopf bifurcation is sub-critical. Numerical simulation coincides with the theoretical analyze result.

Keywords:Hematological Cell Model, Subcritical Hopf Bifurcation, Delay

造血系统的下临界Hopf分岔

马苏奇

中国农业大学数学系，北京

收稿日期：2017年9月1日；录用日期：2017年9月15日；发布日期：2017年9月21日

摘 要

本文讨论了造血干细胞引发周期噬中性粒子血液病的两房模型。数学上，该模型是一多时滞的时滞微分系统的复杂非线性模型。变化时滞可带来系统的Hopf分岔行为，通过计算范式，分析得出该分岔点处发生了下临界Hopf分岔。数值模拟与理论分析结果一致。

关键词 :造血系统，下临界Hopf分岔，时滞

Copyright © 2017 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

在系统生态学动力学行为研究中，造血干细胞模型一直是人们非常关注的一个话题。在数学上，人们致力于应用不同的数学模型来描述干细胞增殖分化的时空动力系统。造血干细胞经分化得到人类必须的三类血液细胞：白细胞、红细胞和血小板。造血功能失调会带来多种血液学疾病，典型的例子有长期发烧、“打摆子”，和慢性的嗜中性粒子白血球减小症。病人呈现出周期性的嗜中性粒子低水平振荡现象，振荡周期经生物学实验证实为11天~19天。应用非线性动力学模型来刻画血液细胞的动力学振荡现象可以帮助理解其复杂动力学机制。

近十几年来，Mackey等人给出了建立血液学疾病机理的基本数学框架 [1] - [7] 。从建模起，描述血液学疾病机理的数学模型是一含有多个滞量的时滞微分系统(DDEs)。造血系统由激发分化和增殖的机制以及成熟机制组成，由于细胞的分化和成熟都需要时间，因此在系统中含有多个时滞。一般地，细胞需要经历的成熟期为2到6天，变化不同参数组合，得到振荡的周期解 [8] 。医学上，人们应用集落刺激因子管理嗜中性粒子减少血液病 [3] [4] [5] 。集落刺激因子管理有助于降低血液细胞的程序性凋亡率和增加增殖分化所需时间。国内Lei等人发展了基于集落刺激因子管理嗜中性粒子的血液病振荡数学模型，一般地，他们给出的数学模型如下 [9] [10] ，

$\begin{array}{l}\frac{\text{d}Q}{\text{d}t}=-\left(\beta \left(Q\right)+k_N\left(N\right)+k_{\delta }\right)Q+2{\text{e}}^{-{\gamma }_s{\tau }_s}\beta \left(Q_{{\tau }_s}\right)Q_{{\tau }_s}\\ \frac{\text{d}N}{\text{d}t}=-{\gamma }_{N}N+{\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}{k}_N\left(N_{{\tau }_N}\right)Q_{{\tau }_N}\end{array}$ (1)

其中Holling函数

$\begin{array}{l}{k}_N\left(N\right)=\frac{f_0{\theta }_1^m}{{\theta }_1^m+N^m},\beta \left(Q\right)=\frac{k_0{\theta }_2^n}{{\theta }_2^n+Q^n}\\ k_N\left(N_{{\tau }_N}\right)=\frac{f_0{\theta }_1^m}{{\theta }_1^m+N{\left(t-{\tau }_N\right)}^m},\beta \left(Q\right)=\frac{k_0{\theta }_2^n}{{\theta }_2^n+Q{\left(t-{\tau }_s\right)}^n}\end{array}$

可以看出，模型(1)中干细胞以速率 $\beta$ 进入其增殖相，并以速率 $k_N$ 进入嗜中性粒子房。其它参数的生理学意义如下，

(1) 造血过程由激发分化和成熟机制组成，具体地，造血干细胞经由细胞分裂分化成为全部血液细胞包括白细胞、红细胞和血小板，并能自我更新， ${\gamma }_s$ 表示凋亡率，增殖相需要的天数是。

(2) 与嗜中性粒细胞不同的是，血液干细胞分化以速率 $k_{\delta }$ 分化成为其它血液细胞；

(3) 嗜中性粒细胞耦合时滞为 ${\tau }_N$ ，由两部分和组成，分别是其前驱细胞的增殖相需要的天数 ${\tau }_{NP}$ 和成熟相需要的天数 ${\tau }_{NM}$ 。 ${\eta }_{NP}$ 是增殖速率， ${\gamma }_0$ 是成熟相中的死亡率。因而，嗜中性粒细胞的组合相中的放大系数为

$A_N={\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}$

注意到速率 ${\gamma }_s,{\eta }_{NP},{\gamma }_0$ 和成熟时间 ${\tau }_{NM}$ 依赖于粒细胞集落刺激因子的管理。考虑到粒细胞集落刺激因子在医学上的成功使用，学者们在周期噬中性粒子的演化模型中，在增殖系数中用时滞来减弱系数。其估计的参数值列出如下 $k_0=2.0$ ， $f_0=0.4$ ， $k_{\delta }=0.01$ ， ${\theta }_1=0.36\times {10}^8细胞/千克$ ，， ${\gamma }_N=2.4$ ， ${\tau }_{NP}=5$ ， ${\eta }_{NP}=2.5379$ ， ${\gamma }_0=0.27$ ， ${\tau }_{NM}=6.6909$ ， ${\tau }_N=11.6909$ 。选取 ${\gamma }_s$ 和 ${\tau }_s$ 为分岔参数。文 [11] 分析了系统(1)的长周期解，以及倍周期分岔和多周期解共存等复杂非线性现象，并证实了Haurie生理学实验的结论 [1] [2] [3] [4] [5] ，即造血模型的周期解存在着一个很广泛的周期，例从11天到89天等。

本文首先分析了系统的Hopf分岔，其次通过范式分析，证实了Hopf分岔现象是下临界的Hopf分岔现象，并给出了下临界Hopf分岔的数值模拟。

2. Hopf分岔分析

设 $\left(Q_0,N_0\right)$ 是系统(1)的平衡点，计算方程在平衡点的特征值，如果方程(1)的特征方程有一对纯虚根 $\pm i{\omega }_0$ ，并满足横截条件，则方程(1)发生了Hopf分岔，且经由Hopf分岔，系统(1)出现稳定或不稳定的周期解。作坐标平移，系统的三次Taylor截断展开式为

$\begin{array}{l}{x}^{\prime }=a_{11}x+a_{12}y+b_{1}x\left(t-{\tau }_s\right)+f\left(x,x\left(t-{\tau }_s\right),y\right)\\ y^{\prime }=-{\gamma }_{N}y+b_{2}x\left(t-{\tau }_N\right)+b_{3}y\left(t-{\tau }_N\right)+g\left(x\left(t-{\tau }_N\right),y,y\left(t-{\tau }_N\right)\right)\end{array}$

其中

$\begin{array}{l}{a}_{11}=-\beta \left(Q^{*}\right)-k_N\left(N^{*}\right)-k_{\delta }-{\beta }^{\prime }\left(Q^{*}\right)Q^{*},a_{12}=-{k^{\prime }}_N\left(N^{*}\right)Q^{*},\\ b_1=2{\text{e}}^{-r_s{\tau }_s}\left(\beta \left(Q^{*}\right)+{\beta }^{\prime }\left(Q^{*}\right)Q^{*}\right),b_2={\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}k\left(N^{*}\right),\\ b_3={\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k}^{\prime }\left(N^{*}\right)Q^{*},\end{array}$

和

$\begin{array}{c}f\left(x,x\left(t-{\tau }_s\right),y\right)=-{\beta }^{\prime }\left(Q^{*}\right)x^2-\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)Q^{*}{x}^2-{k^{\prime }}_N\left(N^{*}\right)xy-\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}{y}^2\\ +2{\text{e}}^{-{\gamma }_s{\tau }_s}{\beta }^{\prime }\left(Q^{*}\right)x^2\left(t-{\tau }_s\right)+{\text{e}}^{-{\gamma }_s{\tau }_s}{\beta }^{″}\left(Q^{*}\right)Q^{*}{x}^2\left(t-{\tau }_s\right)\\ +{\text{e}}^{-{\gamma }_s{\tau }_s}{\beta }^{″}\left(Q^{*}\right)x^3\left(t-{\tau }_s\right)-\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)x^3-\frac{1}{6}{\beta }^{‴}\left(Q^{*}\right)Q^{*}{x}^3-\frac{1}{2}{k^{″}}_N\left(N^{*}\right)x{y}^2\\ -\frac{1}{6}{k^{‴}}_N\left(N^{*}\right)Q^{*}{y}^3+{\text{e}}^{-r_s{\tau }_s}{\beta }^{″}\left(Q^{*}\right)x^3\left(t-{\tau }_s\right)+\frac{1}{3}{\text{e}}^{-r_s{\tau }_s}{\beta }^{‴}\left(Q^{*}\right)Q^{*}{x}^3\left(t-{\tau }_s\right)\\ +o\left(\Vert x^3\Vert ,\Vert x^3(t-{\tau }_s)\Vert ,\Vert y^3\Vert \right)\end{array}$

$\begin{array}{l}g\left(x\left(t-{\tau }_N\right),y,y\left(t-{\tau }_N\right)\right)\\ ={\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{\prime }}_N\left(N^{*}\right)x\left(t-{\tau }_N\right)y\left(t-{\tau }_N\right)\\ +\frac{1}{2}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{″}}_N\left(N^{*}\right)Q^{*}{y}^2\left(t-{\tau }_N\right)+\frac{1}{2}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{″}}_N\left(N^{*}\right)x\left(t-{\tau }_N\right)y^2\left(t-{\tau }_N\right)\\ +\frac{1}{6}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{‴}}_N\left(N^{*}\right)Q^{*}{y}^3\left(t-{\tau }_N\right)+\frac{1}{6}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{‴}}_N\left(N^{*}\right)x\left(t-{\tau }_N\right)y^3\left(t-{\tau }_N\right)\end{array}$

我们引入向量 $m={\left(x,y\right)}^{\text{T}}$ ，设参数向量 $\mu =\left({\gamma }_s,{\tau }_s\right)$ ，那么方程(1)可以写成其向量形式

$\begin{array}{c}{m}^{\prime }\left(t\right)=L\left(\mu \right)m\left(t\right)+R_1\left(\mu \right)m\left(t-{\tau }_s\right)+R_2\left(\mu \right)m\left(t-{\tau }_N\right)\\ +F\left(m\left(t\right),m\left(t-{\tau }_s\right),m\left(t-{\tau }_N\right),\mu \right)\end{array}$ (2)

其中

$L=\left(\begin{array}{cc}a{}_{11}& a_{12}\\ 0& -{\gamma }_N\end{array}\right)$ , $R_1=\left(\begin{array}{cc}{b}_1& 0\\ 0& 0\end{array}\right)$ , $R_2=\left(\begin{array}{cc}0& 0\\ b_2& b_3\end{array}\right)$

和

$\begin{array}{l}F\left(m\left(t\right),m\left(t-{\tau }_s\right),m\left(t-{\tau }_N\right),\mu \right)\\ =\frac{1}{2}B\left(m\left(t\right),m\left(t\right)\right)+\frac{1}{2}\tilde{B}\left(m\left(t-{\tau }_s\right),m\left(t-{\tau }_s\right)\right)\\ +\frac{1}{2}\widehat{B}\left(m\left(t-{\tau }_N\right),m\left(t-{\tau }_N\right)\right)+\frac{1}{6}C\left(m\left(t\right),m\left(t\right),m\left(t\right)\right)\\ +\frac{1}{6}\tilde{C}\left(m\left(t-{\tau }_s\right),m\left(t-{\tau }_s\right),m\left(t-{\tau }_s\right)\right)+\frac{1}{6}\widehat{C}\left(m\left(t-{\tau }_N\right),m\left(t-{\tau }_N\right),m\left(t-{\tau }_N\right)\right)\end{array}$ (3)

设 $C$ 表示连续可微函数 $u:\left[-{\tau }_N,0\right]\to R^2$ 构成的巴拿赫空间，定义时间平移 $u_t\left(\theta \right)=m\left(t+\theta \right)$ , $\theta \in \left[-{\tau }_N,0\right]$ ，那么方程(2)的算子微分方程形式为

${m^{\prime }}_t=A{m}_t+F\left(m_t\right)$ (4)

这里，设在Hopf分岔临界值处的线性算子 $A$ 为

$Au\left(\theta \right)=\{\begin{array}{l}\frac{\text{d}}{\text{d}\theta }u\left(\theta \right),\theta \in \left[-{\tau }_N,0\right),\\ Lu\left(0\right)+R_{1}u\left(-{\tau }_s\right)+R_{2}u\left(-{\tau }_N\right),\theta =0,\end{array}$ (5)

非线性算子 $H$ 可被写作(6)

$H\left(u\right)\left(\theta \right)=\{\begin{array}{l}0,\theta \in \left[-{\tau }_N,0\right),\\ F\left(u\left(0\right),u\left(-{\tau }_s\right),u\left(-{\tau }_N\right)\right),\theta =0,\end{array}$ (6)

其中 $u\in \sum$ 。

设连续可微函数 $v:\left[0,{\tau }_N\right]\to R^2$ 构成空间 $C^{*}$ ，在上定义线性算子 $A$ 的伴随线性算子为

$A^{*}v\left(s\right)=\{\begin{array}{l}-\frac{\text{d}}{\text{d}s}v\left(s\right),s\in \left(0,{\tau }_N\right)\\ L^{\text{T}}v\left(0\right)+R_1^{\text{T}}v\left({\tau }_s\right)+R_2^{\text{T}}v\left({\tau }_N\right),s=0\end{array}$

定义双线性形式 $(,):C^{*}\times C\to R$ 为

$\langle v,u\rangle ={\overline{v}}^{\text{T}}\left(0\right)u\left(0\right)+{\displaystyle {\int }_{-{\tau }_s}^0{\overline{v}}^{\text{T}}\left(\xi +{\tau }_s\right)R_{1}u\left(\xi \right)\text{d}\xi }+{\displaystyle {\int }_{-{\tau }_N}^0{\overline{v}}^{\text{T}}\left(\xi +{\tau }_N\right)R_{2}u\left(\xi \right)\text{d}\xi }$ (7)

容易证明，有

$\langle v,Au\rangle =\langle A^{*}v,u\rangle$

成立。

在临界Hopf分岔值 $u_0=\left({\gamma }_s,{\tau }_s\right)$ 处，有一对纯虚特征值 $\Lambda =\left\{-i{\omega }_0,i{\omega }_0\right\}$ ，且其余特征值都具有负实部。因而 $C$ 可被分解为

$C=P_{\Lambda }\oplus Q_{\Lambda }$ ，其中 $P_{\Lambda }$ 为相应于 $\Lambda$ 的特征空间，而 $Q_{\Lambda }$ 是其补空间。计算得 $A\left(0\right)$ 相应于特征值 $i{\omega }_0$ 的特征向量为 $q\left(\theta \right)=\left(\begin{array}{c}{\gamma }_N+i{\omega }_0-b_3{\text{e}}^{-i{\omega }_0{\tau }_N}\\ b_2{\text{e}}^{-i{\omega }_0{\tau }_N}\end{array}\right){\text{e}}^{i{\omega }_0\theta }$ ， $-{\tau }_N\le \theta \le 0$ ， $A^{*}\left(0\right)$ 相应于 $-i{\omega }_0$ 的特征向量为 $q^{*}\left(s\right)=l\left(\begin{array}{c}-b_2{\text{e}}^{i{\omega }_0{\tau }_N}\\ a_{11}+i{\omega }_0+b_1{\text{e}}^{i{\omega }_0{\tau }_s}\end{array}\right){\text{e}}^{i{\omega }_{0}s}$ ， $0\le s\le {\tau }_N$ ,满足 $\langle q^{*},q\rangle =1$ ，进一步可证明 $\langle q^{*},\overline{q}\rangle =0$ 。

注意到

$\begin{array}{l}Aq\left(\theta \right)=i{\omega }_{0}q\left(\theta \right),-{\tau }_N\le \theta \le 0,\\ A^{*}{q}^{*}\left(s\right)=-i{\omega }_0{q}^{*}\left(s\right),0\le s\le {\tau }_N\end{array}$ (8)

每个 $u\in C$ 可以被写作

$u_t=zq+\overline{z}\overline{q}+w$ (9)

使得

$z=\langle q^{*},m_t\rangle ,w\in Q_{\Lambda }$ (10)

对(9)两边微分，注意到(8)式成立，因而可得

$\begin{array}{c}{z}^{\prime }\left(t\right)=\langle q^{*},{m^{\prime }}_t\rangle =\langle q^{*},Azq\rangle +\langle q^{*},H\left(zq+\overline{z}\overline{q}+w\left(t\right)\right)\rangle \\ =i{\omega }_{0}z+{\overline{q}}^{*}\left(0\right)F^{*}\left(z,\overline{z}\right)\end{array}$ (11)

其中

$F^{*}\left(z,\overline{z}\right)=F\left(2\mathrm{Re}\left\{zq\left(\theta \right)\right\}+w\left(t\right)\right)$ (12)

3. 二维中心流形

将 $w=v\left(z,\overline{z}\right)$ 展开成级数

$w=v\left(z,\overline{z}\right)=\left(\begin{array}{c}{\displaystyle \sum _{2\le i+j}\frac{v_{ij}^{\left(1\right)}}{i!j!}\left(z,\overline{z}\right)}\\ {\displaystyle \sum _{2\le i+j}\frac{v_{ij}^{\left(2\right)}}{i!j!}\left(z,\overline{z}\right)}\end{array}\right)$ (13)

其中系数 $\left\{v_{ij}\right\}$ 依赖于 $\theta$ ，令 $G\left(z,\overline{z}\right)={\overline{q}}^{*}\left(0\right)F^{*}\left(z,\overline{z}\right)$ ，把方程(11)写作

$z^{\prime }\left(t\right)=i{\omega }_{0}z+G\left(z,\overline{z}\right)$ (13)

其中 $G\left(z,\overline{z}\right)$ 有展式

$G\left(z,\overline{z}\right)={\displaystyle \sum _{2\le i+j}\frac{G_{ij}}{i!j!}\left(z,\overline{z}\right)}$ (14)

其中 $G_{ij}$ 为常数，计算见附录。

从方程(9)和(11)可得

$\frac{\text{d}v}{\text{d}t}=\{\begin{array}{l}Av-2\mathrm{Re}\left\{{\overline{q}}^{*}\left(0\right)\widehat{F}q\left(\theta \right)\right\},-{\tau }_N\le \theta <0,\\ Av-2\mathrm{Re}\left\{{\overline{q}}^{*}\left(0\right)\widehat{F}q\left(\theta \right)\right\}+\widehat{F},\theta =0,\end{array}$ (15)

方程(14)的第一个方程可以写成

$\frac{\text{d}v}{\text{d}t}=Av\left(z,\overline{z}\right)+R\left(z,\overline{z}\right)$ (16)

其中

$R\left(z,\overline{z}\right)=-2\mathrm{Re}\left\{{\overline{q}}^{*}\left(0\right)\widehat{F}q\left(\theta \right)\right\}=\left(\begin{array}{c}\frac{R_{20}^{\left(1\right)}}{2}{z}^2+R_{11}^{\left(1\right)}z\overline{z}+\frac{R_{02}^{\left(1\right)}}{2}{\overline{z}}^2+\cdots \\ \frac{R_{20}^{\left(2\right)}}{2}{z}^2+R_{11}^{\left(2\right)}z\overline{z}+\frac{R_{02}^{\left(2\right)}}{2}{\overline{z}}^2+\cdots \end{array}\right)$

其中系数 $R_{ij}^{\left(1,2\right)}=R_{ij}^{\left(1,2\right)}\left(\theta \right)$ 。

由于

$\begin{array}{l}-2\mathrm{Re}\left\{{\overline{q}}^{*}\left(0\right)\widehat{F}q\left(\theta \right)\right\}\\ =-G\left(z,\overline{z}\right)q\left(\theta \right)-\overline{G}\left(z,\overline{z}\right)\overline{q}\left(\theta \right)\\ =-\left(\frac{G_{20}}{2}{z}^2+G_{11}z\overline{z}+\frac{G_{02}}{2}{\overline{z}}^2\right)q\left(\theta \right)-\left(\frac{{\overline{G}}_{20}}{2}{\overline{z}}^2+{\overline{G}}_{11}z\overline{z}+\frac{{\overline{G}}_{02}}{2}{z}^2\right)\overline{q}\left(\theta \right)\end{array}$

所以

$R_{20}^{\left(1\right)}=-G_{20}{q}_1{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{02}{\overline{q}}_1{\text{e}}^{-i{\omega }_0\theta }$ , $R_{20}^{\left(2\right)}=-G_{20}{q}_2{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{02}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }$

$R_{11}^{\left(1\right)}=-G_{11}{q}_1{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{11}{\overline{q}}_1{\text{e}}^{-i{\omega }_0\theta }$ , $R_{11}^{\left(2\right)}=-G_{11}{q}_2{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{11}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }$

$R_{02}^{\left(1\right)}=-G_{02}{q}_1{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{20}{\overline{q}}_1{\text{e}}^{-i{\omega }_0\theta }$ , $R_{02}^{\left(2\right)}=-G_{02}{q}_2{\text{e}}^{i{\omega }_0\theta }-{\overline{G}}_{20}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }$

为了决定分岔解的行为，通过下式计算 $w=v\left(z,\overline{z}\right)$ 。

$\begin{array}{l}A{v}_{20}^{\left(1,2\right)}=2i{\omega }_0{v}_{20}^{\left(1,2\right)}-R_{20}^{\left(1,2\right)},\\ A{v}_{11}^{\left(1,2\right)}=-R_{11}^{\left(1,2\right)},\\ A{v}_{02}^{\left(1,2\right)}=-2i{\omega }_0{v}_{02}^{\left(1,2\right)}-R_{02}^{\left(1,2\right)}\end{array}$ (17)

注意到算子 $A$ 的定义，(17)式即

$\begin{array}{l}\{\begin{array}{l}\frac{\text{d}{v}_{20}^{\left(1\right)}}{\text{d}t}=2i{\omega }_0{v}_{20}^{\left(1\right)}+G_{20}{q}_1{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{02}{\overline{q}}_1{\text{e}}^{-i{\omega }_0\theta }\\ \frac{\text{d}{v}_{20}^{\left(2\right)}}{\text{d}t}=2i{\omega }_0{v}_{20}^{\left(2\right)}+G_{20}{q}_2{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{02}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }\end{array}\\ \{\begin{array}{l}\frac{\text{d}{v}_{11}^{\left(1\right)}}{\text{d}t}=G_{11}{q}_1{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{11}{\overline{q}}_1{\text{e}}^{-i{\omega }_0\theta }\\ \frac{\text{d}{v}_{11}^{\left(2\right)}}{\text{d}t}=G_{11}{q}_2{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{11}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }\end{array}\\ \{\begin{array}{l}\frac{\text{d}{v}_{02}^{\left(1\right)}}{\text{d}t}=-2i{\omega }_0{v}_{02}^{\left(1\right)}+G_{02}{q}_1{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{20}{\overline{\sout{q}}}_1{\text{e}}^{-i{\omega }_0\theta }\\ \frac{\text{d}{v}_{02}^{\left(2\right)}}{\text{d}t}=-2i{\omega }_0{v}_{02}^{\left(1\right)}+G_{02}{q}_2{\text{e}}^{i{\omega }_0\theta }+{\overline{G}}_{20}{\overline{q}}_2{\text{e}}^{-i{\omega }_0\theta }\end{array}\end{array}$ ， $-{\tau }_N\le \theta <0$ (18)

相似地，由方程(15)的第二个方程得到初值条件

$\{\begin{array}{l}{a}_{11}{v_{20}^{\left(1\right)}|}_{\theta =0}+a_{12}{v_{20}^{\left(2\right)}|}_{\theta =0}-2i{\omega }_0{v_{20}^{\left(1\right)}|}_{\theta =0}+b_1{v_{20}^{\left(1\right)}|}_{\theta =-{\tau }_s}+R_{20}^{\left(1\right)}+2f_0=0\\ -{\gamma }_N{v_{20}^{\left(2\right)}|}_{\theta =0}-2i{\omega }_0{v_{20}^{\left(2\right)}|}_{\theta =0}+b_2{v_{20}^{\left(1\right)}|}_{\theta =-{\tau }_N}+b_3{v_{20}^{\left(2\right)}|}_{\theta =-{\tau }_N}+R_{20}^{\left(2\right)}+2g_0=0\end{array}$

$\{\begin{array}{l}{a}_{11}{v_{11}^{\left(1\right)}|}_{\theta =0}+a_{12}{v_{11}^{(2)}|}_{\theta =0}+b_1{v_{11}^{\left(1\right)}|}_{\theta =-{\tau }_s}+R_{11}^{\left(1\right)}+f_1=0\\ -{\gamma }_N{v_{11}^{\left(2\right)}|}_{\theta =0}+b_2{v_{11}^{\left(1\right)}|}_{\theta =-{\tau }_N}+b_3{v_{11}^{\left(2\right)}|}_{\theta =-{\tau }_N}+R_{11}^{\left(2\right)}+g_1=0\end{array}$ (19)

$\{\begin{array}{l}{a}_{11}{v_{02}^{\left(1\right)}|}_{\theta =0}+a_{12}{v_{02}^{\left(2\right)}|}_{\theta =0}+2i{\omega }_0{v_{02}^{\left(1\right)}|}_{\theta =0}+b_1{v_{02}^{\left(1\right)}|}_{\theta =-{\tau }_s}+R_{02}^{\left(1\right)}+2f_2=0\\ -{\gamma }_N{v_{02}^{\left(2\right)}|}_{\theta =0}+2i{\omega }_0{v_{02}^{\left(2\right)}|}_{\theta =0}+b_2{v_{02}^{\left(1\right)}|}_{\theta =-{\tau }_N}+b_3{v_{02}^{\left(2\right)}|}_{\theta =-{\tau }_N}+R_{02}^{\left(2\right)}+2g_2=0\end{array}$

对方程(18)积分在初值条件(19)下积分，可得出

$\{\begin{array}{l}{v}_{20}^{\left(1\right)}=E_1{\text{e}}^{2i{\omega }_0\theta }+\frac{i{G}_{20}{q}_1}{{\omega }_0}{\text{e}}^{i{\omega }_0\theta }+\frac{i{\overline{G}}_{02}{\overline{q}}_1}{3{\omega }_0}{\text{e}}^{-i{\omega }_0\theta }\\ v_{20}^{\left(2\right)}=E_2{\text{e}}^{2i{\omega }_0\theta }+\frac{i{G}_{20}{q}_2}{{\omega }_0}{\text{e}}^{i{\omega }_0\theta }+\frac{i{\overline{G}}_{02}{\overline{q}}_2}{3{\omega }_0}{\text{e}}^{-i{\omega }_0\theta }\end{array}$ (20.1)

$\{\begin{array}{l}{v}_{11}^{\left(1\right)}=E_3-\frac{i{G}_{11}{q}_1}{{\omega }_0}{\text{e}}^{i\omega {}_0\theta }+\frac{i{\overline{G}}_{11}{\overline{q}}_1}{{\omega }_0}{\text{e}}^{-i{\omega }_0\theta }\\ v_{11}^{\left(2\right)}=E_4-\frac{i{G}_{11}{q}_2}{{\omega }_0}{\text{e}}^{i\omega {}_0\theta }+\frac{i{\overline{G}}_{11}{\overline{q}}_2}{{\omega }_0}{\text{e}}^{-i{\omega }_0\theta }\end{array}$ (20.2)

(20.3)

其中系数 $E_1,E_2,E_3,E_4,E_5,E_6$ 由方程初值条件(19)决定。

4. 下临界Hopf分岔

当特征值穿越虚轴时，系统(1)的平衡点处发生了Hopf分岔，一般的Hopf分岔分为超临界的Hopf分岔和下临界的Hopf分岔。若当参数值小于临界值时分岔出不稳定的周期解，Hopf分岔被称为是下临界的；若当参数值大于临界值时分岔出的周期解是稳定的，Hopf分岔被称为是超临界的。对平衡点作分析，方程(1)的线性化方程的特征方程为

$\Delta \left(\lambda ,{\gamma }_s,{\tau }_s\right)=\left|\begin{array}{cc}{a}_{11}-\lambda +b_1{\text{e}}^{-\lambda {\tau }_s}& a_{12}\\ b_2{\text{e}}^{-\lambda {\tau }_N}& -{\gamma }_N-\lambda +b_3{\text{e}}^{-\lambda {\tau }_N}\end{array}\right|=0$ (21)

当有一对共轭复特征值穿越虚轴时，方程(1)的平衡点的稳定性发生改变，这时方程(21)确定了一对

图1. 下临界hopf点附近的周期解。红线表示不稳定周期解，蓝线表示稳定的周期解

Hopf分岔的临界值 ${\mu }_0=\left({\gamma }_s,{\tau }_s\right)=\left(0.0969,4.9789\right)$ 。并由方程(21)确定 $\mathrm{Re}\left\{{\frac{\text{d}\lambda }{\text{d}\tau }|}_{\tau ={\tau }_s}\right\}\ne 0$ ，即系统在Hopf临界值处满足横截性。

由特征方程(21)可以计算 $\frac{\text{d}\lambda }{\text{d}{\tau }_s}>0$ ，以上计算得到大于零的第一Lyapunov数，又由于 $\mu =-\frac{l_1\left(0\right)}{\mathrm{Re}\left\{{\frac{\text{d}\lambda }{\text{d}\tau }|}_{\tau ={\tau }_s}\right\}}<0$ ，从而Hopf分岔是下临界的。数值模拟得到该Hopf分岔点附近的周期解，如图1所示。

5. 总结

本文应用Hopf分岔原理研究造血系统中噬中性粒子数周期性振荡现象发生的机制。模型中引入集落刺激因子管理，有助于降低血液细胞的程序性凋亡率和增加增殖分化所需时间，在数学上，使得造血系统成为含时滞依赖物理参数的多时滞系统。变化时滞得到系统的Hopf分岔点，应用范式理论进一步得到该分岔是下临界Hopf分岔现象。连续延拓得到系统共存的稳定周期解和不稳定周期解。

文章引用

马苏奇. 造血系统的下临界Hopf分岔
Subcritical Hopf Bifurcation of a Hematopological System[J]. 应用数学进展, 2017, 06(06): 740-751. http://dx.doi.org/10.12677/AAM.2017.66089

参考文献 (References)