# Complex Dynamics of Discrete SEIS Models with Simple Demography

Authors: Hui Cao; Yicang Zhou; Baojun Song
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

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


---

## Body

## 1. Introduction

Pure demographic (or ecological) discrete models can have very complicated dynamics. It is well known that the logistic modelxn+1=rxn(1-xn/K), Ricker's model xn+1=xner-xn/K [1–3], and Hassell model xn+1=λxn(1+axn)-b [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, xn+1=Λ+μxn, 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 Section2. 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 toall 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. LetS(t) the number of susceptible individuals at time t, E(t) be the number of latent individuals, and I(t) the number of infectious individuals. N(t)=S(t)+E(t)+I(t) 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 arek1, k2, and k3 with k1+k2+k3=1. Incorporating the exogenous reinfections and treatments into the models in [20, 21], our model framework takes the following form:(2.1)S(t+1)=Λ+pS(t)G1(I(t)N(t))+pH(E(t))+pγI(t),E(t+1)=pS(t)(1-G1(I(t)N(t)))+pk1E(t)G2(I(t)N(t))+pk2(1-α)E(t)+p(k3E(t)-H(E(t))),I(t+1)=pk1E(t)(1-G2(I(t)N(t)))+pαk2E(t)+p(1-γ)I(t),
where Λ is the recruitment rate into the population, p 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-G1 is the probability that susceptible individuals become infected. 1-G2 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 G1 and G2 that are considered as functions of the prevalence I(t)/N(t). exp(-βjI(t)/N(t)) and 1-(βjI(t)/N(t)) have been used for the probabilities of not becoming infected [20, 21]. We choose Gj(I(t)/N(t))=1-(βjI(t)/N(t)),j=1,2, where βj characterizes the disease transmission probability.In model (2.1), H(E(t)) is the successful treatment rate of the latent individuals and k3E(t)-H(E(t)) is the failure rate of the treatment. We consider two different types of treatment strategies for the latent individuals:(2.2)H(E(t))=mk3E(t)orH(E(t))=mk3{E(t),E(t)≤E0,E0,E(t)>E0,
where m is the conditional probability that a latent individual is treated successfully given that the individual receives the treatment and E0 is the cutting point of treatment determined by policy makers. Both m and E0 should rely on the capability of finding latent individuals and budgetary issues.We summarize some important terms in model (2.1) here. From time t to time t+1, the number of the susceptibles who get infected and become latent individuals is pS(t)(1-G1(I(t)/N(t))); the number of the latent who get re-infected and become infectious is pk1E(t)(1-G2(I(t)/N(t))); the number of the latent who naturally progress into the infectious compartment is pk2αE(t); the number of the latent who get recovered by successful treatment and become susceptible is pH(E(t)); the number of the infectious who get recovered and become susceptible is pγI(t). The number of the latent individuals at time t that remain as latent ones at time t+1 is pk1E(t)G2(I(t)/N(t))+pk2(1-α)E(t)+p(k3E(t)-H(E(t))). Similarly, the number of infectious individuals who do not change their epidemiological status from t to t+1 is counted by p(1-γ)I(t).


## 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,H(E(t))=mk3E(t), and model (2.1) can be rewritten as follows:(3.1)S(t+1)=Λ+pS(t)(1-β1I(t)N(t))+pmk3E(t)+pγI(t),E(t+1)=pS(t)β1I(t)N(t)+pk1E(t)(1-β2I(t)N(t))+p(1-α)k2E(t)+p(1-m)k3E(t),I(t+1)=pk1E(t)β2I(t)N(t)+pαk2E(t)+p(1-γ)I(t).
On the other hand, in the case of limited resources, the treatment function of latent individuals is a piecewise function,
(3.2)H(E(t))={mk3E(t),E(t)≤E0,m0,E(t)>E0,
where m0=mk3E0. Model (2.1) takes the following form:(3.3)S(t+1)=Λ+pS(t)(1-β1I(t)N(t))+pk3E(t)-T(E(t))+pγI(t),E(t+1)=pS(t)β1I(t)N(t)+pk1E(t)(1-β2I(t)N(t))+p(1-α)k2E(t)+T(E(t)),I(t+1)=pk1E(t)β2I(t)N(t)+pαk2E(t)+p(1-γ)I(t),
where(3.4)T(E(t))={p(1-m)k3E(t),ifE(t)≤E0,p(k3E(t)-m0),ifE(t)>E0.The disease-free equilibrium,P0(Λ/(1-p),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(3.5)F=(pS(t)β1I(t)N(t)pk1E(t)β2I(t)N(t)0),V=(p(1-k1β2I(t)N(t)-αk2-mk3)E(t)pαk2E(t)+p(1-γ)I(t)Λ+pS(t)(1-β1I(t)N(t))+pmk3E(t)+pγI(t)).
At the disease-free equilibrium, the straight forward calculation yields F=(0pβ100) and V=(p-pαk2-pmk30pαk2p(1-γ)). Therefore, the basic reproductive number of model (2.1) is(3.6)R0=p2β1αk2(1-p(1-γ))(1-p+pαk2+pmk3).
Rearrange terms in the expression of R0 as(3.7)R0=(pβ11-p(1-γ))(pαk21-p+pαk2+pmk3).
Each term in R0 has clear epidemiological interpretation. 1/(1-p(1-γ)) is the average infection period. pβ1/(1-p(1-γ)) is average new cases generated by a typical infectious member in the entire infection period. pαk2/(1-p+pαk2+pmk3) is the proportion that latent individuals become infectious by “natural" progression.It is known that the disease-free equilibriumP0 of (2.1) is locally asymptotically stable if R0<1 and unstable if R0>1 (see [34, Theorem  2.1]). We summarize stability results in terms of R0<1 as in the following theorem.Theorem 3.1.
The disease-free equilibriumP0 of model (2.1) is locally asymptotically stable if R0<1 and unstable if R0>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 N(t)=S(t)+E(t)+I(t):(3.8)N(t+1)=Λ+pN(t).
One can see that N*=Λ/(1-p) is a global attractor for (3.8). Using limiting equations and the limt→∞N(t)=N* [35], we, respectively, reduce the three-dimensional systems (3.1) and (3.3) into two-dimensional ones(3.9)E(t+1)=p(N*-E(t)-I(t))β1I(t)N*+pk1E(t)(1-β2I(t)N*)+p(1-α)k2E(t)+p(1-m)k3E(t),I(t+1)=pk1E(t)β2I(t)N*+pαk2E(t)+p(1-γ)I(t),(3.10)E(t+1)=p(N*-E(t)-I(t))β1I(t)N*+pk1E(t)(1-β2I(t)N*)+p(1-α)k2E(t)+T(E(t)),I(t+1)=pk1E(t)β2I(t)N*+pαk2E(t)+p(1-γ)I(t).
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,k1=0, H(E)=mk3E, and model (3.9) is reduced to(4.1)E(t+1)=p(N*-E(t)-I(t))β1I(t)N*+p(1-α)k2E(t)+p(1-m)k3E(t),I(t+1)=pαk2E(t)+p(1-γ)I(t).
The disease-free equilibrium of (4.1) is P10(E10,I10)=(0,0) and the endemic equilibrium P1*(E1*,I1*) exists if R0>1, where E1*=((1-p(1-γ))N*/(1-p(1-γ)+pk2α))(1-(1/R0)) and I1*=pk2αE1*/(1-p(1-γ)).


#### 4.1.1. Global Stability ofP10

The global stability of the disease-free equilibriumP10(E10,I10)=(0,0) of (4.1) is given in Theorem 4.1.Theorem 4.1.
IfR0<1, then the disease-free equilibrium P10 of (4.1) is globally asymptotically stable.Proof.
DefineF:[0,N*]×[0,N*]→[0,N*]×[0,N*] by
(4.2)F(E,I)=p((N*-E-I)β1IN*+(1-α)k2E+(1-m)k3E,αk2E+(1-γ)I).
Obviously, F is the mapping derived by system (4.1), and (0,0) is a fixed point of F. Linear function V(E,I)=E+(pβ1/(1-p(1-γ)))I on [0,N*]×[0,N*] is continuous and positive definite with respect to (0,0). Therefore, V is a Lyapunov function on the domain of F. For any (E,I)∈[0,N*]×[0,N*],
(4.3)ΔV(E,I)=p(N*-E-I)β1IN*+p(1-α)k2E+p(1-m)k3E+pβ11-p(1-γ)(pαk2E+p(1-γ)I)-E-pβ11-p(1-γ)I≤pβ1I-(1-p(1-α)k2-p(1-m)k3)E+p2β1αk2E1-p(1-γ)+pβ1I=(1-p(1-α)k2-p(1-m)k3)(R0-1)E≤0.
Hence, if R0<1, then ΔV(E,I)<0 holds for (E,I)∈(0,N*]×(0,N*]. It follows from Theorem 4.22 in [36] that P10 is globally asymptotically stable.Theorem4.1 ensures that the disease-free equilibrium P0(Λ/(1-p),0,0) of model (3.1) is globally stable as R0<1 and k1=0.


#### 4.1.2. Stability ofP1*

The stability of the endemic equilibrium is given in Theorem4.2.Theorem 4.2.
If1<R0<1+σ, then the endemic equilibrium P1* of (4.1) is asymptotically stable, where
(4.4)σ=2p(1-αk2-mk3+(1-γ))(1-p(1-γ)+pαk2)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0.Proof.
The linearization matrix of (4.1) at the endemic equilibrium P1* is
(4.5)J1=(-pβ1N*I1*+p(1-α)k2+p(1-m)k3pβ1N*(N*-E1*-2I1*)pαk2p(1-γ)).
The corresponding characteristic equation is
(4.6)f1(λ)=λ2-b1λ+b2=0,
where
(4.7)b1=p-pαk2-pmk3+p(1-γ)-(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2,b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3)+(pαk2-p(1-γ))(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
We claim |b2|<1. When pαk2=p(1-γ), we have b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3), and |b2|<1 holds. When pαk2≠p(1-γ), because of
(4.8)2p(1-αk2-mk3+(1-γ))(1+p(1-γ)-pαk2)<min⁡{2-p+pαk2+pmk3-p(1-γ)pαk2-p(1-γ),p(1-αk2-mk3+(1-γ))p(1-γ)-pαk2},
we obtain that |b2|<1 holds if 1<R0<1+σ. Furthermore, we have
(4.9)f1(1)=1-b1+b2=(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1),(-1)2f1(-1)=1+b1+b2=2(p-pαk2-pmk3+p(1-γ))-(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
It is clear that f1(1)>0 holds if R0>1 and
(4.10)(-1)2f1(-1)=σ(1+σ-R0)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0
holds if R0<1+σ. Therefore, if 1<R0<1+σ, we have f1(1)>0, (-1)2f1(-1)>0, and |b2|<1. It follows from the Jury criteria that the endemic equilibrium P1*(E1*,I1*) of (4.1) is asymptotically stable.It follows from Theorem4.2 that when exogenous re-infection is not included (k1β2=0) the endemic equilibrium P1*(Λ/(1-p))-E1*-I1*,E1*,I1*) of model (3.1) is asymptotically stable if 1<R0<1+σ. Since the Jury criteria is the necessary and sufficient condition for the local stability of the equilibrium, the endemic equilibrium P1* of (3.1) becomes unstable if R0≥1+σ. We summarize Theorems 4.1 and 4.2 for model (3.1) without exogenous reinfection as follows.  (1)
The disease-free equilibriumP0(Λ/(1-p),0,0) is globally stable as R0<1.(2)
WhenR0>1, P0(Λ/(1-p),0,0) becomes unstable and endemic equilibrium P1* is born. The endemic P1* is locally asymptotically stable if 1<R0<1+σ.(3)
WhenR0≥1+σ,P1* becomes unstable.Theorems4.1 and 4.2 state that a forward bifurcation takes place at R0=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 β1 while R0 depends on β1 linearly. For example, given a set of the parameter values, p=0.994, Λ=10, α=0.03, γ=0.1, m=0.1, k1=0, k2=0.8, and k3=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<R0<1+σ, as shown in Figure 1(a) for β1=0.245. When β1=10.972015, we have R0=49.63168 and 1+σ=49.63168. The endemic equilibrium of model (4.1) loses the stability when β1 is greater than 10.972015, and there is a 2-period solution when β1∈[10.972015,11.109). For example, if we take β1=11, then R0=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 k1β2=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 S* 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-G1(I(t)/N(t))=β1I(t)/N(t) means that β1 should not be bigger, and we should chose β1 such that 0≤1-G1(I(t)/N(t))≤1. Otherwise, not only does 1-G1(I(t)/N(t)) lose its meaning but also the number of the susceptible individuals has negative values.The stability and the attraction of the endemic equilibrium of (3.1): (a) the endemic equilibrium and its stability and (b) the attraction area.
(a)(b)


#### 4.1.3. Ultimate Boundary

We now estimate an ultimate boundary to the model (4.1) when R0>1.Theorem 4.3.
If0<β1≤1 and 1<R0<1+(1/(1-p(1-γ))), then the solutions of model (4.1) are bounded by
(4.11)lim⁡⁡sup⁡t→∞(αk2E(t)+(1-p+pαk2+pmk3)I(t))≤(1-p+pαk2+pmk3)N*(1-1R0),lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1-1R0),limsup⁡t→∞I(t)≤N*(1-1R0).Proof.
Let us consider the value ofV(E(t),I(t)) for any (E(t),I(t))∈Ω, where
(4.12)Ω={(E(t),I(t))∣E(t)≥0,I(t)≥0,E(t)+I(t)≤N*}.
Define
(4.13)V(E,I)=pαk2E+(1-p+pαk2+pmk3)I,V0=(1-p+pαk2+pmk3)N*(1-1R0).
Then,
(4.14)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=q(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t)))I(t)≤q(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t),
where q=(1-p(1-γ)(1-p+pαk2+pmk3)). Therefore, ΔV≤0 holds if V≥V0.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)≤V0 for a t0≥0, then we claim that V(E(t),I(t))≤V0 for all t≥t0. In fact, let V(E(t0),I(t0))=V0-ɛ with some ɛ≥0. Then, the inequality V(E(t0),I(t0))=V0-ɛ≤V0 yields
(4.15)I(t0)≤V01-p+pαk2+pmk3=N*(1-1R0).
The inequality in (4.14) and the condition in Theorem 4.3 as well as the inequality for I(t0) imply
(4.16)ΔV=V(E(t0+1),I(t0+1))-V(E(t0),I(t0))≤(1-p(1-γ))R0N*(V0-V(E(t0),I(t0)))I(t0),≤(1-p(1-γ))R0(V0-V(E(t0),I(t0)))(1-1R0)≤(1-p(1-γ))(R0-1)ɛ≤ɛ.
The inequality in (4.16) leads to
(4.17)V(E(t0+1),I(t0+1))≤V(E(t0),I(t0))+ɛ=V0-ɛ+ɛ=V0.
Then, our claim is followed from the mathematical induction.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)>V0 for a t0≥0, then there are two cases that we have to deal with. Case 1.
There exists an integert1>t0 such that V(E(t1),I(t1)≤V0.Case 2.
The inequalityV(E(t),I(t))>V0 holds for all t≥t0.For Case 1 we can prove that V(E(t),I(t))≤V0 for all t≥t1 by the similar argument adopted for the case V(E(t0),I(t0))≤V0. The boundary of V(E(t),I(t))≤V0 for t>t1 gives rise to limsupt→∞V(E(t),I(t))≤V0. For Case 2, the fact that ΔV(E(t),I(t))≤0 for V(E(t),I(t))>V0 leads to that the V(E(t),I(t)) is a monotonic decreasing sequence. Hence, the limit limt→∞V(E(t),I(t)) exits. We claim that limt→∞V(E(t),I(t))≤V0. Otherwise, there exists V*>V0, such that limt→∞V(E(t),I(t))=V*>V0. It follows that there exists a subsequence tk, such that limk→∞E(tk)=E*, limk→∞I(tk)=I*, and limk→∞V(E(tk),I(tk))=V(E*,I*)=V*>V0. The fact that ΔV(E*,I*)<0 and the continuity of ΔV(E,I) imply that there exists an η>0, such that ΔV(E,I)<ΔV(E*,I*)/2<0 for any (E,I) satisfying ∥(E,I)-(E*,I*)∥<η. Consequently, there exists a large integer K such that ∥(E(tk),I(tk))-(E*,I*)∥<η for k>K, |V(E(tk),I(tk))-V(E*,I*)|<|ΔV(E*,I*)|/3, and ΔV(E(tk),I(tk))<ΔV(E*,I*)/2. Then we have
(4.18)V(E(tk+1),I(tk+1))=V(E(tk),I(tk))+∑j=tktk+1-1ΔV(E(j),I(j))<V(E(tk),I(tk))+ΔV(E(tk),I(tk))≤V*+|ΔV(E*,I*)|3+ΔV(E*,I*)2.=V*+ΔV(E*,I*)6.
Inequality V(E(tk+1),I(tk+1))≤V*+ΔV(E*,I*)/6 implies that limk→∞V(E(tk),I(tk))≤V*+ΔV(E*,I*)/6, which is a contradiction. Therefore, limk→∞V(E(tk),I(tk))≤V0, which in turn leads to
(4.19)lim⁡⁡sup⁡t→∞V(E(t),I(t))≤(1-p+pαk2+pmk3)N*(1-1R0)
by the definition of V0. The proof of the theorem is completed.The requirement forβ1 in Theorem 4.3 is technically strict. The interval for β1 should be much longer than the unit interval given in Theorem 4.3. Theorem 4.3 holds as long as 1-β1I(t)/N*≥0; that is, the range of β1 must guarantee that the solutions of model (4.1) are all positive.If0<β1≤1 and k1=0, limk→∞N(t)=N* and Theorem 4.3 implies that the solutions of model (3.1) are ultimately boundaries as 1<R0<1+(1/(1-p(1-γ))).Theorem4.3 provides an estimation on the ultimate boundary based on the combination of the latent and the infectious when R0>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 R0>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.Figure1(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 β1=0.3. For that set of parameter values, R0=1.35704>1, the endemic equilibrium P1*(S1*,E1*,I1*)=(1228.16000,357.57414, 80.93252), and N*=1666.66667. The ultimate boundary domain in EI-plane is given in Figure 1(b). The largest triangle in Figure 1(b) is the feasible domain E+I≤N*=1666.66667. The domain under the line V(E,I)=pαk2E+(1-p+pαk2+pmk3)I=V0 is the domain of the ultimate boundary. Three solutions of model (3.1) with initial values: E(0)=35,I(0)=200, E(0)=200,I(0)=1000, and E(0)=700,I(0)=50, are also displayed in Figure 1(b). Two solutions under the ultimate boundary V(E,I)=V0 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

Whenk1>0, the last equation of (3.9) leads to E2*=(1-p(1-γ))I2*/((pk1β2I2*/N*)+pk2α). Substituting the expression E2* into the first equation of (3.9) yields (use I instead of I2* for convenience)(4.20)a2I2+a1I+a0=0,
where (4.21)a2=p2β1β2k1(N*)2>0,a1=p2β1k2αN*+p(β1+β2k1)(1-p(1-γ))N*-p2β1β2k1N*,≜-(R1-1)-(1-R0)(1+1-p(1-γ)-pk1β2pαk2),a0=(1-p(1-γ))(1-p+pk2α+pk3m)-p2β1k2α≜-c0(R0-1),
where(4.22)R1=pk1β2pαk2+1-p(1-γ)×1-p+pmk31-p+αk2+pmk3,c0=(1-p(1-γ))(1-p+pk2α+pk3m)>0.
The type of bifurcation at R0=1 is totally determined by R1. Near R0=1, a1<0 holds if and only if R1>1. Hence, at R0=1, if R1>1, bifurcation is forward, while if R1<1, bifurcation is backward. R1 gives a measure for exogenous reinfections. The first term in R1, 1/(pαk2+1-p(1-γ)), is the average re-infection period, which is a little bit shorter than the regular average infection period 1/(1-p(1-γ)). The second term in R1 gives the proportion of the latent population. These analysis result in the following theorem that characterizes the feature of bifurcation at R0=1.Theorem 4.4.
Consider model (3.9). If R1<1, then the bifurcation at R0=1 is of forward type. If R1>1, the bifurcation at R0=1 is of backward type.Looking at the expression ofR1, one can see that R1>1 necessarily requires k1β2≠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 p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, and  m=0.1. Figure 2(a) gives three domains in β1β2-plane with the unique disease-free equilibrium (DFE), two positive endemicequilibria(E2±*,I2±*), and one endemic equilibrium (E2+*,I2+*). Taking the above parameter values and fixing β2=0.5, we get the backward bifurcation curve of model (3.9) (see Figure 2(b)).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).
(a)(b)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 P2±*(E2±*,I2±*) is(4.23)J2±=(p-pk2α-pk3m-p(β1+β2k1)I2±*N*pβ1-p(β1+β2k1)E2±*N*-2pβ1I2±*N*pβ2k1I2±*N*+pαk2pβ2k1E2±*N*+p(1-γ)).
We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1,andβ2=0.5 and change the value of β1. When β1∈[0.31831317735,0.41278978), we have R0∈[0.7711266,1), and there are two positive equilibrium points, which are (E2+*,I2+*) and (E2-*,I2-*), respectively. Fixing β1=0.35, we have (E2+*,I2+*)=(601.51345,342.43485) and (E2-*,I2-*)=(329.51923,52.38933). Accordingly, we have ρ(J2+)=0.979777<1 and ρ(J2-)=1.00415>1. Therefore, (E2+*,I2+*) is local stable as R0<1, while (E2-*,I2-*) is unstable as R0<1. When β1>0.41278978, we have R0>1, and there is only one positive equilibrium point (E2+*,I2+*). Fixing β1=1, we have (E2+*,I2+*)=(657.40092,796.73372), and ρ(J2+)=0.88069<1. ρ(J2+) will change when β1 becomes bigger. The numerical simulation shows that ρ(J2+)=1 and R0=9.91134 as β1=4.0913. Therefore, we obtain that (E2+*,I2+*) is locally stable as 1<R0<9.91134, and, when R0≥9.91134, (E2+*,I2+*) 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 inEI-plane when R0<1. Here, we use the same parameter values as for the analysis of the endemic equilibrium bifurcation when R0<1. There are three equilibrium points, the disease-free equilibrium (0,0), the endemic equilibria (E2+*,I2+*), and (E2-*,I2-*), 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 (E2+*,I2+*), and the attraction domain of the disease-free equilibrium (0,0). Figure 3 shows that (E2-*,I2-*) 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 (E2+*,I2+*). It appears that the attraction domain of (E2+*,I2+*) 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 R0<1.Figure 3
The attraction separatrix for those two equilibria of model (3.1).


#### 4.2.2. Model Behavior for LargeR0

This subsection is mainly concerned with the behavior of the model whenR0 is very large. Numerical approach is used to further explore the model. We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, β2=0.5, and let β1 change. Figure 4(a) shows the existence of the stable endemic equilibrium of when R0<1. In fact, model (3.9) has two endemic equilibrium points when 0.77113<R0<1, the small endemic equilibrium is unstable, and the large one is locally stable. If 1<R0<8.6485, model (3.9) has only one endemic equilibrium, which is locally stable. When 8.6485≤R0<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≤R0<9.9566, and a stable periodic solution with period 4 appears. When R0 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 R0 increases.Period doubling bifurcation of model (3.9). Parameter values are p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, andβ2=0.5.
(a)(b)Figure4 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 R0∈[0,10]. However, the number of the latent class may be negative when R0≥8.6485. R0=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 R0 becomes larger. The numerical simulation seems to hint that the endemic equilibrium may be stable if β1≤1 and β2≤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.
If1<R0≤1+((pβ1αk2-(1-p(1-γ))(1-p+pmk3)k1β2)/(1-p(1-γ))(pβ1αk2+(1-p+pmk3)k1β2)), pβ1αk2>(1-p(1-γ))(1-p+pmk3)k1β2, and 0<β1,β2≤1 hold, then the solutions of model (3.9) are ultimately bounded by
(4.24)lim⁡⁡sup⁡t→∞V(E(t),I(t))=lim⁡⁡sup⁡t→∞pαk2E(t)+(1-p+pαk2+pmk3)I(t)≤V0,whereV0≜(1-p+pαk2+pmk3)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).The inequalitylimsupt→∞V(E(t),I(t))≤V0 implies that(4.25)lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0),lim⁡⁡sup⁡t→∞I(t)≤N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).
Introducing V(E(t),I(t))=pαk2E(t)+(1-p+pαk2+pmk3)I(t) and then examining the difference of V yields that(4.26)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t))+(1-p+pmk3)k1β2R0E(t)pβ1αk2N*)I(t)≤(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3+(1-p+pmk3)k1β2R0pβ1αk2)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t).
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

