Periodic Oscillation in Neutrophil Models with Time Delays

Suqi Ma

Mathematical Department, College of Science, China Agricultural University, Beijing, China.

DOI: 10.4236/ijmnta.2017.64011                

Abstract

To understand dynamical characters of neutrophil dynamical behavior, the sensitivity of delay factors which has effects on system dynamic behavior is ubiquitous due to system’s highly nonlinearity. Here we prove that delay supports a subcritical Hopf bifurcation, underlying a feedback mechanism during stem cells proliferation process while changing its coefficient of amplification. The given cell model reproduces a bistable dynamic regime of blood cells and hysteresis. Applying multiple scale method, oscillation motion near Hopf point is discussed. The stability limit of steady state to be abruptly periodic solution is detected.

Keywords

Bautin Bifurcation, Time Delay, Neutrophil Dynamics

1. Introduction

Within its scope of dynamical discussion and biomedical discipline, stem cells are of great interest to give insight in understanding research work in its tissue organization. As a mathematical work, people devotes to use different methods to investigate both its spatial and temporal dynamics based on the understanding of stem cell proliferation and differentiation. Mathematical model which originates from hematopoietic stem cell compartment to periodic neutrophil blood diseases [1] [2] [3] [4] , also the production and regulation of blood cells are governed by delay differential equations. Account for the variable success of granulocyte-colony stimulating factor, reduce of the amplification coefficient in neutrophil line for the treatment of cyclical neutropenia in reducing oscillation is proposed [5] . Recently, physiologically and mathematically, the dependence of neutrophil response on the period of simulated chemotherapy and the secondary response to G-CSF administration is reported [6] [7] [8] [9] .

To support its function on tissue coupling, a clear link between the stem cell compartment and the differentiated mature cell lineages are motivated by delayed items to develop into fully hematopoietic system. The differentiation of hematopoietic stem cells take place in $G_0$ phase. Quiescent phase HSCs (hematopoietic stem cells) can enter into its proliferative phase with the assumption of undergoing mitosis during time ${\tau }_s$ . After a cell division, the total duration of both proliferative phase and maturation phase of neutropenia are experienced to release into circulation through the body. In general, considering tissue system composed of stem cells and cyclical neutrophils population, the mathematical model is formed by two equations with time delays which beyond illustration with schematic picture tools as shown in Figure 1.

The simple version related to DDEs (delay differential equations) to govern the system dynamics with above description have been reported as the following [6] [8] [9] :

$\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}}^{-r_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)

where holling function

Figure 1. A cartoon representation of hematopoietic system response to stem cell population and neutrophil population.

$\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_{{\tau }_s}\right)=\frac{k_0{\theta }_2^n}{{\theta }_2^n+Q{\left(t-{\tau }_s\right)}^n}\end{array}$ (2)

HSCs (hematopoietic stem cells) enter into the proliferative phase with rate $\beta$ (days⁻¹), and differentiate into the committed neutrophil compartment at a rate $k_N$ (days⁻¹). The physiological meaning of other parameters and the corresponding values are explained as the following:

1) The hematopoiesis process consists of mechanism of triggering differentiation and maturation. HSCs produce differentiated cells via cell division to form into all kinds of cell types(white cells, red blood cells and platelets) and capable of self-renewal. $r_s$ denotes the apoptosis rate (days⁻¹) and the duration of the proliferation phase is taken to be ${\tau }_s$ (days).

2) Different to the committed neutrophil line, cells in HSC proliferative phase are assumed to enter into the combined megaryocyte/erythrocyte lines at a rate $k_{\delta }$ .

3) The time delay considered in neutrophil coupling is ${\tau }_N$ , which denotes the total duration of the proliferation and maturation phases of a neutrophil precursor, that is, ${\tau }_N={\tau }_{NP}+{\tau }_{NM}$ . The proliferation rate is assumed as ${\eta }_{NP}$ whilst the death rate during maturation phase is ${\gamma }_0$ . Therefore, the amplification coefficient in neutrophil combined phases is listed as formula

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

Notice that the rates $r_s\mathrm{,}{\eta }_{NP}\mathrm{,}{\gamma }_0$ and maturation time ${\tau }_{NM}$ are dependent on the effects of chemotherapy and G-CSF administration. The estimated parameter values are listed in Table 1.

From point view of dynamics, stable equilibrium solution can enter into instability regime to bifurcate into oscillation mode periodically via super-critical Hopf/sub-critical Hopf bifurcation [10] - [15] . Commonly, people use method of linearization model near the steady states to get the exponential polynomial characteristic equation to analyze its asymptotic stability via roots attribution to such polynomials [4] [5] . Mathematically, the difficulty resides in the presence of two independent delays and the fact that some coefficients in the model equations depend upon these delays [16] [17] . Therefore, there is very

Table 1. The values of parameters used in system (1).

few studies on this topic theoretically. With the fast development of computation methods in all kinds of mathematical questions, authors detect the different dynamic regime corresponding to stable steady states/periodic oscillating modes/bistable periodic oscillating modes and phase-locked resonance phenomena, etc. [5] [18] [19] [20] [21] .

The observation of bistability regime of steady state and periodic solution is found in Model (1). To understand the dynamics of the hematopoietic system, we take regard the amplification coefficient in stem cells as bifurcation parameter. Hopf bifurcation is tracked at its unique critical value. By means of the perturbation procedure such as the method of multiple scales, a codimension-2 bifurcation of periodic motion is theoretically discussed. Finally, varying maturation time delay simultaneously, dynamical bifurcations respond to periodic solutions is discovered and verified near condimension-2 bifurcation point by applying numerical simulation methods.

2. Linear Stability Analysis

Assume system (1) has an equilibrium solution $E\left(Q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{N}^{\mathrm{<mo>\; *\; </mo>}}\right)$ which is assumed to be positive. Notice coefficients as $\beta \mathrm{,}{k}_N$ are expressed by Hilling function which has obvious different analytical characters when varying $n\mathrm{,}m$ . Hereafter $n=2,m=1$ are assumed. By simple calculation, it is proved the formula

$Q^{*}=\frac{{\gamma }_N{N}^{*}}{{\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}{k}_N\left(N^{*}\right)}$ (4)

