Successive Vaccination and Difference in Immunity of a Delay SIR Model with a General Incidence Rate

Abstract

A delay SIR epidemic model with difference in immunity and successive vaccination is proposed to understand their effects on the disease spread. From theorems, it is obtained that the basic reproduction number governs the dynamic behavior of the system. The existence and stability of the possible equilibria are examined in terms of a certain threshold condition about the basic reproduction number. By use of new computational techniques for delay differential equations, we prove that the system is permanent. Our results indicate that the recovery rate and the vaccination rate are two factors for the dynamic behavior of the system. Numerical simulations are carried out to investigate the influence of the key parameters on the spread of the disease, to support the analytical conclusion, and to illustrate possible behavioral scenarios of the model.

1. Introduction

The current threat of some new type diseases has raised our awareness that curbing the spread of some emerging and reemerging human diseases is of public health importance such as H1N1. This emerging disease, which was first reported in Mexico, spread very quickly, due to the travel of infected persons by airplanes, trains, and buses to some other regions. It continued to spread around the world and caused about 5000 deaths. In recent years, many mathematical models have been developed for the transmission dynamics of infectious diseases such as SARS, HIV/AIDS, measles, and smallpox ([1–6], to name a few, and the references therein). These models have provided understanding of the underlying mechanisms which influence the spread of diseases and suggested some control strategies. Moreover, to our knowledge, the first effective control strategy for the elimination of infectious diseases is obtaining immunity. It has been reported that “People’s immunity to A/H1N1 flu virus is greater than previously thought after access vaccines. The WHO is working to give more nations access vaccines to fight the H1N1 flu pandemic.”

Besides these above studies, many authors formulated and analyzed SIR epidemic models for the control of diseases [7–20]. In particular, some authors have studied the effects of vaccination on the spread of diseases [7–11]; others have studied the effects of treatment on the spread of diseases [12–15]. Gao et al. have proposed an epidemic model with density-dependent birth pulses and seasonal prevention [16]. Recently, some works have investigated permanent and temporary immunity [17–20]. However, in these SIR models, an unrealistic assumption is that all the rest of infected individuals acquire immunity besides death. Measles encephalitis in adults in [21, 22] shows that there is difference in immunity of infected individuals. That is, some infected individuals can acquire immunity after recovery, but some do not acquire immunity and can be infected once more. At the same time, vaccination is an important strategy for the elimination of infectious diseases [7–11].

Vaccinations have many types; impulsive vaccination and successive vaccination are two main policies. Successive vaccination is that people have been vaccinated at birth to protect themselves from disease; the studies can be found in [23, 24]. Makinde in [23] studied a SIR model for the transmission dynamics of a childhood disease in the presence of a preventive vaccine and analyzed the vaccination reproductive number for disease control and eradication qualitatively. Impulsive vaccination (only at fixed time sequence we execute effectively the vaccination for the disease) is an important and effective strategy for the elimination of infectious diseases and has been studied in the literature. For example, see [10, 17, 25–27]. In above-mentioned papers, authors almost considered the vaccination of susceptible population. But, in fact, under a certain situation, the vaccine treatment also should be considered for the newborns of the susceptible, the exposed, and the removed. We find that there are few studies on the aspect of the vaccination of newborns. In this case, successive vaccination seems more reasonable than impulsive vaccination. Therefore, in this paper, a SIR model with difference in immunity and successive vaccination is considered.

As far as disease transmission is concerned, the incidence rate, defined as the rate of new infection, plays a very important role in modelling infectious diseases. Bilinear incidence rate 𝛽𝑆𝐼 in [28, 29] and standard incidence rates 𝜆𝑆𝐼/𝑁 in [30, 31] have often been used in epidemic models. However, it is unreasonable to consider the bilinear incidence rate (based on the law of mass action) as the number of susceptibles is large, owing to the number of susceptibles with which every infective contact within a certain time is limited. Standard incidence rate may be a good approximation if the number of available partners is large enough but it is not possible to make more contacts when the population 𝑁 is small. Combine the two previous approaches by assuming that if the number of available partners 𝑁 is low, the number of actual per capita partners is proportional to 𝑁, whereas if the number of available partners is large, there is a saturation effect which makes the number of actual partners constant. Considering this case, a saturation incidence rate of type 𝑓(𝐼)𝑆 with 𝑓(𝐼)=𝑘𝐼/(1+𝛼𝐼) is being proposed in [32]. More general incidence rate used in the literature is the one for which 𝑓(𝐼)=𝑘𝐼^𝑙/(1+𝛼𝐼^ℎ) [33, 34], where 𝐼^𝑙 measures the infection force of the disease and 𝑓(𝐼)=1/(1+𝛼𝐼^ℎ) measures the inhibitory effect caused by behavioral changes. Note that if 𝑓(𝐼) is decreasing when 𝐼 is large, this may be interpreted as the fact that susceptibles tend to reduce their social contacts if the perceived number of infectives increases over a psychologically significant value. The above saturation incidence rate depends also on the size of the infectives 𝐼 termed as infectives-dependent. Particular examples of susceptibles-dependent incidence rate are 𝑓(𝑆)=𝑘𝑆/(1+𝛼𝑆) [35]. Very general incidence rates which are not linear in 𝑆 are also used in Derrick and van den Driessche [36] (𝑓(𝑆,𝐼,𝑁)=𝐼Φ(𝑆,𝐼,𝑁), where 𝑁=𝑆+𝐼), Korobeinikov and Maini [37] (𝑓(𝑆,𝐼)=ℎ₁(𝐼)ℎ₂(𝑆)), and Moghadas and Alexander [38] (𝑓(𝑆,𝐼)=𝛽(1+𝑔(𝐼,𝑣))𝐼𝑆).

Nie et al. [19] and Ji et al. [39] respectively considered a delay SIR epidemic model with nonlinear incidence rate and density-dependent birth and death rates. Motivated by the main idea described in [6, 39], in this paper, we consider a delay SIR model with difference in immunity and successive vaccination and an abstract incidence rate. The main difference between our study and those described in [6, 39] is the difference in immunity and successive vaccination and an abstract incidence rate. An abstract incidence rate of type 𝑓(𝐼)𝑆 is employed to model the spread of the disease which is propagated through the infective individuals, under a few biologically feasible assumptions upon 𝑓(𝐼).