LetP3*(E3*,I3*) be the positive equilibrium of model (3.10). Then the second equation of model (3.10) leads to I3*=pαk2E3*/(1-p(1-γ)-(pβ2k1E3*/N*)), where E3*<N*(1-p(1-r))/pβ2k1.We already have a clear picture about the equilibrium in the regionE(t)≤E0. Now we only seek for equilibria in the region E(t)>E0 with T(E(t))=p(k3E(t)-m0). Similarly, E3* satisfies (4.27)σ3E3+σ2E2+σ1E+σ0=0,
where
(4.28)σ3=p2β2k1(pβ1αk2-β2k1(1-p))(N*)2,σ2=2pβ2k1(1-p+pαk2)N*(1-p(1-γ))+p2αk2(β2k1+β1)N*(1-p(1-γ))-p3β1αk2(β2k1+αk2)N*-p3m0β22k12(N*)2,σ1=(1-p(1-γ))(2p2m0β2k1N*+p2β1αk2-(1-p(1-γ))(1-p+pαk2)),σ0=-pm0(1-p(1-γ))2.From the continuity argument, we can imposeE0 to be a root to (4.27). Obviously, R0 and R1 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<0 relieves our work a little bit. If σ3<0, then (4.27) has exact two positive roots. If σ3>0 and σ1<0, then (4.27) has exact one positive root. If σ3>0 and σ1>0, then (4.27) has exact one or three positive roots. For all situations, E0 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 valuesp=0.994, Λ=10, k1=0.2, k2=0.4, k3=0.4, α=0.03, γ=0.0001, m=0.005, and β2=0.375.As can be seen in Figure5, a backward bifurcation occurs at R0=1 because R1=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 R0<1 and when saturated treatment rate is used. Figure 6 shows the tristable situation on the phase plane.Double backward bifurcations of model (3.10). Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply Backward Bifurcation(b)
A part of (a)Phase portrait of model (3.10): a case of multiple equilibria.
(a)(b)If we chooseβ2=0.26, thenR1=1.1458>1. The bifurcation at R0=1 remains backward. However, the second backward bifurcation shifts to the right, happening at a place of R0>1. One can observe bistability not only for R0<1 but also for R0>1. Figure 7(a) illustrates this bistable dynamics. Further decreasing of β2 can produce a new bifurcation, which will be studied in the  next subsection.(a) Double backward bifurcations for model (3.10). The upper one occurs when R0>1. (b) Forward-backward bifurcation for model (3.10). The upper one occurs when R0<1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply backward bifurcation(b)
Forward-backward bifurcation


