key: cord-0011265-xlx73lhe
authors: Darlai, Rachadawan; Moore, Elvin J.; Koonprasert, Sanoe
title: Andronov–Hopf and Neimark–Sacker bifurcations in time-delay differential equations and difference equations with applications to models for diseases and animal populations
date: 2020-04-29
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-02646-5
sha: 8ceaaa87a76db5f305006a242a690c597b302a5c
doc_id: 11265
cord_uid: xlx73lhe

In many areas, researchers might think that a differential equation model is required, but one might be forced to use an approximate difference equation model if data is only available at discrete points in time. In this paper, a detailed comparison is given of the behavior of continuous and discrete models for two representative time-delay models, namely a model for HIV and an extended logistic growth model. For each model, there are seven different time-delay versions because there are seven different positions to include time delays. For the seven different time-delay versions of each model, proofs are given of necessary and sufficient conditions for the existence and stability of equilibrium points and for the existence of Andronov–Hopf bifurcations in the differential equations and Neimark–Sacker bifurcations in the difference equations. We show that only five of the seven time-delay versions have bifurcations and that all bifurcation versions have supercritical limit cycles with one having a repelling cycle and four having attracting cycles. Numerical simulations are used to illustrate the analytical results and to show that critical times for Neimark–Sacker bifurcations are less than critical times for Andronov–Hopf bifurcations but converge to them as the time step of the discretization tends to zero.

In recent years, time-delay differential equation and difference equation models have been studied by many authors (see, e.g., [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] ) as they are useful tools for modeling a wide variety of systems in areas including traditional areas such as physics and engineering and newer areas such as disease transmission, medical research, optimal drug treatment, bioeconomics, agriculture, finance, insurance, and environmental protection. In many of these time-delay models, bifurcations occur as values of parameters are changed. For example, Andronov-Hopf bifurcations can occur in differential equation systems and Neimark-Sacker bifurcations can occur in difference equation systems at critical values of time delays. In some cases, the changes in system parameters can also lead to chaotic-type solutions (see, e.g. [17, 18] ).

In many applications, the preferred model is a differential equation model. However, as the growth rate of diseases (such as HIV) or other kinds of populations (such as fish or animal populations) can be a slow process or the collection of data may often only be carried out at regular intervals such as a month or a year, it is often only possible to construct difference equation models. One method that is often used to construct a difference equation model is to use a first-order Euler method to approximate the differential equation model [1, 2, [19] [20] [21] [22] . This method is also often used to solve Itó stochastic differential equations as, for example, in the Euler-Mayurama method [1, 23] .

As stated above, an important property of many time-delay differential equation models is that they have Andronov-Hopf bifurcations and an important property of many time-delay difference equation models is that they have Neimark-Sacker bifurcations [24] . Because difference equation models are often used to approximate differential equation models, we believe that it is important to study the relationship between the bifurcations in the differential equations and the approximating difference equations.

In this paper, we analyze the bifurcation properties of one-dimensional time-delay differential equation models and the corresponding discrete models obtained by the forward Euler approximation. As examples, we consider time-delay versions of two systems. The first system is a model discussed by Roberts and Saha [25] and Ding et al. [26] for transmission of HIV in a human population. This model includes the effects of vertical HIV transmission from mother to baby, the effects of births and deaths and of treatment by antivirals. The second system is an extended logistic growth model (ELM) that has been applied to forecasting short lifecycle products and services by Trappey and Wu [27] , and to growth of single-species populations by He et al. [1] and Sakanoue [28] .

For both the HIV and the ELM models, there are seven different versions of time-delay models that can be created from the original models. For each model, we first investigate the properties of the seven different versions of differential equation models and prove conditions for the existence and stability of equilibrium points and for the existence of Andronov-Hopf bifurcations at critical values of the time delays. We then investigate the properties of discretized versions of the models and prove conditions for the existence and stability of equilibrium points and for the existence of Neimark-Sacker bifurcations at critical values of the time delays. We show that both differential equation and difference equation models may have bifurcations from both disease-free and endemic equilibrium points. For the disease-free equilibrium, four of the seven different versions satisfy the bifurcation conditions, and for the endemic equilibrium five of the seven versions have bifurcations. For the bifurcations from the endemic equilibrium points, we show that for each version with bifurcations the critical delay times for Neimark-Sacker bifurcations are less than the critical times for Andronov-Hopf bifurcations but converge to them as the time step of the discretization in the discretized model tends to zero. Numerical simulations are presented for a range of parameter values to illustrate the analytical results. Applications of the results are given to the modeling of HIV transmission and to the modeling of populations.

In this section, we describe the differential equation models for HIV and extended logistic growth (ELM) that we study in this paper.

Many researchers have developed mathematical models to study the transmission of HIV at either the CD4+ T-cell level (see, e.g., [3, 4, 7, 29-33]) or at the population level (see, e.g., [25, 26, 34, 35] ). In the present research, we study a model for transmission of HIV at the population level originally proposed by Roberts and Saha [25] as a general epidemic model, and later developed by Ding et al. [26] as a stochastic differential equation model for the progression to AIDS in a population infected with HIV. The nonlinear, logistictype, differential equation model discussed by Ding et al. is as follows (see also [35] ):

where x(t) is the proportion of the total population that is infected by HIV at time t, p (0 < p < 1) is the vertical transmission probability (the fraction of babies born with HIV infection), B is the birth rate for the population, β is the transmission rate on contact between an infected and an uninfected individual, C is the contact rate between infected and uninfected individuals, and α is the increase of the death rate due to the HIV infection. The development of antiretroviral therapy using reverse transcriptase inhibitors (RTI) and protease inhibitors (PI) has been shown to be an effective method of controlling the spread of HIV by depressing the level of virus in an HIV+ person below a detectable level [29] [30] [31] [32] [33] 36] and to effectively stop transmission of HIV from an HIV+ person to an uninfected person [37] [38] [39] [40] [41] [42] . In the present model, we assume that the effects of both the RTI and the PI can be included in the model in (1) as factors reducing the value of β (the rate of infection on contact) and p (the vertical transmission probability). Following Darlai et al. [35] , we assume that

