Complex Dynamics of Discrete SEIS Models with Simple Demography

Abstract

We investigate bifurcations and dynamical behaviors of discrete SEIS models with exogenous reinfections and a variety of treatment strategies. Bifurcations identified from the models include period doubling, backward, forward-backward, and multiple backward bifurcations. Multiple attractors, such as bistability and tristability, are observed. We also estimate the ultimate boundary of the infected regardless of initial status. Our rigorously mathematical analysis together with numerical simulations show that epidemiological factors alone can generate complex dynamics, though demographic factors only support simple equilibrium dynamics. Our model analysis supports and urges to treat a fixed percentage of exposed individuals.

1. Introduction

Pure demographic (or ecological) discrete models can have very complicated dynamics. It is well known that the logistic model 𝑥_𝑛+1=𝑟𝑥_𝑛(1−𝑥_𝑛/𝐾), Ricker's model 𝑥_𝑛+1=𝑥_𝑛𝑒^{𝑟−𝑥𝑛/𝐾} [1–3], and Hassell model 𝑥_𝑛+1=𝜆𝑥_𝑛(1+𝑎𝑥_𝑛)^−𝑏 [4] all have complex dynamics through the mechanics of period-doubling bifurcation. When adding epidemiological effects to these demographic models to describe the dynamics of an infectious disease, it is not surprising at all that the complex dynamics still retain [5, 6]. Yet an interesting question is that if adding epidemiological effects to a simple demographic model that has only trivial dynamics, for example, 𝑥_𝑛+1=Λ+𝜇𝑥_𝑛, is it possible to observe complex dynamics? Our work in this paper finds a quite positive answer by bifurcation analysis. It is shown that the chaotic dynamics is possible.

The bifurcation approach has been extensively utilized in theoretical epidemiology. The use of equilibrium bifurcation to investigate the dynamical behavior of continuous epidemic models has been successful, for instance, the work in [7–15]. Meanwhile, discrete epidemic models have gained more popularity [5, 6, 16–25], since epidemiological data are usually collected at discrete times and it becomes easier to compare data with models. We will use the bifurcation approach to analyze discrete SEIS epidemic models in this paper.

The appearance of multiple attractors enhances uncertainty of outcomes because of the lack of information to the initial epidemiological status, even though the domain of attractions of each attractor is fully depicted. In this case, it would be necessary to estimate the prevalence regardless of the initial data.

The management measures to completely wipe out an infectious disease may be too costly, far beyond the capability of any society. A practical and feasible way is to have an infectious disease under control, that is, to keep the number of the infectious individuals under a certain level. Thus, estimating an upper bound of the prevalence appears necessary. A large pool of excellent researches have focused on the lower bound of prevalence in the application of permanence theory [26–29], while there are relatively few studies on the upper boundary of infections when a disease persists, especially when multiple attractors exist. In this paper, we also estimate ultimate boundaries, which are given in terms of the basic reproductive number if the disease persists. We believe, in some sense, that an ultimate boundary estimation for an infectious disease is more significant than a basic reproductive number, especially when multiple attractors appear.

Our mathematical models are presented in Section 2. In Section 3, we compute the reproductive number and simplify our models. In Section 4, we analyze the different bifurcations of the models under different treatment strategies and with or without considering exogenous reinfections. Numerical simulations are carried out to test and confirm each type of bifurcations. The ultimate boundaries of the latent and the infectious are obtained under some sufficient conditions. Section 5 includes some concluding remarks and discussions.

2. The Model

Our models could be built within the context of the dynamical and epidemiological study of tuberculosis, where both the exogenous reinfections and the different treatment strategies play an important role. Existing studies have shown that the exogenous re-infection of the latent individuals increases the risk of developing into the active disease [20, 30, 31]. The effect of the treatment for the latent individuals on the dynamics of disease has been investigated using continuous models [32, 33]. Here, we will investigate the treatment effect of the latent individuals on the dynamics of the disease using discrete epidemic models. The models will incorporate two distinct features: exogenous re-infection and flexible treatment strategies.

In modeling the treatment for the latent and the infectious individuals, we should consider both the prevalence of the disease and availability of a budget. The treatment rate for infectious ones is rather easy; we simply take a linear function, as has been typically used. However, a practice of treating the latent individuals, like inactive cases of tuberculosis, is highly theoretic, mainly because it is hard to find asymptomatic individuals without mass-screening. When the pool of the latent individuals is small, it is reasonable to assume that the treatment rate is proportional to the number of the latent individuals (a linear function of the number of the latent individuals). However, if the pool is big, such as the case of tuberculosis in India, China, Eastern Europe, or South Sahara, taking a linear function as a treatment rate would be inappropriate due to the limited budgets. We cannot provide treatment to all individuals with latent tuberculosis, since the people living with latent tuberculosis are approximately one-third of the total population. The other obvious reason is our inability of identifying asymptomatic individuals. We have to design the treatment strategy within our reach of case finding and budgets. It then becomes more plausible to assume a constant treatment rate when the pool of latent individuals is large, as has been used in [13, 15], for instance. Overall, treatment strategies should be designed by considering both the prevalence of the disease and available budget. Part of our goal in this paper is to compare these strategies.

The total population is epidemiologically divided into classes of susceptible, latent, and infectious. Let 𝑆(𝑡) the number of susceptible individuals at time 𝑡, 𝐸(𝑡) be the number of latent individuals, and 𝐼(𝑡) the number of infectious individuals. 𝑁(𝑡)=𝑆(𝑡)+𝐸(𝑡)+𝐼(𝑡) is the total population size. We assume that recovered individuals by successful treatment do not acquire immunity and will become member of susceptible compartment. Hence our model is of SEIS type.

Susceptible individuals may get infected and enter the latent compartment. The survived latent individuals will experience one of the three mutually exclusive events to leave the latent compartment, getting exogenous re-infection through contacting with infectious individuals, or naturally progressing into infectious class, or receiving the treatment. These three events happen randomly. The respective probabilities for exogenous re-infection, natural progression, and receiving the treatment are 𝑘₁, 𝑘₂, and 𝑘₃ with 𝑘₁+𝑘₂+𝑘₃=1. Incorporating the exogenous reinfections and treatments into the models in [20, 21], our model framework takes the following form:𝑆(𝑡+1)=Λ+𝑝𝑆(𝑡)𝐺₁𝐼(𝑡)𝑁(𝑡)(𝐸(𝑡(𝑡),𝐸()(𝑡+𝑝𝐻))+𝑝𝛾𝐼𝑡+1=𝑝𝑆)1−𝐺₁𝐼(𝑡)𝑁(𝑡)+𝑝𝑘₁𝐸(𝑡)𝐺₂𝐼(𝑡)𝑁(𝑡)+𝑝𝑘₂(1−𝛼)𝐸(𝑡)+𝑝𝑘₃𝐸(𝑡)(𝐸(𝑡,𝐼()−𝐻))𝑡+1=𝑝𝑘₁𝐸(𝑡)1−𝐺₂𝐼(𝑡)𝑁(𝑡)+𝑝𝛼𝑘₂𝐸(𝑡)+𝑝1−𝛾𝐼(𝑡), (2.1) where Λ is the recruitment rate into the population, 𝑝 is the survival probability, and 𝛾 is the probability that an infectious individual recovers successfully. 𝛼 is the conditional probability that a latent individual becomes infectious successfully given that the natural progression happens. 1−𝐺₁ is the probability that susceptible individuals become infected. 1−𝐺₂ is the conditional probability that a latent individual gets re-infected successfully given that the exogenous re-infection happens. The probabilities of not becoming infected are 𝐺₁ and 𝐺₂ that are considered as functions of the prevalence 𝐼(𝑡)/𝑁(𝑡). exp(−𝛽_𝑗𝐼(𝑡)/𝑁(𝑡)) and 1−(𝛽_𝑗𝐼(𝑡)/𝑁(𝑡)) have been used for the probabilities of not becoming infected [20, 21]. We choose 𝐺_𝑗(𝐼(𝑡)/𝑁(𝑡))=1−(𝛽_𝑗𝐼(𝑡)/𝑁(𝑡)),𝑗=1,2, where 𝛽_𝑗 characterizes the disease transmission probability.