#### 4.3.3. Forward-Backward Bifurcation

If we setp=0.994, Λ=8, k1=0.3, k2=0.3, k3=0.4,α=0.03,γ=0.001, m=0.01, m0=0.2, and β2=0.1, then the bifurcation at R0=1 is forward because R1=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 R0<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 R0<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 R0<1.A less important forward-backward bifurcation is the case where the whole backward bifurcation curve exists only in the region ofR0>1, thus only making endemic equilibrium bigger. Figure 8 shows this kind of bifurcation, where β2=0.18 in order to make R1=0.7933<1 and the rest of parameters are the same as in the previous subsection.Forward-backward bifurcation for model (3.10). The upper one occurs when R0>1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Forward-backward bifurcation(b)
A part of (a)


## 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,k1=0, H(E)=mk3E, and model (3.9) is reduced to(4.1)E(t+1)=p(N*-E(t)-I(t))β1I(t)N*+p(1-α)k2E(t)+p(1-m)k3E(t),I(t+1)=pαk2E(t)+p(1-γ)I(t).
The disease-free equilibrium of (4.1) is P10(E10,I10)=(0,0) and the endemic equilibrium P1*(E1*,I1*) exists if R0>1, where E1*=((1-p(1-γ))N*/(1-p(1-γ)+pk2α))(1-(1/R0)) and I1*=pk2αE1*/(1-p(1-γ)).