where n av is an antiretroviral therapy factor (0 ≤ n av < 1) and β 0 and p 0 are, respectively, the infection rate of a susceptible person and the vertical transmission probability in the absence of antiretroviral therapy. For simplicity, we rewrite (1) as

where δ = (1p)B > 0 and ε = βCα > 0. Equation (3) is the well-known logistic equation model with an added death rate term -δx(t).

As stated previously, extended logistic growth models have been applied to forecasting short lifecycle products and services by Trappey and Wu [27] , and to growth of singlespecies populations by Sakanoue [28] . 

In this paper, we study the extended logistic growth model discussed by He et al. [1] :

where x(t) is the population at time t, and r is the natural death rate of the population at low population levels. The parameters β, γ , K are positive parameters, where β is a natural birth rate at low population levels, K is the "carrying capacity" of the system, and the factor (x(t)/K) γ is the rate at which the birth rate decreases as x(t) → K . This reduction in the birth rate could be due, for example, to factors such as overcrowding or limits in the available food supply.

In this paper, we extend the models in Eqs. (3) and (4) by introducing time delays into the models. As shown in Table 1 , there are seven time-delay versions of the two differential equation models. We use code "n" for no time delay and "d" for time delay. Then, for example, "dnd" means time delays in positions 1 and 3 and no time delay in position 2.

Darlai et al. [35] have studied the Andronov-Hopf and Neimark-Sacker bifurcations in HIV3:nnd and HIV7:ddd, and He et al. [1] have studied the Neimark-Sacker bifurcations in ELM6:ndd.

Using the first-order Euler method, we discretize the independent variable t with a step size h, and replace a first-order system of differential equations of the form

where t n = t 0 + nh, n = 0, 1, 2, . . . , and w n = x(t n ). The HIV and ELM difference equation models for zero time delays can then be written in the form HIV: w n+1 = w nδhw n + εhw n (1w n ), ELM: w n+1 = w nrhw n + βhw n [1 -(w n /K) γ ].

Note As in all difference equation approximations to differential equations, it is necessary to check that the approximation is a numerically stable approximation to the original 

differential equation. It is well known (see, e.g., [43] ) that the forward Euler approximation can be numerically unstable unless an upper limit is placed on the step size h. It is therefore necessary to check that a Neimark-Sacker bifurcation corresponds to an Andronov-Hopf bifurcation of the original differential equation model and that it is not a bifurcation arising from the instability of the forward-Euler method.

To obtain the first-order Euler approximations to the seven HIV and ELM time-delay differential equations in Table 1 , we assume that the time delay τ is divided into m equal intervals. Then the step size is h = h m τ , where h m = 1/m, and we obtain the seven difference equations in Table 2 . It can be seen that each of these Euler approximation equations are difference equations of order m + 1.

For each model, the equilibrium points x * for the differential equations (3) and (4), and the seven time-delay differential equation versions in Table 1 are the same and are obtained by setting dx dt = 0. For each model, there is a trivial equilibrium point x * = 0 and an endemic equilibrium point which exists only if x * > 0. The equilibrium points for the models are

For both models, the endemic equilbrium exists only if R 0 > 1.

The conditions for local asymptotic stability of the equilibrium points can be obtained by using a standard approach, such as the next-generation method [44] , or by checking the 

eigenvalues of the linearized system at the equilibrium points [45] . In this case, it is easy to check the eigenvalues of the linearized equations about the equilibrium points.

We can obtain the linearized versions of the delay equations in Table 1 by defining perturbations y(t) = x(t)x * and y(tτ ) = x(tτ )x * . The linearized versions can then be written in the standard form

where the ρ and η values for the linearized equations at the disease-free and endemic equilibrium points are shown in Table 3 . As usual, we assume a trial solution y(t) = e λt . The characteristic equation from the trial solution is then

