Delay differential equations are close to reality and reveal complex life phenomena. Due to the increase of time delay, the topological structure of differential equations may change qualitatively, and some properties of the system, such as equilibrium state, stability and periodic phenomenon will change suddenly, that is the so-called bifurcation phenomenon. In the past half century, many scholars began to pay attention to the research of bifurcation theory of delay differential equations and obtained lots of meaningful achievements [1-5].

In practical application, it is very important to study the bifurcation problem by numerical method that the continuous time model is usually discretized for the purpose of experiment or calculation. If both the discrete-time model and the continuous model exhibit similar dynamic behaviors, such as the steady-state stability behavior and persistence of the solution, the boundedness, chaos and bifurcation, then we say it is dynamically consistent[6-8]. In order to reproduce the dynamic behavior accurately, some dynamic consistent numerical methods are needed. In 2015, Jiang et al.[9] studied the Hopf bifurcation for a kind of discrete Gause-type predator-prey system with time delay, which is obtained by Euler method. They gave some parameter conditions for the existence of a unique positive fixed point, and obtained the stability result of the positive fixed point. In 2020, Euler approximation is implemented to obtain discrete version of Schnakenberg model by Din and Haider [10]. They also proved that discrete-time system via Euler approximation undergoes Hopf bifurcation as well as period-doubling bifurcation was also examined at its unique positive steady-state. Moreover, they proposed a nonstandard finite difference method (NSFDM) for Schnakenberg model and proved that NSFDM could preserve the corresponding dynamic behavior. Using time delay $ \tau $ as a parameter, Ding et al.[11] studied the dynamics of Mackey-Glass system by using trapezoidal method, and proved that with the increase of delay, the positive equilibrium point lost stability and Hopf bifurcation occured.

Compared with the Euler method, the NSFDM becomes an effective tool for nonlinear dynamic system based on its good dynamic consistency and good accuracy [12, 13]. Compared with the trapezoidal method, the NSFDM is less of computation. In this paper, we propose a NSFDM for Mackey-Glass system such that it can preserve its dynamic properties.

For the following nonlinear delay differential equations

$ \begin{equation} \dot{p}(t) = \frac{\beta \theta^n}{\theta^n + p^n(t-\tau)} - \gamma p(t), \quad t \geq 0 \end{equation} $

(1.1)

which was described by Mackey and Glass [14] as physiological control systems in 1977, where $ \beta, \theta, n $ and $ \gamma $ are all positive constants, $ p(t) $ is the density of mature cells in blood circulation, $ \tau $ is the time delay from immature cells in bone marrow to maturation and release in circulating blood. For more information on this model, the interested reader can refer to [15]. Throughout this paper, we suppose that

$ \begin{equation} \frac{\beta}{\theta}>\gamma. \end{equation} $

(1.2)

If the condition

$ \begin{equation} \frac{\gamma n \tau k^{n}}{1+k^{n}} e^{\gamma \tau}>1 / e \end{equation} $

(1.3)

is satisfied, then every positive solution of (1.1) oscillates about its positive equilibrium point[16]. Whether (1.1) is sustained for oscillations and stability arouses our great interest. Symptoms of chronic granulocytic leukemia (CGL) can be described by this model. For normal adults, the circulating granulocyte density is either stable or there is a small vibration. This vibration period is generally between 14 and 24 days. The vibration period is about 21 days, which is a change from health to sub-health. In the 30-70 days of the cycle, the density of granulocytes will appear large vibration, that is unhealthy phenomenon, the disease of CGL. In this paper, we focus on the stability and Hopf bifurcation of discrete scheme of (1.1). The effect of time delay on the dynamic behavior of the system is studied, and the conditions for the generation of Hopf bifurcation are also given.

Under transformation $ p(t) = \theta x(t) $, (1.1) changes into

$ \begin{equation} \dot{x}(t) = \frac{a}{ 1+x^n(t-\tau)} - \gamma x(t), \end{equation} $

(2.1)

let $ u(t) = x(\tau t) $, then (2.1) becomes

$ \begin{equation} \dot{u}(t) = \frac{\tau a}{1+u^{n}(t-1)}-\tau \gamma u(t), \end{equation} $

(2.2)