### 4.1.1. Global Stability ofP10

The global stability of the disease-free equilibriumP10(E10,I10)=(0,0) of (4.1) is given in Theorem 4.1.Theorem 4.1.
IfR0<1, then the disease-free equilibrium P10 of (4.1) is globally asymptotically stable.Proof.
DefineF:[0,N*]×[0,N*]→[0,N*]×[0,N*] by
(4.2)F(E,I)=p((N*-E-I)β1IN*+(1-α)k2E+(1-m)k3E,αk2E+(1-γ)I).
Obviously, F is the mapping derived by system (4.1), and (0,0) is a fixed point of F. Linear function V(E,I)=E+(pβ1/(1-p(1-γ)))I on [0,N*]×[0,N*] is continuous and positive definite with respect to (0,0). Therefore, V is a Lyapunov function on the domain of F. For any (E,I)∈[0,N*]×[0,N*],
(4.3)ΔV(E,I)=p(N*-E-I)β1IN*+p(1-α)k2E+p(1-m)k3E+pβ11-p(1-γ)(pαk2E+p(1-γ)I)-E-pβ11-p(1-γ)I≤pβ1I-(1-p(1-α)k2-p(1-m)k3)E+p2β1αk2E1-p(1-γ)+pβ1I=(1-p(1-α)k2-p(1-m)k3)(R0-1)E≤0.
Hence, if R0<1, then ΔV(E,I)<0 holds for (E,I)∈(0,N*]×(0,N*]. It follows from Theorem 4.22 in [36] that P10 is globally asymptotically stable.Theorem4.1 ensures that the disease-free equilibrium P0(Λ/(1-p),0,0) of model (3.1) is globally stable as R0<1 and k1=0.


### 4.1.2. Stability ofP1*

The stability of the endemic equilibrium is given in Theorem4.2.Theorem 4.2.
If1<R0<1+σ, then the endemic equilibrium P1* of (4.1) is asymptotically stable, where
(4.4)σ=2p(1-αk2-mk3+(1-γ))(1-p(1-γ)+pαk2)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0.Proof.
The linearization matrix of (4.1) at the endemic equilibrium P1* is
(4.5)J1=(-pβ1N*I1*+p(1-α)k2+p(1-m)k3pβ1N*(N*-E1*-2I1*)pαk2p(1-γ)).
The corresponding characteristic equation is
(4.6)f1(λ)=λ2-b1λ+b2=0,
where
(4.7)b1=p-pαk2-pmk3+p(1-γ)-(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2,b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3)+(pαk2-p(1-γ))(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
We claim |b2|<1. When pαk2=p(1-γ), we have b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3), and |b2|<1 holds. When pαk2≠p(1-γ), because of
(4.8)2p(1-αk2-mk3+(1-γ))(1+p(1-γ)-pαk2)<min⁡{2-p+pαk2+pmk3-p(1-γ)pαk2-p(1-γ),p(1-αk2-mk3+(1-γ))p(1-γ)-pαk2},
we obtain that |b2|<1 holds if 1<R0<1+σ. Furthermore, we have
(4.9)f1(1)=1-b1+b2=(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1),(-1)2f1(-1)=1+b1+b2=2(p-pαk2-pmk3+p(1-γ))-(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
It is clear that f1(1)>0 holds if R0>1 and
(4.10)(-1)2f1(-1)=σ(1+σ-R0)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0
holds if R0<1+σ. Therefore, if 1<R0<1+σ, we have f1(1)>0, (-1)2f1(-1)>0, and |b2|<1. It follows from the Jury criteria that the endemic equilibrium P1*(E1*,I1*) of (4.1) is asymptotically stable.It follows from Theorem4.2 that when exogenous re-infection is not included (k1β2=0) the endemic equilibrium P1*(Λ/(1-p))-E1*-I1*,E1*,I1*) of model (3.1) is asymptotically stable if 1<R0<1+σ. Since the Jury criteria is the necessary and sufficient condition for the local stability of the equilibrium, the endemic equilibrium P1* of (3.1) becomes unstable if R0≥1+σ. We summarize Theorems 4.1 and 4.2 for model (3.1) without exogenous reinfection as follows.  (1)
The disease-free equilibriumP0(Λ/(1-p),0,0) is globally stable as R0<1.(2)
WhenR0>1, P0(Λ/(1-p),0,0) becomes unstable and endemic equilibrium P1* is born. The endemic P1* is locally asymptotically stable if 1<R0<1+σ.(3)
WhenR0≥1+σ,P1* becomes unstable.Theorems4.1 and 4.2 state that a forward bifurcation takes place at R0=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 β1 while R0 depends on β1 linearly. For example, given a set of the parameter values, p=0.994, Λ=10, α=0.03, γ=0.1, m=0.1, k1=0, k2=0.8, and k3=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<R0<1+σ, as shown in Figure 1(a) for β1=0.245. When β1=10.972015, we have R0=49.63168 and 1+σ=49.63168. The endemic equilibrium of model (4.1) loses the stability when β1 is greater than 10.972015, and there is a 2-period solution when β1∈[10.972015,11.109). For example, if we take β1=11, then R0=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 k1β2=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 S* 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-G1(I(t)/N(t))=β1I(t)/N(t) means that β1 should not be bigger, and we should chose β1 such that 0≤1-G1(I(t)/N(t))≤1. Otherwise, not only does 1-G1(I(t)/N(t)) lose its meaning but also the number of the susceptible individuals has negative values.The stability and the attraction of the endemic equilibrium of (3.1): (a) the endemic equilibrium and its stability and (b) the attraction area.
(a)(b)


### 4.1.3. Ultimate Boundary

