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

$\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

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.

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)

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)