In model (2.1), 𝐻(𝐸(𝑡)) is the successful treatment rate of the latent individuals and 𝑘₃𝐸(𝑡)−𝐻(𝐸(𝑡)) is the failure rate of the treatment. We consider two different types of treatment strategies for the latent individuals:𝐻(𝐸(𝑡))=𝑚𝑘₃𝐸(𝑡)or𝐻(𝐸(𝑡))=𝑚𝑘₃𝐸(𝑡),𝐸(𝑡)≤𝐸⁰,𝐸⁰,𝐸(𝑡)>𝐸⁰, (2.2) where 𝑚 is the conditional probability that a latent individual is treated successfully given that the individual receives the treatment and 𝐸⁰ is the cutting point of treatment determined by policy makers. Both 𝑚 and 𝐸⁰ should rely on the capability of finding latent individuals and budgetary issues.

We summarize some important terms in model (2.1) here. From time 𝑡 to time 𝑡+1, the number of the susceptibles who get infected and become latent individuals is 𝑝𝑆(𝑡)(1−𝐺₁(𝐼(𝑡)/𝑁(𝑡))); the number of the latent who get re-infected and become infectious is 𝑝𝑘₁𝐸(𝑡)(1−𝐺₂(𝐼(𝑡)/𝑁(𝑡))); the number of the latent who naturally progress into the infectious compartment is 𝑝𝑘₂𝛼𝐸(𝑡); the number of the latent who get recovered by successful treatment and become susceptible is 𝑝𝐻(𝐸(𝑡)); the number of the infectious who get recovered and become susceptible is 𝑝𝛾𝐼(𝑡). The number of the latent individuals at time 𝑡 that remain as latent ones at time 𝑡+1 is 𝑝𝑘₁𝐸(𝑡)𝐺₂(𝐼(𝑡)/𝑁(𝑡))+𝑝𝑘₂(1−𝛼)𝐸(𝑡)+𝑝(𝑘₃𝐸(𝑡)−𝐻(𝐸(𝑡))). Similarly, the number of infectious individuals who do not change their epidemiological status from 𝑡 to 𝑡+1 is counted by 𝑝(1−𝛾)𝐼(𝑡).

3. Basic Reproductive Number and Reduced Model

Firstly we explicitly write our models for two specific treatment regimens. In the case of proportional treatment rate of the latent individuals, the treatment function is linear, 𝐻(𝐸(𝑡))=𝑚𝑘₃𝐸(𝑡), and model (2.1) can be rewritten as follows:𝑆()(𝑡)𝛽𝑡+1=Λ+𝑝𝑆1−¹𝐼(𝑡)𝑁(𝑡)+𝑝𝑚𝑘₃𝐸(𝑡)(𝑡),𝐸()(𝑡)𝛽+𝑝𝛾𝐼𝑡+1=𝑝𝑆¹𝐼(𝑡)𝑁(𝑡)+𝑝𝑘₁𝐸(𝑡)𝛽1−²𝐼(𝑡)𝑁(𝑡)+𝑝(1−𝛼)𝑘₂𝐸(𝑡)+𝑝(1−𝑚)𝑘₃𝐸(𝑡),𝐼(𝑡+1)=𝑝𝑘₁𝐸(𝑡)𝛽₂𝐼(𝑡)𝑁(𝑡)+𝑝𝛼𝑘₂𝐸(𝑡)+𝑝1−𝛾𝐼(𝑡). (3.1) On the other hand, in the case of limited resources, the treatment function of latent individuals is a piecewise function, 𝐻(𝐸(𝑡))=𝑚𝑘₃𝐸(𝑡),𝐸(𝑡)≤𝐸⁰,𝑚₀,𝐸(𝑡)>𝐸⁰, (3.2) where 𝑚₀=𝑚𝑘₃𝐸⁰. Model (2.1) takes the following form:𝑆()(𝑡)𝛽𝑡+1=Λ+𝑝𝑆1−¹𝐼(𝑡)𝑁(𝑡)+𝑝𝑘₃𝐸(𝑡)(𝐸(𝑡(𝑡),𝐸()(𝑡)𝛽−𝑇))+𝑝𝛾𝐼𝑡+1=𝑝𝑆¹𝐼(𝑡)𝑁(𝑡)+𝑝𝑘₁𝐸(𝑡)𝛽1−²𝐼(𝑡)𝑁(𝑡)+𝑝(1−𝛼)𝑘₂𝐸(𝑡)(𝐸(𝑡,𝐼()+𝑇))𝑡+1=𝑝𝑘₁𝐸(𝑡)𝛽₂𝐼(𝑡)𝑁(𝑡)+𝑝𝛼𝑘₂𝐸(𝑡)+𝑝1−𝛾𝐼(𝑡), (3.3) where𝑇(𝐸(𝑡=𝑝()𝑘))1−𝑚₃𝐸(𝑡),if𝐸(𝑡)≤𝐸⁰,𝑝𝑘₃𝐸(𝑡)−𝑚₀,if𝐸(𝑡)>𝐸⁰. (3.4)

The disease-free equilibrium, 𝑃₀(Λ/(1−𝑝),0,0), and the basic reproductive number of models (3.1) and (3.3) are identical. Applying the approach in [34] to our model, we have(𝑡)𝛽ℱ=𝑝𝑆¹𝐼(𝑡)𝑁(𝑡)𝑝𝑘₁𝐸(𝑡)𝛽₂𝐼(𝑡)𝑁(𝑡)0𝑝,𝒱=1−𝑘₁𝛽₂𝐼(𝑡)𝑁(𝑡)−𝛼𝑘₂−𝑚𝑘₃𝐸(𝑡)𝑝𝛼𝑘₂𝐸(𝑡)𝐼(𝑡)(𝑡)𝛽+𝑝1−𝛾Λ+𝑝𝑆1−¹𝐼(𝑡)𝑁(𝑡)+𝑝𝑚𝑘₃𝐸(𝑡)+𝑝𝛾𝐼(𝑡). (3.5) At the disease-free equilibrium, the straight forward calculation yields 𝐹=^0𝑝𝛽1₀₀ and 𝑉=^{𝑝−𝑝𝛼𝑘2−𝑝𝑚𝑘30}_{𝑝𝛼𝑘2𝑝(1−𝛾)}. Therefore, the basic reproductive number of model (2.1) is𝑅₀=𝑝²𝛽₁𝛼𝑘₂1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃. (3.6) Rearrange terms in the expression of 𝑅₀ as𝑅₀=𝑝𝛽₁1−𝑝1−𝛾𝑝𝛼𝑘₂1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃. (3.7) Each term in 𝑅₀ has clear epidemiological interpretation. 1/(1−𝑝(1−𝛾)) is the average infection period. 𝑝𝛽₁/(1−𝑝(1−𝛾)) is average new cases generated by a typical infectious member in the entire infection period. 𝑝𝛼𝑘₂/(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃) is the proportion that latent individuals become infectious by “natural" progression.

It is known that the disease-free equilibrium 𝑃₀ of (2.1) is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1 (see [34, Theorem  2.1]). We summarize stability results in terms of 𝑅₀<1 as in the following theorem.

Theorem 3.1. The disease-free equilibrium 𝑃₀ of model (2.1) is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