Then, if the real parts of all eigenvalues λ of (9) are negative, the general solution y(t) → 0 as t → ∞ and the equilibrium point x * is locally asymptotically stable. Zero time delay. For both disease-free and endemic equilibrium points, the condition for local asymptotic stability is λ = ρη < 0. For the HIV model, the stability condition for the disease-free equilibrium becomes ρη = εδ < 0 and for the endemic equilibrium it becomes ρη = δε < 0. Therefore, the HIV disease-free equilibrium is locally asymptotically stable if εδ < 0, and the basic reproductive number is then R 0 = ε δ < 1. Similarly, the ELM disease-free equilibrium is locally asymptotically stable if βr < 0, and the basic reproductive number is R 0 = β r < 1. Note that, for zero time delay, the disease-free equilibrium is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1, and that the endemic equilibrium exists only if R 0 > 1. 

We first look for possible solutions of the characteristic equation (9) for the endemic equilibrium that satisfy condition (C1) by looking for a purely imaginary solution λ c = iφ c with φ c > 0 for τ = τ c > 0.

If these conditions are satisfied, then the corresponding φ c and τ c values are given by

Note: Since complex solutions must exist in complex conjugate pairs, we will assume that φ c = η 2ρ 2 > 0.

Proof Substituting λ = iφ into (9), separating real and imaginary parts, and solving for φ, we obtain

where φ = η 2ρ 2 is real and nonzero if and only if |ρ| < |η|. Then, solving for τ in (11), we obtain the expressions for τ c in (10), and from the second expression for τ c in (10), we see that a real positive solution for τ c exists if and only if η > 0. The minimum value of τ c satisfying Eq. (10) also satisfies condition (C1) and therefore it is a possible Andronov-Hopf bifurcation point.

Then, using the values of ρ and η for the HIV and ELM delay equations given in Table 3 , we find that, as shown in Table 4 , possible Andronov-Hopf bifurcations can occur from disease-free equilibrium states for four of the time-delay versions and from endemic equilibrium states for five of the versions. However, we have found that bifurcations from the disease-free states give limit cycles that contain negative population values. For these "limit cycle" regions, the extra condition that the state variable cannot be negative must be added to the mathematical model.

In the remainder of this paper, we will only consider bifurcations from endemic equilibrium states. Proof Let λ = μ + iφ, where μ and φ are real, and we can assume that φ ≥ 0. Note: An example of the variation of μ and φ as τ is increased that is used in this lemma to prove condition (C2) is shown in Fig. 1 . The real and imaginary parts of (9) are μ = ρηe -μτ cos(φτ ), φ = ηe -μτ sin(φτ ), and the real and imaginary parts of the derivative of the characteristic equation (9) with respect to τ can be written as:

Clearly, μ and φ are continuous functions of τ , but the derivatives can be discontinuous if = 0. We consider three cases:

Further, since μ < 0 and = 1 > 0, we have dμ dτ < 0 and dφ dτ = 0 at τ = 0. As τ increases, μ will decrease and remain negative and φ will remain zero. Therefore, the eigenvalue will be real and negative.

The derivatives dμ dτ and dφ dτ will therefore be discontinuous at τ d . Then, since φ = 0 is actually an unstable solution of (13), the discontinuity at τ d will make φ become nonzero. 3. 1τ ηe -μτ < 0. As τ increases, φ becomes non-negative. If φ > 0, then, for τ in the positive neighborhood of τ d , we have μ < ρη < 0, and therefore for small φ, dφ dτ > 0 if φ > 0 and dφ dτ < 0 if φ < 0. Figure 1 shows the case of φ > 0. In both cases, μ < 0 and therefore dμ dτ > 0 and μ increases monotonically until μ = 0. At this point, we see that the Andronov-Hopf bifurcation condition (C1) will be satisfied and we can identify the τ value as the critical value τ c and the φ value as the critical value φ c defined in Lemma 1. Since μ < 0 for τ < τ c , we have proved condition (C2) for τ c . Finally, substituting μ = 0, φ = φ c and τ = τ c into (12) and (14), we have

Since ηρ > 0, the denominator is always positive and therefore condition (C3) is satisfied. Also, since dμ dτ | τ =τ c > 0, the bifurcation will occur as τ increases to τ c . The proof is complete.

Combining the results of Lemmas 1 and 2, we obtain the following theorem.

Theorem 2 Necessary and sufficient conditions for the existence of an Andronov-Hopf bifurcation in the endemic equilibrium points in equation (8) are η > 0 and |ρ| < η. Then the Andronov-Hopf bifurcation exists at a critical time delay:

where φ c = η 2ρ 2 . Further, the bifurcation occurs as τ increases through τ c .

The equilibrium points w * for the difference equation models are obtained by setting w n+1 = w n = w n-m = w * in equation (6) . For each model, we obtain a trivial equilibrium point w * = 0 and an endemic equilibrium point given by