We now estimate an ultimate boundary to the model (4.1) when R0>1.Theorem 4.3.
If0<β1≤1 and 1<R0<1+(1/(1-p(1-γ))), then the solutions of model (4.1) are bounded by
(4.11)lim⁡⁡sup⁡t→∞(αk2E(t)+(1-p+pαk2+pmk3)I(t))≤(1-p+pαk2+pmk3)N*(1-1R0),lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1-1R0),limsup⁡t→∞I(t)≤N*(1-1R0).Proof.
Let us consider the value ofV(E(t),I(t)) for any (E(t),I(t))∈Ω, where
(4.12)Ω={(E(t),I(t))∣E(t)≥0,I(t)≥0,E(t)+I(t)≤N*}.
Define
(4.13)V(E,I)=pαk2E+(1-p+pαk2+pmk3)I,V0=(1-p+pαk2+pmk3)N*(1-1R0).
Then,
(4.14)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=q(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t)))I(t)≤q(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t),
where q=(1-p(1-γ)(1-p+pαk2+pmk3)). Therefore, ΔV≤0 holds if V≥V0.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)≤V0 for a t0≥0, then we claim that V(E(t),I(t))≤V0 for all t≥t0. In fact, let V(E(t0),I(t0))=V0-ɛ with some ɛ≥0. Then, the inequality V(E(t0),I(t0))=V0-ɛ≤V0 yields
(4.15)I(t0)≤V01-p+pαk2+pmk3=N*(1-1R0).
The inequality in (4.14) and the condition in Theorem 4.3 as well as the inequality for I(t0) imply
(4.16)ΔV=V(E(t0+1),I(t0+1))-V(E(t0),I(t0))≤(1-p(1-γ))R0N*(V0-V(E(t0),I(t0)))I(t0),≤(1-p(1-γ))R0(V0-V(E(t0),I(t0)))(1-1R0)≤(1-p(1-γ))(R0-1)ɛ≤ɛ.
The inequality in (4.16) leads to
(4.17)V(E(t0+1),I(t0+1))≤V(E(t0),I(t0))+ɛ=V0-ɛ+ɛ=V0.
Then, our claim is followed from the mathematical induction.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)>V0 for a t0≥0, then there are two cases that we have to deal with. Case 1.
There exists an integert1>t0 such that V(E(t1),I(t1)≤V0.Case 2.
The inequalityV(E(t),I(t))>V0 holds for all t≥t0.For Case 1 we can prove that V(E(t),I(t))≤V0 for all t≥t1 by the similar argument adopted for the case V(E(t0),I(t0))≤V0. The boundary of V(E(t),I(t))≤V0 for t>t1 gives rise to limsupt→∞V(E(t),I(t))≤V0. For Case 2, the fact that ΔV(E(t),I(t))≤0 for V(E(t),I(t))>V0 leads to that the V(E(t),I(t)) is a monotonic decreasing sequence. Hence, the limit limt→∞V(E(t),I(t)) exits. We claim that limt→∞V(E(t),I(t))≤V0. Otherwise, there exists V*>V0, such that limt→∞V(E(t),I(t))=V*>V0. It follows that there exists a subsequence tk, such that limk→∞E(tk)=E*, limk→∞I(tk)=I*, and limk→∞V(E(tk),I(tk))=V(E*,I*)=V*>V0. The fact that ΔV(E*,I*)<0 and the continuity of ΔV(E,I) imply that there exists an η>0, such that ΔV(E,I)<ΔV(E*,I*)/2<0 for any (E,I) satisfying ∥(E,I)-(E*,I*)∥<η. Consequently, there exists a large integer K such that ∥(E(tk),I(tk))-(E*,I*)∥<η for k>K, |V(E(tk),I(tk))-V(E*,I*)|<|ΔV(E*,I*)|/3, and ΔV(E(tk),I(tk))<ΔV(E*,I*)/2. Then we have
(4.18)V(E(tk+1),I(tk+1))=V(E(tk),I(tk))+∑j=tktk+1-1ΔV(E(j),I(j))<V(E(tk),I(tk))+ΔV(E(tk),I(tk))≤V*+|ΔV(E*,I*)|3+ΔV(E*,I*)2.=V*+ΔV(E*,I*)6.
Inequality V(E(tk+1),I(tk+1))≤V*+ΔV(E*,I*)/6 implies that limk→∞V(E(tk),I(tk))≤V*+ΔV(E*,I*)/6, which is a contradiction. Therefore, limk→∞V(E(tk),I(tk))≤V0, which in turn leads to
(4.19)lim⁡⁡sup⁡t→∞V(E(t),I(t))≤(1-p+pαk2+pmk3)N*(1-1R0)
by the definition of V0. The proof of the theorem is completed.The requirement forβ1 in Theorem 4.3 is technically strict. The interval for β1 should be much longer than the unit interval given in Theorem 4.3. Theorem 4.3 holds as long as 1-β1I(t)/N*≥0; that is, the range of β1 must guarantee that the solutions of model (4.1) are all positive.If0<β1≤1 and k1=0, limk→∞N(t)=N* and Theorem 4.3 implies that the solutions of model (3.1) are ultimately boundaries as 1<R0<1+(1/(1-p(1-γ))).Theorem4.3 provides an estimation on the ultimate boundary based on the combination of the latent and the infectious when R0>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 R0>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.Figure1(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 β1=0.3. For that set of parameter values, R0=1.35704>1, the endemic equilibrium P1*(S1*,E1*,I1*)=(1228.16000,357.57414, 80.93252), and N*=1666.66667. The ultimate boundary domain in EI-plane is given in Figure 1(b). The largest triangle in Figure 1(b) is the feasible domain E+I≤N*=1666.66667. The domain under the line V(E,I)=pαk2E+(1-p+pαk2+pmk3)I=V0 is the domain of the ultimate boundary. Three solutions of model (3.1) with initial values: E(0)=35,I(0)=200, E(0)=200,I(0)=1000, and E(0)=700,I(0)=50, are also displayed in Figure 1(b). Two solutions under the ultimate boundary V(E,I)=V0 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.1.1. Global Stability ofP10

The global stability of the disease-free equilibriumP10(E10,I10)=(0,0) of (4.1) is given in Theorem 4.1.Theorem 4.1.
IfR0<1, then the disease-free equilibrium P10 of (4.1) is globally asymptotically stable.Proof.
DefineF:[0,N*]×[0,N*]→[0,N*]×[0,N*] by
(4.2)F(E,I)=p((N*-E-I)β1IN*+(1-α)k2E+(1-m)k3E,αk2E+(1-γ)I).
Obviously, F is the mapping derived by system (4.1), and (0,0) is a fixed point of F. Linear function V(E,I)=E+(pβ1/(1-p(1-γ)))I on [0,N*]×[0,N*] is continuous and positive definite with respect to (0,0). Therefore, V is a Lyapunov function on the domain of F. For any (E,I)∈[0,N*]×[0,N*],
(4.3)ΔV(E,I)=p(N*-E-I)β1IN*+p(1-α)k2E+p(1-m)k3E+pβ11-p(1-γ)(pαk2E+p(1-γ)I)-E-pβ11-p(1-γ)I≤pβ1I-(1-p(1-α)k2-p(1-m)k3)E+p2β1αk2E1-p(1-γ)+pβ1I=(1-p(1-α)k2-p(1-m)k3)(R0-1)E≤0.
Hence, if R0<1, then ΔV(E,I)<0 holds for (E,I)∈(0,N*]×(0,N*]. It follows from Theorem 4.22 in [36] that P10 is globally asymptotically stable.Theorem4.1 ensures that the disease-free equilibrium P0(Λ/(1-p),0,0) of model (3.1) is globally stable as R0<1 and k1=0.


## 4.1.2. Stability ofP1*

The stability of the endemic equilibrium is given in Theorem4.2.Theorem 4.2.
If1<R0<1+σ, then the endemic equilibrium P1* of (4.1) is asymptotically stable, where
(4.4)σ=2p(1-αk2-mk3+(1-γ))(1-p(1-γ)+pαk2)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0.Proof.
The linearization matrix of (4.1) at the endemic equilibrium P1* is
(4.5)J1=(-pβ1N*I1*+p(1-α)k2+p(1-m)k3pβ1N*(N*-E1*-2I1*)pαk2p(1-γ)).
The corresponding characteristic equation is
(4.6)f1(λ)=λ2-b1λ+b2=0,
where
(4.7)b1=p-pαk2-pmk3+p(1-γ)-(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2,b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3)+(pαk2-p(1-γ))(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
We claim |b2|<1. When pαk2=p(1-γ), we have b2=p(1-γ)(p-pαk2-pmk3)-(1-p(1-γ))(1-p+pαk2+pmk3), and |b2|<1 holds. When pαk2≠p(1-γ), because of
(4.8)2p(1-αk2-mk3+(1-γ))(1+p(1-γ)-pαk2)<min⁡{2-p+pαk2+pmk3-p(1-γ)pαk2-p(1-γ),p(1-αk2-mk3+(1-γ))p(1-γ)-pαk2},
we obtain that |b2|<1 holds if 1<R0<1+σ. Furthermore, we have
(4.9)f1(1)=1-b1+b2=(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1),(-1)2f1(-1)=1+b1+b2=2(p-pαk2-pmk3+p(1-γ))-(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)(R0-1)1-p(1-γ)+pαk2.
It is clear that f1(1)>0 holds if R0>1 and
(4.10)(-1)2f1(-1)=σ(1+σ-R0)(1+p(1-γ)-pαk2)(1-p(1-γ))(1-p+pαk2+pmk3)>0
holds if R0<1+σ. Therefore, if 1<R0<1+σ, we have f1(1)>0, (-1)2f1(-1)>0, and |b2|<1. It follows from the Jury criteria that the endemic equilibrium P1*(E1*,I1*) of (4.1) is asymptotically stable.It follows from Theorem4.2 that when exogenous re-infection is not included (k1β2=0) the endemic equilibrium P1*(Λ/(1-p))-E1*-I1*,E1*,I1*) of model (3.1) is asymptotically stable if 1<R0<1+σ. Since the Jury criteria is the necessary and sufficient condition for the local stability of the equilibrium, the endemic equilibrium P1* of (3.1) becomes unstable if R0≥1+σ. We summarize Theorems 4.1 and 4.2 for model (3.1) without exogenous reinfection as follows.  (1)
The disease-free equilibriumP0(Λ/(1-p),0,0) is globally stable as R0<1.(2)
WhenR0>1, P0(Λ/(1-p),0,0) becomes unstable and endemic equilibrium P1* is born. The endemic P1* is locally asymptotically stable if 1<R0<1+σ.(3)
WhenR0≥1+σ,P1* becomes unstable.Theorems4.1 and 4.2 state that a forward bifurcation takes place at R0=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 β1 while R0 depends on β1 linearly. For example, given a set of the parameter values, p=0.994, Λ=10, α=0.03, γ=0.1, m=0.1, k1=0, k2=0.8, and k3=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<R0<1+σ, as shown in Figure 1(a) for β1=0.245. When β1=10.972015, we have R0=49.63168 and 1+σ=49.63168. The endemic equilibrium of model (4.1) loses the stability when β1 is greater than 10.972015, and there is a 2-period solution when β1∈[10.972015,11.109). For example, if we take β1=11, then R0=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 k1β2=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 S* 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-G1(I(t)/N(t))=β1I(t)/N(t) means that β1 should not be bigger, and we should chose β1 such that 0≤1-G1(I(t)/N(t))≤1. Otherwise, not only does 1-G1(I(t)/N(t)) lose its meaning but also the number of the susceptible individuals has negative values.The stability and the attraction of the endemic equilibrium of (3.1): (a) the endemic equilibrium and its stability and (b) the attraction area.
(a)(b)


## 4.1.3. Ultimate Boundary

We now estimate an ultimate boundary to the model (4.1) when R0>1.Theorem 4.3.
If0<β1≤1 and 1<R0<1+(1/(1-p(1-γ))), then the solutions of model (4.1) are bounded by
(4.11)lim⁡⁡sup⁡t→∞(αk2E(t)+(1-p+pαk2+pmk3)I(t))≤(1-p+pαk2+pmk3)N*(1-1R0),lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1-1R0),limsup⁡t→∞I(t)≤N*(1-1R0).Proof.
Let us consider the value ofV(E(t),I(t)) for any (E(t),I(t))∈Ω, where
(4.12)Ω={(E(t),I(t))∣E(t)≥0,I(t)≥0,E(t)+I(t)≤N*}.
Define
(4.13)V(E,I)=pαk2E+(1-p+pαk2+pmk3)I,V0=(1-p+pαk2+pmk3)N*(1-1R0).
Then,
(4.14)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=q(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t)))I(t)≤q(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t),
where q=(1-p(1-γ)(1-p+pαk2+pmk3)). Therefore, ΔV≤0 holds if V≥V0.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)≤V0 for a t0≥0, then we claim that V(E(t),I(t))≤V0 for all t≥t0. In fact, let V(E(t0),I(t0))=V0-ɛ with some ɛ≥0. Then, the inequality V(E(t0),I(t0))=V0-ɛ≤V0 yields
(4.15)I(t0)≤V01-p+pαk2+pmk3=N*(1-1R0).
The inequality in (4.14) and the condition in Theorem 4.3 as well as the inequality for I(t0) imply
(4.16)ΔV=V(E(t0+1),I(t0+1))-V(E(t0),I(t0))≤(1-p(1-γ))R0N*(V0-V(E(t0),I(t0)))I(t0),≤(1-p(1-γ))R0(V0-V(E(t0),I(t0)))(1-1R0)≤(1-p(1-γ))(R0-1)ɛ≤ɛ.
The inequality in (4.16) leads to
(4.17)V(E(t0+1),I(t0+1))≤V(E(t0),I(t0))+ɛ=V0-ɛ+ɛ=V0.
Then, our claim is followed from the mathematical induction.
IfV(E(t0),I(t0))=pαk2E(t0)+(1-p+pαk2+pmk3)I(t0)>V0 for a t0≥0, then there are two cases that we have to deal with. Case 1.
There exists an integert1>t0 such that V(E(t1),I(t1)≤V0.Case 2.
The inequalityV(E(t),I(t))>V0 holds for all t≥t0.For Case 1 we can prove that V(E(t),I(t))≤V0 for all t≥t1 by the similar argument adopted for the case V(E(t0),I(t0))≤V0. The boundary of V(E(t),I(t))≤V0 for t>t1 gives rise to limsupt→∞V(E(t),I(t))≤V0. For Case 2, the fact that ΔV(E(t),I(t))≤0 for V(E(t),I(t))>V0 leads to that the V(E(t),I(t)) is a monotonic decreasing sequence. Hence, the limit limt→∞V(E(t),I(t)) exits. We claim that limt→∞V(E(t),I(t))≤V0. Otherwise, there exists V*>V0, such that limt→∞V(E(t),I(t))=V*>V0. It follows that there exists a subsequence tk, such that limk→∞E(tk)=E*, limk→∞I(tk)=I*, and limk→∞V(E(tk),I(tk))=V(E*,I*)=V*>V0. The fact that ΔV(E*,I*)<0 and the continuity of ΔV(E,I) imply that there exists an η>0, such that ΔV(E,I)<ΔV(E*,I*)/2<0 for any (E,I) satisfying ∥(E,I)-(E*,I*)∥<η. Consequently, there exists a large integer K such that ∥(E(tk),I(tk))-(E*,I*)∥<η for k>K, |V(E(tk),I(tk))-V(E*,I*)|<|ΔV(E*,I*)|/3, and ΔV(E(tk),I(tk))<ΔV(E*,I*)/2. Then we have
(4.18)V(E(tk+1),I(tk+1))=V(E(tk),I(tk))+∑j=tktk+1-1ΔV(E(j),I(j))<V(E(tk),I(tk))+ΔV(E(tk),I(tk))≤V*+|ΔV(E*,I*)|3+ΔV(E*,I*)2.=V*+ΔV(E*,I*)6.
Inequality V(E(tk+1),I(tk+1))≤V*+ΔV(E*,I*)/6 implies that limk→∞V(E(tk),I(tk))≤V*+ΔV(E*,I*)/6, which is a contradiction. Therefore, limk→∞V(E(tk),I(tk))≤V0, which in turn leads to
(4.19)lim⁡⁡sup⁡t→∞V(E(t),I(t))≤(1-p+pαk2+pmk3)N*(1-1R0)
by the definition of V0. The proof of the theorem is completed.The requirement forβ1 in Theorem 4.3 is technically strict. The interval for β1 should be much longer than the unit interval given in Theorem 4.3. Theorem 4.3 holds as long as 1-β1I(t)/N*≥0; that is, the range of β1 must guarantee that the solutions of model (4.1) are all positive.If0<β1≤1 and k1=0, limk→∞N(t)=N* and Theorem 4.3 implies that the solutions of model (3.1) are ultimately boundaries as 1<R0<1+(1/(1-p(1-γ))).Theorem4.3 provides an estimation on the ultimate boundary based on the combination of the latent and the infectious when R0>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 R0>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.Figure1(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 β1=0.3. For that set of parameter values, R0=1.35704>1, the endemic equilibrium P1*(S1*,E1*,I1*)=(1228.16000,357.57414, 80.93252), and N*=1666.66667. The ultimate boundary domain in EI-plane is given in Figure 1(b). The largest triangle in Figure 1(b) is the feasible domain E+I≤N*=1666.66667. The domain under the line V(E,I)=pαk2E+(1-p+pαk2+pmk3)I=V0 is the domain of the ultimate boundary. Three solutions of model (3.1) with initial values: E(0)=35,I(0)=200, E(0)=200,I(0)=1000, and E(0)=700,I(0)=50, are also displayed in Figure 1(b). Two solutions under the ultimate boundary V(E,I)=V0 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

Whenk1>0, the last equation of (3.9) leads to E2*=(1-p(1-γ))I2*/((pk1β2I2*/N*)+pk2α). Substituting the expression E2* into the first equation of (3.9) yields (use I instead of I2* for convenience)(4.20)a2I2+a1I+a0=0,
where (4.21)a2=p2β1β2k1(N*)2>0,a1=p2β1k2αN*+p(β1+β2k1)(1-p(1-γ))N*-p2β1β2k1N*,≜-(R1-1)-(1-R0)(1+1-p(1-γ)-pk1β2pαk2),a0=(1-p(1-γ))(1-p+pk2α+pk3m)-p2β1k2α≜-c0(R0-1),
where(4.22)R1=pk1β2pαk2+1-p(1-γ)×1-p+pmk31-p+αk2+pmk3,c0=(1-p(1-γ))(1-p+pk2α+pk3m)>0.
The type of bifurcation at R0=1 is totally determined by R1. Near R0=1, a1<0 holds if and only if R1>1. Hence, at R0=1, if R1>1, bifurcation is forward, while if R1<1, bifurcation is backward. R1 gives a measure for exogenous reinfections. The first term in R1, 1/(pαk2+1-p(1-γ)), is the average re-infection period, which is a little bit shorter than the regular average infection period 1/(1-p(1-γ)). The second term in R1 gives the proportion of the latent population. These analysis result in the following theorem that characterizes the feature of bifurcation at R0=1.Theorem 4.4.
Consider model (3.9). If R1<1, then the bifurcation at R0=1 is of forward type. If R1>1, the bifurcation at R0=1 is of backward type.Looking at the expression ofR1, one can see that R1>1 necessarily requires k1β2≠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 p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, and  m=0.1. Figure 2(a) gives three domains in β1β2-plane with the unique disease-free equilibrium (DFE), two positive endemicequilibria(E2±*,I2±*), and one endemic equilibrium (E2+*,I2+*). Taking the above parameter values and fixing β2=0.5, we get the backward bifurcation curve of model (3.9) (see Figure 2(b)).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).
(a)(b)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 P2±*(E2±*,I2±*) is(4.23)J2±=(p-pk2α-pk3m-p(β1+β2k1)I2±*N*pβ1-p(β1+β2k1)E2±*N*-2pβ1I2±*N*pβ2k1I2±*N*+pαk2pβ2k1E2±*N*+p(1-γ)).
We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1,andβ2=0.5 and change the value of β1. When β1∈[0.31831317735,0.41278978), we have R0∈[0.7711266,1), and there are two positive equilibrium points, which are (E2+*,I2+*) and (E2-*,I2-*), respectively. Fixing β1=0.35, we have (E2+*,I2+*)=(601.51345,342.43485) and (E2-*,I2-*)=(329.51923,52.38933). Accordingly, we have ρ(J2+)=0.979777<1 and ρ(J2-)=1.00415>1. Therefore, (E2+*,I2+*) is local stable as R0<1, while (E2-*,I2-*) is unstable as R0<1. When β1>0.41278978, we have R0>1, and there is only one positive equilibrium point (E2+*,I2+*). Fixing β1=1, we have (E2+*,I2+*)=(657.40092,796.73372), and ρ(J2+)=0.88069<1. ρ(J2+) will change when β1 becomes bigger. The numerical simulation shows that ρ(J2+)=1 and R0=9.91134 as β1=4.0913. Therefore, we obtain that (E2+*,I2+*) is locally stable as 1<R0<9.91134, and, when R0≥9.91134, (E2+*,I2+*) 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 inEI-plane when R0<1. Here, we use the same parameter values as for the analysis of the endemic equilibrium bifurcation when R0<1. There are three equilibrium points, the disease-free equilibrium (0,0), the endemic equilibria (E2+*,I2+*), and (E2-*,I2-*), 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 (E2+*,I2+*), and the attraction domain of the disease-free equilibrium (0,0). Figure 3 shows that (E2-*,I2-*) 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 (E2+*,I2+*). It appears that the attraction domain of (E2+*,I2+*) 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 R0<1.Figure 3
The attraction separatrix for those two equilibria of model (3.1).


### 4.2.2. Model Behavior for LargeR0

This subsection is mainly concerned with the behavior of the model whenR0 is very large. Numerical approach is used to further explore the model. We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, β2=0.5, and let β1 change. Figure 4(a) shows the existence of the stable endemic equilibrium of when R0<1. In fact, model (3.9) has two endemic equilibrium points when 0.77113<R0<1, the small endemic equilibrium is unstable, and the large one is locally stable. If 1<R0<8.6485, model (3.9) has only one endemic equilibrium, which is locally stable. When 8.6485≤R0<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≤R0<9.9566, and a stable periodic solution with period 4 appears. When R0 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 R0 increases.Period doubling bifurcation of model (3.9). Parameter values are p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, andβ2=0.5.
(a)(b)Figure4 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 R0∈[0,10]. However, the number of the latent class may be negative when R0≥8.6485. R0=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 R0 becomes larger. The numerical simulation seems to hint that the endemic equilibrium may be stable if β1≤1 and β2≤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.
If1<R0≤1+((pβ1αk2-(1-p(1-γ))(1-p+pmk3)k1β2)/(1-p(1-γ))(pβ1αk2+(1-p+pmk3)k1β2)), pβ1αk2>(1-p(1-γ))(1-p+pmk3)k1β2, and 0<β1,β2≤1 hold, then the solutions of model (3.9) are ultimately bounded by
(4.24)lim⁡⁡sup⁡t→∞V(E(t),I(t))=lim⁡⁡sup⁡t→∞pαk2E(t)+(1-p+pαk2+pmk3)I(t)≤V0,whereV0≜(1-p+pαk2+pmk3)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).The inequalitylimsupt→∞V(E(t),I(t))≤V0 implies that(4.25)lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0),lim⁡⁡sup⁡t→∞I(t)≤N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).
Introducing V(E(t),I(t))=pαk2E(t)+(1-p+pαk2+pmk3)I(t) and then examining the difference of V yields that(4.26)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t))+(1-p+pmk3)k1β2R0E(t)pβ1αk2N*)I(t)≤(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3+(1-p+pmk3)k1β2R0pβ1αk2)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t).
The rest of the proof is similar to that of Theorem 4.4.