In view of above facts, we will formulate a mathematical model in Section 2. We provide the region of biologically feasible solutions in Section 3. Then, we study the existence and stability of the steady states in the next section, analyze the permanence result in Section 5, and give some numerical simulations in Section 6. Lastly, we end the paper with a brief discussion of our results in Section 6.

2. Model Formulation and Invariant Region

In this section, we will present a delay SIR epidemic model with a general nonlinear incidence rate. The total population 𝑁(𝑡) is divided into three subclasses, namely, the susceptibles 𝑆(𝑡), the infectives 𝐼(𝑡), and the recovered individuals 𝑅(𝑡). Based on the SIR model in [12, 39], we considered following system: 𝑑𝑆=𝑑𝑡𝑏−𝑎𝑟𝑁𝐾𝑁−𝛽𝑒^−𝑑1𝜏(𝐼(−()𝑆𝑓𝑡−𝜏))𝑑+1−𝑎𝑟𝑁𝐾𝑆−𝜃𝑆+𝜇₁𝐼,𝑑𝐼𝑑𝑡=𝛽𝑒^−𝑑1𝜏𝑓(𝐼(𝑆−()𝑡−𝜏))𝑑+1−𝑎𝑟𝑁𝐾𝐼−𝜇₁𝐼−𝑒^−𝑑1𝜔𝛿𝐼(𝑡−𝜔),𝑑𝑅𝑑𝑡=𝑒^−𝑑1𝜔()−()𝛿𝐼𝑡−𝜔𝑑+1−𝑎𝑟𝑁𝐾𝑁(0)𝑅+𝜃𝑆,=𝑆₀(0)(𝜃)[−>0,𝐼=𝜑≥0,∀𝜃∈𝜏,0],𝑅(0)=𝑅₀≥0, (1) where 𝜏=max{𝜏,𝜔} and 𝜑∈𝐶([−𝜏,0],𝑅). We give the following useful assumptions.(1)There are no disease induced deaths, and all the newborns are susceptible.(2)𝑓(𝐼) is the nonlinear incidence rate satisfying the following assumptions: 𝑓(0)=0,𝑓^′(𝐼)>0,𝑓^′′(𝐼)<0,lim_𝑡→∞𝑓(𝐼)=𝑐<+∞. (2)(3)The force of infection at any time 𝑡 is dominated by 𝛽𝑒^−𝑑1𝜏𝑆(𝑡)𝑓(𝐼(𝑡−𝜏)), where 𝜏 is incubation period and 0<𝑒^−𝑑1𝜏≤1 represents the survival probability of individuals in the population after time 𝜏 [20]. It is also assumed that 𝑑₁≤𝑑 in [−𝜏,0], where 𝑑 is the death rate and 𝑑₁ is the death rate in the time interval [−𝜏,0].(4)The parameters 𝑎 are a convex combination constant, 𝑟=𝑏−𝑑>0 is the intrinsic growth rate (𝑏 is the birth rate), and 𝐾>0 is the carrying capacity of the population. The term (𝑏−(𝑎𝑟𝑁(𝑡)/𝐾)) has a density-dependent per capita birth rate and the term (𝑑+((1−𝑎)𝑟𝑁(𝑡)/𝐾)) has a density-dependent per capita death rate [39].(5)For 0<𝑎<1, the birth and death rates are consistent with the limited resources associated with density dependence. The birth rate is density independent when 𝑎=0 and the death rate is density independent when 𝑎=1. Thus, the spread of the disease (animals such as rodents, etc.) is assumed to be governed by the following system of logistic equations with time delay.(6)The total population is assumed to be large enough to be adequately described by a deterministic model and is divided into compartments based on the disease status [40].(7)The successive vaccination rate 𝜃 is positive. The positive constant 𝜇₁ is the recovery rate of the infectious individuals from compartment 𝐼 to 𝑆. The parameters 𝛽 are the effective per capita contact rate constant of infected individuals. The parameters 𝛿 are the recovery rate of infected individuals.(8)Models are formulated as functional differential and/or integral equations when time delay is included [40]. Ours follows the former with the assumption that the 𝐼-equation satisfies a certain integral condition [41].

Models with multiple delays are not common, but few authors have in the past considered these-Beretta et al. [42]-to name but a few. Since 𝑁(𝑡)=𝑆(𝑡)+𝐼(𝑡)+𝑅(𝑡), thus the governing equation (1) can be rewritten as 𝑑𝑁𝑁𝑑𝑡=𝑟1−𝐾𝑁,𝑑𝐼𝑑𝑡=𝛽𝑒^−𝑑1𝜏()𝑓(𝐼(−()𝑁−𝐼−𝑅𝑡−𝜏))𝑑+1−𝑎𝑟𝑁𝐾𝐼−𝜇₁𝐼−𝑒^−𝑑1𝜔𝛿𝐼(𝑡−𝜔),𝑑𝑅𝑑𝑡=𝑒^−𝑑1𝜔()−()𝛿𝐼𝑡−𝜔𝑑+1−𝑎𝑟𝑁𝐾𝑅+𝜃(𝑁−𝐼−𝑅). (3) Let 𝜏=max(𝜔,𝜏). Then (3) satisfies the following initial conditions 𝑁(0)=𝑆₀(0)(𝜃)[−>0,𝐼=𝜑≥0,∀𝜃∈𝜏,0],𝑅(0)=𝑅₀≥0. (4)

In this paper, we will consider two different delays 𝜏, 𝜔 which are important parameters on the dynamic behavior. So, the present study is continuation of the previous work 𝜏=𝜔 by Naresh et al. [43].

Lemma 1. All solutions of the model system (3) starting in 𝑅³₊ are bounded and eventually enter the compact attracting set Φ=(𝑆,𝐼,𝑅)∈𝑅³₊∶𝑆(𝑡)+𝐼(𝑡)+𝑅(𝑡)=𝑁(𝑡)≤𝐾. (5)

