Dynamics on semi-discrete Mackey-Glass model

1.   Introduction

Red blood cells (RBCs) are derived from stem cells in the bone marrow, and their main function is to distribute oxygen in the human body. In a healthy human body, RBCs tend to a stable state, and if their numbers suddenly decrease, it will stimulate the stem cells in the bone marrow to produce significant excitement, thereby generating new cells that return to a stable state. There must be a period of time between the moment RBCs decline and the moment they recover. Under normal circumstances, the process of forming new cells takes approximately 5 to 7 days ^[1].

Currently, some diseases are known to cause the number of RBCs to gradually oscillate, which is represented by a periodic function. If RBCs are insufficient for a period of time, the human body will be unable to recover and die. To check the health status of the body, it can be achieved by detecting RBCs. In 1976, Lasota et al. ^[24] proposed the Lasota hematopoietic model, which perfectly described the changes in RBCs in the human body. Later, Mackey and Glass ^[1] combined the respiratory system to further characterize RBCs in the blood.

In 1976, Mackey and Glass introduced the following hetmatopoietic system

where $\delta, p\in(0, +\infty)$, $k, m\in \mathbb{N}$. $\xi(t)$ denotes the population of mature cells at time $t$, and the cells disappear from cycle with the velocity $\delta$, and the delay $\tau > 0$ denotes the time of generating new RBCs.

EI Sheikh et al. ^[4] provided some sufficient and necessary conditions for the oscillation of all positive solutions when $k = 0$, while Gopalsamy ^[17] provided some sufficient and necessary conditions for the oscillation of all positive solutions and global attractiveness when $k = 0$ or $1$. Through numerical simulations of the system, Hela and Sternberg ^[18] found that the degree of chaos increases with the increase of time delay $\tau$. System (1.1) with $k = 1$ attracted many authors' attentions ^{[3,12,28,29,36]}. While for the case $k = 0$ and $k > 1$, there were less researchers concerned. Therefore, the paper deals with system (1.1) with $k = 0$ and $k > 1$.

Hematopoietic models are usually divided into two categories: Continuous models described by differential equations, and discrete models obtained by discretizing continuous systems. These mathematical methods include the Euler difference method, Adams Bashforth method ^[49], piecewise-constant approximation method ^[22], the conformable derivatives ^[21], non-standard finite difference scheme (NSFDS) of Mickens-type ^[32], and semi discrete method ^[17,38,39]. The main methods are analyzed by using the Euler difference method and semi discrete method.

In recent years, the advantages of discrete models have become increasingly apparent in various fields. To begin in scientific research, data is collected over a period of time, such as year, month, day, etc. Therefore, considering a discrete hematopoietic model is more in line with practical significance ^[38,39].

Moreover, due to the fact that medical tests using RBC counting are often conducted at regular intervals, and changes in the number of RBCs are difficult to observe continuously over a long period of time, the discrete case is considered to be more appropriate to describle the real world sometimes, and it is obtained in the literatures that discrete systems have richer dynamical behaviors. Various of discrete-time models have been studied in detail recently ^{[9,19,25,31,41,48,50]}. There are a series of results on the qualitativeness and bifurcations of discrete-time models like flip bifurcation and Neimark-Sacker bifurcation ^{[8,10,26,31,34,41,47]}.

In 2002, Tamas and Gabor ^[13] proposed a semi-discretization method to study the delayed system, which is a simple and effective way to deal with the delayed term. Then, by using this method some new semi-discrete models have been proposed and analyzed. For example, in ^[39], the authors proposed a semi-discrete hematopoiesis model and discussed its dynamical behavior which included the stabilities of the fixed points and the Neimark-Sacker bifurcation. For other studies for semi-discrete models, one may see ^[9,20] and the references therein. Moreover, Yao and Li ^[40] provided bifurcation difference induced by different discrete methods in a discrete predator-prey model, which shows that it is meaningful to investigate the model with different discretization methods.

When analyzing the semi-discrete Mackey-Glass model, we need to consider the situation where multiple parameters change simultaneously. Complicated bifurcations likely occur when more than one systemic parameter is varied at the same time ^[6]. Codimension-2 bifurcations in the discrete-time systems will occur when the dimension of the center manifold is changed by the approximation of extra multipliers to the unit circle, or some of the nondegeneracy conditions for the one-parameter bifurcations are violated ^[44,45]. Kuznetzov derived explicit formulas for the normal form coefficients to verify the nondegeneracy of codim-2 bifurcations of fixed points with 1 or 2 critical multipliers ^[45]. Li and He considered 1:2 and 1:4 resonances in a discrete-time Hindmarsh-Rose model ^[6], and Wu and Zhao derived codimension-two bifurcations with 1:2 resonance in a discrete-time predator-prey model ^[11]. Rana and Uddin analytically showed a flip-Neimark Sacker bifurcation of the discretized L$\ddot{\text{u}}$ system ^[35]. More types of complicated bifurcations of discrete-time or continuous-time systems have been studied in ^{[2,5,6,7,11,14,15,27,30,33,35,42,43,44,45,46]}. Moreover, from the perspective of numerical convergence, discrete systems can reflect the dynamical behaviors of the original continuous system and generate richer dynamical phenomena ^{[9,25,38,39,41,48,50]}.

In this paper, the semi-discrete hetmatopoietic Mackey-Glass model is analyzed. We discuss the dynamical behaviors of system (1.1), including the existence of fixed points, local stability, and codim 1 and 2 bifurcations. We have found that with the parameters varying, there exist invariant curves when the system is undergoing Neimark-Sacker bifurcation, period-4 saddle points, and heteroclinic cycle composed by the separatrices of them when undergoing 1:4 resonance. The existence of Li-Yorke chaos is proved in this case. In the case of $2 \leq k \leq m$, the system undergoes saddle-node bifurcation. In order to illustrate the correctness of theoretical analysis, numerical simulations are provided and some other complex dynamical behaviors are reviewed.