For the HIV model, the endemic equilibrium exists only if ε > δ, and for the extended logistic model, the endemic equilibrium exists only if β > r. The equilibrium points in (17) are also equilibrium points for the differential equation models.

To obtain the basic reproductive numbers R 0 for local stability of the disease-free equilibrium points, we check the eigenvalues of the linearized system at the disease-free equilibrium. We let y n = w nw * and assume that y n is small. Then the linearized equations are:

Therefore, the HIV disease-free equilibrium is locally asymptotically

Similarly, the ELM disease-free equilibrium is locally asymptotically stable if R 0 = |1h m (rβ)| < 1, i.e., if R 0 = β r < 1. Note that in both cases, the disease-free equilibrium is locally asymptotically stable if R 0 < 1, whereas the endemic equilibrium exists only if R 0 > 1. Also, as noted in section (4), the stability conditions for the numerical stability of the forward Euler approximation put an upper limit on the step size h = h m τ in the difference equations (18).

We obtain the conditions for Neimark-Sacker bifurcations of the endemic equilibrium solutions of the nonlinear HIV and ELM equations (18) from the linearized equations about the endemic equilibrium points.

To obtain the linearized, time-delayed versions of the nonlinear difference equations, we define y n+1 = w n+1w * , y n = w nw * and y n-m = w n-mw * . Then we obtain

where ρ and η are constants and h m = 1/m. The values of these constants for the seven linearized versions of the discrete HIV and ELM equations are the same as the values given in Table 3 for the seven linearized versions of the HIV and ELM differential equations in Table 1 .

Then, assuming a trial solution of the form y n = λ n for (19), we obtain the characteristic equation, which it is convenient to write in the alternative forms

From the real and imaginary parts of the characteristic equations in (20)- (22), we obtain the conditions mr m+1 cos (m + 1)ω = (m + ρτ )r m cos(mω) + ητ ,

mr m+1 sin (m + 1)ω = (m + ρτ )r m sin(mω),

mr m+1 cos(ω) = (m + ρτ )r m + ητ cos(mω),

From bifurcation theory [24] , the conditions for the existence of a Neimark-Sacker bifurcation are as follows. 

,

Note: Since any complex eigenvalues must occur in complex conjugate pairs, existence of an eigenvalue with ω = ω c > 0 implies existence of an eigenvalue with ω = -ω c < 0.

Proof Taking the real part of the characteristic equation P 3 (λ) = 0 in (22) for λ = e ±iω , we obtain the condition in (29) . Then, taking the imaginary part of P 1 (λ) = 0, and the real and imaginary parts of P 2 (λ) = 0, for λ = e ±iω , we obtain the formulas for τ 1 (ω), τ 2 (ω), τ 3 (ω) in (30) . Therefore, if conditions 1 and 2 of the lemma are satisfied, then λ c = e iω c is a solution of the characteristic equations in (20)-(22) for a real positive value of the time delay.

We now prove the following lemma.

Necessary and sufficient conditions for the critical value τ c > 0 defined in Lemma 3 to exist are: |ρ| < η and η > 0.

Proof If λ = r for real r, then the characteristic equation P 2 (λ) = 0 in (21) is m(r -1) = τ (ρηr -m ). Then, differentiating this equation with respect to τ , we have

For τ = 0, a real solution of the characteristic equation m(r -1) = τ (ρηr -m ) is r = 1 and, at r = 1, the derivative dr dτ = ρ-η m . Therefore, if ρη ≥ 0, then a real eigenvalue will remain on or move outside the unit circle for τ > 0, and hence τ c cannot be a Neimark-Sacker bifurcation point. Therefore, ρη < 0 is a necessary condition for real eigenvalues to lie inside the unit circle in the positive neighborhood of τ = 0.

We next prove that necessary and sufficient conditions for τ c > 0 to exist are that η 2ρ 2 > 0. From the formulas for τ 2 (ω c ) and τ 3 (ω c ) in (30) and the condition that 0 < ω c < π , it can be seen that necessary conditions for τ 2 (ω c ) and τ 3 (ω c ) to be positive are that sin(ω c ) > 0 and sin( 1 2 ω c ) > 0. Therefore, necessary conditions for τ c = τ 2 (ω c ) = τ 3 (ω c ) > 0 to exist are η cos(mω c ) > ρ, η sin(mω c ) > 0, and therefore η 2 > ρ 2 .

Next, if ρη < 0 and η 2ρ 2 > 0, then ρ + η > 0. Then, since ρη < 0 and ρ + η > 0, it is necessary that η > 0 and |ρ| < η. We now prove that the conditions are sufficient. In Eq. (29), the function R(ω) is a continuous function of ω for 0 ≤ ω ≤ π . We have R(0) = ρη < 0 and R( 2π 2m+1 ) = ρ cos( π 2m+1 ) + η > 0. Therefore, from the intermediate value theorem of elementary calculus, R(ω) = 0 for at least one ω ∈ (0, 2π 2m+1 ). Therefore, from Lemma 3, τ = τ c exists and satisfies condition (C1) of the Neimark-Sacker theorem. The proof is complete We now prove condition (C2). Note: An example of the variation of r and ω as τ is increased that is used in this lemma to prove condition (C2) is shown in Fig. 2 .