## 4.2.1. Backward Bifurcation

Whenk1>0, the last equation of (3.9) leads to E2*=(1-p(1-γ))I2*/((pk1β2I2*/N*)+pk2α). Substituting the expression E2* into the first equation of (3.9) yields (use I instead of I2* for convenience)(4.20)a2I2+a1I+a0=0,
where (4.21)a2=p2β1β2k1(N*)2>0,a1=p2β1k2αN*+p(β1+β2k1)(1-p(1-γ))N*-p2β1β2k1N*,≜-(R1-1)-(1-R0)(1+1-p(1-γ)-pk1β2pαk2),a0=(1-p(1-γ))(1-p+pk2α+pk3m)-p2β1k2α≜-c0(R0-1),
where(4.22)R1=pk1β2pαk2+1-p(1-γ)×1-p+pmk31-p+αk2+pmk3,c0=(1-p(1-γ))(1-p+pk2α+pk3m)>0.
The type of bifurcation at R0=1 is totally determined by R1. Near R0=1, a1<0 holds if and only if R1>1. Hence, at R0=1, if R1>1, bifurcation is forward, while if R1<1, bifurcation is backward. R1 gives a measure for exogenous reinfections. The first term in R1, 1/(pαk2+1-p(1-γ)), is the average re-infection period, which is a little bit shorter than the regular average infection period 1/(1-p(1-γ)). The second term in R1 gives the proportion of the latent population. These analysis result in the following theorem that characterizes the feature of bifurcation at R0=1.Theorem 4.4.
Consider model (3.9). If R1<1, then the bifurcation at R0=1 is of forward type. If R1>1, the bifurcation at R0=1 is of backward type.Looking at the expression ofR1, one can see that R1>1 necessarily requires k1β2≠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 p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, and  m=0.1. Figure 2(a) gives three domains in β1β2-plane with the unique disease-free equilibrium (DFE), two positive endemicequilibria(E2±*,I2±*), and one endemic equilibrium (E2+*,I2+*). Taking the above parameter values and fixing β2=0.5, we get the backward bifurcation curve of model (3.9) (see Figure 2(b)).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).
(a)(b)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 P2±*(E2±*,I2±*) is(4.23)J2±=(p-pk2α-pk3m-p(β1+β2k1)I2±*N*pβ1-p(β1+β2k1)E2±*N*-2pβ1I2±*N*pβ2k1I2±*N*+pαk2pβ2k1E2±*N*+p(1-γ)).
We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1,andβ2=0.5 and change the value of β1. When β1∈[0.31831317735,0.41278978), we have R0∈[0.7711266,1), and there are two positive equilibrium points, which are (E2+*,I2+*) and (E2-*,I2-*), respectively. Fixing β1=0.35, we have (E2+*,I2+*)=(601.51345,342.43485) and (E2-*,I2-*)=(329.51923,52.38933). Accordingly, we have ρ(J2+)=0.979777<1 and ρ(J2-)=1.00415>1. Therefore, (E2+*,I2+*) is local stable as R0<1, while (E2-*,I2-*) is unstable as R0<1. When β1>0.41278978, we have R0>1, and there is only one positive equilibrium point (E2+*,I2+*). Fixing β1=1, we have (E2+*,I2+*)=(657.40092,796.73372), and ρ(J2+)=0.88069<1. ρ(J2+) will change when β1 becomes bigger. The numerical simulation shows that ρ(J2+)=1 and R0=9.91134 as β1=4.0913. Therefore, we obtain that (E2+*,I2+*) is locally stable as 1<R0<9.91134, and, when R0≥9.91134, (E2+*,I2+*) 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 inEI-plane when R0<1. Here, we use the same parameter values as for the analysis of the endemic equilibrium bifurcation when R0<1. There are three equilibrium points, the disease-free equilibrium (0,0), the endemic equilibria (E2+*,I2+*), and (E2-*,I2-*), 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 (E2+*,I2+*), and the attraction domain of the disease-free equilibrium (0,0). Figure 3 shows that (E2-*,I2-*) 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 (E2+*,I2+*). It appears that the attraction domain of (E2+*,I2+*) 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 R0<1.Figure 3
The attraction separatrix for those two equilibria of model (3.1).