The rest of paper is organized as follows: In Section 2, we discussed the existences and stabilities of fixed points of the M-G model with $k = 0$ and theoretically verified the existence of Neimark-Sacker bifurcations and 1:4 resonances. In Section 3, these theoretical results are supported by the numerical simulations. Then, we consider the case of $1 < k \leq m$ in Section 4. We analyzed the existences and stabilities of fixed points and verified the existence of saddle-node bifurcations, Neimark-Sacker bifurcations, and 1:4 resonances, which are more complex than the case in Section 2. In Section 5, the numerical simulations are presented to support the theoretical results and show the impact of negative feedback, mixed feedback, and positive feedback depending on the parameter $k$. Finally, a brief conclusion is organized in Section 6 which provides more directions for future works.

2.   M-G model with $k=0$

In this section we discuss system (1.1) with $k = 0$ as follows:

where $\delta, \tau, p \in (0, +\infty)$ and $m \in \mathbb{N}.$ Next, we establish the semi-discrete model for system (2.1). For the sake of simplicity, intoducing the transformations $s = \frac{t}{h}$ and $N(t) = N(s\tau) = \eta(s)$, (2.1) is then changed to

For simplicity, resetting $\delta = \delta\tau$, $p = p\tau$ in the above equation, one gets

Through this transformation, the reduction of the time-delay parameter $\tau$ in the system is completed, which makes the problem simple.

Let $[s]$ denote the integer that is not bigger than $s$, and consider the semi-discrete model of (2.3)

It is easy to see that the right side of the Eq (2.4) is constant on the interval $[n, n+1)$. Obviously, the following conclusions about Eq (2.4) are true.

Lemma 2.1. The solution $\eta(s)$ of Eq (2.4) satisfies

1. $\eta(s)$ is continuous on $[0, +\infty)$,

2. $\frac{\operatorname{d}\eta}{\operatorname{d}s}$ exists on $\bigcup\limits_{s = 0}^{+\infty}(s, s+1)$,

3. Equation (2.4) is true on every interval $[m, m+1)$ for $m = 0, 1, 2, \dots$.

For any $s\in[n, n+1)$, we have $[s] = n$. Substituting it into Eq (2.4) we get

and then integrate the Eq (2.5) from $n$ to $s$ for any $s\in[n, n+1)$ where $n\in\{ 0, 1, 2, 3, \dots\}$. The following difference equation can be found for $s \in [n, n+1)$

Let $s \to (n+1)^{-}$ in the above equation, and then the semi-discrete system with no delay for system (2.1) will be obtained as follows:

which is a second-order difference equation. Take the normal transformation

and we may achieve the discrete planar dynmaical system of (2.1) as follows:

where $\delta > 0$, $p > 0$, $x_0, y_0 \in (0, +\infty)$. We shall discuss the dynamics on system (2.6) in detail.

2.1.   Local dynamics for the fixed points

At the beginning of the discussion, we consider the existence of the fixed point of (2.1).

Theorem 2.1. The conclusions of the fixed point for system (2.6) are true:

${({{\rm{i}}})}$ there is a unique fixed point $E_*(x_*, y_*)$ of system (2.6) satisfying

${({\rm{ii}})}$ $x_*$ is strictly increasing at $p$.

Proof. It is easy to see that the fixed point of system (2.6) satisfies the following equations

Next, we show the existence and uniqueness of the fixed point. Set

so we may study its positive zero points.

Notice that

and

Hence, there is a unique zero point of $F(x, p)$ in $D_x$, and $0 < x_* < \frac{p}{\delta}$.

$({\rm{ii}})$ Since $F_x^\prime(x, p) > 0$, the implicit function theorem shows that there exists a function $x_* = f(p)$ in the field $D_{xp} = \{(x, p) \big| 0 < x_* < \frac{p}{\delta}, p > 0\}$ and

which means $x_*$ is strictly increasing in $p$. Next, we deal with the stability of the fixed point of system (2.6). The Jacobian matrix of (2.6) at $E_*$ is

then $\operatorname{Tr}(J_{E_*}) = 1-\delta$, $\operatorname{Det}(J_{E_*}) = \frac{\delta m x^m_*}{1 + x^m_*}$, and the corresponding characteristic equation is

hence,

Lemma 2.2. ^[39] Consider $F(\lambda) = \lambda^2+B\lambda+C$ with two constant real parameters $B$ and $C$. Suppose $\lambda_1$ and $\lambda_2$ are two zero points of $F(\lambda)$. Then, the zero points of the equation are entirely determined by both $F(-1)$ and the parameters.

Case I. If $F(1) > 0$, then

(1) $|\lambda_1| < 1$ and $|\lambda_2| < 1$ if, and only if, $F(-1) > 0$ and $C < 1$.

(2) $|\lambda_1| > 1$ and $|\lambda_2| > 1$ if, and only if, $F(-1) > 0$ and $C > 1$.

(3) $|\lambda_1| < 1$ and $|\lambda_2| > 1$ if, and only if, $F(-1) < 0$.

(4) There exist $\lambda_1 = -1$ and $\lambda_2\neq-1$ if, and only if, $F(-1) = 0$ and $B\neq2$.

(5) $\lambda_1 = \lambda_2 = -1$ if, and only if, $F(-1) = 0$ and $B = 2$.

(6) There exist a pair of conjugate complex roots with $|\lambda_1| = |\lambda_2| = 1$ if, and only if, $-2 < B < 2$ and $C = 1$.

Case II. If $F(1) = 0$, i.e., $\lambda_1 = 1$, then $|\lambda_2| > 1(resp., < 1)$ if, and only if, $|C| > 1$ $(resp., < 1).$

Case III. If $F(1) < 0$, then $F(\lambda) = 0$ has one root $\lambda_1 \in(1, +\infty)$. Then, the following statements about the other root $\lambda_2$ hold.

(1) $\lambda_2 < ( = )-1$ if, and only if, $F(-1) < ( = )0$;

(2) $-1 < \lambda_2 < 1$ if, and only if, $F(-1) > 0$.

Therefore, the following conclusions can be obtained.

Theorem 2.2. Let $p_{_{NS}} = \frac{\delta^2 m}{\delta m - 1} \big(\frac{1}{\delta m-1}\big)^{\frac{1}{m}}$, then

${\rm{(i)}}$the fixed point $E_*$ of system (2.6) is a sink for $\delta m \le 1$ or $\delta m > 1$ and $p < p_{_{NS}}$;

