International Journal of Analysis and Applications

Volume 17, Number 5 (2019), 809-820

URL: https://doi.org/10.28924/2291-8639

DOI: 10.28924/2291-8639-17-2019-809

BIFURCATION AND STABILITY ANALYSIS OF A DISCRETE TIME SIR EPIDEMIC

MODEL WITH VACCINATION

ÖZLEM AK GÜMÜŞ1,∗, A. GEORGE MARIA SELVAM2 AND D. ABRAHAM VIANNY2

1Adıyaman University, Faculty of Arts and Sciences, Department of Mathematics, 02040, Adiyaman,

Turkey

2Sacred Heart College (Autonomous), Department of Mathematics, Tirupattur-635 601, Vellore Dt., Tamil

Nadu, India

∗Corresponding author: akgumus@adiyaman.edu.tr

Abstract. In this paper, we study the qualitative behavior of a discrete-time epidemic model with vacci-

nation. Analysis of the model shows forth that the Disease Free Equilibrium (DFE) point is asymptotically

stable if the basic reproduction number R0 is less than one, while the Endemic Equilibrium (EE) point

is asymptotically stable if the basic reproduction number R0 is greater than one. The results are rein-

forced with numerical simulations and enhanced with graphical representations like time trajectories, phase

portraits and bifurcation diagrams for different sets of parameter values.

1. Introduction

Mathematical models defining biological events has an important place in the study of population dynam-

ics. Most of the biological occurrences in nature are illustrated by discrete time, which point to, that there

are particular time instants at which the basic events in the system can occur, and it is not essential that at

these discrete time instants only a exclusive event happens. The most realistic approach to non-overlapping

generations like fish or insect populations, is created with discrete time system ( [6], [9], [15], [16], [17], [18]).

Received 2019-05-26; accepted 2019-07-08; published 2019-09-02.

2010 Mathematics Subject Classification. 39A10, 39A28, 39A30.

Key words and phrases. difference equations; epidemic model; bifurcation; stability theory.

c©2019 Authors retain the copyrights

of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License.

809

https://doi.org/10.28924/2291-8639
https://doi.org/10.28924/2291-8639-17-2019-809


Int. J. Anal. Appl. 17 (5) (2019) 810

One of the other famous examples identified by these systems are epidemic models [10]. Modeling an

outbreak that is progressing through the population allows us to examine the consequences of ways of

preventing or controlling the disease. Although the work of Kermack and McKendrick have the basic

foundation for epidemics, the first attempt in explaining, predicting or modeling of epidemics dates back to

over a century is made by Hamer (1906), Ross (1911). These early models operate on the principle where

individuals can be classified by their epidemiological status which are susceptible to the infection, infected

and recovered (immune) ( [12], [13], [14]).

Discrete time models are more suitable than continuous time models to examine infectious diseases due to

many reasons. Statistical data on diseases are collected at a specific time. In this case, the appropriate model

defining the disease will be the discrete time model [19]. On the other hand, the studies on discrete time

models obtained from the continuous time model by using nonstandard discretization technique are more

suitable to avoid mathematical complexity with regularity of solutions ( [20], [21]). Furthermore, although

the continuous-time logistic equation has only equilibrium dynamics, the well known discrete logistic equation

which is discrete counterpart of it exhibits period doubling bifurcation to chaos ( [22], [23]).

2. SIR Epidemic Model with Vaccination

The mathematical modelling of infectious diseases has a significant role in the studies of dynamical system.

Because studies on the dynamics of these models help us to control diseases like swine flu, bird flue and AIDS.

SIR models are suitable to define the transmission of infectious diseases with lifelong immunity such as chicken

pox, measles, smallpox, mumps and SARS. The SIR model is one of the simplest and most fundamental

of all epidemiological models and in these models with a single epidemic, births and deaths are ignored,

and so, only two transitions are possible: infection (moving individuals from the susceptible to the infected

class) and recovery (moving individuals from the infected to the recovered class). The assumptions in this

model are that the per capita rate that a given susceptible individual becomes infected is proportional to

the prevalence of infection in the population and that infected individuals recovers at a constant rate. The