here $ a = \beta / \theta $. In classical FDM, the first derivative $ d u / d t $ is replaced by $ (u(t+h)-u(t))/h $, where $ h $ is the step size, however, $ d u / d t $ is replaced by $ (u(t+h)-u(t))/\phi(h) $ in NSFDM, where $ \phi(h) $ called denominator function, is the continuous function of step size $ h $, which satisfies $ \phi(h) = h+O\left(h^{2}\right), 0<\phi(h)<1, h \rightarrow 0 $.

In this paper, we select the denominator function of NSFDM which is $ \phi(h) = (1-e^{-a \tau h})/(a \tau) $, where $ h = 1 / m $ is the step size. Using the NSFDM to (2.2), we obtain the following difference scheme

$ \begin{equation} u_{k+1} = \big(1-\frac{\left(1-e^{-a \tau h}\right) \gamma}{a}\big) u_{k}+\frac{1-e^{-a \tau h}}{1+u_{k-m}^{n}}, \end{equation} $

(2.3)

here $ u_{k} $ and $ u_{k-m} $ are the approximate values to $ u(t_k) $ and $ u\left(t_{k}-\tau\right) $, respectively.

It is obvious that (2.2) and (2.3) have the same fixed point $ u_{*} $, which satisfies $ \gamma u^{n+1}+\gamma u-a = 0 $. Let $ F(x) = \gamma x^{n+1}+\gamma x-a $, then $ F^{\prime}(x) = (\mathrm{n}+1) \gamma x^{n}+\gamma>0 $ for $ x \geq 0 $. Therefore, (2.2) has a unique positive fixed point $ u_{*} $. At the same time, $ a/(1+u_{*}^{n})-\gamma u_{*} = 0 $ implies that

$ \begin{equation} u_{*}<\frac{a}{\gamma}. \end{equation} $

(2.4)

Set $ y_{k} = u_{k}-u_{*} $, then $ y_{k} $ satisfies

$ \begin{equation} \ y_{k+1} = \big(1-\frac{\left(1-e^{-a \tau h}\right) \gamma}{a}\big) y_{k}+\frac{1-e^{-a \tau h}}{1+\left(y_{k-m}+u_{*}\right)^{n}}-\frac{1-e^{-a \tau h}}{1+u_{*}^{n}}. \end{equation} $

(2.5)

Let $ Y_{n} = \left(y_{n}, y_{n-1}, \cdots, y_{n-m}\right)^{T} $, we introduce a map $ Y_{n+1} = F\left(Y_{n}, \tau\right) $, where $ F = (F_0, F_1, \cdots, F_m)^T $ and

$ \begin{equation} F_{k} = \left\{\begin{array}{l} y_{k+1} = \big(1-\frac{\left(1-e^{-a \tau h}\right) \gamma}{a}\big) y_{k}+\frac{1-e^{-a \tau h}}{1+\left(y_{k-m}+u_{*}\right)^{n}}-\frac{1-e^{-a \tau h}}{1+u_{*}^{n}}, \quad k = 0, \\ y_{n-k+1}, \quad 1\leq k \leq m. \end{array}\right. \end{equation} $

(2.6)

Obviously, the origin is a fixed point of $ Y_{n+1} = F\left(Y_{n}, \tau\right) $, and its linear part is $ Y_{k+1} = A Y_{k} $, where

$ \begin{equation} A = \begin{bmatrix} a_m & 0 & \cdots & 0 & 0 & a_0 \\ 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 & 0 \\ 0 & 0 & \cdots & 0 & 1 & 0 \end{bmatrix}, \end{equation} $

(2.7)

and

$ \begin{equation} a_m = 1-\frac{(1-e^{-a \tau h})\gamma}{a}, \quad a_0 = -\frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2}. \end{equation} $

(2.8)

Therefore, the characteristic equation of A is

$ \begin{equation} \lambda^{m+1} - \left (1-\frac{(1-e^{-a \tau h})\gamma}{a}\right)\lambda^m + \frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2} = 0. \end{equation} $

(2.9)

Lemma 2.1  For sufficiently small $ \tau>0 $, all roots of (2.9) are less than one.