${\rm{(ii)}}$ the fixed point $E_*$ of system (2.6) is a source for $\delta m > 1$ and $p > p_{_{NS}}$;

${\rm{(iii)}}$ system (2.6) may undergo Neimark-Sacker bifurcation at the fixed point $E_*$ for $\delta m > 1$ and $p = p_{_{NS}}$.

Proof. Notice that $\operatorname{Det}(J_{E_*}) = \frac{\delta m x^m_*}{1 + x^m_*}$; it is easy to derive that $\operatorname{Det}(J_{E_*}) < 1$ holds for $\delta m \le 1$.

Meanwhile, for $\delta m > 1$ and $p = p_{_{NS}}$, we have $x_* = \big(\frac{1}{\delta m -1}\big)^{\frac{1}{m}}$, which leads to $\operatorname{Det}(J_{E_*}) = 1$. Since $F(-1) > 0$ and $F(1) > 0$, by Lemma 1, one may know that $\lambda_{1, 2}$ are a pair of conjugate complex roots; moreover, $|\lambda_{1, 2}| = 1$. Hence, system (2.6) probably undergoes Neimark-Sacker bifurcation at $E_*$ (we will give the theoretical proof in the next subsection).

Furthermore, we may see that $\operatorname{Det}(J_{E_*})$ is strictly increasing in $x_*$ in the field $D_x$. Combining the results in Theorem 2.1 that $x_*$ is strictly increasing in $p$, which in turn leads to the fact $\operatorname{Det}(J_{E_*})$ is strictly increasing in $p$, so for $\delta m > 1$ and $p < p_{_{NS}}$, one gets $\operatorname{Det}(J_{E_*}) < 1$ and $|\lambda_{1, 2}| < 1$, which implies the fixed point $E_*$ of system (2.6) is a sink.

Whereas for $\delta m > 1$ and $p > p_{_{NS}}$, we have $\operatorname{Det}(J_{E_*}) > 1$ and $|\lambda_{1, 2}| > 1$, which shows that the fixed point $E_*$ of system (2.6) is a source.

The proof is then completed.

2.2.   Bifurcation analysis for M-G model with $k = 0$

2.2.1.   Neimark-Sacker bifurcation

In this sbusection, we shall show the existence of Neimark-Sacker bifurcation and the stability of the bifurcated invariant curve.

To start, we choose $p$ as a bifurcation parameter, giving a perturbation $p_*$ of parameter $p_{NS}$, then we obtain the following pertubed system:

where $|p_* | \ll 1$.

Next, taking $u = x - x_*, v = y - y_*$, then the fixed point $E_*$ can be moved to the origin $O(0, 0)$ and system (2.7) becomes

The characteristic function of the linearization of system (2.8) at its fixed point $O$ is:

where $a(p_*) = 1-\delta$, $b(p_*) = \frac{m (p_{_{NS}}+p_*)x^{m-1}_*}{(1 + x^m_*)^2} = \frac{\delta m x^m_*(p)}{1+x^m_*(p)}$. It is easy to see

is a pair of conjugate complex roots of $F(\lambda) = 0$. As we know, to guarantee the existence of Neimark-Sacker bifurcation, the following conditions must be fufilled:

Notice that $a(p_*)\big |_{p_* = 0} = 1-\delta$ and $b(p_*) = \frac{m (p_{_{NS}}+p_*)x^{m-1}_*}{(1 + x^m_*)^2} = \frac{\delta m x^m_*(p)}{1+x^m_*(p)}$, so

which means the condition (${\mathbf{C.2}}$) is true. Furthermore, $|\lambda_{1, 2}| = \sqrt{b(p_*)} = \sqrt{\operatorname{Det}(J_{E_*})}$,

which guarantees that the condition $({\mathbf{C.1}})$ holds. Therefore, system (2.6) undergoes Neimark-Sacker Bifrucation as $p_*$ varies in the neighborhood of $p_{NS}$.

In the following, we shall discuss the stability of the invariant curve by three steps.

The first step: Expand the right side of system (2.8) as Taylor series at $(u, v) = (0, 0)$

where $\rho = \sqrt{\|u\|^2+\|v\|^2}$,

The second step: Take the matrix

and it is not difficult to see

introduce the invertable transformation $(u, v)^T = T(X, Y)^T$, and the normal form of system (2.9) is arrived as:

where $F(X, Y) = \frac{2b_{20}}{\sqrt{3+2\delta-\delta^2}}Y^2 + \frac{2b_{30}}{\sqrt{3+2\delta-\delta^2}}Y^3,$ $G(X, Y) = 0, \rho = \sqrt{|X|^2+|Y|^2}.$ By computation, we get

The last step is to compute the first Lyapunuov coefficient by

where

Therefore, by computation and simplification, one may achieve

Next, we discuss the sign of $a^*$. Denote $h(m) = -(\delta^3m-\delta^3+\delta^2m^2-2\delta^2-2\delta m +2)$,

Obviously, $h^\prime(m) < 0$ for $\delta m > 1$. Therefore, $h(m)\leq h(2) < 0$ for $\delta m > 1$ (in fact, $\delta m > 1$ and $m\in\mathbb{N}^+$ means $m\geq 2$). This leads to $a^* < 0$ directly. Hence, the following theorem can be obtained.

Theorem 2.3. For $\delta m > 1$, the first Lyapunov coefficient $a^* < 0$ is always true, which means when $p$ is varying in the sufficient small neighborhood of $P_{NS}$, system (2.6) undergoes Neimark-Sacker bifurcation at the fixed point $E_*$; moreover, the bifurcated invariant curve is stable.

2.2.2.   Resonance 1:4 for the fixed point $E_{*}$

Now we choose $p$ and $\delta$ as the bifurcation parameters. Since $\lambda_{1, 2} = \pm \operatorname{i}$ when $p = p_{SN}$ and $\delta = \delta_4$ $= 1$, give a perturbation $p_*$ and $\delta_{*}$ of parameter $p_{NS}$ and $\delta_4 = 1$. In order to ensure that the value of $p_{NS}$ is meaningful, $m$ is not equal to 1 in this case. Then, we obtain the following pertubed system:

where $p = p_{NS} + p_*$, $\delta = \delta_4 +\delta_*$, $|p_* |, |\delta_*| \ll 1$.

Next, taking $u = x - x_*, v = y - y_*$, then the fixed point $E_*$ can be moved to the origin $O(0, 0)$ and system (2.13) becomes:

The Jacobian of system (2.14) at its fixed point $O$ is:

System (2.14) is expanded as the Taylor series at $(u, v) = (0, 0)$

where $\rho = \sqrt{\|u\|^2+\|v\|^2}$,

If $p_* = \delta_* = 0$, the Jacobian will be equivalent to

The characteristic value of the linearization of system (2.14) at its fixed point $O$ is

and the eigenvectors are

The map (2.15) can be written as

where

To analyze the bifurcation, we take any vector $(u, v)^T \in \mathbb{R}^2$ in the form $(u, v)^T = z q_1 + \bar{z} q_2$, where $z$ is a complex variable which transforms the map (2.16) to the complex form

where

with

Take an invertible parameter-dependent change of complex coordinate

with

Then, eliminate all the quadratic terms in (2.17), and transform it to

with

Then, by setting

we take another invertible parameter-dependent change of complex coordinate

which changes the map (2.19) into

where $g_{kl}$ are given by (2.18) and

When $m \neq 5$ ($D \neq 0$), the fourth iterate of $\Gamma$ allows for approximation by a complex flow, and the bifurcation of 1:4 resonance will be determined by

and

We have

By $m \in \mathbb{N}_+$, we have

When dealing with strong resonances, we will repeatedly use the approximation of maps near their fixed points by shifts along the orbits of certain systems of autonomous ordinary differential equations ^[6,44].

Lemma 2.3. ^[44] The fourth iterate of the map (2.20) can be represented in the form

where $\varphi^t$ is the flow of a planar system

Theorem 2.4. If $m \neq 1, 2, 3, 5, 7$ (means that $p_{NS}$ is meaningful, $\operatorname{Re} \tilde{A} \neq 0$ and $\operatorname{Im} \tilde{A} \neq 0$), and the parameters $(p, \delta)$ vary in a small neighborhood of $(p_{SN}, 1)$, the system (2.6) undergoes a bifurcation of 1:2 resonance around the fixed point $E_{*}$, and admits the following bifurcation curves:

● There is a Neimark-Sacker bifurcation curve at the fixed point $E_{*}$ of the system (2.6).

● There exist four saddle fixed points $S_k$ when $|\tilde{A}| > 1$, which are the corresponding period-4 points of (2.21). As the parameters $(p, \delta)$ in (2.6) varying in the neighborhood of $(p_{NS}, 1)$, these nontrivial periodic points appear or disappear in the neighborhood of $E_*$.

● There is a square heteroclinic cycle around $E_{*}$ composed by the separatrices of $S_k$, which is stable from the inside.

3.   Numerical simulation for M-G model with $k=0$

3.1.   Neimark-Sacker bifurcation

Select $\delta = 0.2$, $m = 10$, and the initial point $E_0(1.5, 1.5)$. By the theoretical analysis in Section 2.1, one has $p_{NS} = 0.4$, and at this case the fixed point of system (2.6) is $E_*(1, 1)$, $a^* = -0.5687$. See Figures 1–3 by numerical simulation.

Figures 1(a)–(c) show that the fixed point $E_*$ is a sink when $p < 0.4$; meanwhile for $p = 0.41$, a stable invariant curve occurs, see the Figures 1(d), (f); it can be also shown in Figure 2 which is consistent with the results in Theorem 2.3.

The Figure 3 shows the maximum Lyapunov exponent (MLE) varying with the parameter $p$, which shows the MLE is always less than zero when $p < 0.4$; meanwhile for $\gamma > 0.4$, the MLE will fluctuate below 0, which also means the oribit from the initial $(1.5, 1.5)$ will be periodic. Also the MLE is not greater than 0, which implies the chaos does not occour for system (2.6) at this time.

3.2.   1:4 resonance

Taking $m = 4$ as an example, the same applies to other situations.

Let $\delta$ vary form the neighborhood of $\delta_4 = 1$, while $p_{NS} = 1.0131$, $x_*(p_{NS}) = 11^{-1/4}$ and let $p$ vary in $[p_{NS}-0.02, p_{NS}+0.2]$. Thus, it follows from Theorem 4.6 that the fixed point $E_*$ of system (2.6) is a 1:4 resonance point.

Figures 4(a), (b) show that there exists the invariant cycle around the fixed point when $\delta$ varies near $\delta = 1$. A nondegenerate Neimark-Sacker bifurcation curve when $\delta = 0.8$ is illustrated in Figure 4(c), and four period points when $\delta = 0.9$ and $\delta = 1.01$ are given in Figures 4(d), (e). Figures 4(b), (d) show that a complex heteroclinic curve will occur near the Neimark-Sacker bifurcation curve or period points when $\delta$ varies in the neighborhood of $\delta = 0.9$ or $1.01$. A possible chaotic phenomenon is shown in Figure 4f, and in fact, we have verified it through numerical simulations (see Figures 6(a) and 7(b)).

Figure 5 shows the phase portraits of system (2.6) corresponding to Figures 4(a)–(e). Since the system undergoes Neimark-Sacker bifurcation, the fixed point $E_*$ changes from local asymptotically stable to unstable and there occurs a stable closed invariant curve around it, which is shown in Figures 5(a), (b)

From Figures 5(d)–(f), we observe that a heteroclinic curve consists of the separatrices of four saddles; specially, Figure 5f shows the collision between the saddle points and the invariant circle. By selecting $\delta = 1.685$, $p = 2.01311$ (see Figures 6(a), (b)), we calculate that the max Lyapunov exponent is $L = 0.0715457 > 0$, and the system has 3-period points, which implies chaos ^[37].

The 3-dimensional bifurcation diagrams of map (2.6) in $(p, \delta, y)$ and the maximium Lyapunov exponents corresponding to $\delta$ and $p$ are shown in Figure 7(a). The maximum Lyapunov exponents corresponding to $p$ and $\delta$ are calculated and plotted in Figure 7(b) that confirms the existence of the period orbits and chaotic regions near the 1:4 resonance point $E_*$ in the parametric space.