Do transformation of phase variables $x\to Q-Q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}y\to N-N^{\mathrm{<mo>\; *\; </mo>}}$ , Equation (1) can be rewritten and truncated into its three order form as

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

where

$\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}$ (6)

and

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

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

We introduce variable vector z by assuming $z=\left(\begin{array}{c}x\\ y\end{array}\right)$ , and parameter vector $\mu$ as $\mu =\left({\tau }_s,r_s\right)$ . Equation (1) is written as its vector form

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

where

$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)$ (9)

and

$F\left(z\left(t\right),z\left(t-{\tau }_s\right),z\left(t-{\tau }_N\right),\mu \right)=\left(\begin{array}{c}f\left(x,x_{{\tau }_s},y\right)\\ g\left(x_{{\tau }_N},y,y_{{\tau }_N}\right)\end{array}\right)$ (10)

The hematopoietic model (1) has two delays in its dynamical description. As a common point view, we use the related linear system to analyze its asymptotic stability.

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

The corresponding characteristic equation is derived from Equation (11) as

$\Delta \left(\lambda ,{\tau }_s,{\gamma }_0\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$ (12)

At critical values of parameters, the onset of instability of linear system (11) may lead to oscillation solutions of nonlinear system. Hopf bifurcation can be tracked via the change of direction of characteristic root distribution which determined by Equation (12). Herein we solve possible values of Hopf bifurcation by software DDE-Biftool. The imaginary roots in complex plane are plotted in Figure 2(a) which shows Hopf bifurcation occurs since a pair of pure imaginary roots with zero real part appears in characteristic equation. If denote critical value of time delay as ${\tau }_s={\tau }_{sc}$ , the direction of imaginary roots transverse in complex plane is shown in Figure 2(b). As time delay ${\tau }_s$ decrease through ${\tau }_{sc}$ , the head imaginary root(which is complex) transverse from right

plane to left plane and $\Re \frac{\text{d}\lambda \left({\mu }_{cr}\right)}{\text{d}{\tau }_s}\ne 0$ is satisfied. Therefore, periodic

solutions may appear near the Hopf point.

(a) (b)

Figure 2. The characteristic roots of Equation (12) are plotted in complex plane: (a) While ${\tau }_s=2.9310$ ; (b) Varying time delay, head imaginary roots are observed to transverse through the vertical axis from right plane to left plane.

Figure 3 depicts a branch diagram of periodic solution emanating from the Hopf point as a function of ${\tau }_s$ and $r_s$ , respectively. Variation of solutions along the branch is characterized by their maximal and minimal values over the period for each computed point on the branch. As ${\tau }_s\mathrm{,}{r}_s$ grows from its Hopf point value, an instable periodic solution emanates which designates Hopf point as a subcritical Hopf point which brings bi-stability of system and hysteresis.

3. Multiple Scales Analysis

We consider the related vector form (8) in this section to discuss periodic solutions of system (1). Variable vector z depends on control vector parameter $\mu$ and has subcritical Hopf bifurcation point at critical value ${\mu }_c=\left(r_{sc}\mathrm{,}{\tau }_{sc}\right)$ . In the following, the multiple scale method is used to investigate system dynamical behavior around subcritical Hopf point.

Without loss of generality, set characteristic matrix $M+iN$ to satisfy

$M+iN=i\omega I-L\left({\mu }_c\right)-R_1\left({\mu }_c\right){\text{e}}^{-i\omega {\tau }_s}-R_2\left({\mu }_c\right){\text{e}}^{-i\omega {\tau }_N}$ (13)

then the corresponding critical eigenvectors are governed by

$Mp=Nq,Mq=-NP$ (14)

and

$rM=sN,sM=-rN$ (15)

with the assumption $r+is=\left(r_1+i{s}_1,r_2+i{s}_2\right),p+iq={\left(p_1+i{q}_1,p_2+i{q}_2\right)}^{\text{T}}$ .

Subsequently, we introduce two time scales $T_0=t$ , $T_1=ϵt$ , $T_2={ϵ}^{2}t$ in Equation(8) to rewrite it into the following formula

$\frac{\text{d}z}{\text{d}t}=\frac{\partial z}{\partial T_0}+ϵ\frac{\partial z}{\partial T_1}+ϵ\frac{\partial z}{\partial T_2}$ (16)

(a) (b)

Figure 3. Branch of periodic solutions is seen to emanate from its subcritical Hopf point to form hysteresis which give rise to a stable and an unstable periodic solutions simultaneously: (a) Subcritical Hopf bifurcation occurs at critical value ${\tau }_s=2.9310$ ; (b) Subcritical Hopf bifurcation occurs at critical value $r_s=0.1982$ .

According to multiple scale analysis, the mono-parameter family of solutions of Equation (16) may expand into its two scale form

$z\left(t\right)=z\left(T_0,T_1,T_2\right)=z_0\left(T_0,T_1,T_2\right)+ϵz_1\left(T_0,T_1,T_2\right)+{ϵ}^2{z}_2\left(T_0,T_1,T_2\right)$ (17)

and correspondingly,

$\begin{array}{c}z\left(t-{\tau }_{sc}\right)=z\left(T_0-{\tau }_{sc},T_1-{\tau }_{sc},T_2-{\tau }_{sc}\right)\\ =z_0\left(T_0-{\tau }_{sc},T_1-{\tau }_{sc},T_2-{\tau }_{sc}\right)\\ +ϵz_1\left(T_0-{\tau }_{sc},T_1-{\tau }_{sc},T_2-{\tau }_{sc}\right)\\ +{ϵ}^2{z}_2\left(T_0-{\tau }_{sc},T_1-{\tau }_{sc},T_2-{\tau }_{sc}\right)\end{array}$

$\begin{array}{c}z\left(t-{\tau }_N\right)=z\left(T_0-{\tau }_N,T_1-{\tau }_N,T_2-{\tau }_N\right)\\ =z_0\left(T_0-{\tau }_N,T_1-{\tau }_N,T_2-{\tau }_N\right)\\ +ϵz_1\left(T_0-{\tau }_N,T_1-{\tau }_N,T_2-{\tau }_N\right)\\ +{ϵ}^2{z}_2\left(T_0-{\tau }_N,T_1-{\tau }_N,T_2-{\tau }_N\right).\end{array}$ (18)

where $z_i\left(T_0,T_1,T_2\right)=\left(\begin{array}{c}{z}_{i1}\left(T_0,T_1,T_2\right)\\ z_{i2}\left(T_0,T_1,T_2\right)\end{array}\right)$ for $i=0,1,2$ . We set ${\tau }_s={\tau }_{sc}+{ϵ}^2{\tau }_{ϵ}$ , the above delay terms can both expand into their Taylor’s series as

$\begin{array}{c}z\left(t-{\tau }_s\right)=z_0\left(T_0-{\tau }_{sc},T_1,T_2\right)-{ϵ}^2{\tau }_{ϵ}\frac{\partial z_0}{\partial T_0}\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ -ϵ\left({\tau }_{sc}\frac{\partial z_0}{\partial T_1}\left(T_0-{\tau }_s,T_1,T_2\right)-z_1\left(T_0-{\tau }_{sc},T_1,T_2\right)\right)\\ -{ϵ}^2{\tau }_{sc}\left(\frac{\partial z_0}{\partial T_2}\left(T_0-{\tau }_{sc},T_1,T_2\right)+\frac{\partial z_1}{\partial T_1}\left(T_0-{\tau }_{sc},T_1,T_2\right)\right)\\ +{ϵ}^2{z}_2\left(T_0-{\tau }_{sc},T_1,T_2\right)+\frac{1}{2}{ϵ}^2{\tau }_{sc}^2\frac{{\partial }^2{z}_0}{\partial T_1^2}\left(T_0-{\tau }_{sc},T_1,T_2\right)\end{array}$

$\begin{array}{c}z\left(t-{\tau }_N\right)=z_0\left(T_0-{\tau }_N,T_1,T_2\right)\\ -ϵ\left({\tau }_N\frac{\partial z_0}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)-z_1\left(T_0-{\tau }_N,T_1,T_2\right)\right)\\ -{ϵ}^2{\tau }_N\left(\frac{\partial z_0}{\partial T_2}\left(T_0-{\tau }_N,T_1,T_2\right)+\frac{\partial z_1}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)\right)\\ +{ϵ}^2{z}_2\left(T_0-{\tau }_N,T_1-{\tau }_N,T_2-{\tau }_N\right)+\frac{1}{2}{ϵ}^2{\tau }_N^2\frac{{\partial }^2{z}_0}{\partial T_1^2}\left(T_0-{\tau }_N,T_1,T_2\right)\end{array}$ (19)