Proof Let λ = re iω , where r ≥ 0 and ω = 0 are real. Also, since any complex eigenvalues must occur in complex conjugate pairs, we can assume that 0 < ω < π .

Then, differentiating Eqs. (23) and (24) with respect to τ and simplifying using (23)-(26), we can write the derivatives for r and ω in the forms dr dτ = 1 r m+1ητ ρr mη cos(ω) + ρητ r m cos(ω) -cos(mω) + ηr m+1 cos(ω) -cos (m + 1)ω + η 2 τ 1 -cos(ω) = 1 τ r m+1ητ ρr mη + 2ητ r m sin 2 1 2 mω + 2 ηρr m r m+1 sin 2 1 2 (m + 1)ω + sin 2 1 2 ω , 

As for the Andronov-Hopf bifurcation, we consider three cases. 1. r m+1ητ > 0. Then, at τ = 0, we have from (20) λ = r = 1, and ω = 0, and from (34) that dω dτ = 0. Hence, in the positive neighborhood of τ = 0, we have ω = 0 and then = m r m+1ητ

Therefore, as shown in Fig. 2 , for r m+1ητ > 0, the modulus r decreases and the angle ω remains zero. 2. r m+1ητ = 0. At τ = τ d = r m+1 η , the denominator = 0 and, as shown in Fig. 2 , the derivatives dr dτ and dω dτ are discontinuous. Then, since ω = 0 is an unstable solution of (34), ω becomes nonzero (positive in Fig. 2 ) in the positive neighborhood of τ d . 3. r m+1ητ < 0. In this case, it can be seen from Eq. (33) that dr dτ > 0, and therefore r increases.

We now show that r m+1ητ < 0 for all τ ≤ τ c . From (26), we have

Also, from Lemma 4, we have the condition that, for ω = 0, the critical value of ω c ∈ (0, 2π 2m+1 ). It is straightforward to prove that sin(mω) m sin(ω) < 1 for ω ∈ 0, 2π 2m + 1 .

Therefore r m+1ητ < 0 for τ d < τ ≤ τ c , and hence r is a continuous, monotonically increasing function of τ for τ d ≤ τ ≤ τ c . Therefore τ = τ c is the minimum value of τ for which r = 1. The proof is complete.

We now prove condition (C3).

If the necessary and sufficient conditions stated in Lemma 4 are satisfied, then the critical τ c > 0 defined in Lemma 3 also satisfies condition (C3) of the Neimark-Sacker bifurcation theorem.

Proof As shown in Lemma 5, dr dτ > 0 for τ d < τ ≤ τ c , and therefore condition (C3) is satisfied.

We now prove condition (C4). Proof Since 0 < ω c < 2π 2m+1 and m > 1, we have 0 < kω c < 2π for k = 1, 2, 3, 4 and therefore e ikω c = 1.

Then, combining Lemmas 3-7, we have the following theorem. Table 4 for the seven versions of the HIV and ELM differential equations in Table 1 .

In this section, we find the direction and stability of Neimark-Sacker bifurcations by studying condition (C5) of the Neimark-Sacker theorem for the five versions (dnn, nnd, dnd, ndd, ddd) of the time-delay difference equations that show bifurcations. To prove condition (C5), we follow the method given in Li [19] and Kuznetsov [24] to compute the function called c 1 (τ ) for the critical value τ = τ c . The main ideas in the derivation are to replace the single nonlinear difference equations of order m + 1 in the perturbations y n = w nw * by a system of m + 1 first-order difference equations and then to examine the terms up to third order in y n in the Taylor series expansions of the right-hand sides of these systems of first-order nonlinear difference equations. The critical function c 1 (τ c ) can then be computed from these terms.

The terms up to third order of the nonlinear difference equations for the five versions of the HIV and ELM models are shown in Table 5 .

We now introduce a vector Y n = (y n , y n-1 , y n-2 , . . . , y n-m ) T , n ≥ m and write the system of first-order difference equations up to third order in the Taylor series expansion in the form

In (39), A is the Jacobian of the system of m+1 first-order difference equations and B and C are the second-and third-order terms in the system. We will first discuss the Jacobian and eigenvectors and then give the formulas for the second-and third-order terms. Finally, we will derive the formula for the critical coefficient c 1 (τ c ) for the five HIV and ELM versions that have Neimark-Sacker bifurcations.