4.   M-G model with $1 < k\leq m$

In this section, we study the Makey-Glass model with $1 < k\leq m$.

We shall see in this section, that the model with $1 < k\leq m$ exhibits more complex dynamics (saddle-node bifurcation and chaos) than the model with $k = 0$ and $k = 1$.

where $\delta, \tau, p \in (0, +\infty)$, $m, k \in \mathbb{N^+}, 1 < k \le m.$ With the same way one can derive the semi-discrete system of system (4.1):

Taking the transfomations $\begin{cases} x_n = \eta(n-1), \\ y_n = \eta(n), \end{cases}$ we can obtain the planar dynamical system

In the following, we discuss the dynamics of system (4.3). It is shown that the origin $O(0, 0)$ is always the fixed point of system (4.3), which is a trivial fixed point. With the increasing of parameter $p$, the number and the stability of the fixed point will change.

When the parameter $p$ meets the first critical value $p_0$, there will exist a positive fixed point $E_*$, and the system may undergo saddle-node bifurcation at $E_*$ for $p = p_0$. Meanwhile, for $p > p_0$, there are two fixed points $E_{1, 2}$ bifurcated from $E_*$.

4.1.   The local dynamics of the fixed points

Suppose $E(x, y)$ is the fixed point of system (4.3), and it must satisfy the following equations

Denote $f(x) = \delta x^{m+1} - px^k +\delta x$. Obvilusly if $\bar{x}$ is the root of $f(x) = 0$, then $(\bar{x}, \bar{x})$ is also the fixed point of system (4.3). It is easy to see that $x = 0$ is always a root of $f(x) = 0$, then $E_0(0, 0)$ is the trivial fixed point of (4.3).

To find the nonzero roots of $f(x) = 0$, take $g(x) = \frac{f(x)}{x} = \delta x^{m} - px^{k-1} +\delta$, and we may discuss the positive solutions of $g(x) = 0$. Notice that

Since $x^{m-k+1}$ is strictly increasing at $x$, $g(x)$ will get to its minimum value at $x = \sqrt[m-k+1]{\frac{p(k-1)}{m\delta}} \triangleq x_0$. Let $\rho_1 = \frac{p}{\delta}, \rho_2 = \frac{k-1}{m}$, then we have $\rho_1 > 1, \rho_2 < 1$; hence, $x_0 = (\rho_1\rho_2)^{\frac{1}{m-k+1}}$, and the minimum value of $g(x)$ is

Theorem 4.1. Let $p_{_{SN}} = \delta \rho_2^{-\rho_2}(1-\rho_2)^{\rho_2 - 1}$, $x_0 = (\rho_1 \rho_2)^{\frac{1}{m-k+1}}$, then the following conclusions about the fixed points of system (4.3) are true:

$({{\rm{i}}})$ when $p < p_{_{SN}}$, there is only one trivial fixed point $E_0(0, 0)$ of system (4.3);

$({\rm{ii}})$ when $p = p_{_{SN}}$, there are two fixed points of system (4.3), which are the trivial fixed point $E_0(0, 0)$ and the positive fixed point $E_*(x_*, y_*)$, $x_* = y_* = x_0$;

$({\rm{iii}})$ when $p > p_{_{SN}}$, there are three fixed points of system (4.3), which are $E_0(0, 0)$, $E_1(x_1, y_1), E_2(x_2, y_2)$, where $x_i = y_i$ statisfies $\delta x_i^{m+1} - px_i^k +\delta x_i = 0$, $i = 1, 2$.

Proof. It is not difficult to see that when $p = p_{_{SN}}$, $\rho_1 = \rho_2^{-\rho_2}(1-\rho_2)^{\rho_2 - 1}$, the minimum value of $g(x)$ satisfies

which shows there is only one solution $x_0$ of $g(x) = 0$.

Whereas for $p < p_{_{SN}}$, we have

that is to say, $g(x_0)$ is strictly decreasing at $p$ for $\rho_1 > 1$, $0 < \rho_2 < 1$. Hence, $g(x_0) > 0$ for $p < p_{_{SN}}$, which means there is no positive solution of $g(x) = 0$; as to $p > p_{_{SN}}$, we have $g(x_0) < 0$, combined with the fact

which means there are only two solutions $x_1, x_2$ for $g(x) = 0$ and $0 < x_1 < x_0 < x_2$.

Notice that $E_0(0, 0)$ is always the fixed point of system (4.3), hence, we can finish the proof.

It is easy to see from the above theorem that system (4.3) probably undergoes saddle-node bifurcation. In the next subsection we shall verify the existence of saddle-node bifurcation theoretically.

Theorem 4.2. Let $p > p_{_{SN}}$, $E_1(x_1, y_1)$, and $E_2(x_2, y_2)$ be the fixed points of system (4.3), then $x_1, x_2$, and the parameter $p$ satisfy that $x_i$ is strictly increasing with respect to $p$ in $(p_{_{SN}}, +\infty)$, $i = 1, 2.$

Proof. Denote $D_2 = \{(x_i, p)\ \big|\ x_i\in(x_0, +\infty), p\in(p_{_{SN}}, +\infty)\}$, and $x_i$ satisfies $g(x_i, p) = 0$. Notice that $g^\prime_p(x_i, p) = -x_i^{k-1}$, $g_{x_i}^\prime(x_i, p) = \delta mx_i^{m-1}-p(k-1)x_i^{k-2}$ are continuous in $D_2$, and

The implicit theorem shows that there exists a function $x_i = x_i(p)$ satisfying $g(x_i(p), p) = 0$ when $p > p_{_{SN}}$, and

then the proof is finished.

Next, we study the stabilities of the fixed points.

To begin, we discuss the stability of the fixed point $E_0(0, 0)$.

The Jocabian Matrix of system (4.3) at $E_0(0, 0)$ is

and the corresponding characteristic equation is

hence, $\lambda_1 = 0, \lambda_2 = 1-\delta < 1$, which means the fixed point is a sink.