Since we do not consider the disease-induced death rate, the total population size is governed by an extremely simple equation. In fact, by adding all equations in model (2.1), we obtain the following equation for the total population 𝑁(𝑡)=𝑆(𝑡)+𝐸(𝑡)+𝐼(𝑡):𝑁(𝑡+1)=Λ+𝑝𝑁(𝑡). (3.8) One can see that 𝑁^∗=Λ/(1−𝑝) is a global attractor for (3.8). Using limiting equations and the lim_𝑡→∞𝑁(𝑡)=𝑁^∗ [35], we, respectively, reduce the three-dimensional systems (3.1) and (3.3) into two-dimensional ones𝐸(𝑡+1)=𝑝𝑁^∗(𝑡)(𝑡)𝛽−𝐸−𝐼¹𝐼(𝑡)𝑁^∗+𝑝𝑘₁𝐸(𝑡)𝛽1−²𝐼(𝑡)𝑁^∗+𝑝(1−𝛼)𝑘₂𝐸(𝑡)+𝑝(1−𝑚)𝑘₃𝐸(𝑡),𝐼(𝑡+1)=𝑝𝑘₁𝐸(𝑡)𝛽₂𝐼(𝑡)𝑁^∗+𝑝𝛼𝑘₂𝐸(𝑡)𝐼(𝑡),𝐸()𝑁+𝑝1−𝛾𝑡+1=𝑝^∗(𝑡)(𝑡)𝛽−𝐸−𝐼¹𝐼(𝑡)𝑁^∗+𝑝𝑘₁𝐸(𝑡)𝛽1−²𝐼(𝑡)𝑁^∗+𝑝(1−𝛼)𝑘₂𝐸(𝑡)(𝐸(𝑡,𝐼()+𝑇))𝑡+1=𝑝𝑘₁𝐸(𝑡)𝛽₂𝐼(𝑡)𝑁^∗+𝑝𝛼𝑘₂𝐸(𝑡)+𝑝1−𝛾𝐼(𝑡). (3.9)(3.10) In the rest of the paper, we analyze the dynamic behavior of the limiting system (3.9) and (3.10) under different treatment strategies.

4. Analysis

In this section, we study the bifurcation and stability of equilibrium of model (3.9) and (3.10). Numerical simulations are also presented to demonstrate the theoretical results.

4.1. Without Exogenous Reinfections

The SEIS model is the simplest one if there is no exogenous re-infection and the treatment rate takes a linear form. For this case, 𝑘₁=0, 𝐻(𝐸)=𝑚𝑘₃𝐸, and model (3.9) is reduced to𝐸(𝑡+1)=𝑝𝑁^∗(𝑡)(𝑡)𝛽−𝐸−𝐼¹𝐼(𝑡)𝑁^∗+𝑝(1−𝛼)𝑘₂𝐸(𝑡)+𝑝(1−𝑚)𝑘₃𝐸(𝑡),𝐼(𝑡+1)=𝑝𝛼𝑘₂𝐸(𝑡)+𝑝1−𝛾𝐼(𝑡). (4.1) The disease-free equilibrium of (4.1) is 𝑃⁰₁(𝐸⁰₁,𝐼⁰₁)=(0,0) and the endemic equilibrium 𝑃^∗₁(𝐸^∗₁,𝐼^∗₁) exists if 𝑅₀>1, where 𝐸^∗₁=((1−𝑝(1−𝛾))𝑁^∗/(1−𝑝(1−𝛾)+𝑝𝑘₂𝛼))(1−(1/𝑅₀)) and 𝐼^∗₁=𝑝𝑘₂𝛼𝐸^∗₁/(1−𝑝(1−𝛾)).

4.1.1. Global Stability of 𝑃⁰₁

The global stability of the disease-free equilibrium 𝑃⁰₁(𝐸⁰₁,𝐼⁰₁)=(0,0) of (4.1) is given in Theorem 4.1.

Theorem 4.1. If 𝑅₀<1, then the disease-free equilibrium 𝑃⁰₁ of (4.1) is globally asymptotically stable.

Proof. Define 𝐹∶[0,𝑁^∗]×[0,𝑁^∗]→[0,𝑁^∗]×[0,𝑁^∗] by 𝐹(𝐸,𝐼)=𝑝𝑁^∗𝛽−𝐸−𝐼¹𝐼𝑁^∗+(1−𝛼)𝑘₂𝐸+(1−𝑚)𝑘₃𝐸,𝛼𝑘₂𝐸+1−𝛾𝐼. (4.2) Obviously, 𝐹 is the mapping derived by system (4.1), and (0,0) is a fixed point of 𝐹. Linear function 𝑉(𝐸,𝐼)=𝐸+(𝑝𝛽₁/(1−𝑝(1−𝛾)))𝐼 on [0,𝑁^∗]×[0,𝑁^∗] is continuous and positive definite with respect to (0,0). Therefore, 𝑉 is a Lyapunov function on the domain of 𝐹. For any (𝐸,𝐼)∈[0,𝑁^∗]×[0,𝑁^∗], Δ𝑉(𝐸,𝐼)=𝑝𝑁^∗𝛽−𝐸−𝐼¹𝐼𝑁^∗+𝑝(1−𝛼)𝑘₂𝐸+𝑝(1−𝑚)𝑘₃𝐸+𝑝𝛽₁1−𝑝1−𝛾𝑝𝛼𝑘₂𝐸+𝑝1−𝛾𝐼−𝐸−𝑝𝛽₁𝐼1−𝑝1−𝛾≤𝑝𝛽₁𝐼−1−𝑝(1−𝛼)𝑘₂−𝑝(1−𝑚)𝑘₃𝑝𝐸+²𝛽₁𝛼𝑘₂𝐸1−𝑝1−𝛾+𝑝𝛽₁𝐼=1−𝑝(1−𝛼)𝑘₂−𝑝(1−𝑚)𝑘₃𝑅₀−1𝐸≤0. (4.3) Hence, if 𝑅₀<1, then Δ𝑉(𝐸,𝐼)<0 holds for (𝐸,𝐼)∈(0,𝑁^∗]×(0,𝑁^∗]. It follows from Theorem 4.22 in [36] that 𝑃⁰₁ is globally asymptotically stable.

Theorem 4.1 ensures that the disease-free equilibrium 𝑃₀(Λ/(1−𝑝),0,0) of model (3.1) is globally stable as 𝑅₀<1 and 𝑘₁=0.

4.1.2. Stability of 𝑃^∗₁

The stability of the endemic equilibrium is given in Theorem 4.2.

Theorem 4.2. If 1<𝑅₀<1+𝜎, then the endemic equilibrium 𝑃^∗₁ of (4.1) is asymptotically stable, where 𝜎=2𝑝1−𝛼𝑘₂−𝑚𝑘₃+1−𝛾1−𝑝1−𝛾+𝑝𝛼𝑘₂1+𝑝1−𝛾−𝑝𝛼𝑘₂1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃>0. (4.4)