## 4.2.2. Model Behavior for LargeR0

This subsection is mainly concerned with the behavior of the model whenR0 is very large. Numerical approach is used to further explore the model. We fix parameter values p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, β2=0.5, and let β1 change. Figure 4(a) shows the existence of the stable endemic equilibrium of when R0<1. In fact, model (3.9) has two endemic equilibrium points when 0.77113<R0<1, the small endemic equilibrium is unstable, and the large one is locally stable. If 1<R0<8.6485, model (3.9) has only one endemic equilibrium, which is locally stable. When 8.6485≤R0<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≤R0<9.9566, and a stable periodic solution with period 4 appears. When R0 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 R0 increases.Period doubling bifurcation of model (3.9). Parameter values are p=0.994, Λ=10, α=0.03, γ=0.1, k1=0.5, k2=0.3, k3=0.2, m=0.1, andβ2=0.5.
(a)(b)Figure4 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 R0∈[0,10]. However, the number of the latent class may be negative when R0≥8.6485. R0=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 R0 becomes larger. The numerical simulation seems to hint that the endemic equilibrium may be stable if β1≤1 and β2≤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.
If1<R0≤1+((pβ1αk2-(1-p(1-γ))(1-p+pmk3)k1β2)/(1-p(1-γ))(pβ1αk2+(1-p+pmk3)k1β2)), pβ1αk2>(1-p(1-γ))(1-p+pmk3)k1β2, and 0<β1,β2≤1 hold, then the solutions of model (3.9) are ultimately bounded by
(4.24)lim⁡⁡sup⁡t→∞V(E(t),I(t))=lim⁡⁡sup⁡t→∞pαk2E(t)+(1-p+pαk2+pmk3)I(t)≤V0,whereV0≜(1-p+pαk2+pmk3)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).The inequalitylimsupt→∞V(E(t),I(t))≤V0 implies that(4.25)lim⁡⁡sup⁡t→∞E(t)≤(1+1-p+pmk3pαk2)N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0),lim⁡⁡sup⁡t→∞I(t)≤N*(1+(1-p+pmk3)k1β2pβ1αk2-1R0).
Introducing V(E(t),I(t))=pαk2E(t)+(1-p+pαk2+pmk3)I(t) and then examining the difference of V yields that(4.26)ΔV=V(E(t+1),I(t+1))-V(E(t),I(t))=(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*1-p+pαk2+pmk31-p+pαk2+pmk3(E(t)+I(t))+(1-p+pmk3)k1β2R0E(t)pβ1αk2N*)I(t)≤(1-p(1-γ)(1-p+pαk2+pmk3))×(R0-1-R0N*V(E(t),I(t))1-p+pαk2+pmk3+(1-p+pmk3)k1β2R0pβ1αk2)I(t),=(1-p(1-γ))R0N*(V0-V(E(t),I(t)))I(t).
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