Next, we study the stability of $E_*(x_*, y_*)$ for $p = p_{_{SN}}$. The Jacobian matrix of (4.3) at $E_*(x_*, y_*)$ is

so $\operatorname{Tr}(J_{E_*}) = 1-\delta$, $\operatorname{Det}(J_{E_*}) = -\frac{\delta(k-m)x^m_*+\delta k}{1 + x^m_*}$, and the corresponding characteristic equation is

It is not difficult to see

and $C = \operatorname{Det}(J_{E_*}) = -\delta$, which means $|C| < 1$; Lemma 2.2 implies that $\lambda_1 = 1$ and $|\lambda_2| < 1$, that is to say, $E_*$ is a nonhyperbolic fixed point.

Last, we discuss the stabilities of the fixed points $E_1(x_1, y_1)$ and $E_2(x_2, y_2)$ for $p > p_{_{SN}}$.

and the corresponding characteristic equation is

where $B = (\delta-1)$, $C = -\frac{\delta(k-m)x^m_2 + \delta k}{1 + x^m_2}$,

It is easy to see that $G(x_1)$ is increasing at $x_1$, while $x_1$ is decreasing at $p$; hence,

Combined with Theorem 4.2, we know that there is an eigenvalue $\lambda_1 > 1$.

Notice $F(-1) = 2-2\delta +G(x_1)$, which is decreasing at $p$, and

So when $2-\delta-\delta k < 0$, we have

such that $F(-1) = 0$. At this case, $x_1 = \Big(\frac{\delta k + \delta - 2}{\delta(m-k)+(2-\delta)}\Big)^{\frac{1}{m}}$, $\lambda_2 = -1$. Moreover, when $p_{_{SN}} < p < p_{change}$, $F(-1) > 0$, which implies $|\lambda_2| < 1$; while for $p > p_{change}$, one has $F(-1) < 0$, and then $|\lambda_2| > 1$.

Theorem 4.3. The conclusions of the fixed point $E_1$ for system (4.3) are true:

${({{\rm{i}}})}$ If $\delta k + \delta - 2 \le 0$, then the fixed point $E_1$ is a saddle;

${({\rm{ii}})}$ if $\delta k + \delta - 2 > 0$, then

$(a)$ when $p_{_{SN}} < p < p_{change}$, $E_1$ is a saddle;

$(b)$ when $p > p_{change}$, $E_1$ is a source, where $p_{change}$ is given below.

The theorem shows that $p_{change}$ is a critical value. When the parameter $p$ varies in the neighborhood, the stability of the fixed point will be changed. So, the system may undergo bifurcation at $E_1$, which will be discussed in the next section.

With similar arguments, we may know that the characteristic equation of the Jacobian matrix for system (4.3) at $E_2(x_2, y_2)$ is:

where $B = (\delta-1)$, $C = -\frac{\delta(k-m)x^m_2 + \delta k}{1 + x^m_2}$, denote

It is easy to verify that $G(x_2)$ is increasing at $x_2$, and by Theorem 4.2, $x_2$ is also increasing at $p$ for $p > p_{_{SN}}$. As a result, $G(x_2(p))$ is increasing at $p$ for $p > p_{_{SN}}$ which leads to

Notice that $F(-1) = 2-2\delta + F(1)$, then $F(-1) > 0$ holds. As we know, $C(p) = G(x_2(p))-\delta$ is increasing at $p$, moreover,

so there exists a point

for $\delta(m-k) > 1$, where $x_2 = \Big(\frac{\delta k+1}{\delta m - \delta k - 1}\Big)^{\frac{1}{m}}$, such that $C(p_{_{NS}}) = 1$. According to Lemma 2.2, we know that the Jacobian matrix of system (4.3) at $E_1$ has a pair of conjugate complex roots $\lambda_{1, 2}$ with $|\lambda_{1, 2}| = 1$. When $p_{_{SN}} < p < p_{_{NS}}$, $C < 1$, $|\lambda_{1, 2}| < 1$; while for $p > p_{_{NS}}$, $C > 1$, $|\lambda_{1, 2}| > 1$. Hence, we get the following results.

Theorem 4.4. The fixed point $E_2$ of system (4.3) has the following properties:

(ⅰ) in the case of $\delta (m-k) \le 1$, the fixed point $E_2$ is a sink;

(ⅱ) in the case of $\delta (m-k) > 1$,

$(a)$ if $p_{_{SN}} < p < p_{_{NS}}$, then the fixed point $E_2$ is a sink;

$(b)$ if $p = p_{_{NS}}$, then system (4.3) perhaps undergoes Neimark-Sacker bifurcation at $E_2$;

$(c)$ if $p > p_{_{NS}}$, then the fixed point $E_2$ is a source.

4.2.   Bifurcation analysis for M-G model with $1 < k\leq m$

4.2.1.   Saddle-node bifurcation for the fixed point $E_{*}$

Theorem 4.2 shows that system (4.3) may undergo saddle-node bifurcation for $p = p_{_{SN}}$, and Theorem 4.4 shows system (4.3) may undergo Neimark-Sacker bifurcation at $E_2$ for $p = p_{_{NS}}$. In this section, we shall verify the truth theoretically.

To start, according to the analysis in the above section we know that for $p = p_{_{SN}}$, the eigenvalues of Jacobian matrix for system (4.3) at $E_*(x_*, y_*)$ (where $x_* = y_* = (\rho_1\rho_2)^{\frac{1}{m-k+1}}$) are $\lambda_1 = 1$, $\lambda_2 = -\delta$. To verify the existence of saddle-node bifurcation, the following conditions must be fulfilled:

where $W$ is the eigenvector of $J^T_{E_*}$ (namely, the transverse of $J_{E_*}$) corresponding to the eigenvalue $\lambda = 1$, $V$ is the eigencector of $J_{E_*}$ corresponding to the eigenvalue $\lambda = 1$, and

It is easy to compute

therefore

which means ${\mathbf{(SN.1)}}$ is true.

$D^2G(x, y; p)$ is a double linear function, and

where

Hence, we have

where

