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 , Ricker's model [13], and Hassell model [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, , 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 [715]. Meanwhile, discrete epidemic models have gained more popularity [5, 6, 1625], 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 [2629], 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 . Incorporating the exogenous reinfections and treatments into the models in [20, 21], our model framework takes the following form: 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. is the probability that susceptible individuals become infected. 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 . and have been used for the probabilities of not becoming infected [20, 21]. We choose , 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: 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 , the number of the susceptibles who get infected and become latent individuals is ; the number of the latent who get re-infected and become infectious is ; 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 is . Similarly, the number of infectious individuals who do not change their epidemiological status from to is counted by .

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: On the other hand, in the case of limited resources, the treatment function of latent individuals is a piecewise function, where . Model (2.1) takes the following form: where

The disease-free equilibrium, , and the basic reproductive number of models (3.1) and (3.3) are identical. Applying the approach in [34] to our model, we have At the disease-free equilibrium, the straight forward calculation yields and . Therefore, the basic reproductive number of model (2.1) is Rearrange terms in the expression of as Each term in has clear epidemiological interpretation. is the average infection period. is average new cases generated by a typical infectious member in the entire infection period. 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 and unstable if (see [34, Theorem  2.1]). We summarize stability results in terms of as in the following theorem.

Theorem 3.1. The disease-free equilibrium of model (2.1) is locally asymptotically stable if and unstable if .

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 : One can see that is a global attractor for (3.8). Using limiting equations and the [35], we, respectively, reduce the three-dimensional systems (3.1) and (3.3) into two-dimensional ones 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, , , and model (3.9) is reduced to The disease-free equilibrium of (4.1) is and the endemic equilibrium exists if , where and .

4.1.1. Global Stability of

The global stability of the disease-free equilibrium of (4.1) is given in Theorem 4.1.

Theorem 4.1. If , then the disease-free equilibrium of (4.1) is globally asymptotically stable.

Proof. Define by Obviously, is the mapping derived by system (4.1), and is a fixed point of . Linear function on is continuous and positive definite with respect to . Therefore, is a Lyapunov function on the domain of . For any , Hence, if , then holds for . It follows from Theorem 4.22 in [36] that is globally asymptotically stable.

Theorem 4.1 ensures that the disease-free equilibrium of model (3.1) is globally stable as and .

4.1.2. Stability of

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

Theorem 4.2. If , then the endemic equilibrium of (4.1) is asymptotically stable, where

Proof. The linearization matrix of (4.1) at the endemic equilibrium is The corresponding characteristic equation is where We claim . When , we have , and holds. When , because of we obtain that holds if . Furthermore, we have It is clear that holds if and holds if . Therefore, if , we have , , and . 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 () the endemic equilibrium of model (3.1) is asymptotically stable if . 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 . We summarize Theorems 4.1 and 4.2 for model (3.1) without exogenous reinfection as follows. The disease-free equilibrium is globally stable as .When , becomes unstable and endemic equilibrium is born. The endemic is locally asymptotically stable if .When becomes unstable.

Theorems 4.1 and 4.2 state that a forward bifurcation takes place at . 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, , , , , , , , and , 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 , as shown in Figure 1(a) for . When , we have and . The endemic equilibrium of model (4.1) loses the stability when is greater than 10.972015, and there is a 2-period solution when . For example, if we take , then , and . The 2-period solution is . 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 , and the numerical computation shows that the corresponding periodic solution of model (3.1) is , . 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, means that should not be bigger, and we should chose such that . Otherwise, not only does lose its meaning but also the number of the susceptible individuals has negative values.

4.1.3. Ultimate Boundary

We now estimate an ultimate boundary to the model (4.1) when .

Theorem 4.3. If and , then the solutions of model (4.1) are bounded by

Proof. Let us consider the value of for any , where Define Then, where . Therefore, holds if .
If for a , then we claim that for all . In fact, let with some . Then, the inequality yields The inequality in (4.14) and the condition in Theorem 4.3 as well as the inequality for imply The inequality in (4.16) leads to Then, our claim is followed from the mathematical induction.
If for a , 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 . For Case 2, the fact that for leads to that the is a monotonic decreasing sequence. Hence, the limit exits. We claim that . Otherwise, there exists , such that . It follows that there exists a subsequence , such that , , and . The fact that and the continuity of imply that there exists an , such that for any satisfying . Consequently, there exists a large integer such that for , , and . Then we have Inequality implies that , which is a contradiction. Therefore, , which in turn leads to 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 ; that is, the range of must guarantee that the solutions of model (4.1) are all positive.