Proof  When $ \tau = 0 $, (2.9) is equivalent to $ \lambda^{m+1} - \lambda^m = 0 $. It has an $ m $-fold root zero and a single root $ \lambda = 1 $. Considering the root $ \lambda(\tau) $ of (2.9), making $ \lambda(0) = 1 $, this root depends continuously on $ \tau $ and (2.9) is differentiable about $ \tau $, then we have

$ \begin{equation} \frac{d\lambda}{d\tau} = -\frac{e^{-a \tau h} h \gamma \left(n(a-\gamma u_{*}) + a\lambda^m\right)}{\lambda^{m-1}\left(m\gamma(1-e^{-a \tau h}) + a(\lambda(m+1)-m)\right)} \end{equation} $

(2.10)

and

$ \begin{equation} \frac{d\bar{\lambda}}{d\tau} = -\frac{e^{-a \tau h} h \gamma \left(n(a-\gamma u_{*}) + a\bar{\lambda}^m\right)}{\bar{\lambda}^{m-1}\left(m\gamma(1-e^{-a \tau h}) + a(\bar{\lambda}(m+1)-m)\right)} . \end{equation} $

(2.11)

Since

$ \begin{equation} \frac{d |\lambda|^2}{d \tau} = \lambda \frac{d \bar{\lambda}}{d \tau} + \bar{\lambda}\frac{d \lambda}{d \tau}, \end{equation} $

(2.12)

from (2.4) and (2.12), we obtain that

$ \begin{equation} \left. \frac{d |\lambda|^2}{d \tau} \right|_{\tau = 0, \lambda = 1} = -\frac{2 h \gamma \left(n(a-\gamma u_{*}) + a\right)}{a} < 0. \end{equation} $

(2.13)

So $ \lambda $ can not go through unit circle. Therefore, for sufficiently small $ \tau>0 $, all characteristic roots of (2.9) are within the unit circle.

Suppose $ e^{i \omega} $ is a root on the unit circle, when $ \omega \in (0, \pi] $, $ e^{i \omega} $ is the root of (2.9), so we have

$ \begin{equation} e^{i \omega} + \frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2} e^{- i m \omega} = 1-\frac{(1-e^{-a \tau h})\gamma}{a}. \end{equation} $

(2.14)

Separating the real part and the imaginary part from (2.14), we give

$ \begin{equation} \cos\omega - \frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2} \cos m \omega = 1-\frac{(1-e^{-a \tau h})\gamma}{a} \end{equation} $

(2.15)

and

$ \begin{equation} \sin\omega - \frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2} \sin m \omega = 0. \end{equation} $

(2.16)

So we obtain

$ \begin{equation} \cos\omega = 1+\frac{\gamma^2(1-e^{-a \tau h})^2 \left(a+n(a-\gamma u_{*})\right)\left(a-n(a-\gamma u_{*})\right)}{2a^3 \left(a-\gamma(1-e^{-a \tau h})\right)}. \end{equation} $

(2.17)

If $ n(a-\gamma u_{*})/a<1 $, then $ \cos\omega>1 $, which is a conflict, so the lemma is proved.

Lemma 2.2  Supposing that $ n(a-\gamma u_{*})/a<1 $, (2.9) has no modules of roots more than one.

From (2.16), we get that

$ \begin{equation} \frac{(1-e^{-a \tau h})(a - \gamma u_{*})n \gamma}{a^2} = \frac{\sin \omega}{\sin m \omega} \end{equation} $

(2.18)

is positive, $ \sin\omega $ has the same symbol as $ \sin m \omega $. So there exists a real sequence $ \omega_i $ which satisfies

$ \begin{equation} \omega_i \in \left(\frac{2i \pi}{m}, \frac{(2i + 1)\pi}{m}\right), i = 0, 1, 2, \cdots, \left[\frac{m-1}{2}\right], \end{equation} $

(2.19)

where $ [.] $ is the greatest integer function.

Let $ \lambda_i(\tau) = r_i(\tau) e^{i \omega_i(\tau)} $ is the root of (2.9), $ \tau = \tau_i $ satisfyies $ r_i(\tau_i) = 1 $ and $ \omega_i(\tau_i) = \omega_i $, we have the following lemma.

Lemma 2.3  Suppose that $ n(a-\gamma u_{*})/a>1 $, then