Proof. The linearization matrix of (4.1) at the endemic equilibrium 𝑃^∗₁ is 𝐽₁=−𝑝𝛽₁𝑁^∗𝐼^∗₁+𝑝(1−𝛼)𝑘₂+𝑝(1−𝑚)𝑘₃𝑝𝛽₁𝑁^∗𝑁^∗−𝐸^∗₁−2𝐼^∗₁𝑝𝛼𝑘₂𝑝1−𝛾. (4.5) The corresponding characteristic equation is 𝑓₁(𝜆)=𝜆²−𝑏₁𝜆+𝑏₂=0, (4.6) where 𝑏₁=𝑝−𝑝𝛼𝑘₂−𝑝𝑚𝑘₃+𝑝1−𝛾−1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑅₀−11−𝑝1−𝛾+𝑝𝛼𝑘₂,𝑏₂=𝑝1−𝛾𝑝−𝑝𝛼𝑘₂−𝑝𝑚𝑘₃−1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃+𝑝𝛼𝑘₂−𝑝1−𝛾1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑅₀−11−𝑝1−𝛾+𝑝𝛼𝑘₂. (4.7) We claim |𝑏₂|<1. When 𝑝𝛼𝑘₂=𝑝(1−𝛾), we have 𝑏₂=𝑝(1−𝛾)(𝑝−𝑝𝛼𝑘₂−𝑝𝑚𝑘₃)−(1−𝑝(1−𝛾))(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃), and |𝑏₂|<1 holds. When 𝑝𝛼𝑘₂≠𝑝(1−𝛾), because of 2𝑝1−𝛼𝑘₂−𝑚𝑘₃+1−𝛾1+𝑝1−𝛾−𝑝𝛼𝑘₂<min2−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃−𝑝1−𝛾𝑝𝛼𝑘₂,𝑝−𝑝1−𝛾1−𝛼𝑘₂−𝑚𝑘₃+1−𝛾𝑝1−𝛾−𝑝𝛼𝑘₂, (4.8) we obtain that |𝑏₂|<1 holds if 1<𝑅₀<1+𝜎. Furthermore, we have 𝑓₁(1)=1−𝑏₁+𝑏₂=1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑅₀,()−1−1²𝑓₁(−1)=1+𝑏₁+𝑏₂=2𝑝−𝑝𝛼𝑘₂−𝑝𝑚𝑘₃−+𝑝1−𝛾1+𝑝1−𝛾−𝑝𝛼𝑘₂1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑅₀−11−𝑝1−𝛾+𝑝𝛼𝑘₂. (4.9) It is clear that 𝑓₁(1)>0 holds if 𝑅₀>1 and (−1)²𝑓₁()=𝜎−11+𝜎−𝑅₀1+𝑝1−𝛾−𝑝𝛼𝑘₂1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃>0 (4.10) holds if 𝑅₀<1+𝜎. Therefore, if 1<𝑅₀<1+𝜎, we have 𝑓₁(1)>0, (−1)²𝑓₁(−1)>0, and |𝑏₂|<1. It follows from the Jury criteria that the endemic equilibrium 𝑃^∗₁(𝐸^∗₁,𝐼^∗₁) of (4.1) is asymptotically stable.

It follows from Theorem 4.2 that when exogenous re-infection is not included (𝑘₁𝛽₂=0) the endemic equilibrium 𝑃^∗₁(Λ/(1−𝑝))−𝐸^∗₁−𝐼^∗₁,𝐸^∗₁,𝐼^∗₁) of model (3.1) is asymptotically stable if 1<𝑅₀<1+𝜎. Since the Jury criteria is the necessary and sufficient condition for the local stability of the equilibrium, the endemic equilibrium 𝑃^∗₁ of (3.1) becomes unstable if 𝑅₀≥1+𝜎. We summarize Theorems 4.1 and 4.2 for model (3.1) without exogenous reinfection as follows. (1)The disease-free equilibrium 𝑃₀(Λ/(1−𝑝),0,0) is globally stable as 𝑅₀<1.(2)When 𝑅₀>1, 𝑃₀(Λ/(1−𝑝),0,0) becomes unstable and endemic equilibrium 𝑃^∗₁ is born. The endemic 𝑃^∗₁ is locally asymptotically stable if 1<𝑅₀<1+𝜎.(3)When 𝑅₀≥1+𝜎,𝑃^∗₁ becomes unstable.

Theorems 4.1 and 4.2 state that a forward bifurcation takes place at 𝑅₀=1. This is consistent with continuous dynamical models for tuberculosis when exogenous re-infection is not considered [30, 31].

Although the expression of 𝜎 is complicated, the condition in Theorem 4.2 is easy to check since 𝜎 is independent of 𝛽₁ while 𝑅₀ depends on 𝛽₁ linearly. For example, given a set of the parameter values, 𝑝=0.994, Λ=10, 𝛼=0.03, 𝛾=0.1, 𝑚=0.1, 𝑘₁=0, 𝑘₂=0.8, and 𝑘₃=0.2, one can determine if the endemic equilibrium is stable or not. Numerical simulations suggest that the endemic equilibrium of model (3.1) may be globally asymptotically stable when 1<𝑅₀<1+𝜎, as shown in Figure 1(a) for 𝛽₁=0.245. When 𝛽₁=10.972015, we have 𝑅₀=49.63168 and 1+𝜎=49.63168. The endemic equilibrium of model (4.1) loses the stability when 𝛽₁ is greater than 10.972015, and there is a 2-period solution when 𝛽₁∈[10.972015,11.109). For example, if we take 𝛽₁=11, then 𝑅₀=49.75827, and 1+𝜎=49.63168. The 2-period solution is ((838.46721,189.77679),(1591.60886,360.24119)). Although model (4.1) has a period solution of period 2 if the endemic equilibrium of model (4.1) is unstable, the periodic solution may not be realistic since model (4.1) is the limiting system of model (3.1) as 𝑘₁𝛽₂=0, and the numerical computation shows that the corresponding periodic solution of model (3.1) is ((838.46721,189.77679,638.42267), (1591.60886,360.24119,−285.18339)). Obviously, one of 𝑆^∗ is negative, which hints that the period solution of model (3.1) can be ruled out if we confine all solutions as positive. In fact, 1−𝐺₁(𝐼(𝑡)/𝑁(𝑡))=𝛽₁𝐼(𝑡)/𝑁(𝑡) means that 𝛽₁ should not be bigger, and we should chose 𝛽₁ such that 0≤1−𝐺₁(𝐼(𝑡)/𝑁(𝑡))≤1. Otherwise, not only does 1−𝐺₁(𝐼(𝑡)/𝑁(𝑡)) lose its meaning but also the number of the susceptible individuals has negative values.

(a)

(b)

Figure 1 

The stability and the attraction of the endemic equilibrium of (3.1): (a) the endemic equilibrium and its stability and (b) the attraction area.

4.1.3. Ultimate Boundary

We now estimate an ultimate boundary to the model (4.1) when 𝑅₀>1.

Theorem 4.3. If 0<𝛽₁≤1 and 1<𝑅₀<1+(1/(1−𝑝(1−𝛾))), then the solutions of model (4.1) are bounded by limsup_𝑡→∞𝛼𝑘₂𝐸(𝑡)+1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝐼(𝑡)≤1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑁^∗11−𝑅₀,limsup_𝑡→∞𝐸(𝑡)≤1+1−𝑝+𝑝𝑚𝑘₃𝑝𝛼𝑘₂𝑁^∗11−𝑅₀,limsup_𝑡→∞𝐼(𝑡)≤𝑁^∗11−𝑅₀. (4.11)