Lemma 2. Let the initial data be 𝑁(0)=𝑆₀>0, 𝐼(0)=𝐼₀(𝑢)≥0, for all 𝑢∈[−𝜏,0], with 𝐼₀(0)>0, 𝑅(0)=𝑅₀≥0. Then, the solution (𝑆(𝑡),𝐼(𝑡),𝑅(𝑡)) of the model remains positive for all time 𝑡>0.

Lemma 3 (see [44]). For the characteristic equation in the form 𝑝(𝜆)+𝑞(𝜆)𝑒^−𝑟𝜆=0, where 𝑝 and 𝑞 are polynomials with real coefficients and 𝑟>0 is the delay, suppose(a)𝑝(𝜆)≠0, 𝑅(𝜆)>0;(b)|𝑞(𝑖𝑦)|<|𝑝(𝑖𝑦)|; 0≤𝑦<∞;(c)lim_{|𝜆|→∞,𝑅(𝜆)≥0}|𝑞(𝜆)/𝑝(𝜆)|=0.
Then 𝑅(𝜆)<0 for every root 𝜆 and all 𝑟>0.

3. Equilibrium and Stability Analysis

In this section, we focus on the existence and local stability of equilibria. Let the right-hand side of equalities in model (3) be zero. Then, there are two equilibria; namely,(i)𝐸₀=(𝐾,0,𝑝), 𝑝=𝐾𝜃/[𝑑+(1−𝑎)𝑟+𝜃]), disease-free equilibrium;(ii)𝐸^∗=(𝑁^∗,𝐼^∗,𝑅^∗), endemic equilibrium,

where the values of 𝑁^∗, 𝐼^∗, and 𝑅^∗ are given in Section 3.2.

3.1. Community Matrix

Firstly, after computing the Jacobian or community matrix of model (3) at point (𝑁,𝐼,𝑅), the characteristic equation is given by𝑟−2𝑟𝐾𝑁−𝜆00𝛽𝑒^−𝑑1𝜏𝑓(𝐼)−(1−𝑎)𝑟𝐾𝐼𝑚𝑒^{−(𝑑1+𝜆)𝜏}−𝛿𝑒^{−(𝑑1+𝜆)𝜔}−𝑛−𝜆𝛽𝑒^−𝑑1𝜏𝑓(𝐼)()𝑟𝜃−1−𝑎𝐾𝑅𝛿𝑒^{−(𝑑1+𝜆)𝜔}()𝑟−𝜃−𝑑+1−𝑎𝐾𝑁+𝜃−𝜆=0, (6)where 𝑚=𝛽(𝑁−𝐼−𝑅)𝑓^′(𝐼), 𝑛=𝛽𝑒^−𝑑1𝜏𝑓(𝐼)+𝑑+((1−𝑎)𝑟/𝐾)𝑁+𝜇₁.

Now, we analyze the equilibria stability of system (3). Computing the Jacobian of system (3) evaluated at 𝐸₀, one gets the following matrix𝐽𝐸₀=−𝑟−𝜆000𝑎₂₂()𝜃−𝜆0𝑑+𝜃𝑑+(1−𝑎)𝑟+𝜃𝛿𝑒^{−(𝑑1+𝜆)𝜏}[()],−𝜃−𝑑+1−𝑎𝑟+𝜃−𝜆 (7)where 𝑎₂₂=𝛽𝐾−𝑝𝑓^′(0)𝑒^{−(𝑑1+𝜆)𝜏}−𝛿𝑒^{−(𝑑1+𝜆)𝜔}−𝑑+(1−𝑎)𝑟+𝜇₁. (8) Denote 𝐴=𝛽𝐾−𝑝𝑓^′(0),𝐶=𝑑+(1−𝑎)𝑟+𝜇₁; (9) then 𝑎₂₂=𝐴𝑒^{−(𝑑1+𝜆)𝜏}−𝛿𝑒^{−(𝑑1+𝜆)𝜔}−𝐶. (10) Denote ℎ(𝜆)=𝐴𝑒^{−(𝑑1+𝜆)𝜏}−𝛿𝑒^{−(𝑑1+𝜆)𝜔}−𝐶−𝜆. (11)

The eigenvalues of the system (3) about the steady state 𝐸₀ are 𝜆₁=−𝑟, ℎ(𝜆)=0 and 𝜆₃=−[𝑑+(1−𝑎)𝑟+𝜃]. All the parameters of the model are assumed to be nonnegative. Therefore, 𝜆₁ and 𝜆₃ are negative. Next, we discuss the roots of ℎ(𝜆)=0 in five cases.

Case 1. For 𝜏=𝜔≠0, from the second equation of the system (3), we can get the following.

Proposition 4. For 𝜏=𝜔>0, 𝑅(𝜆)<0 for every root 𝜆 of ℎ(𝜆)=0 when (𝐴−𝛿)𝑒^−𝑑1𝜏<𝐶. (12)

Proof. From the above analysis, 𝜆₂ satisfies the following characteristic equation: 𝑔(𝜆)=(𝐴−𝛿)𝑒^{−(𝑑1+𝜆)𝜏}−𝐶−𝜆=0. (13)(1)Clearly, 𝜆=0 is not a root of 𝑔(𝜆)=0.(2)From the fact that 𝑔(0)<0, 𝑔^′(𝜆)<0 for 𝜆>0, it is obtained that 𝑔(𝜆)=0 has no positive real root.(3)It is sufficient to show that 𝑔(𝜆)=0 does not admit a purely imaginary root. In fact, if 𝜆=𝑖𝑣(𝑣>0) is a root of (𝑔(𝜆)=0), then by separating the real part, one gets (𝐴−𝛿)𝑒^−𝑑1𝜏cos(𝑣𝜏)=𝐶. (14)Together with the condition of Proposition 4, we have cos(𝑣𝜏)>1. (15) This is impossible.(4)It is easy to show that 𝑔(𝜆)=0 has no imaginary root whose real part is positive. Otherwise, there is an imaginary root 𝜆=𝑢+𝑖𝑣 with 𝑢>0. Without any loss of generality, we consider 𝑣>0. Then, we take the real and imaginary parts of 𝑔(𝜆)=0; namely, (𝐴−𝛿)𝑒^{−(𝑑1+𝑢)𝜏}cos(𝑣𝜏)=𝐶+𝑢. (16)Combined with (9), we have 𝐶>(𝐴−𝛿)𝑒^−𝑑1𝜏>(𝐴−𝛿)𝑒^{−(𝑑1+𝑢)𝜏}cos(𝑣𝜏)=𝐶+𝑢. (17) This is a contradiction which implies that all eigenvalues roots of 𝑔(𝜆) have negative real parts. Therefore, the disease-free equilibrium of the system (3) is locally asymptotically stable when (9) holds. The proof is completed.