namely, ${\mathbf{(SN.2)}}$ is fulfilled. So, system (4.3) undergoes saddle-node bifurcation when $p$ varies in the neiborhood of $p_{_{SN}}$, and the semi-stable fixed point $E_*$ is bifurcated to a sink $E_2$ and a saddle $E_1$.

4.2.2.   Neimark-Sacker bifurcation for the fixed point $E_{2}$

Next, we shall illustrate that system (4.3) undergoes Neimark-Sacker bifurcation. Similar to the proof in Section 2.2, to guarantee the existence of Neimark-Sacker bifurcation, we need to verify the transversal condition ${\mathbf{(C.1)}}$ and the nondegenerate condition ${\mathbf{(C.2)}}$. By virtue of arguements in Section 4.2, we know that $\operatorname{Det}(J_{E_2}) = 1$ when

and the Jacobian matrix of the perturbed system to system (4.3) at $E_2$ is

where $x_2 = \big(\frac{\delta k + 1}{\delta m -\delta k - 1}\big)^{\frac{1}{m}}$. It is easy to compute the eigenvalues

which satisfy the nondegereate condition ${\mathbf{(C.2)}}$.

Notice that

and $|\lambda_{1, 2}| = \sqrt{b(p_*)}.$ Combining with Theorem 4.2, we have

which verifies the transversal condition ${\mathbf{(C.1)}}$. As a result, there will be an invariant curve bifurcated from the fixed point $E_2$ of system (4.3).

What we need to do is to determine the stability of the invariant curve bifurcated from the fixed point $E_2$ of system (4.3). With the similar arguments in Section 3.2, we need to compute the first Lyapunov coefficient

where

As a result, we obtain the following theorem.

Theorem 4.5. There are no $k$ and $m$ satisfying $a^* = 0$, so $a^* \neq 0$, moreover,

${\rm{(i)}}$ if $a^* < 0$, then the fixed point $E_2$ of system (4.3) is a sink for $p < p_{_{NS}}$, while for $p > p_{_{NS}}$, there exists a stable invariant curve for system (4.3) and the fixed point $E_2$ is a source;

${\rm{(ii)}}$ if $a^* > 0$, then the fixed point $E_2$ of system (4.3) is a source for $p > p_{_{NS}}$, while for $p < p_{_{NS}}$, there exists an unstable invariant curve for system (4.3) and the fixed point $E_2$ is a sink.

4.2.3.   Resonance 1:4 for the fixed point $E_{2}$

We still choose $p$ and $\delta$ as the bifurcation parameters of (4.3), Now, consider the case of a perturbation $p_*$ and $\delta_{*}$ of parameter $p_{NS}$ and $\delta_4 = 1$. In order to ensure that the value of $p_{NS}$ is meaningful, take $m-k \neq 1$ in this case. Then, we obtain the following pertubed system:

where $p = p_{NS} + p_*$, $\delta = \delta_4 +\delta_*$, $|p_* |, |\delta_*| \ll 1$.

Next, taking $u = x - x_2, v = y - y_2$, then the fixed point $E_2$ can be moved to the origin $O(0, 0)$ and system (4.4) becomes

The Jacobian of system (4.5) at its fixed point $O$ is:

Expand the right of system (4.5) as Taylor series at $(u, v) = (0, 0)$

where $\rho = \sqrt{\|u\|^2+\|v\|^2}$,

Substituting $p_* = \delta_* = 0$ and $\delta_4 = 1$ into (4.6), we have

with

and

The characteristic value of the linearization of system (4.3) at its fixed point $O$ is

and the corresponding eigenvectors are

Use the same discussion in Section 2.2.2. We take any vector $(u, v)^T \in \mathbb{R}^2$ in the form $(u, v)^T = z q_1 + \bar{z} q_2$, where $z$ is a complex variable which transforms the map (4.7) to the complex form

where

with

Take the invertible parameter-dependent change of complex coordinate the same as Section 2.2.3:

by setting

and

by setting

with $\lambda_1 = \operatorname{i}$. The map (4.8) is changed into

where $g_{kl}$ are given by (4.9) and

which is determined by the coefficients in (4.7).

Then, we have

When $1 < k \leq m$, there is no $k \in \mathbb{N}_+$ and $m \in \mathbb{N}_+$ satisfying $C = 0$ or $D = 0$, so we have that $C \neq 0$ and $D \neq 0$ holds. Then, the fourth iterate of $\Gamma$ allows for approximation by a complex flow, and the bifurcation of 1:4 resonance will be determined by

and by $1 < k\leq m$, $k, m \in \mathbb{N}_+$, and $m-k \neq 1$, we have

Lemma 4.1. ^[44] The fourth iterate of the map (4.10) can be represented in the form

where $\varphi^t$ is the flow of a planar system

Theorem 4.6. Let $m - k \neq 1$. When the parameters $(p, \delta)$ vary in a small neighborhood of $(p_{NS}, 1)$, the system (4.3) undergoes a bifurcation of 1:4 resonance around the fixed point $E_{2}$ and admits the following bifurcation curves:

● There is a Neimark-Sacker bifurcation curve at the fixed point $E_{*}$ of system (4.3).

● Besides $E_1$, there exist four saddle fixed points $S_k$ of (4.12) around $E_2$ when $|\tilde{A}| < 1$, which are the corresponding period-4 points of (4.3). As the parameters $(p, \delta)$ vary in the neighborhood of $(p_{NS}, 1)$, these nontrivial periodic points appear or disappear in the neighborhood of $E_2$.

● There is a square heteroclinic cycle around $E_{2}$ composed by the separatrices of $S_k$ which is stable from the inside.

5.   Numerical simulation for $M-G$ model with $1 < k\leq m$

5.1.   Codim 1 bifurcations

To illustrate the theoretical results, we present some numerical simulations by virtue of the Python program.

Choose parameters $\delta = 0.2$, $k = 5$, $m = 10$. At this case, $p_{_{SN}} = 0.39$. Let $p$ vary in a neighborhood of $p_{_{SN}} = 0.39$, and we obtain the following figure.