fundamental parameter that governs the behavior of the epidemic is the basic reproductive ratio, R0 which

is defined as the average number of secondary cases produced by a single infectious individual in a totally

susceptible population.

Vaccination can be included in a epidemic model by assuming a proportion of susceptible individuals

are vaccinated during each time interval ( [1], [3], [4], [5]). Vaccination operates by reducing the pool of

susceptible individuals, and when this is reduced sufficiently, an infectious disease cannot spread within the

population. Most importantly, it is not necessary to vaccinate everyone to prevent an epidemic, immunizing

someone not only protects that person but confers some protection to the population in general [11].



Int. J. Anal. Appl. 17 (5) (2019) 811

3. The Discrete Time System

The author of [2] has presented the dynamics of the SIR epidemic model which is as follows:

St+1 =St −
a

N
ItSt + β(Rt + It)

It+1 =
a

N
ItSt + (1 −β −γ)It

Rt+1 =(1 −β)Rt + γIt + pSt

(3.1)

where a > 0, 0 < β < 1 and 0 < γ < 1.

In this paper, we focus on the dynamics of a SIR epidemic model by including vaccination to the model

as presented in [2]. The general SIR epidemic model is of the following form [1]:

St+1 =(1 −p)St −
a

N
ItSt + β(Rt + It)

It+1 =
a

N
ItSt + (1 −β −γ)It

Rt+1 =(1 −β)Rt + γIt + pSt

such that the initial conditions S0,I0 and R0 which are positive real numbers with (S0 +I0 +R0 = N). Here

0 < p + a < 1 and 0 < β + γ < 1. Also, β is the probability of birth, γ is the probability of recovery, p is the

proportion of vaccinated, a is the contact rate and N is the total population size. Moreover, we have the

following equivalent two dimensional system using the relation St + It + Rt = N.

St+1 =(1 −p)St −
a

N
ItSt + β(N −St)

It+1 =
a

N
ItSt + (1 −β −γ)It

(3.2)

where p, a, β and γ have positive values.

4. Stability of Equilibrium points and Numerical Simulations

In this section, we consider the discrete-time system (3.2). Foremost, we discuss the existence of equi-

librium points for (3.2), and then study the stability of the equilibrium points by using the characteristic

polynomial or the eigenvalues of the Jacobian matrix evaluated at each of the fixed points.

Lemma 4.1. [7] Let Q(x) = x2 + Bx + C. Suppose that Q(1) > 0, x1 and x2 are two roots of Q(x) = 0.

Then

(i) |x1| < 1 and |x2| < 1 if and only if Q(−1) > 0 and C < 1;

(ii) |x1| < 1 and |x2| > 1 (or |x1| > 1 and |x2| < 1) if and only if Q(−1) < 0;

(iii) |x1| > 1 and |x2| > 1 if and only if Q(−1) > 0 and C > 1;

(iv) x1 = −1 and |x2| 6= 1 if and only if Q(−1) = 0 and B 6= 0, 2;

(v) x1 and x2 are complex and |x1| = |x2| = 1 if and only if B2 − 4C < 0 and C = 1.



Int. J. Anal. Appl. 17 (5) (2019) 812

Lemma 4.2. [8]The characteristic polynomial

Q(x) = x2 + Bx + C

where B=-(Trace of the Jacobian matrix) and C= Determinant of the Jacobian matrix has all its roots

inside the unit open disk (|x| < 1) if and only if (i) Q(1) > 0 and Q(−1) > 0. (ii) D+1 = 1 + C > 0 and

D−1 = 1 −C > 0

Now, we will investigate the equilibrium points and then analyze the stability of these equilibrium points.

For analyzing the local stability of equilibrium points (S∗,I∗), we give the following theorems.

Theorem 4.1. The model (3.2) has two equilibrium points, P0 =
(
βN
β+p

, 0
)

and P1 =
(

(β+γ)N
a

,
N(aβ−(β+γ)(p+β))

a(β+γ)

)
.

Proof. When we examine the following equilibrium points (S∗,I∗) of the model (3.2), we easily obtain the