Case 2. For 𝜏≠0,𝜔=0, by the same way as in Case 1, one gets the following.

Proposition 5. For all 𝜏≠0,𝜔=0, 𝑅(𝜆)<0 for every root 𝜆 of ℎ(𝜆)=0 when 𝐴𝑒^−𝑑1𝜏−𝛿<𝐶. (18)

Case 3. For 𝜔≠0,𝜏=0, one gets the following.

Proposition 6. For 𝜔≠0,𝜏=0, 𝑅(𝜆)<0 for every root 𝜆 of ℎ(𝜆)=0 when 𝛿𝑒^−𝑑1𝜔<|𝐴−𝐶|. (19)

Proof. By the fact that ℎ(𝜆)=0 is equivalent to 𝑝(𝜆)+𝑞(𝜆)𝑒^−𝜆𝜔=0 with 𝑝(𝜆)=𝐴−𝐶−𝜆, 𝑞(𝜆)=−𝛿𝑒^−𝑑1𝜔.(i)Suppose 𝜆=𝑢+𝑖𝑣(𝑢>0). Then, 𝑝(𝜆)=𝐴−𝐶−𝑢−𝑖𝑣≠0.(ii)By |𝑞(𝑖𝑣)|=𝛿𝑒^−𝑑1𝜔, |𝑝(𝑖𝑣)|=|𝐴−𝐶−𝑖𝑣|=(𝐴−𝑐)²+𝑣², together with the condition of Proposition 6, we know that |𝑞(𝑖𝑣)|<|𝑝(𝑖𝑣)|.(iii)Suppose 𝜆=𝑢+𝑖𝑣,(𝑢>0). Then, lim_{𝑢2+𝑣2→+∞}𝑞(𝜆)𝑝(𝜆)=𝛿𝑒^−𝑑1𝜔,lim_{𝑢2+𝑣2→+∞}1(𝐴−𝐶−𝑢)²+𝑣²=0. (20)Then, using Lemma 3, we have 𝑅(𝜆)<0 for all 𝜔.

Case 4. For 𝜔>𝜏>0, let 𝜀=𝜔−𝜏. Then, 𝜔=𝜏+𝜀. For fixed 𝜏, ℎ(𝜆)=𝐴𝑒^{−(𝑑1+𝜆)𝜏}−𝛿𝑒^{−(𝑑1+𝜆)(𝜏+𝜀)}−𝐶−𝜆. (21) Let 𝜆=𝑢+𝑖𝑣(𝑢>0). Then, we take the real and imaginary parts of ℎ(𝜆)=0; namely, 𝐴𝑒^{−(𝑑1+𝑢)𝜏}cos(𝑣𝜏)−𝛿𝑒^{−(𝑑1+𝑢)(𝜏+𝜀)}cos(𝑣(𝜏+𝜀))=𝐶+𝑢,−𝐴𝑒^{−(𝑑1+𝑢)𝜏}sin(𝑣𝜏)−𝛿𝑒^{−(𝑑1+𝑢)(𝜏+𝜀)}sin(𝑣(𝜏+𝜀))=𝑣. (22) Sum of squares of the above equalities is 𝐴²𝑒^{−2(𝑑1+𝑢)𝜏}+𝛿²𝑒^{−2(𝑑1+𝑢)(𝜏+𝜀)}−2𝐴𝛿𝑒^{−(𝑑1+𝑢)(2𝜏+𝜀)}cos(𝑣𝜀)−(𝐶+𝑢)²−𝑣²=0. (23) Then, we have 𝜕𝑢𝜕𝜀_𝜀=0=𝛿𝑢+𝑑₁(𝐴−𝛿)(𝐴−𝛿)²𝜏+(𝐶+𝑢)𝑒^{2𝜏(𝑢+𝛿)}. (24) Obviously, 𝜕𝑢/𝜕𝜀|_𝜀=0<0 when 𝐴−𝛿<0. Combined with (𝐴−𝛿)𝑒^−𝑑1𝜏<𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔, we have the following.

Proposition 7. For 𝜔>𝜏>0, 𝑅(𝜆)<0 for every root 𝜆 of ℎ(𝜆)=0 when 𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔<𝐶,𝐴−𝛿<0,0<𝜔−𝜏≪1. (25)

Case 5. For 0<𝜏<𝜔, in the same way as in Case 4, we have the following.

Proposition 8. For 0<𝜏<𝜔, 𝑅(𝜆)<0 for every root 𝜆 of ℎ(𝜆)=0 when 𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔<𝐶,𝐴−𝛿>0,0<𝜏−𝜔≪1. (26)

From what has been discussed above, we get the following.

Theorem 9. The disease-free equilibrium of the system (3) is locally asymptotically stable if one of the following conditions holds.(a)𝜏=𝜔≠0, (𝐴−𝛿)𝑒^−𝑑1𝜏<𝐶.(b)𝜏≠0, 𝜔=0, 𝐴𝑒^−𝑑1𝜏−𝛿<𝐶.(c)𝜏=0, 𝜔≠0, 𝛿𝑒^−𝑑1𝜔<|𝐴−𝐶|.(d)𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔<𝐶, 𝐴−𝛿<0, 0<𝜔−𝜏≪1.(e)𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔<𝐶, 𝐴−𝛿>0, 0<𝜏−𝜔≪1.

3.2. Existence of Endemic Equilibrium

Thus, by Theorem 9, we may define the basic reproduction number as 𝑅₀=𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔𝐶. (27)

This threshold 𝑅₀ defines the average number of secondary infections generated by a typical infectious individual in a completely susceptible population in a steady demographic state.