We substitute Equation (17) and Equation (19) into the vector differential equation (16), and equating separately coefficients of ${ϵ}^0$ to obtain the perturbation equation

$\begin{array}{c}\frac{\partial z_0}{\partial T_0}\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)=L\left({\mu }_c\right)z_0\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)+R_1{\mu }_c{z}_0\left(T_0-{\tau }_s\mathrm{,}{T}_1\mathrm{,}{T}_2\right)\\ +R_2\left({\mu }_c\right)z_0\left(T-{\tau }_N\mathrm{,}{T}_1\mathrm{,}{T}_2\right)\end{array}$ (20)

Therefore, we have

$\frac{\partial z_{01}}{\partial T_0}\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)=a_{11}{z}_{01}\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)+a_{12}{z}_{02}\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)$ (21)

$\begin{array}{c}\frac{\partial z_{02}}{\partial T_0}\left(T_0,T_1,T_2\right)=-{\gamma }_N{z}_{02}\left(T_0,T_1,T_2\right)+b_2{z}_{01}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +b_3{z}_{02}\left(T_0-{\tau }_N,T_1,T_2\right)\end{array}$ (22)

The sequence of perturbation equations can be obtained by equating separately coefficients of ${ϵ}^1\mathrm{,}{ϵ}^2$ too.

$\begin{array}{l}{ϵ}^1:\frac{\partial z_{01}}{\partial T_1}\left(T_0,T_1,T_2\right)+\frac{\partial z_{11}}{\partial T_0}\left(T_0,T_1,T_2\right)\\ =a_{11}{z}_{11}\left(T_0,T_1,T_2\right)+a_{12}{z}_{12}\left(T_0,T_1,T_2\right)-b_1{\tau }_{sc}\frac{\partial z_{01}}{\partial T_1}\left(T_0-{\tau }_s,T_1,T_2\right)\\ +b_1{z}_{11}\left(T_0-{\tau }_{sc},T_1,T_2\right)+2{\text{e}}^{-r_s{\tau }_{sc}}{\beta }^{\prime }\left(Q^{*}\right)z_{01}^2\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ -\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)Q^{*}{z}_{01}^2\left(T_0,T_1,T_2\right)-\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}{z}_{02}^2\left(T_0,T_1,T_2\right)\\ +{\text{e}}^{-r_{sc}{\tau }_{sc}}{\beta }^{″}\left(Q^{*}\right)Q^{*}{z}_{01}^2\left(T_0-{\tau }_{sc},T_1,T_2\right)-{\beta }^{\prime }\left(Q^{*}\right)z_{01}^2\left(T_0,T_1,T_2\right)\\ -{k^{\prime }}_N\left(N^{*}\right)z_{01}\left(T_0,T_1,T_2\right)z_{02}\left(T_0,T_1,T_2\right),\end{array}$ (23)

$\begin{array}{l}{ϵ}^1:\frac{\partial z_{02}}{\partial T_1}\left(T_0,T_1,T_2\right)+\frac{\partial z_{12}}{\partial T_0}\left(T_0,T_1,T_2\right)\\ =-{\gamma }_N{z}_{12}\left(T_0,T_1,T_2\right)-b_2{\tau }_N\frac{\partial z_{01}}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)+b_2{z}_{11}\left(T_0-{\tau }_N,T_1,T_2\right)\\ -b_3{\tau }_N\frac{\partial z_{02}}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)+b_3{z}_{12}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{\prime }}_N\left(N^{*}\right)z_{01}\left(T_0-{\tau }_N,T_1,T_2\right)z_{02}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +\frac{1}{2}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{″}}_N\left(N^{*}\right)Q^{*}{z}_{02}^2\left(T_0-{\tau }_N,T_1,T_2\right),\end{array}$ (24)