If and , and Theorem 4.3 implies that the solutions of model (3.1) are ultimately boundaries as .

Theorem 4.3 provides an estimation on the ultimate boundary based on the combination of the latent and the infectious when . 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 . 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 . For that set of parameter values, , the endemic equilibrium , , and . The ultimate boundary domain in -plane is given in Figure 1(b). The largest triangle in Figure 1(b) is the feasible domain . The domain under the line is the domain of the ultimate boundary. Three solutions of model (3.1) with initial values: , , and , 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 , the last equation of (3.9) leads to . Substituting the expression into the first equation of (3.9) yields (use instead of for convenience) where where The type of bifurcation at is totally determined by . Near , holds if and only if . Hence, at , if , bifurcation is forward, while if , bifurcation is backward. gives a measure for exogenous reinfections. The first term in , , is the average re-infection period, which is a little bit shorter than the regular average infection period . 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 .

Theorem 4.4. Consider model (3.9). If , then the bifurcation at is of forward type. If , the bifurcation at is of backward type.

Looking at the expression of , one can see that necessarily requires . 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 , , , , , , , and . Figure 2(a) gives three domains in -plane with the unique disease-free equilibrium (DFE), two positive endemicequilibria, and one endemic equilibrium . Taking the above parameter values and fixing , we get the backward bifurcation curve of model (3.9) (see Figure 2(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 is We fix parameter values , , , , , , , and change the value of . When , we have , and there are two positive equilibrium points, which are and , respectively. Fixing , we have and . Accordingly, we have and . Therefore, is local stable as , while is unstable as . When , we have , and there is only one positive equilibrium point . Fixing , we have , and . will change when becomes bigger. The numerical simulation shows that and as . Therefore, we obtain that is locally stable as , and, when , 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 . Here, we use the same parameter values as for the analysis of the endemic equilibrium bifurcation when . There are three equilibrium points, the disease-free equilibrium , 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 . 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 , 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 . This fact implies that we should be careful in assessing a disease control program, since the attraction domain of is small, the disease may not be eliminated even though .

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 , , , , , , , , , and let change. Figure 4(a) shows the existence of the stable endemic equilibrium of when . In fact, model (3.9) has two endemic equilibrium points when , the small endemic equilibrium is unstable, and the large one is locally stable. If , model (3.9) has only one endemic equilibrium, which is locally stable. When 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 , 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.

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 . However, the number of the latent class may be negative when . 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 and , 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 , , and hold, then the solutions of model (3.9) are ultimately bounded by

The inequality implies that Introducing and then examining the difference of yields that 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 , where .

We already have a clear picture about the equilibrium in the region . Now we only seek for equilibria in the region with . Similarly, satisfies where

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 relieves our work a little bit. If , then (4.27) has exact two positive roots. If and , then (4.27) has exact one positive root. If and , 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 , , , , , , , , and .

As can be seen in Figure 5, a backward bifurcation occurs at because . 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 and when saturated treatment rate is used. Figure 6 shows the tristable situation on the phase plane.

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

4.3.3. Forward-Backward Bifurcation

If we set , , , , , , , and , then the bifurcation at is forward because , 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 , 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 . 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 .

A less important forward-backward bifurcation is the case where the whole backward bifurcation curve exists only in the region of , thus only making endemic equilibrium bigger. Figure 8 shows this kind of bifurcation, where in order to make and the rest of parameters are the same as in the previous subsection.

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 , bistability for both and , and period doubling for a reduced model.

We now know that controls the appearance of the backward bifurcation at , creating an endemic equilibrium when , an example of bistable scenario. The treatment regime does not participate in the backward bifurcation . 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 is forward, meaning that makes the disease die out, and if the treatment practice is changed through reducing the treatment effort, then an endemic is born when . 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.


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