The Jacobian matrix A = J(τ ) for the linearized system of m + 1 first-order equations corresponding to the nonlinear difference equations in (19) is given by Table 3) Version HIV ψ = h m τ dnn y n+1 = y n + εψy n -δψy n-m -εψ(y n + w * ) 2 second order y n+1 = y n + ρψy n -ηψy n-m -2εψy 2 n nnd y n+1 = y n -δψy n + εψy n-m -εψ(y n + w * )(y n-m + w * ) second order y n+1 = y n -ηψy n-m -εψ(y n y n-m + y n-m y n ) dnd y n+1 = y n -εψy n -δψy n-m -εψ(y n + w * )(y n-m + w * ) second order y n+1 = y n + ρψy n -ηψy n-m -εψ(y n y n-m + y n-m y n ) ndd y n+1 = y n -δψy n + εψy n-m -εψ(y n-m + w * ) 2 second order y n+1 = y n + ρψy n -ηψy n-m -2εψy 2 n-m ddd y n+1 = y n + (ε -δ)ψy n-m -εψ(y n-m + w * ) 2 second order y n+1 = y n -ηψy n-m -2εψy 2

γ -2 γ dnn y n+1 = y n + βψy n -rψy n-m -αψ β K γ (y n + w * ) γ +1 third order y n+1 = y n + ρψy n -ηψy n-m -(γ 2 + γ )νψy 2 n -(γ 3 -γ )ξ ψy 3 n nnd y n+1 = y n + ψΛy n -αψ β K γ (y n + w * )(y n-m + w * ) γ third order y n+1 = y n -ηψy n-m -νψ[(γ 2 -γ )y 2 n-m + γ (y n y n-m + y n-m y n )] -ξψ[(γ 3 -3γ 2 + 2γ )y 3 n-m + (γ 2 -γ )(y n y 2 n-m + y n-m y n y n-m + y 2 n-m y n )] dnd y n+1 = y n + αψ(1 + β)y n -rψy n-m αψ β K γ (y n + w * )(y n-m + w * ) γ third order y n+1 = y n + ρψy n -ηψy n-m νψ[(γ 2 -γ )y 2 n-m + γ (y n y n-m + y n-m y n )] -ξψ[(γ 3 -3γ 2 + 2γ )y 3 n-m + (γ 2 -γ )(y n y 2 n-m + y n-m y n y n-m + y 2 n-m y n )] ndd y n+1 = y n -rψy n + αψ(1 + β)y n-m -αψ β K γ (y n-m + w * ) γ +1 third order y n+1 = y n + ρψy n -ηψy n-m -(γ 2 + γ )νψy 2 n-m -(γ 3 -γ )ξ ψy 3 n-m ddd y n+1 = y n + ψΛy n-m -αψ β K γ (y n-m + w * ) γ +1 third order y n+1 = y n -ηψy n-m -(γ 2 + γ )νψy 2 n-m -(γ 3 -γ )ξ ψy 3

Since the Jacobian matrix is the companion matrix [43] of the characteristic equation (20) 

where we have chosen q m = 1 as a normalization constant. Following Li [19] , we also introduce the adjoint eigenvector r e iω c = D 1, σ e imω c , σ e i(m-1)ω c , . . . , σ e i2ω c , σ e iω c T ,

of the transposed Jacobian matrix J(τ c ) T . In (42) , σ = -ηh m τ c , and D = (e -imω c + mσ e iω c ) -1 is a normalization constant chosen so that the inner product of the eigenvector (41) and adjoint eigenvector (42) satisfy q * , q = 1. Table 6 Table of b 0 (Y n , Y n ) and c 0 (Y n , Y n , Y n ) formulas for HIV and ELM models b 0 (Y n , Y n ) HIV

dnn -2εh m τ y n y n -(γ 2 + γ )νh m τ y n y n nnd, dnd -εh m τ (y n y n-m + y n-m y n ) -νh m τ [(γ 2 -γ )y n-m y n-m + γ (y n y n-m + y n-m y n )] ndd, ddd -2εh m τ y n-m y n-m -νh m τ (γ 2 + γ )y n-m y n-m

dnn -(γ 3 -γ )ξ h m τ y n y n y n nnd, dnd -ξ h m τ [(γ 3 -3γ 2 + 2γ )y n-m y n-m y n-m + (γ 2 -γ )(y n y n-m y n-m + y n-m y n y n-m + y n y n-m y n-m )]

ndd, ddd -ξ h m τ (γ 3 -γ )y n-m y n-m y n-m

The formulas for the second-and third-order terms in Eq. (39) are given by

where the b 0 and c 0 terms are second-and third-order partial derivatives, respectively, of the right-hand sides of the system of m + 1 first-order nonlinear difference equations corresponding to the nonlinear difference equations in Table 5 . The values of the b 0 and c 0 terms for the five versions of the HIV and ELM models are shown in Table 6 . Following the algorithms in Li [19] and Kuznetsov [24] , we concentrate on the expression for the critical coefficient c 1 (τ c ),

where g 02 = q * , B(q, q) = Db 0 (q, q), , q) ,

p e 2iω c -D e 2iω ce iω c q -D e 2iω ce -iω c q , g 21 = q * , B(q, ω 20 ) + 2 q * , B(q, ω 11 ) + q * , C(q, q, q)

= Db 0 (q, ω 20 ) + 2Db 0 (q, ω 11 ) + Dc 0 (q, q, q), 