$\begin{array}{l}{ϵ}^2:\frac{\partial z_{01}}{\partial T_2}\left(T_0,T_1,T_2\right)+\frac{\partial z_{11}}{\partial T_1}\left(T_0,T_1,T_2\right)+\frac{\partial z_{21}}{\partial T_0}\left(T_0,T_1,T_2\right)\\ =a_{11}{z}_{21}\left(T_0,T_1,T_2\right)+a_{12}{z}_{22}\left(T_0,T_1,T_2\right)-b_1{\tau }_{sc}\frac{\partial z_{01}}{\partial T_2}\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ +b_1{\tau }_{sc}\frac{\partial z_{11}}{\partial T_1}\left(T_0-{\tau }_{sc},T_1,T_2\right)+b_1{z}_{21}\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ +\frac{1}{2}{b}_1{\tau }_{sc}^2\frac{{\partial }^2{z}_{01}}{\partial T_1^2}\left(T_0-{\tau }_{sc},T_1,T_2\right)-b_1{\tau }_{ϵ}\frac{\partial z_{01}}{\partial T_0}\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ +{\beta }^{″}\left(Q^{*}\right){\text{e}}^{-r_{sc}{\tau }_{sc}}{z}_{01}^3\left(T_0-{\tau }_{sc},T_1,T_2\right)-\frac{1}{6}{\beta }^{‴}\left(Q^{*}\right)Q^{*}{Z}_{01}^3\left(T_0,T_1,T_2\right)\end{array}$

$\begin{array}{l}-\frac{1}{6}{k^{‴}}_N\left(N^{*}\right)Q^{*}{z}_{02}^3\left(T_0,T_1,T_2\right)+\frac{1}{3}{\text{e}}^{-r_{sc}{\tau }_{sc}}{\beta }^{‴}\left(Q^{*}\right)Q^{*}{z}_{01}^3\left(T_0-{\tau }_{sc},T_1,T_2\right)\\ -\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)z_{01}^3\left(T_0,T_1,T_2\right)-\frac{1}{2}{k^{″}}_N\left(N^{*}\right)z_{01}\left(T_0,T_1,T_2\right)z_{02}^2\left(T_0,T_1,T_2\right)\\ -{\beta }^{″}\left(Q^{*}\right)q^{*}{z}_{01}\left(T_0,T_1,T_2\right)z_{11}\left(T_0,T_1,T_2\right)-{k^{″}}_N\left(N^{*}\right)Q^{*}{z}_{02}\left(T_0,T_1,T_2\right)z_{12}\left(T_0,T_1,T_2\right)\\ +2{\text{e}}^{-r_{sc}{\tau }_{sc}}{\beta }^{″}\left(Q^{*}\right)Q^{*}{z}_{01}\left(T_0-{\tau }_s,T_1,T_2\right)z_{11}\left(T_0-{\tau }_s,T_1,T_2\right)\\ -{k^{\prime }}_N\left(N^{*}\right)\left(z_{11}\left(T_0,T_1,T_2\right)z_{02}\left(T_0,T_1,T_2\right)+z_{01}\left(T_0,T_1,T_2\right)z_{12}\left(T_0,T_1,T_2\right)\right),\end{array}$ (25)

$\begin{array}{l}{ϵ}^2:\frac{\partial z_{02}}{\partial T_2}\left(T_0,T_1,T_2\right)+\frac{\partial z_{12}}{\partial T_1}\left(T_0,T_1,T_2\right)+\frac{\partial z_{22}}{\partial T_0}\left(T_0,T_1,T_2\right)\\ =-{\gamma }_N{z}_{22}\left(T_0,T_1,T_2\right)-b_2{\tau }_N\frac{\partial z_{01}}{\partial T_2}\left(T_0-{\tau }_N,T_1,T_2\right)+b_2{\tau }_N\frac{\partial z_{11}}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +b_2{z}_{21}\left(T_0-{\tau }_N,T_1,T_2\right)+\frac{1}{2}{b}_2{\tau }_N^2\frac{{\partial }^2{z}_{01}}{\partial T_1^2}\left(T_0-{\tau }_N,T_1,T_2\right)\\ -b_3{\tau }_N\frac{\partial z_{02}}{\partial T_2}\left(T_0-{\tau }_N,T_1,T_2\right)+b_3{\tau }_N\frac{\partial z_{12}}{\partial T_1}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +b_3{z}_{22}\left(T_0-{\tau }_N,T_1,T_2\right)+\frac{1}{2}{b}_3{\tau }_N^2\frac{{\partial }^2{z}_{02}}{\partial T_1^2}\left(T_0-{\tau }_N,T_1,T_2\right)\end{array}$

$\begin{array}{l}+{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{″}}_N\left(N^{*}\right)z_{01}\left(T_0-{\tau }_N,T_1,T_2\right)z_{02}^2\left(T_0-{\tau }_N,T_1,T_2\right)\\ +\frac{1}{6}{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{‴}}_N\left(N^{*}\right)Q^{*}{z}_{02}^3\left(T_0-{\tau }_N,T_1,T_2\right)\\ +{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{\prime }}_N\left(N^{*}\right)z_{11}\left(T_0-{\tau }_N,T_1,T_2\right)z_{02}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{\prime }}_N\left(N^{*}\right)z_{01}\left(T_0-{\tau }_N,T_1,T_2\right)z_{12}\left(T_0-{\tau }_N,T_1,T_2\right)\\ +{\text{e}}^{{\tau }_{NP}{\eta }_{NP}-{\gamma }_0{\tau }_{NM}}{k^{″}}_N\left(N^{*}\right)Q^{*}{z}_{02}\left(T-{\tau }_N,T_1,T_2\right)z_{12}\left(T_0-{\tau }_N,T_1,T_2\right)\end{array}$ (26)

Further, we suppose Equation (8) has the following general solution