$ \begin{equation} \left.\frac{dr_{i}^2(\tau)}{d\tau}\right|_{\tau = \tau_i, \omega = \omega_i}>0. \end{equation} $

(2.20)

Proof  From (2.8) and (2.9), we have

$ \begin{equation} \bar{\lambda}^m = \frac{\lambda-a_m}{a_0} = -\frac{a^2 (\lambda -1) + a(1-e^{-a \tau h})\gamma}{(1- e^{-a \tau h})(a- \gamma u_{*})n \gamma}. \end{equation} $

(2.21)

Then by (2.10) and (2.21), we obtain

$ \begin{equation} \begin{aligned} &\bar{\lambda}\frac{d\lambda}{d\tau}\\ = & -\frac{\bar{\lambda} e^{-a \tau h }h \gamma\left(n(a-\gamma u_{*})+a\lambda^m\right)}{\lambda^{m-1}\left(m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)\right)} = -\frac{\bar{\lambda}^m \lambda \bar{\lambda} e^{-a \tau h }h \gamma\left(n(a-\gamma u_{*})+a\lambda^m\right)}{\bar{\lambda}^m \lambda^{m}\left(m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)\right)}\\ = & -\frac{ e^{-a \tau h }h \gamma\left(n(a-\gamma u_{*})\bar{\lambda}^m+a\right)}{ m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)} = -\frac{ e^{-a \tau h }h \gamma\left(n(a-\gamma u_{*})\bar{\lambda}^m+a\right)}{ m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)}\\ = & -\frac{ e^{-a \tau h }h \gamma\left(n(a-\gamma u_{*})(-\frac{a^2 (\lambda -1) + a(1-e^{-a \tau h})\gamma}{(1- e^{-a \tau h})(a- \gamma u_{*})n \gamma})+a\right)}{ m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)} = \frac{\frac{e^{-a \tau h} h a^2 (\lambda -1)}{1- e^{-a \tau h}}}{m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)}. \end{aligned} \end{equation} $

(2.22)

From (1.2) we get

$ \begin{equation} \begin{aligned} &\frac{d r_{i}^2}{d \tau}\\ = & \lambda \frac{d \bar{\lambda}}{d \tau} + \bar{\lambda} \left. \frac{d \lambda}{d \tau} \right|_{\tau = \tau_i, \omega = \omega_i} = 2\operatorname{Re} \big(\bar{\lambda}\left. \frac{d \lambda}{d \tau}\right|_{\tau = \tau_i, \omega = \omega_i}\big) = 2\operatorname{Re}\big( \frac{\frac{e^{-a \tau h} h a^2 (\lambda -1)}{1- e^{-a \tau h}}}{m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right)} \big)\\ = & 2\operatorname{Re}\big(\frac{e^{-a \tau h} h a^2(\lambda-1)(m \gamma (1-e^{-a \tau h})+a\left(\bar{\lambda}(m+1)-m\right))}{(1-e^{-a \tau h}) \left (m \gamma (1-e^{-a \tau h})+a\left(\lambda(m+1)-m\right) \right) \left(m \gamma (1-e^{-a \tau h})+a\left(\bar{\lambda}(m+1)-m\right) \right)} \big)\\ = & \frac{2 e^{-a \tau h} h a^2 (1-\cos\omega)\left(2a m + a -m\gamma(1-e^{-a\tau h}) \right)}{(1-e^{-a \tau h})\left|m\gamma (1-e^{-a \tau h}) - am + a(m+1)e^{i\omega} \right|^2}\\ = & \frac{2 e^{-a \tau h} h a^2 (1-\cos\omega)\left( m \left(a -\gamma(1-e^{-a\tau h}) \right)+ ma +a \right)}{(1-e^{-a \tau h})\left|m\gamma (1-e^{-a \tau h}) - am + a(m+1)e^{i\omega} \right|^2} > 0. \end{aligned} \end{equation} $

(2.23)

Lemma 2.4

(i) If $ n(a-\gamma u_{*})/a < 1 $, then for any $ \tau>0 $, all roots of (2.9) are in the unit circle.