The formulas for the b 0 , c 0 and inner products required to compute the terms in (46) are shown in Table 7 . The values of the characteristic polynomial P 1 (λ) (Eq. (20) ) required in (46) are P 1 (1) = h m τ (ηρ) and P 1 (e 2iω c ) = e 2i(m+1)ω ce 2imω c + h m τ (ηρe 2imω c ). The values of p(1) and p(λ 2 c ) can be obtained from (48).

In this section, we present results of numerical simulations to illustrate the analytical results obtained in previous sections. We show the following. Table 10 Values of equilibrium points and R 0 of HIV and ELM models from Tables 8 and 9 Parameters Disease-free Endemic 

The parameter values used in the HIV simulations are given in Table 8 and were selected from [46] . The parameter values used in the ELM simulations are shown in Table 9 and were selected from [1] .

As shown previously, the equilibrium populations and values of R 0 for the differential equation and difference equation models are the same. The equilibrium population and R 0 values for the data sets in Tables 8 and 9 are shown in Table 10 . Therefore, for zero time delay for both the differential equation and the difference equation models, the disease-free equilibrium is locally asymptotically stable for HIV set 3 of Table 8 and the endemic equilibria exist and are locally asymptotically stable for HIV sets 1 and 2 of Table 8 and ELM set 1 of Table 9 .

For the parameter values in Tables 8 and 9 , critical delay conditions for Andronov-Hopf bifurcation (Theorem 1) and Neimark-Sacker bifurcation (Theorem 3) of the endemic equilibria for the five HIV and ELM versions are shown in Tables 11 and 12 . It can be seen that the results for the HIV and ELM models are similar. For the HIV model, the Andronov-Hopf and Neimark-Sacker bifurcations occur for the bounded region 1 < R 0 < 3 for version 1, and for unbounded regions R 0 > 1 for versions 3, 5, 7, and for R 0 > 3 for version 6. Similarly, for the ELM model, the Andronov-Hopf and Neimark-Sacker bifurcations occur for the bounded region 1 < R 0 < 2 for version 1, and for unbounded regions R 0 > 1 for versions 3, 5, 7, and for R 0 > 2 for version 6. It can also be seen from these tables, that, for version 1 of both models, the C5 conditions of Theorem 3 correspond to a repelling invariant closed curve for τ > τ c , whereas for versions 3, 5, 6, 7 the C5 conditions correspond to an attracting invariant closed curve for τ > τ c .

Examples of plots of the time dependence of the solutions for time delays τ < τ c are shown in Fig. 3 for the differential equations in Table 1 and in Fig. 4 for the difference Table 2 . Examples of plots of the time dependence of the solutions for time delays τ > τ c are shown in Fig. 5 for the differential equations in Table 1 and in Fig. 6 for the difference equations in Table 2 . The figures show plots for HIV and ELM versions 1 and 6. The plots for HIV and ELM models are qualitatively similar, but slightly different in detail. As noted in Sect. 8.2, the bifurcation occurs as τ increases through τ c and the invariant closed curve is a repelling curve for version 1 and an attracting curve for version 6.

As shown in Tables 4, 11 Table 11 and 12) 

From Tables 11 and 12 , we have for condition (C3) of the Andronov-Hopf theorem that the values for the derivatives of the real part of the eigenvalues at the critical point dμ dτ | τ c are greater than zero for all versions. Therefore, the Andronov-Hopf bifurcations will occur as τ increases through τ c for all versions. It can also be seen from Tables 11 and 12 that case (1) of the condition (C5) of the Neimark-Sacker theorem occurs for version 1 (dnn) of the HIV and ELM discrete models and that case (2) occurs for versions 3 (nnd), 5 (dnd), 6 (ndd) and 7 (ddd). Therefore the Neimark-Sacker bifurcations will occur as τ increases through τ c for all versions, but version 1 will have a repelling limit cycle and versions 3, 5, 6, 7 will have attracting limit cycles.

A comparison of the critical delay values τ c for the Andronov-Hopf and Neimark-Sacker bifurcations are shown in Fig. 8 . It can be seen that the Neimark-Sacker values tend to the Andronov-Hopf values as the value of m increases, i.e., as the step size h = 1/m in the Euler difference equation approximation for the differential equation is reduced. These results suggest that difference equation models based on the Euler approximation can be used to obtain good estimates for critical delay values for models of the type studied in this paper. 

The effects of increasing the antiretroviral therapy factor n av in (2) are shown in Fig. 9 . Figure 9 (a) shows the reduction in the basic reproduction number, Fig. 9 (b) shows the effect on the equilibrium infected population and Fig. 9 (c) shows the effect on the critical Andronov-Hopf bifurcation point. In practise, as stated in the introduction, it is well known (see, e.g., [29, 30] ) that antiretroviral therapy cannot completely eliminate the virus. However, recent studies (see, e.g., [39, 41, 42] ) have suggested that the therapy can reduce the virus sufficiently that HIV transmission from an HIV+ to an uninfected person will not occur.