$z_0\left(T_0,T_1,T_2\right)=\left(\begin{array}{l}{z}_{01}\left(T_0,T_1,T_2\right)\hfill \\ z_{02}\left(T_0,T_1,T_2\right)\hfill \end{array}\right)=B\left(T_1,T_2\right)\exp \left(i{\omega }_0{T}_0\right)\left(p+iq\right)+cc$ (27)

where $cc$ is assumed to be the complex conjugate terms of the preceding terms in equation and frequency ${\omega }_0$ satisfies characteristic Equation (12) by root $\lambda =\pm i{\omega }_0$ .

Substitute $z_0\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)$ into Equation (21) and Equation (22), we have

$\frac{\partial z_1}{\partial T_0}\left(T_0,T_1,T_2\right)-\left(L\left(\mu \right)+R_1{\text{e}}^{i{\omega }_0{\tau }_{sc}}+R_2{\text{e}}^{i{\omega }_0{\tau }_N}\right)z_1\left(T_0,T_1,T_2\right)=S{T}_1+NS{T}_1$ (28)

where

$\begin{array}{l}S{T}_1=-\left(\begin{array}{c}\left(p_1+i{q}_1\right)\left(1+b_1{\tau }_{sc}{\text{e}}^{-i{\omega }_0{\tau }_{sc}}\right)\\ \left(p_1+i{q}_1\right)b_2{\tau }_N{\text{e}}^{-i{\omega }_0{\tau }_N}+\left(p_2+i{q}_2\right)\left(1+b_3{\tau }_N{\text{e}}^{-i{\omega }_0{\tau }_N}\right)\end{array}\right)\\ \frac{\partial B}{\partial T_1}\left(T_1,T_2\right)\exp \left(i{\omega }_0{T}_0\right)\left(p+iq\right)+cc\end{array}$ (29)

and

$S{T}_1=\left(\begin{array}{l}{\gamma }_1\hfill \\ {\gamma }_2\hfill \end{array}\right)B^2{\text{e}}^{2i{\omega }_0{T}_0}+\left(\begin{array}{l}{\gamma }_{g1}\hfill \\ {\gamma }_{g2}\hfill \end{array}\right)B\bar{B}+cc$ (30)

$\begin{array}{l}{\gamma }_1=-\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)Q^{*}{\left(p_1+i{q}_1\right)}^2-\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}{\left(p_2+i{q}_2\right)}^2\\ -{\beta }^{\prime }\left(Q^{*}\right){\left(p_1+i{q}_1\right)}^2+{\beta }^{″}\left(Q^{*}\right)Q^{*}{\left(p_1+i{q}_1\right)}^2{\text{e}}^{-r_s{\tau }_s}{\text{e}}^{2i{\omega }_0{\tau }_s}\\ +2{\beta }^{\prime }\left(Q^{*}\right){\left(p_1+i{q}_1\right)}^2-{k^{\prime }}_N\left(N^{*}\right)\left(p_1+i{q}_1\right)\left(p_2+i{q}_2\right)\\ {\gamma }_{g1}=-{\beta }^{\prime }\left(Q^{*}\right)\left(p_1^2+q_1^2\right)-{k^{\prime }}_N\left(N^{*}\right)\left(p_1+i{q}_1\right)\left(p_2-i{q}_2\right)\\ -\frac{1}{2}{\beta }^{″}\left(Q^{*}\right)Q^{*}\left(p_1^2+q_1^2\right)+{\beta }^{″}\left(Q^{*}\right)Q^{*}\left(p_1^2+q_1^2\right){\text{e}}^{-r_s{\tau }_s}\\ -\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}\left(p_2^2+q_2^2\right)+2{\beta }^{\prime }\left(Q^{*}\right)\left(p_1^2+q_1^2\right){\text{e}}^{-r_s{\tau }_s}\end{array}$

$\begin{array}{l}{\gamma }_2={\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}({k^{\prime }}_N\left(N^{*}\right)\left(p_1+i{q}_1\right)\left(p_2+i{q}_2\right)\underset{}{\overset{}{}}\\ +\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}{\left(p_2+i{q}_2\right)}^2){\text{e}}^{-2i{\omega }_0{\tau }_N},\\ {\gamma }_{g2}={\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}(k{N}^{\prime }\left(N^{*}\right)\left(p_1+i{q}_1\right)\left(p_2-i{q}_2\right)\underset{}{\overset{}{}}\\ +\frac{1}{2}{k^{″}}_N\left(N^{*}\right)Q^{*}\left(p_2^2+q_2^2\right)).\end{array}$ (31)

$S{T}_1$ and $NS{T}_1$ respectively represent the resonance terms and non-resonance terms. It is well known that resonance terms can produce the secular terms and it satisfy the solvability condition

$\left(r+is\right)S{T}_1=0$ (32)

A particular solution of Equation (27) is given by

$z_1^{*}\left(T_0,T_1,T_2\right)=\Phi \left(T_1,T_2\right){\text{e}}^{i{\omega }_0{T}_0}=\left(\begin{array}{l}{\varphi }_1\left(T_1,T_2\right)\hfill \\ {\varphi }_2\left(T_1,T_2\right)\hfill \end{array}\right){\text{e}}^{i{\omega }_0{T}_0}$ (33)

which satisfies

$\left(M+iN\right)z_1^{*}\left(T_0,T_1,T_2\right)=S{T}_1$ (34)

Therefore, solution is calculated as

$z_1\left(T_0,T_1,T_2\right)=\left(\begin{array}{l}{\gamma }_1{\Gamma }_1\hfill \\ {\gamma }_2{\Gamma }_2\hfill \end{array}\right)B^2\left(T_1,T_2\right){\text{e}}^{2i{\omega }_0{T}_0}+\left(\begin{array}{l}{\gamma }_{g1}{\Gamma }_{g1}\hfill \\ {\gamma }_{g2}{\Gamma }_{g2}\hfill \end{array}\right)B\left(T_1,T_2\right)\bar{B}\left(T_1,T_2\right)+cc$ (35)

with