Proof. Let us consider the value of 𝑉(𝐸(𝑡),𝐼(𝑡)) for any (𝐸(𝑡),𝐼(𝑡))∈Ω, where Ω=(𝐸(𝑡),𝐼(𝑡))∣𝐸(𝑡)≥0,𝐼(𝑡)≥0,𝐸(𝑡)+𝐼(𝑡)≤𝑁^∗. (4.12) Define 𝑉(𝐸,𝐼)=𝑝𝛼𝑘₂𝐸+1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑉𝐼,₀=1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑁^∗11−𝑅₀. (4.13) Then, (𝐸()((𝐸(𝑡)(𝑡𝑅Δ𝑉=𝑉𝑡+1,𝐼𝑡+1))−𝑉,𝐼))=𝑞₀𝑅−1−⁰𝑁^∗1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃(𝐸(𝑡)(𝑡𝐼(𝑡)𝑅+𝐼))≤𝑞₀𝑅−1−⁰𝑁^∗𝑉(𝐸(𝑡),𝐼(𝑡))1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝐼(𝑡),=1−𝑝1−𝛾𝑅₀𝑁^∗𝑉₀−𝑉(𝐸(𝑡),𝐼(𝑡))𝐼(𝑡), (4.14) where 𝑞=(1−𝑝(1−𝛾)(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃)). Therefore, Δ𝑉≤0 holds if 𝑉≥𝑉₀.
If 𝑉(𝐸(𝑡₀),𝐼(𝑡₀))=𝑝𝛼𝑘₂𝐸(𝑡₀)+(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃)𝐼(𝑡₀)≤𝑉₀ for a 𝑡₀≥0, then we claim that 𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀ for all 𝑡≥𝑡₀. In fact, let 𝑉(𝐸(𝑡₀),𝐼(𝑡₀))=𝑉₀−𝜀 with some 𝜀≥0. Then, the inequality 𝑉(𝐸(𝑡₀),𝐼(𝑡₀))=𝑉₀−𝜀≤𝑉₀ yields 𝐼𝑡₀≤𝑉₀1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃=𝑁^∗11−𝑅₀. (4.15) The inequality in (4.14) and the condition in Theorem 4.3 as well as the inequality for 𝐼(𝑡₀) imply Δ𝑉=𝑉𝐸𝑡₀+1,𝐼𝑡₀+1−𝑉𝐸𝑡₀,𝐼𝑡₀≤1−𝑝1−𝛾𝑅₀𝑁^∗𝑉₀−𝑉𝐸𝑡₀,𝐼𝑡₀𝐼𝑡₀,≤1−𝑝1−𝛾𝑅₀𝑉₀−𝑉𝐸𝑡₀,𝐼𝑡₀11−𝑅₀≤1−𝑝1−𝛾𝑅₀−1𝜀≤𝜀. (4.16) The inequality in (4.16) leads to 𝑉𝐸𝑡₀+1,𝐼𝑡₀+1≤𝑉𝐸𝑡₀,𝐼𝑡₀+𝜀=𝑉₀−𝜀+𝜀=𝑉₀. (4.17) Then, our claim is followed from the mathematical induction.
If 𝑉(𝐸(𝑡₀),𝐼(𝑡₀))=𝑝𝛼𝑘₂𝐸(𝑡₀)+(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃)𝐼(𝑡₀)>𝑉₀ for a 𝑡₀≥0, then there are two cases that we have to deal with. Case 1. There exists an integer 𝑡₁>𝑡₀ such that 𝑉(𝐸(𝑡₁),𝐼(𝑡₁)≤𝑉₀.Case 2. The inequality 𝑉(𝐸(𝑡),𝐼(𝑡))>𝑉₀ holds for all 𝑡≥𝑡₀. For Case 1 we can prove that 𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀ for all 𝑡≥𝑡₁ by the similar argument adopted for the case 𝑉(𝐸(𝑡₀),𝐼(𝑡₀))≤𝑉₀. The boundary of 𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀ for 𝑡>𝑡₁ gives rise to limsup_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀. For Case 2, the fact that Δ𝑉(𝐸(𝑡),𝐼(𝑡))≤0 for 𝑉(𝐸(𝑡),𝐼(𝑡))>𝑉₀ leads to that the 𝑉(𝐸(𝑡),𝐼(𝑡)) is a monotonic decreasing sequence. Hence, the limit lim_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡)) exits. We claim that lim_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀. Otherwise, there exists 𝑉^∗>𝑉₀, such that lim_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))=𝑉^∗>𝑉₀. It follows that there exists a subsequence 𝑡_𝑘, such that lim_𝑘→∞𝐸(𝑡_𝑘)=𝐸^∗, lim_𝑘→∞𝐼(𝑡_𝑘)=𝐼^∗, and lim_𝑘→∞𝑉(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))=𝑉(𝐸^∗,𝐼^∗)=𝑉^∗>𝑉₀. The fact that Δ𝑉(𝐸^∗,𝐼^∗)<0 and the continuity of Δ𝑉(𝐸,𝐼) imply that there exists an 𝜂>0, such that Δ𝑉(𝐸,𝐼)<Δ𝑉(𝐸^∗,𝐼^∗)/2<0 for any (𝐸,𝐼) satisfying ‖(𝐸,𝐼)−(𝐸^∗,𝐼^∗)‖<𝜂. Consequently, there exists a large integer 𝐾 such that ‖(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))−(𝐸^∗,𝐼^∗)‖<𝜂 for 𝑘>𝐾, |𝑉(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))−𝑉(𝐸^∗,𝐼^∗)|<|Δ𝑉(𝐸^∗,𝐼^∗)|/3, and Δ𝑉(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))<Δ𝑉(𝐸^∗,𝐼^∗)/2. Then we have 𝑉𝐸𝑡_𝑘+1,𝐼𝑡_𝑘+1=𝑉𝐸𝑡_𝑘,𝐼𝑡_𝑘+^{𝑡𝑘+1−1}_{𝑗=𝑡𝑘}𝐸𝑗𝑗𝐸𝑡Δ𝑉,𝐼<𝑉_𝑘,𝐼𝑡_𝑘+Δ𝑉𝐸𝑡_𝑘,𝐼𝑡_𝑘≤𝑉^∗+Δ𝑉𝐸^∗,𝐼^∗3+Δ𝑉𝐸^∗,𝐼^∗2.=𝑉^∗+Δ𝑉𝐸^∗,𝐼^∗6. (4.18) Inequality 𝑉(𝐸(𝑡_𝑘+1),𝐼(𝑡_𝑘+1))≤𝑉^∗+Δ𝑉(𝐸^∗,𝐼^∗)/6 implies that lim_𝑘→∞𝑉(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))≤𝑉^∗+Δ𝑉(𝐸^∗,𝐼^∗)/6, which is a contradiction. Therefore, lim_𝑘→∞𝑉(𝐸(𝑡_𝑘),𝐼(𝑡_𝑘))≤𝑉₀, which in turn leads to limsup_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))≤1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑁^∗11−𝑅₀ (4.19) by the definition of 𝑉₀. The proof of the theorem is completed.

The requirement for 𝛽₁ in Theorem 4.3 is technically strict. The interval for 𝛽₁ should be much longer than the unit interval given in Theorem 4.3. Theorem 4.3 holds as long as 1−𝛽₁𝐼(𝑡)/𝑁^∗≥0; that is, the range of 𝛽₁ must guarantee that the solutions of model (4.1) are all positive.

If 0<𝛽₁≤1 and 𝑘₁=0, lim_𝑘→∞𝑁(𝑡)=𝑁^∗ and Theorem 4.3 implies that the solutions of model (3.1) are ultimately boundaries as 1<𝑅₀<1+(1/(1−𝑝(1−𝛾))).

Theorem 4.3 provides an estimation on the ultimate boundary based on the combination of the latent and the infectious when 𝑅₀>1. We notice that most epidemic researches have focused on the persistence of the diseases by studying the existence and stability of the endemic equilibrium when 𝑅₀>1. Although we cannot obtain the exact quantity for the latent and the infectious, the ultimate boundary in Theorem 4.3 provides a novel estimation.

Figure 1(b) demonstrates the role of Theorem 4.3. The parameter values are the same as the ones we used for Figure 1(a) but with 𝛽₁=0.3. For that set of parameter values, 𝑅₀=1.35704>1, the endemic equilibrium 𝑃^∗₁(𝑆^∗₁,𝐸^∗₁,𝐼^∗₁)=(1228.16000,357.57414, 80.93252), and 𝑁^∗=1666.66667. The ultimate boundary domain in 𝐸𝐼-plane is given in Figure 1(b). The largest triangle in Figure 1(b) is the feasible domain 𝐸+𝐼≤𝑁^∗=1666.66667. The domain under the line 𝑉(𝐸,𝐼)=𝑝𝛼𝑘₂𝐸+(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃)𝐼=𝑉₀ is the domain of the ultimate boundary. Three solutions of model (3.1) with initial values: 𝐸(0)=35,𝐼(0)=200, 𝐸(0)=200,𝐼(0)=1000, and 𝐸(0)=700,𝐼(0)=50, are also displayed in Figure 1(b). Two solutions under the ultimate boundary 𝑉(𝐸,𝐼)=𝑉₀ are kept below the boundary and approach the endemic equilibrium. One solution with the initial value starting outside the ultimate boundary enters the ultimate boundary domain eventually.