An important idea in controlling diseases through vaccination and antiretroviral treatment is that of "herd immunity", i.e., the level of immunity required to make the diseasefree equilibrium stable and the endemic equilibrium negative. From Eqs. (2) and (7), we obtain R 0 = ε δ = β 0 (1n av )Cα 1 -(1n av )p 0 B = 1 for n av = 1 -

For sets 1 and 2 of Table 8 analytically and confirmed numerically that one version for each model has a supercritical repelling limit cycle and the remaining four versions have supercritical attracting limit cycles. We have also shown numerically that the same qualitative behavior occurs for the Andronov-Hopf bifurcations.

We have carried out numerical simulations for a range of biologically reasonable parameter values and obtained results that agree with the analytical results. The numerical results have shown that the equilibrium points and bifurcation behavior of an Euler approximation to a differential equation system have the same qualitative behavior as the differential equation system and that the critical bifurcation points of the difference equation

converge to the critical times of the differential equation as the number of discretization points is increased.

The results of this paper show that a difference equation approximation can be used to study the dynamical behavior of time-delayed differential equation systems and to give reliable information about the existence and properties of bifurcations.

Stability and Neimark-Sacker bifurcation of numerical discretization of delay differential equations

Neimark-Sacker bifurcation in a discrete-time financial system

Stability and Hopf bifurcation of a HIV infection model with CTL-response delay

Bifurcation analysis of a delay differential equation model associated with the induction of long-term memory

Neimark-Sacker bifurcation of a two-dimensional discrete-time predator-prey model

Stability and Hopf bifurcation of a predator-prey model with stage structure and time delay for the prey

Andronov-Hopf bifurcation and sensitivity analysis of a time-delay HIV model with logistic growth and antiretroviral treatment

Stability and Hopf bifurcation of a within-host Chikungunya virus infection model with two delays

On the difference equation x n+1 = ax n-l + bx n-k + f (x n-l , x n-k )

Stability, bifurcation analysis and chaos control in chlorine dioxide-iodine-malonic acid reaction

Delayed loss of stability due to the slow passage through Hopf bifurcations in reaction-diffusion equations

Analysis of bifurcation, chaos and pattern formation in a discrete time and space Gierer Meinhardt system

Neimark-Sacker bifurcation and stability analysis of a discrete-time prey-predator model with Allee effect in prey

Stability, bifurcation analysis and chaos control for a predator-prey system

Chaotic behavior of predator-prey model with group defense and non-linear harvesting in prey

Stability and Hopf bifurcation of a predator-prey model

Bifurcations and chaotic-type solutions in a time-delay model of hematopoietic stem cell growth

A bifurcation path to chaos in a time-delay fisheries predator-prey model with prey consumption by immature and mature predators

Dynamics of a discrete food-limited population model with time delay

Neimark-Sacker bifurcations in the Euler method for a delay differential equation

Neimark-Sacker bifurcation analysis on a numerical discretization of Gause-type predator-prey model with delay

Bifurcations of a two-dimensional discrete time plant-herbivore system

Numerical Solution of SDE Through Computer Experiments

Elements of Applied Bifurcation Theory

The asymptotic behaviour of a logistic epidemic model with stochastic disease transmission

Asymptotic behavior and stability of a stochastic model for AIDS transmission

An evaluation of the extended logistic, simple logistic, and Gompertz models for forecasting short lifecycle products and services

Extended logistic model for growth of single-species populations

Modeling latently infected cell activation: viral and latent reservoir persistence, and viral blips in HIV-infected patients on potent therapy

Modeling HIV persistence, the latent reservoir, and viral blips

Viral dynamics of an HIV model with latent infection incorporating antiretroviral therapy

HIV-1 infection and low steady state viral loads

Asymmetric division of activated latently infected cells may explain the decay kinetics of the HIV-1 latent reservoir and intermittent viral blips

Life expectancy of HIV-positive adults: a review

Andronov-Hopf and Neimark-Sacker bifurcations in time-delay models

Prospects for treatment of latent HIV

Wikipedia: protease inhibitor (pharmacology). Wikipedia

Emedicine: antiretroviral therapy for HIV infection

Unaid: global-AIDS-update2016

Wikipedia: reverse-transcriptase inhibitor

A farewell to condoms

World Health Organization: AIDS. World Health Organization

Numerical Analysis. Cengage Learning

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Introduction to Dynamic Systems: Theory, Models and Applications

Stability analysis of an HIV/AIDS epidemic with treatment

no contribution to the design of the study and collection, analysis and interpretation of any data or in the writing or editing of the manuscript.

A published data set from Cao et al. [46] was used in the simulation of the HIV model and a published data set from He et al. [1] was used in the numerical simulation of the extended logistic model.

The authors declare that they have no competing interests.

The authors made equal contributions to the research. All authors read and approved the final manuscript. 

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 4 February 2020 Accepted: 17 April 2020