equilibrium points of the model (3.2) by using St = St+1 = S
∗ and It = It+1 = I

∗ :

S∗ = (1 −p)S∗ −
a

N
I∗S∗ + β(N −S∗)

I∗ =
a

N
I∗S∗ + (1 −β −γ)I∗

�

Theorem 4.2. Suppose that p + β < 1. The disease-free equilibrium (DFE) point P0 =
(
βN
β+p

, 0
)

of the

system (3.2) is locally asymptotically stable (LAS) if

aβ

(β + p)(β + γ)
< 1. (4.1)

Proof. By considering (3.2), we can get the Jacobian matrix evaluated P0 as

JP0 =


 1 −p−β −aβ(β+p)

0 aβ
(β+p)

+ (1 −β −γ)


 .

The eigenvalues of this matrix are

x1 = 1 −p−β, x2 =
aβ

(p + β)
+ (1 −γ −β).

If β + p < 1, then it is easy to see that x1 = 1 −p−β < 1, and also since β + γ < 1, x2 > 0 is always true.

Consequently, if |x2| = aβ(β+p) + (1 −β −γ) < 1, then we get
aβ

(β+p)(β+γ)
< 1. �

Corollary 4.1. The basic reproductive ratio R0 is referred as
aβ

(β+p)(β+γ)
. This ratio is a threshold parameter.

If R0 < 1, then there exists that the DFE point is LAS.

We consider the initial conditions (S(0),I(0)) = (70, 30) for numerical study.



Int. J. Anal. Appl. 17 (5) (2019) 813

Example 4.1. (a) For the DFE point, we assume the parameter values as N = 100, p = 0.0005, a = 0.1,

β = 0.02, γ = 0.2. The eigenvalues are |x1| = 0.9795 < 1, |x2| = 0.8776 < 1 and R0 = 0.4435 < 1 then the

DFE point P0 = (97.561, 0) of the model (3.2) is LAS (see Figure-1). (b)We take the parameter values as

N = 100,p = 0.05,a = 1.7,β = 1.7,γ = 0.1. The eigenvalues are |x1| = 0.7500 < 1, |x2| = 0.8514 < 1 and

R0 = 0.9175 < 1 then the DFE point P0 = (97.1429, 0) of the model (3.2) is LAS (see Figure-2). Note that

trace JP0 > 0.

Figure 1. Time plots and phase portrait of DFE point P0 with stability R0 < 1

Figure 2. Time plots and phase portrait of DFE point P0 with stability R0 < 1

Theorem 4.3. If 1 < R0 <
2

(p+β)
, then the endemic equilibrium (EE) point P1 =

(
(β+γ)N

a
,
N(aβ−(β+γ)(p+β))

a(β+γ)

)
of the system (3.2) is LAS.



Int. J. Anal. Appl. 17 (5) (2019) 814

Proof. By considering (3.2), we can write the Jacobian matrix evaluated at P1 as

JP1 =


 1 − βaβ+γ −β −γ

βa
β+γ
− (p + β) 1


 , (4.2)

If it is organized as relating to R0, we find

JP1 =


 1 − (β + p)R0 −aβ(β+p)R0

(p + β)(R0 − 1) 1


 ,

The characteristic polynomial of the Jacobian matrix at JP1 is as follows:

Q(x) = x2 − (2 − (β + p)R0)x + 1 − (β + p)R0 + aβ
(

1 −
1

R0

)
. (4.3)

For the stability of the EE point of (3.2), we get

0 < aβ

(
1 −

1

R0

)
< R0(p + β),

or equivalently

1 < R0 <
βa

(β + γ)2(p + β)
+ 1. (4.4)

from Lemma 4.2. Note that

aβ

β + γ
< 2 (4.5)

is always provided. Equivalently, we have

R0 <
2

(p + β)
. (4.6)

If (4.4) and (4.6) are compared, then we get

2

(p + β)
<

βa

(β + γ)2(p + β)
+ 1. (4.7)

Thus the proof is completed. �

Corollary 4.2. If 1 < R0 <
2