In Theorem 9, we have already shown that the system (3) has an infection-free steady state which is locally asymptotically stable under condition 𝑅₀<1. The disease-free equilibrium is unstable when 𝑅₀>1, and the system (3) has a nontrivial endemic equilibrium 𝐸^∗=(𝑁^∗,𝐼^∗,𝑅^∗) when 𝑅₀>1. From (3), 𝑁^∗𝑅=𝐾>0,_∗=𝛿𝑒^−𝑑1𝜔−𝜃()𝐼𝑑+1−𝑎𝑟+𝜃^∗+𝐾𝜃𝑑+(1−𝑎)𝑟+𝜃≐𝑞𝐼^∗+𝑝, (28) where 𝑞=(𝑒^−𝑑1𝜔𝛿−𝜃)/(𝑑+(1−𝑎)𝑟+𝜃). Substituting these values of 𝑁^∗ and 𝑅^∗ in the second equation of (3), we get the following equation for 𝐼: 𝐺(𝐼)=𝛽𝑒^−𝑑1𝜏𝐾−1+𝑞𝐼−𝑝𝑓(𝐼)−𝐶+𝛿𝑒^−𝑑1𝜔𝐼. (29) Obviously, 𝐼=0 is one of the roots of (29) as 𝑓(0)=0. Therefore, to exclude that root, choose 𝐻(𝐼)=𝛽𝑒^−𝑑1𝜏𝑓(𝐼)𝐾−1+𝑞𝐼−𝑝𝐼−𝐶+𝛿𝑒^−𝑑1𝜔. (30)

It can easily be seen that the function 𝐻(𝐼) is negative for large positive 𝐼; that is, 𝐻(𝐾)=−𝛽𝑒^−𝑑1𝜏𝑓(𝑘)𝐾𝑞+𝑝𝐾−𝐶+𝛿𝑒^−𝑑1𝜔<0. (31)

Next, we determine the sign of its derivative 𝐻^′(𝐼)=𝛽𝑒^−𝑑1𝜏𝑓𝐾−𝑝^′(𝐼)𝐼−𝑓(𝐼)𝐼²−𝛽𝑒^−𝑑1𝜏1+𝑞𝑓^′(𝐼). (32) It can easily be seen that 𝐾>𝑝. In addition, from the properties of the function 𝑓(𝐼), in particular from 𝑓(0)=0 and 𝑓^′′(0)<0, it follows that 𝑓(𝐼)−𝑓^′(𝐼)𝐼>0, and consequently 𝐻^′(𝐼)<0 for all 𝐼>0. Therefore, for a positive root of 𝐻(𝐼)=0 to exist, 𝐻(𝐼) has to satisfy 𝐻(0)>0; that is, 𝐻(0)=𝐴𝑒^−𝑑1𝜏−𝛿𝑒^−𝑑1𝜔=−𝐶𝐴𝑒^−𝑑1𝜏−𝑒^−𝑑1𝜔𝛿𝐶𝐶=𝑅−1₀−1𝐶. (33) Hence, one needs the requirement that 𝑅₀>1 to ensure the existence of the endemic equilibrium. From the above analysis, we have the following theorem.

Theorem 10. The system (3) has a nontrivial endemic equilibrium 𝐸^∗=(𝑁^∗,𝐼^∗,𝑅^∗) when 𝑅₀>1.

3.3. Local Stability of the Endemic Equilibrium

In this section, we analyze the local stability of the endemic equilibrium 𝐸^∗ for 𝜏=𝜔. Its characteristic equation is given by−𝑟−𝜆00𝛽𝑒^−𝑑1𝜏𝑓𝐼^∗−(1−𝑎)𝑟𝐼^∗𝐾𝑚𝑒^{−(𝑑1+𝜆)𝜏}−𝑛−𝜆𝛽𝑒^−𝑑1𝜏𝑓𝐼^∗()𝜃−1−𝑎𝑟𝑅^∗𝐾𝛿𝑒^{−(𝑑1+𝜆)𝜏}−𝜃−[𝑑+(1−𝑎)𝑟+𝜃]−𝜆=0, (34)where 𝑚=𝛽(𝐾−𝐼^∗−𝑅^∗)𝑓^′(𝐼^∗)−𝛿, 𝑛=𝛽𝑒^−𝑑1𝜏𝑓(𝐼^∗)+𝑑+(1−𝑎)𝑟+𝜇₁.

The Jacobin matrix leads to the characteristic equation (𝜆+𝑟)𝜆²+𝑚₁𝜆+𝑚₀+𝑛₁𝜆+𝑛₀𝑒^−𝜆𝜏=0, (35) where 𝑚₁=𝛽𝑒^−𝑑1𝜏𝑓𝐼^∗+2𝑑+2(1−𝑎)𝑟+𝜃+𝜇₁𝑚>0,₀=𝛽𝑒^−𝑑1𝜏𝑓𝐼^∗+𝑑+(1−𝑎)𝑟+𝜇₁×(𝑑+(1−𝑎)𝑟+𝜃)+𝑑+(1−𝑎)𝑟+𝜇₁𝑛𝜃>0,₁=−𝛽𝑓^′𝐼^∗𝐾−𝐼^∗−𝑅^∗−𝛿𝑒^−𝑑1𝜏,𝑛₀=𝛽𝑓^′𝐼^∗𝐾−𝐼^∗−𝑅^∗−𝛿[𝑑+(1−𝑎)𝑟+𝜃]×𝑒^−𝑑1𝜏+𝛿𝑓𝐼^∗𝑒^{−2𝑑1𝜏}. (36) Since all the model parameters are assumed to be nonnegative, it follows that one eigenvalue is negative; that is, 𝜆₁=−𝑟. Thus, the stability of 𝐸^∗ depends on the roots of the quasi-polynomial 𝜆²+𝑚₁𝜆+𝑚₀+𝑛₁𝜆+𝑛₀𝑒^−𝜆𝜏=0. (37) We note that 𝑚₁>0 and 𝑚₀>0, whereas 𝑛₁ and 𝑛₀ may be positive or negative. For 𝜏=0, we state the following results that follow directly from (39). The endemic steady state is locally asymptotically stable if the following conditions hold: 𝛽𝑓𝐼^∗+2𝑑+2(1−𝑎)𝑟+𝜃+𝜇₁+𝛿𝛽𝑓^′𝐼^∗>𝛽𝐾−𝐼^∗−𝑅^∗𝑓^′𝐼^∗,𝛽𝑓𝐼^∗+𝑑+(1−𝑎)𝑟+𝜇₁+𝛿−𝛽𝐾−𝐼^∗−𝑅^∗𝑓^′𝐼^∗×[𝑑+(1−𝑎)𝑟+𝜃]+𝑓𝐼^∗𝛿>0. 𝐻₁𝐻₂ The main purpose of this paper is to study the stability behavior of 𝐸^∗ in the case 𝜏≠0. Obviously, 𝑖𝜂(𝜂>0) is the root of (29) if and only if 𝜂 satisfies −𝜂²+𝑚₁𝑖𝜂+𝑚₀=−𝑛₁𝑖𝜂+𝑛₀cos𝜂𝜏−𝑖sin𝜂𝜏. (38) Separating the real and imaginary parts, we have −𝜂²+𝑚₀=−𝑛₀cos𝜂𝜏−𝑛₁𝜂sin𝜂𝜏,𝑚₁𝜂=−𝑛₁𝜂cos𝜂𝜏+𝑛₀sin𝜂𝜏. (39)(40) Eliminating 𝜏 by squaring and adding (39) and (40), we obtain a polynomial in 𝜂 as 𝜂⁴+𝑚²₁−𝑛²₁−2𝑚₀𝜂²+𝑚²₀−𝑛²₀=0. (41) Suppose that the conditions 𝑚²₁>𝑛²₁+2𝑚₀,𝑚²₀>𝑛²₀ 𝐻₃ hold for all 𝜏≥0. Then, the infected steady state of the system (3) is locally asymptotically stable.