4.2. With Exogenous Reinfections and Linear Treatment Rate

4.2.1. Backward Bifurcation

When 𝑘₁>0, the last equation of (3.9) leads to 𝐸^∗₂=(1−𝑝(1−𝛾))𝐼^∗₂/((𝑝𝑘₁𝛽₂𝐼^∗₂/𝑁^∗)+𝑝𝑘₂𝛼). Substituting the expression 𝐸^∗₂ into the first equation of (3.9) yields (use 𝐼 instead of 𝐼^∗₂ for convenience)𝑎₂𝐼²+𝑎₁𝐼+𝑎₀=0, (4.20) where 𝑎₂=𝑝²𝛽₁𝛽₂𝑘₁(𝑁^∗)²𝑎>0,₁=𝑝²𝛽₁𝑘₂𝛼𝑁^∗+𝑝𝛽₁+𝛽₂𝑘₁1−𝑝1−𝛾𝑁^∗−𝑝²𝛽₁𝛽₂𝑘₁𝑁^∗,≜−𝑅₁−1−1−𝑅₀1+1−𝑝1−𝛾−𝑝𝑘₁𝛽₂𝑝𝛼𝑘₂,𝑎₀=1−𝑝1−𝛾1−𝑝+𝑝𝑘₂𝛼+𝑝𝑘₃𝑚−𝑝²𝛽₁𝑘₂𝛼≜−𝑐₀𝑅₀−1, (4.21) where𝑅₁=𝑝𝑘₁𝛽₂𝑝𝛼𝑘₂×+1−𝑝1−𝛾1−𝑝+𝑝𝑚𝑘₃1−𝑝+𝛼𝑘₂+𝑝𝑚𝑘₃,𝑐₀=1−𝑝1−𝛾1−𝑝+𝑝𝑘₂𝛼+𝑝𝑘₃𝑚>0. (4.22) The type of bifurcation at 𝑅₀=1 is totally determined by 𝑅₁. Near 𝑅₀=1, 𝑎₁<0 holds if and only if 𝑅₁>1. Hence, at 𝑅₀=1, if 𝑅₁>1, bifurcation is forward, while if 𝑅₁<1, bifurcation is backward. 𝑅₁ gives a measure for exogenous reinfections. The first term in 𝑅₁, 1/(𝑝𝛼𝑘₂+1−𝑝(1−𝛾)), is the average re-infection period, which is a little bit shorter than the regular average infection period 1/(1−𝑝(1−𝛾)). The second term in 𝑅₁ gives the proportion of the latent population. These analysis result in the following theorem that characterizes the feature of bifurcation at 𝑅₀=1.

Theorem 4.4. Consider model (3.9). If 𝑅₁<1, then the bifurcation at 𝑅₀=1 is of forward type. If 𝑅₁>1, the bifurcation at 𝑅₀=1 is of backward type.

Looking at the expression of 𝑅₁, one can see that 𝑅₁>1 necessarily requires 𝑘₁𝛽₂≠0. Therefore, the exogenous reinfections are the driver behind the occurrence of the backward bifurcation. This is consistent with the continuous epidemic models for tuberculosis [30, 31].

The backward bifurcation of model (3.9) is shown in Figure 2 with following parameter values 𝑝=0.994, Λ=10, 𝛼=0.03, 𝛾=0.1, 𝑘₁=0.5, 𝑘₂=0.3, 𝑘₃=0.2, and 𝑚=0.1. Figure 2(a) gives three domains in 𝛽₁𝛽₂-plane with the unique disease-free equilibrium (DFE), two positive endemicequilibria(𝐸^∗_2±,𝐼^∗_2±), and one endemic equilibrium (𝐸^∗₂₊,𝐼^∗₂₊). Taking the above parameter values and fixing 𝛽₂=0.5, we get the backward bifurcation curve of model (3.9) (see Figure 2(b)).

(a)

(b)

Figure 2 

The backward bifurcation of model (3.9): (a) The domain with different number of endemic equilibria in 𝛽1𝛽2-plane, (b) The bifurcation curve of model (3.9).

We use the numerical simulation to discuss the stability of the endemic equilibria of model (3.9). The linearized matrix of (3.9) at the positive equilibrium points 𝑃^∗_2±(𝐸^∗_2±,𝐼^∗_2±) is𝐽_2±=𝑝−𝑝𝑘₂𝛼−𝑝𝑘₃𝑝𝛽𝑚−¹+𝛽₂𝑘₁𝐼^∗_2±𝑁^∗𝑝𝛽₁−𝑝𝛽₁+𝛽₂𝑘₁𝐸^∗_2±𝑁^∗−2𝑝𝛽₁𝐼^∗_2±𝑁^∗𝑝𝛽₂𝑘₁𝐼^∗_2±𝑁^∗+𝑝𝛼𝑘₂𝑝𝛽₂𝑘₁𝐸^∗_2±𝑁^∗+𝑝1−𝛾. (4.23) We fix parameter values 𝑝=0.994, Λ=10, 𝛼=0.03, 𝛾=0.1, 𝑘₁=0.5, 𝑘₂=0.3, 𝑘₃=0.2, 𝑚=0.1,and𝛽₂=0.5 and change the value of 𝛽₁. When 𝛽₁∈[0.31831317735,0.41278978), we have 𝑅₀∈[0.7711266,1), and there are two positive equilibrium points, which are (𝐸^∗₂₊,𝐼^∗₂₊) and (𝐸^∗₂₋,𝐼^∗₂₋), respectively. Fixing 𝛽₁=0.35, we have (𝐸^∗₂₊,𝐼^∗₂₊)=(601.51345,342.43485) and (𝐸^∗₂₋,𝐼^∗₂₋)=(329.51923,52.38933). Accordingly, we have 𝜌(𝐽₂₊)=0.979777<1 and 𝜌(𝐽₂₋)=1.00415>1. Therefore, (𝐸^∗₂₊,𝐼^∗₂₊) is local stable as 𝑅₀<1, while (𝐸^∗₂₋,𝐼^∗₂₋) is unstable as 𝑅₀<1. When 𝛽₁>0.41278978, we have 𝑅₀>1, and there is only one positive equilibrium point (𝐸^∗₂₊,𝐼^∗₂₊). Fixing 𝛽₁=1, we have (𝐸^∗₂₊,𝐼^∗₂₊)=(657.40092,796.73372), and 𝜌(𝐽₂₊)=0.88069<1. 𝜌(𝐽₂₊) will change when 𝛽₁ becomes bigger. The numerical simulation shows that 𝜌(𝐽₂₊)=1 and 𝑅₀=9.91134 as 𝛽₁=4.0913. Therefore, we obtain that (𝐸^∗₂₊,𝐼^∗₂₊) is locally stable as 1<𝑅₀<9.91134, and, when 𝑅₀≥9.91134, (𝐸^∗₂₊,𝐼^∗₂₊) is unstable.