(p+β)
, then Q(1) > 0, Q(−1) > 0 and C < 1 is always confirmed such that

0 < p + a < 1 and 0 < β + γ < 1.

Proof. From (4.2), the characteristic polynomial is as follows:

Q(x) = x2 +

(
aβ

β + γ
− 2
)
x + 1 −

aβ

β + γ
+ aβ − (β + γ)(β + p). (4.8)

Obviously, Q(1) > 0 is always true, since Q(1) > 0, R0 > 1. Also, we obtain

Q(−1) = 4 + aβ − (β + γ)(β + p) −
2aβ

β + γ
(4.9)

D−1 = 1 −C =
aβ

β + γ
−aβ + (β + γ)(β + p) (4.10)

and

D+1 = 1 + C = 2 −
aβ

β + γ
+ aβ − (β + γ)(β + p). (4.11)



Int. J. Anal. Appl. 17 (5) (2019) 815

Here, we take

C = Q(0) = 1 −
aβ

β + γ
+ aβ − (β + γ)(β + p) (4.12)

such that R0 > 1. Note that whenever Q(−1) > 0, D+1 > 0 is always true. By considering (4.10), we can

write,

aβ −aβ(β + γ) + (β + γ)2(β + p) > 0

such that β + γ < 1. It clear that aβ −aβ(β + γ) > 0. So 1 −C > 0 is always provided. Similarly, we can

write by considering (4.9).

(β + γ)[4 + aβ − (β + γ)(β + p)] − 2aβ > 0

such that aβ − (β + γ)(β + p) > 0. Then, we get Q(−1) > 0. From (4.5) and by the positive state of the EE

point of (3.2), the result is clear. �

Example 4.2. For the EE point, we take the parameter values as (a) N = 100, p = 0.0005, a = 0.6,

β = 0.025, γ = 0.3. Applying the conditions, we get Q(1) = 0.0068 > 0, Q(−1) = 3.9144 > 0, C = 0.9606 < 1

and R0 = 1.81 > 1 and thus the EE point P1 = (54.1667, 3.4423) of the model (3.2) is LAS (see Figure-3).

(b) We take the parameter values as N = 100, p = 0.05, a = 4.2, β = 1.6, γ = 0.1 and applying the

conditions we see that Q(1) = 3.9150 > 0, Q(−1) = 0.0092 > 0, C = 0.9621 < 1 and R0 = 2.3957 > 1 and

so the EE point P1 = (40.4762, 54.8319) of the model (3.2) is LAS (see Figure-4).

Figure 3. Time plots and phase portrait of EE point P1 with stability R0 > 1

5. Bifurcation

In this section, we give the bifurcation diagrams of the susceptible and infected populations of the model

(3.2). The bifurcation diagrams are considered for four cases:

Case (i): Fixing parameters N = 100, β = 0.8, p = 0.0005, γ = 0.1 and varying a.

The bifurcation diagrams of model (3.2) are plotted with contact rate a ∈ (3.0, 4.15) as the bifurcation



Int. J. Anal. Appl. 17 (5) (2019) 816

Figure 4. Time plots of EE point P1 with stability R0 > 1

parameter and the system undergoes periodic doubling or flip bifurcation. When a ∈ (3.0, 3.36) there appears

stability. In the range a ∈ (3.36, 3.8) periodic-2 orbits, for a ∈ (3.8, 3.9) periodic-4 orbits and a ∈ (3.9, 3.95)

periodic-8 orbits occur, leading to chaos for a ∈ (3.95, 4.15). Local amplifications corresponding to figure (5)

for a ∈ [3.75, 4.15] can be seen in figure(6).

Figure 5. Bifurcation diagrams for susceptible and infected populations with a ∈ (3.0, 4.15)

Case (ii): Fixing parameters N = 100, p = 0.0005, β = 0.8, a = 3.5 and varying γ. The bifurcation

diagrams of model (3.2) are plotted with recovery rate γ ∈ (0, 0.4), as the bifurcation parameter. When

γ ∈ (0, 0.023) there appears chaos. In the range γ ∈ (0.023, 0.03) periodic-8 orbits, for γ ∈ (0.03, 0.05)