Figure 8 shows that when $p < p_{_{SN}} = 0.39$, there is only one fixed point $O(0, 0)$ of (4.3) (see Figure 8(a)); when $p = p_{_{SN}} = 0.39$, there are two fixed points $O(0, 0)$ and $E_*$(see Figure 8 (b)); meanwhile for $p > p_{_{SN}} = 0.39$, there are three fixed points $O$, $E_1$, and $E_2$ (see Figures 8(c), (d)). It is in accordance with Theorem 4.1.

Choose the parameters $\delta = 0.433$, $m = 10$, $k = 4$. At this case $p_{_{NS}} = 0.9989$ and $p_{change} = 1.18$. Take the initial value $E_0 = (1.063, 1.063)$, and $p$ varies in the neighborhood of $p_{_{NS}}$. We get the following figures (see Figure 9). Figure 9 shows that when $p < p_{_{NS}} = 0.9989$, the fixed point $E_2$ of system (4.3) is a sink, the orbits from $E_0(1.063, 1.063)$ are all convergent to $E_2$ (see the first four sub-fiures of Figure 9). Meanwhile for $p > p_{_{NS}} = 0.9989$, there will occur a stable invariant curve (see the last four sub-figures of Figure 9), which is in accordance with Theorem 4.2.

To illustrate the existence chaos more clearly, we present the following figures, which show the change of maximum Lypunuov exponent with $p$. The figure exhibits that when $p < 1.32$, the MLE is not greater than 0, while for $p > 1.32$, the MLE is greater than $0$, which means that there exists chaos for the system when $p > 1.32$. The results are in accordance with Figure 10.

5.2.   1:4 resonance

Select $m = 10, k = 4$, and let $\delta$ vary from the neighborhood of $\delta_4 = 1$, while $p_{NS} = 2$, and let $p$ vary in $[p_{NS}-0.02, p_{NS}+0.2]$. By Theorem 4.6, the fixed point $x_2$ is a 1:4 resonance point. The phenomenon is similar to the case of $k = 0$.

Figures 11 and 12 show the complex bifurcation phenomena when $k = 4$, $m = 10$, $\delta = 1.02$ and $p$ vary. The invariant cycle induced by a Neimark-Sacker bifurcation is shown in the Figures 11(a)–(c).

As the parameter $p$ changes, further 4-period saddle points and heteroclinic cycles are generated by undergoing 1:4 resonance, showed in the Figures 12(a)–(c). A similar period-doubling phenomena occurs, and then the system enters chaos. See the Figures 13(a), (b).

It is worth noting that from Figure 14(b) and Figures 15(a)–(d), we can observe that there are period-doubling phenomena from the period-4 orbits occuring by undergoing the bifurcation of 1:4 resonance to chaos, which is different from the case of $k = 0$. From Figure 14(a), there are "period bubbling" phenomena when $\delta$ varies less than the critical value (the combination of period-doubling and inverse period-doubling ^[6,16,23]).

It should be noted that although there is a phenomenon of period doubling, there is no conventional flip bifurcation here, since multipliers of the fixed point are not equal to –1. This does not comply with the conventional flip bifurcation conditions. The phenomenon of doubling the period is a complex bifurcation phenomenon caused by a strong resonance of 1:4, similar phenomena have been mentioned in ^[6].

To describe the bifurcation of 1:4 resonance more clearly, Figure 14(d) is given to show the three-dimensional bifurcation diagrams of map (4.3) in $(p, \delta, y)$ space, respectively. Moreover, by calculating the maximium Lyapunov exponent corresponding to the system (for example, $\delta$, similarly $p$), we see that there exist both positive and negative maximium Lyapunov exponents via the changes of parameters from Figures 16(a), (b), so there exist stable fixed point or stable period windows in the chaotic region ^[6]. In general, the existences of the positive maximium Lyapunov exponents are always considered as characteristics of chaos.

6.   Conclusions

In this paper, the dynamical behaviors of the semi-discrete M-G model with $k = 0$ and $1 < k \leq m$ are analyzed. The local dynamical behaviors of fixed points are theoretically discussd by using qualitative theory, and the possible bifurcations of the system are theoretically revealed, which are verified by numerical simulations, too. We find that when the value of $k$ is different, the system exhibits different dynamic behaviors. When $1 < k < m$, and comparing with $k = 0$, the system produces more complex phenomena.

When $k = 0$ and certain parameter conditions are satisfied, the system will produce Neimark-Sacker bifurcation and 1:4 resonance. By calculating the maximum Lyapunov exponent, we find that the system will generate chaos during the process of Neimark-Sacker bifurcation. The existence of invariant cycles ending heteroclinic loops means that the metabolic level of RBCs will be in a stable state, which is beneficial for human health. The bifurcation of 1:4 resonance proves that, in certain regions, invariant circles may bifurcate from the period-2 orbits, which means that the number of RBCs exhibit periodic oscillations, but overall remain in a stable metabolic state.

Compared to the case of $k = 0$, we find that the Neimark-Saker bifurcation and the 1:4 resonance still exist in the case of $1 < k \leq m$. However, more complex dynamical behaviors will occur. The system with $k \neq 0$ undergoes saddle-node bifurcation before the Neimark-Sacker bifurcation, and there are periodic bubbling phenomena induced by the bifurcation of 1:4 resonance.

The Neimark-Sacker bifurcation and the 1:4 strong resonance of the system induce the existence of chaos, but we are still unclear about the types of chaotic attractors that exist in the system. Theoretical verification and explanation of different types of chaotic attractors will be one of our main works in the next stage. Some symmetric phenomena from the presented phase portraits and the coressponding bifurcation diagrams can be observed, which is a topic worth further exploration in the future.

In addition, we have not yet demonstrated the impact of memory effect and hereditary properties. Whether the memory effect and hereditary properties can be used to improve the adequacy of the model is also a question worth considering.

Author contributions

Yulong Li: Formal analysis, conceptualization, methodology, writing-review and editing, writing-original draft, resources, software; Long Zhou: Formal analysis, writing-review and editing, writing-original draft, software; Fengjie Geng: Supervision, writing-review and editing, funding acquisition, validation, project administration. All the authors have agreed and given their consent for the publication of this research paper.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is partly supported by the Fundamental Research Funds for the Central Universities (#2652023058).

Conflict of interest

The author declares no conflict of interest.