$\begin{array}{l}{\Gamma }_1=-\frac{{\text{e}}^{2i{\omega }_0\left({\tau }_s+{\tau }_N\right)}\left(2i{\omega }_0{\gamma }_1+{\gamma }_N{\gamma }_1+{\gamma }_2{a}_{12}\right)-{\text{e}}^{2i{\omega }_0{\tau }_s}{b}_3{\gamma }_1}{{\gamma }_1{\Delta }_1}\\ {\Gamma }_2=-\frac{{\text{e}}^{2i{\omega }_0{\tau }_s}{\gamma }_1{b}_2+{\text{e}}^{2i{\omega }_0\left({\tau }_s+{\tau }_N\right)}{\gamma }_2\left(2i{\omega }_0-a_{11}\right)-{\text{e}}^{2i{\omega }_0{\tau }_N}{\gamma }_2{b}_1}{{\gamma }_2{\Delta }_1},\\ {\Gamma }_{g1}=\frac{{\gamma }_N{\gamma }_{g1}+a_{12}{\gamma }_{g2}-b_3{\gamma }_{g1}}{{\gamma }_{g1}{\Delta }_2},\\ {\Gamma }_{g2}=-\frac{a_{11}{\gamma }_{g2}+b_1{\gamma }_{g2}-b_2{\gamma }_{g1}}{{\gamma }_{g2}{\Delta }_2}\end{array}$ (36)

$\begin{array}{l}{\Delta }_1=\left(-2i{\omega }_0{\gamma }_N+2i{a}_{11}{\omega }_0+a_{11}{\gamma }_N+4{\omega }_0^2\right){\text{e}}^{2i{\omega }_0\left({\tau }_s+{\tau }_N\right)}\\ +\left(2i{\omega }_0{b}_3-a_{11}{b}_3+a_{12}{b}_2\right){\text{e}}^{2i{\omega }_0{\tau }_s}+\left(b_1{\gamma }_N+2i{b}_1{\omega }_0\right){\text{e}}^{2i{\omega }_0{\tau }_N}-b_1{b}_3\\ {\Delta }_2=-a_{12}{b}_2-\left(a_{11}+b_1\right){\gamma }_N+a_{11}{b}_3+b_1{b}_3\end{array}$ (37)

The above process is repeated in the following analysis. Considering Equation (23) and Equation (24), we have

$\frac{\partial z_2}{\partial T_0}\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)-\left(L\left(\mu \right)+R_1{\text{e}}^{i{\omega }_0{\tau }_s}+R_2{\text{e}}^{i{\omega }_0{\tau }_N}\right)z_2\left(T_0\mathrm{,}{T}_1\mathrm{,}{T}_2\right)=S{T}_2+NS{T}_2$ (38)

where

$\begin{array}{c}S{T}_2=\left(\begin{array}{c}-\left(p_1+i{q}_1\right)\left(1+b_1{\tau }_s{\text{e}}^{-i{\omega }_0{\tau }_s}\right)\\ -\left(p_1+i{q}_1\right)b_2{\tau }_N{\text{e}}^{-i{\omega }_0{\tau }_N}-\left(p_2+i{q}_2\right)\left(1+b_3{\tau }_N{\text{e}}^{-i{\omega }_0{\tau }_N}\right)\end{array}\right)\\ \frac{\partial B}{\partial T_2}\left(T_1,T_2\right)\exp \left(i{\omega }_0{T}_0\right)\\ -i{\omega }_0\left(p_1+i{q}_1\right)\left(\begin{array}{l}{b}_1\hfill \\ 0\hfill \end{array}\right){\tau }_{ϵ}{\text{e}}^{-i{\omega }_0{\tau }_s}B\left(T_1,T_2\right)\exp \left(i{\omega }_0{T}_0\right)\\ +\left(\begin{array}{c}{\Delta }_3\\ {\Delta }_4\end{array}\right)\exp \left(i{\omega }_0{T}_0\right)B^2\left(T_1,T_2\right)\bar{B}\left(T_1,T_2\right)\end{array}$ (39)

with

$\begin{array}{c}{\Delta }_3=(\left(p_1-i{q}_1\right){\left(p_1+i{q}_1\right)}^2\left({\beta }^{‴}\left(Q^{*}\right)Q^{*}+3{\beta }^{″}\left(Q^{*}\right)\right)\\ +2{\beta }^{″}\left(Q^{*}\right)Q^{*}\left({\gamma }_1{\Gamma }_1+2\Re \left\{{\gamma }_{g1}{\Gamma }_{g1}\right\}\right)\left(p_1+i{q}_1\right)){\text{e}}^{-r_s{\tau }_s}{\text{e}}^{i{\omega }_0{\tau }_s}\\ +(\frac{1}{2}\left(p_2+i{q}_2\right)\left(i{q}_2{p}_1-3q_2{q}_1-3p_1{p}_2-i{q}_1{p}_2\right)-2Q^{*}\Re \left\{{\gamma }_{g2}{\Gamma }_{g2}\right\}\left(p_2+i{q}_2\right)\\ \begin{array}{c}\\ \end{array}+Q^{*}{\gamma }_2{\Gamma }_2\left(-p_2+i{q}_2\right)){k^{″}}_N\left(N^{*}\right)+\frac{1}{2}{Q}^{*}\left(i{q}_2-p_2\right){\left(p_2+i{q}_2\right)}^2{k^{‴}}_N\left(N^{*}\right)\\ +\frac{1}{2}{Q}^{*}\left(i{q}_1-p_1\right){\left(p_1+i{q}_1\right)}^2{\beta }^{‴}\left(Q^{*}\right)+(\frac{3}{2}\left(i{q}_1-p_1\right){\left(p_1+i{q}_1\right)}^2\\ \begin{array}{c}\\ \end{array}-Q^{*}{\gamma }_1{\Gamma }_1\left(p_1-i{q}_1\right)-2Q^{*}\Re \left\{{\gamma }_{g1}{\Gamma }_{g1}\right\}\left(p_1+i{q}_1\right)){\beta }^{″}(Q*)\end{array}$