periodic-4 orbits, for γ ∈ (0.05, 0.12) periodic-2 orbits occur which is called as periodic half bifurcation.

Finally for the range γ ∈ (0.12, 0.4) there appears stability (see Figure-7).

Case (iii): Fixing parameters N = 100, β = 0.8, a = 3.5, γ = 0.1 and varying p.

The bifurcation diagrams of model (3.2) are plotted in the particular range of p ∈ (0, 0.2), with proportion



Int. J. Anal. Appl. 17 (5) (2019) 817

Figure 6. Local amplification corresponding to figure (5) for a ∈ (3.75, 4.15)

Figure 7. Bifurcation diagrams for susceptible and infected populations with γ ∈ (0, 0.4)

vaccinated rate as the bifurcation parameter. When p ∈ (0, 0.02) there appears chaos. In the range p ∈

(0.02, 0.056) there appears stability, in the range p ∈ (0.056, 0.127) there appears periodic-2 orbits, in the

range p ∈ (0.127, 0.15) periodic-4 orbits, in the range p ∈ (0.15, 0.16) periodic-8 orbits and in the range

p ∈ (0.16, 0.185) there is chaos (see Figure-8).

Case (iv): Fixing parameters N = 100, p = 0.0005, a = 4.2, γ = 0.1 and varying β. The bifurcation

diagrams of the model (3.2) are plotted in the particular range of β ∈ (1.0, 2.5), with birth rate as the

bifurcation parameter. Local amplification corresponding to figure (9) for β ∈ (1.5, 2.5) can be seen in

Figure (10).



Int. J. Anal. Appl. 17 (5) (2019) 818

Figure 8. Bifurcation diagrams for susceptible and infected populations with p ∈ (0, 0.2)

Figure 9. Bifurcation diagrams for susceptible and infected population with β ∈ (1.0, 2.5)

6. Conclusion

In this paper, we consider an discrete time SIR epidemic model with vaccination and obtained the con-

ditions for the existence of the equilibrium points and discussed the stability of the system at DFE and

EE points. Also the numerical examples ascertain the theoretical findings. Time plots and phase portraits

are presented for the susceptible and infected populations for biological feasible parameters. Bifurcation

diagrams and local amplifications of the same are presented. The discrete model exhibits varied and rich

dynamical behavior.

Estimates on R0 have been obtained to determine the emergence of diseases such as measles, chickenpox

and smallpox [24]. We present the dynamics of the model with the effect of vaccine ( [1], [2]). In Example

4.1-(a) and in Example 4.2-(a), we observe that the diseases free equilibrium is local asymptotic stable since



Int. J. Anal. Appl. 17 (5) (2019) 819

Figure 10. Local amplification corresponding to figure (9) for β ∈ (1.5, 2.5)

R0 < 1 (see Figure-1) and the endemic equilibrium point is local asymptotic stable since R0 > 1 (see Figure-

3) by taking p = 0.0005 and N = 100. Example 4.1-(b) shows that there is a decrease in the number of

susceptible persons even if the vaccination rate increases when the rate of recovery decreases and the rate of

contact increases (see Figure-2). If the rate of contact increases further, Example 4.2-(b) demonstrates an

increase in the number of diseases (see Figure-4). Figure 5 points the bifurcation diagrams for susceptible and

infected populations with changing values of a. In Figure 6, we exhibit the local amplification corresponding

to Figure 5. Figure 7 shows the bifurcation diagrams for susceptible and infected populations with changing

values γ. Figure 8 displays the bifurcation diagrams for susceptible and infected populations with changing

values p. Lastly, for the particular range of β, local amplification corresponding to Figure 9 which shows

bifurcation diagrams are presented in Figure 10. Consequently, the lower contact rate of a has an effect of

reducing the disease ( [1], [24]). Also increasing of the rate of vaccination has a reinforcing effect on the

reduction of the disease as other parameters remain constant.

References

[1] L.J.S. Allen, An Introduction to Mathematical Biology, Pearson, New Jersey, 2007.