After spelling out the stability of the endemic equilibria, we use the numerical method to investigate the attraction domain of the positive equilibrium points in 𝐸𝐼-plane when 𝑅₀<1. Here, we use the same parameter values as for the analysis of the endemic equilibrium bifurcation when 𝑅₀<1. There are three equilibrium points, the disease-free equilibrium (0,0), the endemic equilibria (𝐸^∗₂₊,𝐼^∗₂₊), and (𝐸^∗₂₋,𝐼^∗₂₋), respectively. In Figure 3, the solid red line is the separatrix, which divides the whole domain into two parts, the attraction domain of the larger endemic equilibrium (𝐸^∗₂₊,𝐼^∗₂₊), and the attraction domain of the disease-free equilibrium (0,0). Figure 3 shows that (𝐸^∗₂₋,𝐼^∗₂₋) is saddle, and the ultimate limit of solutions of model (3.1) depend on the initial values. If the initial values below the separatrix, the solutions will tend to the disease-free equilibrium (0,0), while the solutions with initial values above the separatrix will ultimately tend to (𝐸^∗₂₊,𝐼^∗₂₊). It appears that the attraction domain of (𝐸^∗₂₊,𝐼^∗₂₊) is much bigger than that of (0,0). This fact implies that we should be careful in assessing a disease control program, since the attraction domain of (0,0) is small, the disease may not be eliminated even though 𝑅₀<1.

Figure 3 

The attraction separatrix for those two equilibria of model (3.1).

4.2.2. Model Behavior for Large 𝑅₀

This subsection is mainly concerned with the behavior of the model when 𝑅₀ is very large. Numerical approach is used to further explore the model. We fix parameter values 𝑝=0.994, Λ=10, 𝛼=0.03, 𝛾=0.1, 𝑘₁=0.5, 𝑘₂=0.3, 𝑘₃=0.2, 𝑚=0.1, 𝛽₂=0.5, and let 𝛽₁ change. Figure 4(a) shows the existence of the stable endemic equilibrium of when 𝑅₀<1. In fact, model (3.9) has two endemic equilibrium points when 0.77113<𝑅₀<1, the small endemic equilibrium is unstable, and the large one is locally stable. If 1<𝑅₀<8.6485, model (3.9) has only one endemic equilibrium, which is locally stable. When 8.6485≤𝑅₀<9.6417 the endemic equilibrium of model (3.9) is unstable, and there exists a stable periodic solution of model (3.9) with period 2. Furthermore, the periodic solution with period 2 become unstable when 9.6417≤𝑅₀<9.9566, and a stable periodic solution with period 4 appears. When 𝑅₀ becomes larger, the period-4 solution loses its stability and a stable period-8 solution appears (Figure 4(b)). The period-doubling may undergo to chaos as 𝑅₀ increases.

(a)

(b)

Figure 4 

Period doubling bifurcation of model (3.9). Parameter values are 𝑝=0.994, Λ=10, 𝛼=0.03, 𝛾=0.1, 𝑘1=0.5, 𝑘2=0.3, 𝑘3=0.2, 𝑚=0.1, and𝛽2=0.5.

Figure 4 only displays the stable endemic equilibrium or stable periodic solution of model (3.9). The values of the equilibrium or the periodic solution of the infectious individuals are positive for all 𝑅₀∈[0,10]. However, the number of the latent class may be negative when 𝑅₀≥8.6485. 𝑅₀=8.6485 corresponds to the critical value at which the endemic equilibrium loses its stability, and periodic solution with period 2 appears. A reason for model (3.9) undergoing period-doubling solutions may be the appearance of the negative latent solutions when 𝑅₀ becomes larger. The numerical simulation seems to hint that the endemic equilibrium may be stable if 𝛽₁≤1 and 𝛽₂≤1, which gives the positive solution if the initial values are positive. Castillo-Chavez and Yakubu discussed similar problem for a discrete SIS model [19].

4.2.3. Ultimate Boundary

The following theorem gives an ultimate boundary to the model (3.9).

Theorem 4.5. If 1<𝑅₀≤1+((𝑝𝛽₁𝛼𝑘₂−(1−𝑝(1−𝛾))(1−𝑝+𝑝𝑚𝑘₃)𝑘₁𝛽₂)/(1−𝑝(1−𝛾))(𝑝𝛽₁𝛼𝑘₂+(1−𝑝+𝑝𝑚𝑘₃)𝑘₁𝛽₂)), 𝑝𝛽₁𝛼𝑘₂>(1−𝑝(1−𝛾))(1−𝑝+𝑝𝑚𝑘₃)𝑘₁𝛽₂, and 0<𝛽₁,𝛽₂≤1 hold, then the solutions of model (3.9) are ultimately bounded by limsup_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))=limsup_𝑡→∞𝑝𝛼𝑘₂𝐸(𝑡)+1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝐼(𝑡)≤𝑉₀,𝑤ℎ𝑒𝑟𝑒𝑉⁰≜1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃𝑁^∗1+1−𝑝+𝑝𝑚𝑘₃𝑘₁𝛽₂𝑝𝛽₁𝛼𝑘₂−1𝑅₀. (4.24)

The inequality limsup_𝑡→∞𝑉(𝐸(𝑡),𝐼(𝑡))≤𝑉₀ implies thatlimsup_𝑡→∞𝐸(𝑡)≤1+1−𝑝+𝑝𝑚𝑘₃𝑝𝛼𝑘₂𝑁^∗1+1−𝑝+𝑝𝑚𝑘₃𝑘₁𝛽₂𝑝𝛽₁𝛼𝑘₂−1𝑅₀,limsup_𝑡→∞𝐼(𝑡)≤𝑁^∗1+1−𝑝+𝑝𝑚𝑘₃𝑘₁𝛽₂𝑝𝛽₁𝛼𝑘₂−1𝑅₀. (4.25) Introducing 𝑉(𝐸(𝑡),𝐼(𝑡))=𝑝𝛼𝑘₂𝐸(𝑡)+(1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃)𝐼(𝑡) and then examining the difference of 𝑉 yields that(𝐸()((𝐸(𝑡)(𝑡=Δ𝑉=𝑉𝑡+1,𝐼𝑡+1))−𝑉,𝐼))1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃×𝑅₀𝑅−1−⁰𝑁^∗1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃(𝐸(𝑡)+𝐼(𝑡))+1−𝑝+𝑝𝑚𝑘₃𝑘₁𝛽₂𝑅₀𝐸(𝑡)𝑝𝛽₁𝛼𝑘₂𝑁^∗𝐼(𝑡)≤1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃×𝑅₀𝑅−1−⁰𝑁^∗𝑉(𝐸(𝑡),𝐼(𝑡))1−𝑝+𝑝𝛼𝑘₂+𝑝𝑚𝑘₃+1−𝑝+𝑝𝑚𝑘₃𝑘₁𝛽₂𝑅₀𝑝𝛽₁𝛼𝑘₂𝐼(𝑡),=1−𝑝1−𝛾𝑅₀𝑁^∗𝑉⁰−𝑉(𝐸(𝑡),𝐼(𝑡))𝐼(𝑡). (4.26) The rest of the proof is similar to that of Theorem 4.4.

4.3. With Exogenous Reinfections and Saturated Treatment

In this subsection, we analyze model (3.10), that is, the SEIS model with a saturated treatment rate.

4.3.1. Endemic Equilibria

Let 𝑃^∗₃(𝐸^∗₃,𝐼^∗₃) be the positive equilibrium of model (3.10). Then the second equation of model (3.10) leads to 𝐼^∗₃=𝑝𝛼𝑘₂𝐸^∗₃/(1−𝑝(1−𝛾)−(𝑝𝛽₂𝑘₁𝐸^∗₃/𝑁^∗)), where 𝐸^∗₃<𝑁^∗(1−𝑝(1−𝑟))/𝑝𝛽₂𝑘₁.

