Andronov–Hopf and Neimark–Sacker bifurcations in time-delay differential equations and difference equations with applications to models for diseases and animal populations

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.


Introduction
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.

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

HIV models
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).

Extended logistic growth models
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.

Time-delay differential equation models
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.

Time-delay difference equation models
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.

Equilibrium points
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.

Local asymptotic stability
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.

Andronov-Hopf bifurcations of time-delay models
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.

Lemma 1 Necessary and sufficient conditions for existence of purely imaginary solutions
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.
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 .

Equilibrium points, stability and Neimark-Sacker bifurcations of difference equation models 6.1 Equilibrium points and basic reproduction numbers
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).

Conditions for Neimark-Sacker bifurcations
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.

Lemma 4
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.
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 < ω < π .
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.

Lemma 6
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.
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.

Direction and stability of the Neimark-Sacker bifurcations
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 second-and third-order terms
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).

Numerical simulations
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].

Comparison of qualitative behavior of versions
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.

Andronov-Hopf and Neimark-Sacker bifurcations
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)

Comparison of critical values for Andronov-Hopf and Neimark-Sacker bifurcations
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.

Effects of antiretroviral therapy
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.