Theorem 11. For 𝜏=𝜔, if 𝑅₀>1, then the endemic equilibrium of the system (3) is locally asymptotically stable, when conditions 𝐻₁–𝐻₃ hold.

Corollary 12. For 𝜏=𝜔, if 𝜇₁<𝜇^∗₁ or 𝜃<𝜃^∗, then the endemic equilibrium of the system (3) is locally asymptotically stable, when conditions 𝐻₁–𝐻₃ hold.

4. Permanence

In this section, we investigate a permanence result [5]. The following is our main result of this paper. We will give the following result by using some techniques given in [8, 11]. The proof of the permanence with nonlinear incidence is a daunting task. Consequently, for simplicity and mathematical convenience, let us choose a linear incidence rate 𝑓(𝐼)=𝐼. The result holds with the nonlinear incidence, as shown numerically, but the algebraic proof is long and tedious, and the conditions to impose on some of the parameters may be very restrictive. Now, let us firstly give the following theorem.

Theorem 13. If 𝑅₀>1 holds, then the system (1) with 𝜏=𝜔 is permanent; that is, there are positive constants 𝑐_𝑖(𝑖=1,2,3) such that 𝑐₁<lim_𝑡→∞inf𝑆(𝑡)≤lim_𝑡→∞(𝑡)𝑐sup𝑆≤𝐾,₂<lim_𝑡→∞inf𝐼(𝑡)≤lim_𝑡→∞(𝑡)𝑐sup𝐼≤𝐾,₃<lim_𝑡→∞inf𝑅(𝑡)≤lim_𝑡→∞sup𝑅(𝑡)≤𝐾 (42) hold for any solution of (1) with (𝜙₁(𝜃),𝜙₂(𝜃),𝜙₃(𝜃)) in the interior of Φ for all 𝜃∈[−𝜏,0]. In fact, 𝑐_𝑖(𝑖=1,2,3) can be chosen explicitly as 𝑐₁=(𝑏−𝑎𝑟)𝐾𝛽𝐾𝑒^−𝑑1𝜏(),𝑐+𝑑+1−𝑎𝑟+𝜃₂=𝐼^∗𝑒^{−(𝑑+(1−𝑎)𝑟+𝜇1+𝑒1𝜏−𝑑𝛿)},𝑐₃=𝛿𝑒^−𝑑1𝜏𝑐₂+𝜃𝑐₁()𝑟.𝑑+1−𝑎 (43)