We already have a clear picture about the equilibrium in the region 𝐸(𝑡)≤𝐸⁰. Now we only seek for equilibria in the region 𝐸(𝑡)>𝐸⁰ with 𝑇(𝐸(𝑡))=𝑝(𝑘₃𝐸(𝑡)−𝑚₀). Similarly, 𝐸^∗₃ satisfies 𝜎₃𝐸³+𝜎₂𝐸²+𝜎₁𝐸+𝜎₀=0, (4.27) where 𝜎₃=𝑝²𝛽₂𝑘₁𝑝𝛽₁𝛼𝑘₂−𝛽₂𝑘₁1−𝑝(𝑁^∗)²,𝜎₂=2𝑝𝛽₂𝑘₁1−𝑝+𝑝𝛼𝑘₂𝑁^∗+𝑝1−𝑝1−𝛾²𝛼𝑘₂𝛽₂𝑘₁+𝛽₁𝑁^∗−𝑝1−𝑝1−𝛾₃𝛽₁𝛼𝑘₂𝛽₂𝑘₁+𝛼𝑘₂𝑁^∗−𝑝³𝑚₀𝛽²₂𝑘²₁(𝑁^∗)²,𝜎₁=1−𝑝1−𝛾2𝑝²𝑚₀𝛽₂𝑘₁𝑁^∗+𝑝²𝛽₁𝛼𝑘₂−1−𝑝1−𝛾1−𝑝+𝑝𝛼𝑘₂,𝜎₀=−𝑝𝑚₀1−𝑝1−𝛾². (4.28)

From the continuity argument, we can impose 𝐸⁰ to be a root to (4.27). Obviously, 𝑅₀ and 𝑅₁ are not enough to determine all possible configurations of the roots of (4.27). We can directly use these sigmas that have less epidemiological meaning. Notice that 𝜎₀<0 relieves our work a little bit. If 𝜎₃<0, then (4.27) has exact two positive roots. If 𝜎₃>0 and 𝜎₁<0, then (4.27) has exact one positive root. If 𝜎₃>0 and 𝜎₁>0, then (4.27) has exact one or three positive roots. For all situations, 𝐸⁰ is always the smallest positive root.

Next we will explore the dynamics of model (3.10) numerically by carefully selecting parameter values.

4.3.2. Double Backward Bifurcations

A double backward bifurcation is observed for parameter values 𝑝=0.994, Λ=10, 𝑘₁=0.2, 𝑘₂=0.4, 𝑘₃=0.4, 𝛼=0.03, 𝛾=0.0001, 𝑚=0.005, and 𝛽₂=0.375.

As can be seen in Figure 5, a backward bifurcation occurs at 𝑅₀=1 because 𝑅₁=1.6527>1. This backward bifurcation is bifurcated from the disease-free equilibrium. Meanwhile, at an endemic equilibrium, the second backward bifurcation happens, so that we can see that there is a window of five equilibria. Three of them are stable, and two of them are unstable. Indeed, we have found a tristability when 𝑅₀<1 and when saturated treatment rate is used. Figure 6 shows the tristable situation on the phase plane.

(a) Multiply Backward Bifurcation

(b) A part of (a)

Figure 5 

Double backward bifurcations of model (3.10). Dotted lines represent unstable equilibria, and solid lines represent stable ones.

(a)

(b)

Figure 6 

Phase portrait of model (3.10): a case of multiple equilibria.

If we choose 𝛽₂=0.26, then𝑅₁=1.1458>1. The bifurcation at 𝑅₀=1 remains backward. However, the second backward bifurcation shifts to the right, happening at a place of 𝑅₀>1. One can observe bistability not only for 𝑅₀<1 but also for 𝑅₀>1. Figure 7(a) illustrates this bistable dynamics. Further decreasing of 𝛽₂ can produce a new bifurcation, which will be studied in the next subsection.

(a) Multiply backward bifurcation

(b) Forward-backward bifurcation

Figure 7 

(a) Double backward bifurcations for model (3.10). The upper one occurs when 𝑅0>1. (b) Forward-backward bifurcation for model (3.10). The upper one occurs when 𝑅0<1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.

4.3.3. Forward-Backward Bifurcation

If we set 𝑝=0.994, Λ=8, 𝑘₁=0.3, 𝑘₂=0.3, 𝑘₃=0.4,𝛼=0.03,𝛾=0.001, 𝑚=0.01, 𝑚₀=0.2, and 𝛽₂=0.1, then the bifurcation at 𝑅₀=1 is forward because 𝑅₁=0.9835<1, as we can see from Figure 7(b) that a backward bifurcation occurs on the top. The backward bifurcation curve returns all the way back to the region of 𝑅₀<1, generating an endemic in the region. This would be the worst consequence of the saturated treatment strategy because it help to establish the disease when 𝑅₀<1. If a fixed proportion of the latent individuals are treated (linear treatment rate is used), then the backward bifurcation on the top cannot happen, consequently the disease dies out when 𝑅₀<1.

A less important forward-backward bifurcation is the case where the whole backward bifurcation curve exists only in the region of 𝑅₀>1, thus only making endemic equilibrium bigger. Figure 8 shows this kind of bifurcation, where 𝛽₂=0.18 in order to make 𝑅₁=0.7933<1 and the rest of parameters are the same as in the previous subsection.

(a) Forward-backward bifurcation

(b) A part of (a)

Figure 8 

Forward-backward bifurcation for model (3.10). The upper one occurs when 𝑅0>1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.

5. Conclusion and Discussion

In this paper, we analyzed a class of discrete SEIS models and their epidemiological implications. First, the models without exogenous reinfections were considered. We obtained the global stability of the disease-free equilibrium, existence and uniqueness of endemic equilibrium and its stability. We also found the sufficient conditions for the ultimate boundary of the solutions. Then we considered the models with exogenous reinfections. These models undergo a backward bifurcation. Furthermore, we showed that the models with exogenous re-infection and treatment exhibit period-double phenomenon when 𝑅₀ grows larger. The period-double bifurcation was obtained for the reduced model (3.9). One component in the full model may become negative when the reduced model has periodic orbits. Finally we considered the models with exogenous re-infection and saturated treatment. No matter whether there is the exogenous re-infection, the saturated treatment can lead to a backward bifurcation at an endemic, but exogenous re-infection determines the backward bifurcation at the disease-free equilibrium.

From our qualitative and quantitative analysis, we had already found that a combination of exogenous reinfections and treatment regimens is capable of generating complex dynamics. These include tristability (disease-free equilibrium and two endemic equilibria) when 𝑅₀<1, bistability for both 𝑅₀<1 and 𝑅₀>1, and period doubling for a reduced model.

We now know that 𝑘₁𝛽₂ controls the appearance of the backward bifurcation at 𝑅₀=1, creating an endemic equilibrium when 𝑅₀<1, an example of bistable scenario. The treatment regime does not participate in the backward bifurcation 𝑅₀=1. However, the treatment regime is able to create equilibria at a higher level, far away from the disease-free equilibrium, regardless of the values of 𝑅₀. The treatment strategy change actually means stopping the trend of treatment effort and treating fewer individuals. This reduction in treatment efforts leads to the emerging of the larger endemic equilibrium. The worst scenario is that if the bifurcation at 𝑅₀=1 is forward, meaning that 𝑅₀<1 makes the disease die out, and if the treatment practice is changed through reducing the treatment effort, then an endemic is born when 𝑅₀<1. This reduction of treatment effort helps to reestablish the disease at local level, making it impossible to eliminate the disease. A plausible conclusion is the strategy should be to treat as many latent individuals as possible. Our model supports that the strategy that a fixed percentage (proportion) of the latent individuals should be treated consistently to eliminate the disease. For diseases like tuberculosis, failure to take care of the existence of long-lived latent individuals may be a critical reason to the longevity of the disease. Exogenous reinfections intrinsically result in bistability through the occurrence of a backward bifurcation, which is out of our control. However, persistence of the disease resulting from a change of policy should and can be avoided by all means. Obviously, it would be a sad story if reduction of treatment is due to an insufficient health infrastructure or budget.

Acknowledgment

This study was supported by NSFC (10971163) and by IDRC of Canada (104519-010).