[2] Q. Din, Qualitative behavior of a discrete SIR epidemic model, Int. J. Biomath. 9 (6) (2016), 1650092.

[3] A. George Maria Selvam and D. Abraham Vianny, Behavior of a Discrete Fractional Order SIR Epidemic model, Int. J.

Eng. Technol. 7 (2018), 675-680.

[4] A. George Maria Selvam, D. Abraham Vianny and Mary Jacintha, Stability in a fractional order SIR epidemic model of

childhood diseases with discretization, IOP Conf. Ser., J. Phys., Conf. Ser., 1139 (2018), 012009.

[5] A.George Maria Selvam and D.Abraham Vianny, Discrete Fractional Order SIR Epidemic Model of Childhood Diseases

with Constant Vaccination and its Stability, Int. J. Techn. Innov. Modern Eng. Sci., 4 (11) (2018), 405-410.

[6] E.A. Grove and G. Ladas, Periodicities in Nonlinear Difference Equations, Chapman & Hall/CRC Press, Boca Raton,

2004.



Int. J. Anal. Appl. 17 (5) (2019) 820

[7] X. Liu, D. Xiao, Complex dynamic behaviors of a discrete-time predator-prey system. Chaos Solution Fract., 32 (2007),

80-94.

[8] X. Liu, C. Mou, W. Niu, D. Wang, Stability analysis for discrete biological models using algebraic methods. Math. Comput

Sci. 5 (2011), 247-262.

[9] R M. May, Simple mathematical models with very complicated dynamics, Nature, 261 (1976), 459-467.

[10] H. Sedaghat, Nonlinear Difference Equations: Theory with Applications to Social Science Models, Kluwer Academic

Publishers, Dordrecht, Netherlands, 2003.

[11] M. J. Keeling and L. Danon, Mathematical modelling of infectious diseases, Br. Med. Bull. 92 (1) (2009), 33-42.

[12] Kermack, W. O. McKendrick, A. G, A Contribution to the Mathematical Theory of Epidemics, Proc. Royal Soc. Ser A,

115 (772) (1927), 700-721.

[13] W. Hamer, II. Epidemic disease in England. Lancet, 1 (1906), 733-739.

[14] R. Ross, The Prevention of Malaria. 2nd ea. John Murray, London, 1911.

[15] O. Ak Gumus, Global and local stability analysis in a nonlinear discretetime population model, Adv. Difference Equ. 2014

(2014), 299.

[16] H Merdan, O. Ak Gumus, Stability analysis of a general discrete-time population model involving delay and Allee effects.

Appl. Math. Comput. 219 (2012), 1821-1832.

[17] Q Din, ÖA Gümüş and H Khalil, Neimark-Sacker Bifurcation and Chaotic Behavior of a Modified Host-Parasitoid Model,

Z. Naturforsch., A, 72 (1) (2017), 25-37.

[18] H. Merdan, Ö. Ak Gümüş and G. Karahisarli, Global Stability Analysis of a General Scalar Difference Equation, Discon-

tinuity, Nonlinear. Complex. 7 (3) (2018), 225-232.

[19] I. Longili, The generalized discrete-time epidemic model with immunity: a synthesis, Math. Biosci. 82(1986), 19-41.

[20] S. Jang, S. Elaydi, Difference equations from discretization of a continuous epidemic model with immigration of infectives,

Can. Appl. Math. Q. 11 (2003), 93-105.

[21] X. Ma, Y. Zhou and H. Cao, Global stability of the endemic equilibrium of a discrete SIR epidemic model, Adv. Difference

Equ. 2013 (2013), 42.

[22] L. Allen, An Introduction to Deterministic Models in Biology, Prentice-Hall, 2004.

[23] S. Elaydi, Discrete Chaos, Chapman and Hall/CRC, 2007.

[24] May, R. M., Parasitic infections as regulators of animal populations. Amer. Scientist, 71 (1983), 36-45.


	1. Introduction
	2. SIR Epidemic Model with Vaccination
	3. The Discrete Time System
	4. Stability of Equilibrium points and Numerical Simulations
	5. Bifurcation
	6. Conclusion
	References