Proof. Note that 0<𝑁(𝑡)<𝐾 for all 𝑡≥0 and that lim_𝑡→∞𝑁(𝑡)=𝐾. It is easy to see that lim_𝑡→∞inf𝑆(𝑡)≥𝑐₁. In fact, let 𝜖<𝐾 be arbitrary. Choose 𝑇₁>𝜏 so large that 𝑁(𝑡)>𝐾−𝜖 for 𝑡>𝑇₁. We have the following inequality: 𝑆̇(𝑡)>−𝛽𝐾𝑒^−𝑑1𝜏()𝑆(𝑡)+(),+𝑑+1−𝑎𝑟+𝜃𝑏−𝑎𝑟)(𝐾−𝜖 (44) for all 𝑡≥𝑇₁, which implies that liminf_𝑡→∞𝑆(𝑡)≥(𝑏−𝑎𝑟)(𝐾−𝜖)𝛽𝐾𝑒^−𝑑1𝜏().+𝑑+1−𝑎𝑟+𝜃 (45) Note that 𝜖 may be arbitrarily small so that lim_𝑡→∞inf𝑆(𝑡)≥𝑐₁.
Next, we will show lim_𝑡→∞inf𝐼(𝑡)≥𝑐₂. For any 𝜉∶0<𝜉<1, we see the inequality 𝑆^∗<[(𝑏−𝑎𝑟)𝐾+𝜇₁𝐼^∗]/(𝛽𝑒^−𝑑1𝜏𝜉𝐼^∗+𝑑+(1−𝑎)𝑟+𝜃). There exist sufficiently large 𝜌≥1 and sufficiently small 𝜖 such that 𝑆^∗<{[(𝑏−𝑎𝑟)(𝐾−𝜖)+𝜇₁𝐼^∗]/(𝛽𝑒^−𝑑1𝜏𝜉𝐼^∗+𝑑+(1−𝑎)𝑟+𝜃)}(1−𝑒^{−(𝛽𝑒1𝜏−𝑑𝜉𝐼∗+𝑑+(1−𝑎)𝑟+𝜃)𝜌𝜏})≡𝑆^Δ. We show that 𝐼(𝑡₀)>𝑞𝐼^∗ for some 𝑡₀≥𝜌𝜏. In fact, if not, it follows from the first equation of (1) that, for all 𝑡≥𝜌𝜏+𝜏≥𝑇₁+𝜏, 𝑆̇(𝑡)≥(𝑏−𝑎𝑟)(𝐾−𝜖)+𝜇₁𝐼^∗−𝛽𝑒^−𝑑1𝜏𝜉𝐼^∗+𝑑+(1−𝑎)𝑟+𝜃𝑆(𝑡). (46) Hence, for 𝑡≥𝜌𝜏+𝜏, 𝑆(𝑡)≥𝑒^{−(𝛽𝑒1𝜏−𝑑𝜉𝐼∗+𝑑+(1−𝑎)𝑟+𝜃)(𝑡−𝜌𝜏−𝜏)}×𝑆𝜌𝜏+𝜏+(𝑏−𝑎𝑟)(𝐾−𝜖)+𝜇₁𝐼^∗×^𝑡_{𝜌𝜏+𝜏}𝑒^{−(𝛽𝑒1𝜏−𝑑𝜉𝐼∗+𝑑+(1−𝑎)𝑟+𝜃)(𝑡−𝜌𝜏−𝜏)}>)𝑑𝜃(𝑏−𝑎𝑟)(𝐾−𝜖+𝜇₁𝐼^∗𝛽𝑒^−𝑑1𝜏𝜉𝐼^∗()×+𝑑+1−𝑎𝑟+𝜃1−𝑒^{−(𝛽𝑒1𝜏−𝑑𝜉𝐼∗+𝑑+(1−𝑎)𝑟+𝜃)(𝑡−𝜌𝜏−𝜏)}, (47) which gives us, for 𝑡≥2𝜌𝜏+𝜏, 𝑆(𝑡)>𝑆^Δ>𝑆^∗. (48) For 𝑡≥0, we define a positive differentiable function 𝑉(𝑡) as follows: 𝑉(𝑡)(𝑡)+𝛽𝑒=𝐼𝐾−𝑝−𝛿^−𝑑1𝜏𝑅₀^𝑡_𝑡−𝜏𝐼(𝑠)𝑑𝑠. (49) We obtain the inequality, for 𝑡≥2𝜌𝜏+𝜏, 𝑉̇(𝑡)=𝛽𝑒^−𝑑1𝜏𝑆(𝑡)−𝑆^∗𝐼()+()𝑟×𝑁(𝑡)𝑡−𝜏1−𝑎1−𝐾𝐼(𝑡)>𝛽𝑒^−𝑑1𝜏𝑆(𝑡)−𝑆^∗𝐼(𝑡−𝜏)>𝛽𝑒^−𝑑1𝜏𝑆^Δ−𝑆^∗𝐼(𝑡−𝜏). (50)
Let 𝑖=min_{𝜃∈[−𝜏,0]}𝐼(2𝜌𝜏+2𝜏+𝜃). Now, let us show that 𝐼(𝑡)≥𝑖 for all 𝑡≥2𝜌𝜏+𝜏. In fact, if there exists 𝑇₂≥0 such that 𝐼(𝑡)≥𝑖 for 2𝜌𝜏+𝜏≤𝑡≤2𝜌𝜏+2𝜏+𝑇₂, 𝐼(2𝜌𝜏+2𝜏+𝑇₂)=𝑖 and 𝐼(̇2𝜌𝜏+2𝜏+𝑇₂)≤0. Direct calculation using the second equation of (1) and (48) gives 𝐼̇2𝜌𝜏+2𝜏+𝑇₂>𝛽𝑒^−𝑑1𝜏𝑆2𝜌𝜏+2𝜏+𝑇₂−𝑑+(1−𝑎)𝑟+𝜇₁+𝑒^−𝑑1𝜏𝛿𝑖>𝑑+(1−𝑎)𝑟+𝜇₁+𝑒^−𝑑1𝜏𝛿𝑆^Δ𝑆^∗−1𝑖>0. (51) This contradicts the definition of 𝑇₂. Thus, we have shown that 𝐼(𝑡)≥𝑖 for all 𝑡≥2𝜌𝜏+𝜏. Hence, for all 𝑡≥2𝜌𝜏+2𝜏, 𝑉̇(𝑡)>𝛽𝑒^−𝑑1𝜏𝑆^Δ−𝑆^∗𝑖, (52) which implies that 𝑉(𝑡)→+∞ as 𝑡→+∞. This contradicts the boundedness of 𝑉(𝑡). Consequently, 𝐼(𝑡₀)>𝜉𝐼^∗ for some 𝑡₀≥𝜌𝜏.
In the rest, we now need to consider two cases:(i)𝐼(𝑡)≥𝜉𝐼^∗ for all large 𝑡;(ii)𝐼(𝑡) oscillates about 𝜉𝐼^∗ for all large 𝑡.
We now need to show that 𝐼(𝑡)≥𝜉𝑐₂ for large 𝑡. Obviously, it suffices to show that it holds only for case (ii). We suppose that for any large 𝑇 there exists 𝑡₁, 𝑡₂>𝑇 such that 𝐼(𝑡₁)=𝐼(𝑡₂)=𝜉𝐼^∗ and 𝐼(𝑡)<𝜉𝐼^∗ for 𝑡₁<𝑡<𝑡₂. If 𝑡₂−𝑡₁≤𝜏, the second equation (1) gives us 𝐼(̇𝑡)>−(𝑑+(1−𝑎)𝑟+𝜇₁)𝐼(𝑡), which implies that 𝐼(𝑡)>𝐼(𝑡₁)𝑒^{−(𝑑+(1−𝑎)𝑟+𝜇1)(𝑡−𝑡1)} on (𝑡₁,𝑡₂). Thus, 𝐼(𝑡)>𝜉𝑐₂. On the other hand, if 𝑡₂−𝑡₁>𝜏, applying the same manner gives 𝐼(𝑡)≥𝜉𝑐₂ on [𝑡₁,𝑡₁+𝜏], and hence the remaining work is to show 𝐼(𝑡)≥𝜉𝑐₂ on [𝑡₁+𝜏,𝑡₂]. In fact, assuming that there exists 𝑇₃>0 such that 𝐼(𝑡)≥𝜉𝑐₂ on [𝑡₁,𝑡₁+𝜏+𝑇₃], 𝐼(𝑡₁+𝜏+𝑇₃)=𝜉𝑐₂, and 𝐼(̇𝑡₁+𝜏+𝑇₃)≤0, it follows from (1) that 𝐼̇𝑡₁+𝜏+𝑇₃≥𝛽𝑒^−𝑑1𝜏𝑆𝑡₁+𝜏+𝑇₃−𝑑+(1−𝑎)𝑟+𝜇₁+𝑒^−𝑑1𝜏𝛿𝜉𝑐₂>𝑑+(1−𝑎)𝑟+𝜇₁+𝑒^−𝑑1𝜏𝛿𝑆^Δ𝑆^∗−1𝜉𝑐₂>0. (53) This contradicts the definition of 𝑇₃. Hence, 𝐼(𝑡)≥𝜉𝑐₂ on [𝑡₁,𝑡₂]. Consequently, 𝐼(𝑡)≥𝜉𝑐₂ for large 𝑡 in the case (ii). Therefore, lim_𝑡→∞inf𝐼(𝑡)≥𝜉𝑐₂. Note that 𝑞 may be so close to 1 that lim_𝑡→∞inf𝐼(𝑡)≥𝑐₂.
Finally, let us show that lim_𝑡→∞inf𝑅(𝑡)≥(𝛿𝑒^−𝑑1𝜏𝑐₂+𝜃𝑐₁)/(𝑑+(1−𝑎)𝑟). The third equation gives us 𝑅̇(𝑡)≥𝛿𝑒^−𝑑1𝜏[()𝑟]𝑅≥𝐼+𝜃𝑆−𝑑+1−𝑎𝛿𝑒^−𝑑1𝜏𝜉𝑐₂+𝜃𝜉𝑐₁−[𝑑+(1−𝑎)𝑟]𝑅 (54) for large 𝑡. Hence, lim_𝑡→∞inf𝑅(𝑡)≥(𝛿𝑒^−𝑑1𝜏𝜉𝑐₂+𝜃𝜉𝑐₁)/(𝑑+(1−𝑎)𝑟). In a similar manner, we could show lim_𝑡→∞inf𝑅(𝑡)≥𝑐₃. This proves the theorem.

Corollary 14. If 𝜇₁<𝜇^∗₁ and 𝜃<𝜃^∗, then the system (1) with 𝜏=𝜔 is permanent.

5. Numerical Analysis

Since it is important to visualize the dynamical behavior of the model, the model system (3) is integrated numerically with the help of MATLAB 7.0 using the following set of parameters.

(1) Let 𝑟=0.5, 𝑘=8, 𝑑=0.04, 𝑑₁=0.04, 𝛽=1, 𝑎=0.3, 𝛿=0.2, 𝜇₁=0.8, 𝜃=0.01, and 𝜏=5. It is easy to compute that 𝐸₀=(8,0,0.2) and 𝑅₀=0.94<1. In Figure 1, the infective population and recovered population, respectively, are plotted against the total population. We see from the figure that for any initial start the solution curves tend to the equilibrium 𝐸₀. Hence, we infer that the system (3) may be stable about the disease-free equilibrium point 𝐸₀, which satisfies Theorem 9.

Figure 1 

The disease-free equilibrium 𝐸0 is locally asymptotically stable. Variation of infective population 𝐼(𝑡) and recovered population 𝑅(𝑡) with total population 𝑁(𝑡).

(2) Let 𝑟=0.5, 𝑘=8, 𝑑=0.04, 𝑑₁=0.04, 𝛽=1, 𝑎=0.3, 𝛿=0.2, 𝜇₁=0.2, 𝜃=0.02, and 𝜏=5. We get 𝐸^∗=(8,2.23,1.11) and this set of parameter values satisfies the local asymptotic stability conditions of 𝐸^∗. It is easy to verify that 𝑅₀=1.83>1 and all other conditions of Theorem 11 are satisfied. So, we can obtain from Figure 2 that the system (3) is stable at the endemic equilibrium point 𝐸^∗.

Figure 2 

The endemic equilibrium 𝐸∗ is locally asymptotically stable. Variation of infective population 𝐼(𝑡) and recovered population 𝑅(𝑡) with total population 𝑁(𝑡).

(3) The results of numerical simulation are displayed graphically in Figures 3 and 4. In Figure 3, the variation of the infective population and recovered population is shown with time for different values of the removal rate constant from groups 𝐼 to 𝑆, 𝜇₁. It is found that both the infective population and the recovered population decrease as 𝜇₁ increases. Figure 4 depicts the variation of infective population and recovered population, respectively, with time for the different successive vaccination rate, 𝜃. As 𝜃 increases, the infective population decreases whereas the recovered population increases.

Figure 3 

Variation of infective population 𝐼(𝑡) and recovered population 𝑅(𝑡) with time for different values of 𝜇1.

Figure 4 

Variation of infective population 𝐼(𝑡) and recovered population 𝑅(𝑡) with time for different values of 𝜃.

6. Discussion

In this paper, we will consider two different delays which are important parameters on the dynamic behavior. So, the present study is continuation of the previous work by [43]. Furthermore, from biological epidemic point of view, we investigate successive vaccination and difference in immunity in our system. From mathematical point of view, we study the stability of disease-free equilibrium and the existence of endemic equilibrium for different delay and consider the permanence of the system in the new paper.

In Theorems 9, 10, 11, and 13 corresponding to their corollaries, when the effect of the successive vaccination rate and the transfer rate from the infectious group to the susceptible group after treatment is strong, that is, 𝜃>𝜃^∗ and 𝜇₁>𝜇^∗₁, the basic reproduction number 𝑅₀ being unity is a strict threshold for the control of the disease; the disease will be extinct or otherwise will tend to break out and persist. The other results are displayed graphically from our numerical simulation. We show the variation of the infective population and recovered population with time for different values of 𝜇₁. It is found that both the infective population and the recovered population decrease as 𝜇₁ increases. The infective population decreases whereas the recovered population increases as the successive vaccination rate increases, 𝜃, respectively.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work is supported by National Natural Science Foundation of China (11101305).