$\begin{array}{c}{\Delta }_4={\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}{\text{e}}^{-i{\omega }_0{\tau }_N}(\frac{1}{2}\left(p_2+i{q}_2\right)\left(6q_2{q}_1-2i{q}_2{p}_1+6p_1{p}_2+2i{q}_1{p}_2\right)\\ +Q^{*}\left(2\Re \left\{{\gamma }_{g2}{\Gamma }_{g2}\right\}\left(p_2+i{q}_2\right)+\left(p_2-i{q}_2\right){\gamma }_2{\Gamma }_2\right){k^{″}}_N\left(N^{*}\right)\\ +Q^{*}(2\Re \left\{{\gamma }_{g1}{\Gamma }_{g1}\right\}\left(p_2+i{q}_2\right)+2\Re \left\{{\gamma }_{g2}{\Gamma }_{g2}\right\}\left(p_1+i{q}_1\right)+{\gamma }_1{\Gamma }_1\left(p_2-i{q}_2\right)\\ +{\gamma }_2{\Gamma }_2\left(p_1-i{q}_1\right)){k^{\prime }}_N\left(N^{*}\right)+\frac{1}{2}{Q}^{*}\left(p_2+i{q}_2\right)\left(p_2^2+q_2^2\right){k^{‴}}_N\left(N^{*}\right))\end{array}$ (40)

There is also a solvability condition for Equation (37) as

$\left(r+is\right)\cdot S{T}_2=0$ (41)

and one obtain the complex valued norm form as

$\begin{array}{l}{B}^{\prime }=\frac{P_1+i{P}_2}{{\Lambda }_3+i{\Lambda }_4}B+\frac{{\Lambda }_1+i{\Lambda }_2}{{\Lambda }_3+i{\Lambda }_4}{B}^2\bar{B}\hfill \end{array}$ (42)

where

$\begin{array}{l}{P}_1=\Re \left\{i{\omega }_0{b}_1{\tau }_{ϵ}\left(r_1+i{s}_1\right)\exp \left(-i{\omega }_0{\tau }_{sc}\right)\right\},\\ P_2=\Im \left\{i{\omega }_0{b}_1{\tau }_{ϵ}\left(r_1+i{s}_1\right)\exp \left(-i{\omega }_0{\tau }_{sc}\right)\right\},\\ {\Lambda }_1=\Re \left\{-\left(r_1+i{s}_1\right){\Delta }_3-\left(r_2+i{s}_2\right){\Delta }_4\right\},\\ {\Lambda }_2=\Im \left\{-\left(r_1+i{s}_1\right){\Delta }_3-\left(r_2+i{s}_2\right){\Delta }_4\right\},\end{array}$

$\begin{array}{c}{\Lambda }_3=-r_1{p}_1-r_1{p}_1{b}_1{\tau }_s\cos \left({\omega }_0{\tau }_s\right)-b_2{\tau }_N{r}_2{p}_1\cos \left({\omega }_0{\tau }_N\right)-r_2{p}_2{b}_3{\tau }_N\cos \left({\omega }_0{\tau }_N\right)\\ -r_2{p}_2-r_1{q}_1{b}_1{\tau }_s\sin \left({\omega }_0{\tau }_s\right)-s_1{p}_1{b}_1{\tau }_s\sin \left({\omega }_0{\tau }_s\right)+s_1{q}_1{b}_1{\tau }_s\cos \left({\omega }_0{\tau }_s\right)\\ -b_2{\tau }_N{r}_2{q}_1\sin \left({\omega }_0{\tau }_N\right)-b_2{\tau }_N{s}_2{p}_1\sin \left({\omega }_0{\tau }_N\right)+b_2{\tau }_N{s}_2{q}_1\cos \left({\omega }_0{\tau }_N\right)\\ -r_2{q}_2{b}_3{\tau }_N\sin \left({\omega }_0{\tau }_N\right)-s_2{p}_2{b}_3{\tau }_N\sin \left({\omega }_0{\tau }_N\right)+s_2{q}_2{b}_3{\tau }_N\cos \left({\omega }_0{\tau }_N\right)\\ +s_2{q}_2+s_1{q}_1,\end{array}$

$\begin{array}{c}{\Lambda }_4=r_1{p}_1{b}_1{\tau }_s\sin \left({\omega }_0{\tau }_s\right)-r_1{q}_1{b}_1{\tau }_s\cos \left({\omega }_0{\tau }_s\right)-s_1{p}_1{b}_1{\tau }_s\cos \left({\omega }_0{\tau }_s\right)\\ +b_2{\tau }_N{r}_2{p}_1\sin \left({\omega }_0{\tau }_N\right)-b_2{\tau }_N{r}_2{q}_1\cos ({\omega }_0{\tau }_N)-b_2{\tau }_N{s}_2{p}_1\cos \left({\omega }_0{\tau }_N\right)\\ +r_2{p}_2{b}_3{\tau }_N\sin \left({\omega }_0{\tau }_N\right)-r_2{q}_2{b}_3{\tau }_N\cos \left({\omega }_0{\tau }_N\right)-s_2{p}_2{b}_3{\tau }_N\cos \left({\omega }_0{\tau }_N\right)\\ -r_2{q}_2-s_2{p}_2-s_1{q}_1{b}_1{\tau }_s\sin \left({\omega }_0{\tau }_s\right)-b_2{\tau }_N{s}_2{q}_1\sin \left({\omega }_0{\tau }_N\right)\\ -s_2{q}_2{b}_3{\tau }_N\sin \left({\omega }_0{\tau }_N\right)-r_1{q}_1-s_1{p}_1\end{array}$ (43)

Finally, doing the nonlinear transform

$B\left(T_1,T_2\right)=\frac{1}{2}\alpha \left(T_1,T_2\right)\exp \left(i\beta \left(T_1,T_2\right)\right)$ (44)

The first Lyapunov coefficient of the normal form is calculated [11] [13] [14] as

$l_0=\frac{{\Lambda }_1{\Lambda }_3+{\Lambda }_2{\Lambda }_4}{{\Lambda }_3^2+{\Lambda }_4^2}$ (45)

4. Bautin Bifurcation

Bautin bifurcation occurs if $l_0=0$ , as shown in Figure 4, near Bautin point, Super-critical Hopf bifurcation occurs at point 2 and subcritical Hopf bifurcation occurs at Point 1 respectively. As shown in Figure 5, in region (1),

Figure 4. Bautin bifurcation of periodic solutions determined by its first Lyapunov coefficient.

Figure 5. System dynamics near Bautin point.