LetP3*(E3*,I3*) be the positive equilibrium of model (3.10). Then the second equation of model (3.10) leads to I3*=pαk2E3*/(1-p(1-γ)-(pβ2k1E3*/N*)), where E3*<N*(1-p(1-r))/pβ2k1.We already have a clear picture about the equilibrium in the regionE(t)≤E0. Now we only seek for equilibria in the region E(t)>E0 with T(E(t))=p(k3E(t)-m0). Similarly, E3* satisfies (4.27)σ3E3+σ2E2+σ1E+σ0=0,
where
(4.28)σ3=p2β2k1(pβ1αk2-β2k1(1-p))(N*)2,σ2=2pβ2k1(1-p+pαk2)N*(1-p(1-γ))+p2αk2(β2k1+β1)N*(1-p(1-γ))-p3β1αk2(β2k1+αk2)N*-p3m0β22k12(N*)2,σ1=(1-p(1-γ))(2p2m0β2k1N*+p2β1αk2-(1-p(1-γ))(1-p+pαk2)),σ0=-pm0(1-p(1-γ))2.From the continuity argument, we can imposeE0 to be a root to (4.27). Obviously, R0 and R1 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<0 relieves our work a little bit. If σ3<0, then (4.27) has exact two positive roots. If σ3>0 and σ1<0, then (4.27) has exact one positive root. If σ3>0 and σ1>0, then (4.27) has exact one or three positive roots. For all situations, E0 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 valuesp=0.994, Λ=10, k1=0.2, k2=0.4, k3=0.4, α=0.03, γ=0.0001, m=0.005, and β2=0.375.As can be seen in Figure5, a backward bifurcation occurs at R0=1 because R1=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 R0<1 and when saturated treatment rate is used. Figure 6 shows the tristable situation on the phase plane.Double backward bifurcations of model (3.10). Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply Backward Bifurcation(b)
A part of (a)Phase portrait of model (3.10): a case of multiple equilibria.
(a)(b)If we chooseβ2=0.26, thenR1=1.1458>1. The bifurcation at R0=1 remains backward. However, the second backward bifurcation shifts to the right, happening at a place of R0>1. One can observe bistability not only for R0<1 but also for R0>1. Figure 7(a) illustrates this bistable dynamics. Further decreasing of β2 can produce a new bifurcation, which will be studied in the  next subsection.(a) Double backward bifurcations for model (3.10). The upper one occurs when R0>1. (b) Forward-backward bifurcation for model (3.10). The upper one occurs when R0<1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply backward bifurcation(b)
Forward-backward bifurcation


### 4.3.3. Forward-Backward Bifurcation

If we setp=0.994, Λ=8, k1=0.3, k2=0.3, k3=0.4,α=0.03,γ=0.001, m=0.01, m0=0.2, and β2=0.1, then the bifurcation at R0=1 is forward because R1=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 R0<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 R0<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 R0<1.A less important forward-backward bifurcation is the case where the whole backward bifurcation curve exists only in the region ofR0>1, thus only making endemic equilibrium bigger. Figure 8 shows this kind of bifurcation, where β2=0.18 in order to make R1=0.7933<1 and the rest of parameters are the same as in the previous subsection.Forward-backward bifurcation for model (3.10). The upper one occurs when R0>1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Forward-backward bifurcation(b)
A part of (a)


## 4.3.1. Endemic Equilibria

LetP3*(E3*,I3*) be the positive equilibrium of model (3.10). Then the second equation of model (3.10) leads to I3*=pαk2E3*/(1-p(1-γ)-(pβ2k1E3*/N*)), where E3*<N*(1-p(1-r))/pβ2k1.We already have a clear picture about the equilibrium in the regionE(t)≤E0. Now we only seek for equilibria in the region E(t)>E0 with T(E(t))=p(k3E(t)-m0). Similarly, E3* satisfies (4.27)σ3E3+σ2E2+σ1E+σ0=0,
where
(4.28)σ3=p2β2k1(pβ1αk2-β2k1(1-p))(N*)2,σ2=2pβ2k1(1-p+pαk2)N*(1-p(1-γ))+p2αk2(β2k1+β1)N*(1-p(1-γ))-p3β1αk2(β2k1+αk2)N*-p3m0β22k12(N*)2,σ1=(1-p(1-γ))(2p2m0β2k1N*+p2β1αk2-(1-p(1-γ))(1-p+pαk2)),σ0=-pm0(1-p(1-γ))2.From the continuity argument, we can imposeE0 to be a root to (4.27). Obviously, R0 and R1 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<0 relieves our work a little bit. If σ3<0, then (4.27) has exact two positive roots. If σ3>0 and σ1<0, then (4.27) has exact one positive root. If σ3>0 and σ1>0, then (4.27) has exact one or three positive roots. For all situations, E0 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 valuesp=0.994, Λ=10, k1=0.2, k2=0.4, k3=0.4, α=0.03, γ=0.0001, m=0.005, and β2=0.375.As can be seen in Figure5, a backward bifurcation occurs at R0=1 because R1=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 R0<1 and when saturated treatment rate is used. Figure 6 shows the tristable situation on the phase plane.Double backward bifurcations of model (3.10). Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply Backward Bifurcation(b)
A part of (a)Phase portrait of model (3.10): a case of multiple equilibria.
(a)(b)If we chooseβ2=0.26, thenR1=1.1458>1. The bifurcation at R0=1 remains backward. However, the second backward bifurcation shifts to the right, happening at a place of R0>1. One can observe bistability not only for R0<1 but also for R0>1. Figure 7(a) illustrates this bistable dynamics. Further decreasing of β2 can produce a new bifurcation, which will be studied in the  next subsection.(a) Double backward bifurcations for model (3.10). The upper one occurs when R0>1. (b) Forward-backward bifurcation for model (3.10). The upper one occurs when R0<1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Multiply backward bifurcation(b)
Forward-backward bifurcation


## 4.3.3. Forward-Backward Bifurcation

If we setp=0.994, Λ=8, k1=0.3, k2=0.3, k3=0.4,α=0.03,γ=0.001, m=0.01, m0=0.2, and β2=0.1, then the bifurcation at R0=1 is forward because R1=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 R0<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 R0<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 R0<1.A less important forward-backward bifurcation is the case where the whole backward bifurcation curve exists only in the region ofR0>1, thus only making endemic equilibrium bigger. Figure 8 shows this kind of bifurcation, where β2=0.18 in order to make R1=0.7933<1 and the rest of parameters are the same as in the previous subsection.Forward-backward bifurcation for model (3.10). The upper one occurs when R0>1. Dotted lines represent unstable equilibria, and solid lines represent stable ones.
(a)
Forward-backward bifurcation(b)
A part of (a)


## 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 whenR0 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) whenR0<1, bistability for both R0<1 and R0>1, and period doubling for a reduced model.We now know thatk1β2 controls the appearance of the backward bifurcation at R0=1, creating an endemic equilibrium when R0<1, an example of bistable scenario. The treatment regime does not participate in the backward bifurcation R0=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 R0. 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 R0=1 is forward, meaning that R0<1 makes the disease die out, and if the treatment practice is changed through reducing the treatment effort, then an endemic is born when R0<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.



---

*Source: 653937-2011-11-30.xml*