(ii) If $ n(a-\gamma u_{*})/a > 1 $, then (2.9) has a sequence of monotonically increasing numbers column $ \tau_{i}(i = 0, 1 $, $ 2, \cdots, [(m-1) / 2]) $. When $ \tau = \tau_{i} $, (2.9) has a pair of roots $ e^{\pm i\omega_i} $. If $ \tau>\tau_{i} $, the roots of (2.9) are out of unit circle. If $ \tau \in\left[0, \tau_{0}\right) $, the modules of roots in (2.9) are all less than one. If $ \tau \in\left[\tau_{i}, \tau_{i+1}\right) $, then the equation has $ 2(k+1) $ roots that have modulus more than one.

Proof  By Lemma 2.1, Lemma 2.2 and Corollary 2.4 in [17], we can obtain (i).

If $ n(a-\gamma u_{*})/a > 1 $, let $ \tau_{i} $ be the same as that in (2.20), then (2.9) has the root if and only if $ \tau = \tau_{i} $ and $ \omega = \omega_i $. By Lemma 2.1 and Lemma 2.2, the modulus of the root of (2.9) less than one if $ \tau \in [0, \tau_{0}) $. If $ \tau = \tau_{0} $, then the modulus of roots are less than one except $ e^{\pm i\omega_i} $. When $ \tau \geq \tau_{i} $, (2.9) has roots passing through the unit circle from $ e^{\pm i\omega_i} $. In addition, according to Rouche theorem, the number of eigenvalues whose modulus are more than one is obtained naturally.

From Lemma 2.4, the stability of the zero solution can be obtained in the following theorem.

Theorem 2.5

(i) If $ n(a-\gamma u_{*})/a < 1 $, then $ u = u_{*} $ is asymptotically stable for any $ \tau \geq 0 $.

(ii) If $ n(a-\gamma u_{*})/a > 1 $, then $ u = u_{*} $ is asymptotically stable for $ \tau \in\left[0, \tau_{0}\right) $, and unstable for $ \tau>\tau_{0} $.

(iii) For $ n(a-\gamma u_{*})/a > 1 $, (2.5) undergoes a Hopf bifurcation at $ u_{*} $ when $ \tau = \tau_{i}, $ for $ i = 0, 1 $, $ 2, \cdots, [(m-1) / 2] $.

Let $ a = 2, \gamma = 1 $ and $ n = 4 $ in (2.2), it is easy to find that the positive equilibrium point $ u_{*} = 1 $, so the condition $ n(a-\gamma u_{*})/a > 1 $ holds.

In Table 1, we give the absolute errors (AE) and the relative errors (RE) at $ t = 10 $ of the NSFDM with initial value $ u(t) = 1.1 $ and $ \tau = 1 $. From this table we know that the NSFDM has good convergence.

In Table 2, we give the values of $ \tau_k $ for different step size $ h = 1/2, 1/4, 1/8 $ and $ 1/16 $. From Theorem 2.5 and Table 2 we can see that $ \tau_k $ is the bifurcation points. Furthermore, in Figures 1-4, we present the numerical solution and phase diagram of the system discretized by the NSFDM. From Theorem 2.5 we can conclude that the equilibrium is asymptotically stable for $ \tau \in [0, \tau_{0}) $ ($ \tau_{0} \approx 1.0986 $ in Figure 1, $ \tau_{0} \approx 0.7643 $ in Figure 2, $ \tau_{0} \approx 0.6734 $ in Figure 3 and $ \tau_{0} \approx 0.6368 $ in Figure 4), unstable for $ \tau>\tau_{0} $ and an attracting bifurcating periodic solution exists for $ \tau>\tau_{0} $. This is just what Figures 1-4 show intuitively.

It is not difficult to see that the given numerical results illustrate the correctness of the theoretical analysis.

Mackey-Glass system is discretized by NSFDM, the influence of time delay on blood cell density is analyzed. If the time delay exceeds a certain critical value, Hopf bifurcation will occur and result in the density of mature cells produce periodic oscillation. However, if the delay is small enough, the equilibrium is stable.

From a biological point of view, if we can put off the production of immature cells in bone marrow to the normal time of their release in the circulating blood, stabilize the density of mature cells in the blood circulation, the disease will be brought under control. Our analysis results can provide critical insights and guidance for the analysis and design of control schemes from the perspective of dynamics and control theory.