periodic solution with big amplitude appears due to Hopf bifurcation which switches the stability character of the positive steady state. In region (2), the coexistence of both stable and unstable periodic solutions is inspired due to the Bautin bifurcation which further form the bistable dynamic regime and hysteresis. The long time evolution dynamical behavior jumps to oscillating periodic limit cycle and the steady state is stable as shown in region (3). Choose parameters in system (1) with values as shown in Table 1, We calculate lyapunov coefficient which is referred above to be vicinity at Bautin point with codimension 2. The suitable values of parameter $r_s$ and delay ${\tau }_s$ are chosen and the corresponding eigenvector are given as listed in Table 2.

Table 2. Parameters, eigenvalues and eigenvectors at Bautin point.

5. Discussion

In this work, we discussed the subcritical Hopf bifurcation of a hematopoiesis system which models two compartments of cell lines. The bistable regimes and hysteresis which are well known produced by subcritical Hopf bifurcation were detected near codimesion-2 bifurcation point (Bautin point).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Colijn, C. and Mackey, M.C. (2005) A Mathematical Model of Hematopoiesis: I. Periodic Chronic Myelogenous Leukemia. Journal of Theory Biology, 237, 117-132. https://doi.org/10.1016/j.jtbi.2005.03.033
[2]	Colijn, C. and Mackey, M.C. (2005) A Mathematical Model of Hematopoiesis: II. Cyclical Neutropenia. Journal of Theory Biology, 237, 133-146. https://doi.org/10.1016/j.jtbi.2005.03.034
[3]	Colijn, C., Fowler, A. and Mackey, M.C. (2006) High Frequency Spikes in Long Period Blood Cell Oscillations. Journal of Mathematical Biology, 53, 499-519. https://doi.org/10.1007/s00285-006-0027-9
[4]	Adimy, M., Crauste, F. and Ruan, S.G. (2006) Periodic Oscillations in Leukopoiesis Models with Two Delays. Journal of Theory Biology, 242, 288-299. https://doi.org/10.1016/j.jtbi.2006.02.020
[5]	Colijn, C. and Mackey, M.C. (2007) Bifurcation and Bistability in a Model of Hematopoietic Regulation. SIAM Journal of Applied Dynamical Systems, 6, 378-394. https://doi.org/10.1137/050640072
[6]	Zhuge, C.J., Lei, J.Z. and Mackey, M.C. (2012) Neutrophil Dynamics in Response to Chemotherapy and G-SCF. Journal of Theory Biology, 293, 111-120. https://doi.org/10.1016/j.jtbi.2011.10.017
[7]	Brooks, G., Langlois, G.P., Lei, J.Z. and Mackey, M.C. (2012) Neutrophil Dynamics after Chemotherapy and G-CSF: The Role of Pharmacokinetics in Shaping the Response. Journal of Theory Biology, 315, 97-109. https://doi.org/10.1016/j.jtbi.2012.08.028
[8]	Lei, J.Z. and Mackey, M.C. (2014) Understanding and Treating Cytopenia through Mathematical Modeling. In: Corey, S., Kimmel, M. and Leonard, J., Eds., A Systems Biology Approach to Blood. Advances in Experimental Medicine and Biology, Volume 844, Springer, New York, NY, 279-302.
[9]	Foley, C. and Mackey, M.C. (2009) Mathematical Model for G-CSF Administration after Chemotherapy. Journal of Theory Biology, 257, 27-44. https://doi.org/10.1016/j.jtbi.2008.09.043
[10]	Fredrickson-Hemsing, L., Ji, S. and Bozovic, D. (2012) Mode-Locking Dynamics of Hair Cells of the Inner Ear. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 86, 021915. https://doi.org/10.1103/PhysRevE.86.021915
[11]	Strogatz, S.H. (1994) Nonlinear Dynamics and Chaos. Westview Press, Cambridge.
[12]	Shlomovitz, R., Fredrickson-Hemsing, L., Bozovic, D., et al. (2013) Low Frequency Entrainment of Oscillatory Bursts in Hair Cells. Journal of Biophysics, 104, 1661-1669.
[13]	Zhang, S. and Xu, J. (2013) Quasiperiodic Motion Induced by Heterogeneous Delays in a Simplified Internet Congestion Control Model. Nonlinear Analysis: Real World Applications, 14, 661-670.
[14]	Zhang, S. and Xu, J. (2011) Time-Varying Delayed Feedback Control for an Internet Congestion Control Model. Discrete and Continuous Dynamical Systems Series B, 653-665.
[15]	Li, C.G., Chen, G.R., Liao, X.F. and Yu, J.B. (2004) Hopf Bifurcation in an Internet Congestion Control Model. Chaos, Solitons and Fractals, 19, 853-862.
[16]	Hale, J. (2003) Theory of Functional Differential Equations. World Publishing Corporation, Beijing.
[17]	Ma, S.Q., Wang, X.H., Lei, J.H. and Feng, Z.S. (2010) Dynamics of the Delay Hematological Cell Model. International Journal of Biomathematics, 3, 105-125. https://doi.org/10.1142/S1793524510000829
[18]	Bernard, S., Bélair, J. and Mackey, M.C. (2003) Oscillations in Cyclical Neutropenia: New Evidence Based on Mathematical Modeling. Journal of Theory Biology, 223, 293-298.
[19]	Santillan, M., Mahaffy, J.M., Bélair, J. and Mackey, M.C. (2000) Regulation of Platelet Production: The Normal Response to Perturbation and Cyclical Platelet Disease. Journal of Theory Biology, 206, 585-603. https://doi.org/10.1006/jtbi.2000.2149
[20]	Haurie, C., Dale, D.C., Rudnicki, R. and Mackey, M.C. (2000) Modeling Complex Neutrophil Dynamics in the Grey Collie. Journal of Theory Biology, 24, 505-519. https://doi.org/10.1006/jtbi.2000.2034
[21]	Hearm, T., Haurie, C. and Mackey, M.C. (1998) Cyclical Neutropenia and the Peripherial Control of White Blood Cell Production. Journal of Theory Biology, 192, 167-181. https://doi.org/10.1006/jtbi.1997.0589