Stability Analysis of a Fractional Epidemic Model Involving Vaccination Effect

1. Introduction

Compartmental models for infectious diseases separates a population into various classes according to the stages of infection. The reflected rates of transition between compartments are stated by time derivatives of the population sizes in each compartment, and so such type models are formulated by means of differential equations. Mathematical models, in which the rates of transfer between compartments depend on the sizes of compartments in the past or at the moment of transfer, especially require working with differential, integral or integro-differential equations.

The early studies related to the compartmental models in mathematical epidemiology are based on Kermack and McKendrick, [1]. In 1927 they introduced a compartmental epidemic model based on basic transfers between divided groups in a population. To explain the compartments, they divided the population into three groups referred as $S,$I and $R.$ $S\left(t\right)$ denotes the number of individuals who are susceptible against the disease at time t, in other words, who are not yet infected and have not any immunity. $I\left(t\right)$ represents the number of infectious, that is, members of $I\left(t\right)$ are able to spread the disease to susceptibles with effective contact. $R\left(t\right)$ shows the number of recovered individuals, who have immunity against the pathogen and so no there exists any probability of spread of infection via these individuals.

After the model given by Kermack and McKendrick, mathematical epidemiology has been developed with numerous studies providing various contributions to the literature in this field. [2,3,4,5] are just a few of the studies in this area.

Vaccination of susceptibles against the infection is one of the effective control measures in the process related to the struggle with diseases. So far, a lot of studies explaining the effects of vaccination upon the spread of the diseases has been introduced, see [6,7,8,9] and references therein.

In recent years, some fractional-order differential equation models of infectious disease dynamics have been introduced with the Caputo derivation. Some authors have extended classical epidemic models to fractional-order epidemic systems and discussed the stability of equilibrias, [10,11,12].

In this study, it will be followed by the terminology of the $SIR$ model to explain a new model related to disease that confers immunity against reinfection, by adding the vaccinated class, too. This paper reveals and analyzes an epidemic model, considering vaccination effect in the spread of any disease. For analysis of the model, firstly disease-free and endemic equilibrium points have been determined. Then the basic reproduction number related to the model has been obtained by using the next generation matrix method. Also, it is shown that the disease-free and the endemic equilibria are locally and globally asymptotically stable for ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_0>1$, respectively.

Local stabilities of equilibrium points have been researched with the way of analyzing the corresponding characteristic equation. To prove global stabilities of equilibria it have been used LaSalle Invariance Principle associated with Lyapunov function and Dulac Criteria, respectively.

2. Some Basic Mathematical Properties and Model Structure

The definitions of the fractional integral and Caputo fractional derivative are defined as follows, [13].

Definition 1.

Suppose that $\alpha >0,$ $a\in \mathbb{R},$ $t>a$ and x be an integrable function. Then the fractional integral of x of order α is defined as

$$I_{a+}^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_a^t{\left(t-s\right)}^{\alpha -1}x\left(s\right)ds.$$

The Caputo fractional derivative of x of order $\alpha$ is defined as

$${}^{C}D_{a+}^{\alpha }x\left(t\right)={\left({\displaystyle \frac{d}{dt}}\right)}^n{I}_{a+}^{n-\alpha }\left[x\left(t\right)-\sum _{k=0}^{n-1}{\displaystyle \frac{x^{\left(k\right)}\left(a\right)}{k!}}{\left(t-a\right)}^k\right],t>a,$$

where $n=\left[\alpha \right]+1$. When $\alpha$ is an integer number, ${}^{C}D_{a+}^{\alpha }x$ is usual derivative of order $\alpha$ of $x.$ If x is a function of class $C^n,$ its fractional derivative is represented by

$${}^{C}D_{a+}^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\Gamma \left(n-\alpha \right)}}{\int }_a^t{\left(t-s\right)}^{n-\alpha -1}{x}^{\left(n\right)}\left(s\right)ds.$$

Especially we write ${}^{C}D^{\alpha }x$ instead of ${}^{C}D_{0+}^{\alpha }x.$ It is obvious that the Caputo fractional derivative of order $\alpha$ of x

$${}^{C}D_{a+}^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\Gamma \left(1-\alpha \right)}}{\int }_0^t{\displaystyle \frac{x^{\prime }\left(s\right)ds}{{\left(t-s\right)}^{\alpha }}}$$

for $0<\alpha <1.$

On the other hand let us recall some properties of Laplace transforms. Laplace transform of the Caputo fractional derivative ([10]) is given by

$$\mathcal{L}\left\{{}^{C}D^{\alpha }x\left(t\right)\right\}={\lambda }^{\alpha }\mathcal{L}\left\{x\left(t\right)\right\}-\sum _{k=0}^{n-1}{x}^{\left(k\right)}\left(0\right){\lambda }^{\alpha -k-1}.$$

(1)

Also Laplace transform of Mittag–Leffler function defined by the power series

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{z^k}{\Gamma \left(\alpha k+\beta \right)}}$$

hold the following equality, [10]

$$\mathcal{L}\left\{t^{\beta -1}{E}_{\alpha ,\beta }\left(-a{t}^{\alpha }\right)\right\}={\displaystyle \frac{{\lambda }^{\alpha -\beta }}{{\lambda }^{\alpha }+a}}.$$

(2)

Theorem 1.

Let $x^{*}$ be an equilibrium point for ${}^{C}D^{\alpha }x\left(t\right)=f\left(t,x\right)$ and Ω be a domain containing $x^{*}.$ Let $L:\left[0,\infty \right)\times \Omega \to \mathbb{R}$ be a continuously differentiable function such that

$$W_1\left(x\right)\le L\left(t,x\right)\le W_2\left(x\right),$$

$${}^{C}D^{\alpha }L\left(t,x\right)\le -W_3\left(x\right)$$

hold for all $t\ge 0,$ $x\in \Omega$ and $\alpha \in \left(0,1\right),$ where functions $W_i$ are continuous positive definite functions on $\Omega .$ Then the equilibrium point of system ${}^{C}D^{\alpha }x\left(t\right)=f\left(t,x\right)$ is uniformly asymptotically stable, [11,12].

2.1. Model Structure

We introduce a new fractional $SVIR$ mathematical compartmental model expressed by the system of the following nonlinear fractional integro-differential equations with all parameters being positive:

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }S\left(t\right)& =& b-\beta S\left(t\right)I\left(t\right)-qS\left(t\right)-\mu S\left(t\right)\hfill \\ \hfill {}^{C}D^{\alpha }V\left(t\right)& =& qS\left(t\right)-\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\mu V\left(t\right)\hfill \\ \hfill {}^{C}D^{\alpha }I\left(t\right)& =& \beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I\left(t\right)\hfill \\ \hfill {}^{C}D^{\alpha }R\left(t\right)& =& \gamma I\left(t\right)-\mu R\left(t\right),\hfill \end{array}$$

(3)

where $\alpha \in \left(0,1\right]$ and functions $S\left(t\right),$ $V\left(t\right),$ $I\left(t\right)$ and $R\left(t\right)$ represent the numbers of the susceptible, vaccinated, infectious and recovered individuals at time t, respectively. The total population size is $N\left(t\right)$ and $N\left(t\right)=S\left(t\right)+V\left(t\right)+I\left(t\right)+R\left(t\right)$ for all $t\ge 0$. Also all these functions are nonnegative.

In the model, all newborn individuals get involved in the population by entering to S with a constant rate $b.$ $\beta$ denotes the effective contact rate between infectious and susceptible or vaccinated individuals. $\mu$ is the natural death rate in each compartment and $\delta$ is the death rate welded from the outbreak. The rate of vaccinated individuals within susceptibles is represented by q. Also $\gamma$ denotes the transition rate from infectious compartment to compartment R.

Vaccination may not be always completely effective for individuals. Therefore, in this model, we suppose that the vaccine provides temporary immunity effect or permanent effect for some individuals. In this context, we use the distribution parameter $\theta$ for meaning the protection period caused by vaccination. Essentially, the term $\theta$ designates the protection efficacy of the vaccine. That is, $\theta =0$ means that the vaccine is completely ineffective. Also, $0<\theta <\infty$ means that the vaccinated individuals have only partial protection. Otherwise it means that the vaccine is wholly effective against the pathogen.

On the other hand, we assume that the protection provided by the vaccination may vary according to individual. So, we use a distribution function f to reflect this fact to the model. f is a distribution function and $f\left(\theta \right)$ is the ratio of the individuals whose protection period provided by the vaccine is $\theta .$ Classically, it is supposed that $f\left(\theta \right)$ is non-negative for each $\theta$ and f is continuous on ${\mathbb{R}}_{+}$, also f satisfies $\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)d\theta =1$. The term $q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta$ represents the number of surviving individuals at time t who have been vaccinated at time $t-\theta$ and have protection period $\theta$.

Of course, since it is thought that the vaccine may not provide full protection, some of vaccinated individuals may still be susceptible to infection. Therefore, as a result of a sufficiently effective contact with infectious individuals, a vaccinated individual who no longer has any protection enters to the infectious compartment. To reflect this transition to the compartment I, the expression $\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta$ has been used in the model.

3. Analysis of the Model

In this part, for the model, the convenient and positive invariant region, equilibrium points and basic reproduction number have been determined.

3.0.1. Feasible Region

Theorem 2.

For the system (3), the set

$$\Omega =\left\{\left(S,V,I,R\right):S\in BC\left(\left[-\theta ,\infty \right),{\mathbb{R}}_{+}\right);V,I,R\in BC\left(\left[0,\infty \right),{\mathbb{R}}_{+}\right)\mathrm{and}N\left(t\right)\le {\displaystyle \frac{b}{\mu }}\right\}$$

(4)

is positively invariant region in which the solutions of the model is bounded.

Proof. (Following [10]) Adding the four equations in (3), we write

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }N\left(t\right)& =& b-\mu \left(S\left(t\right)+V\left(t\right)+I\left(t\right)+R\left(t\right)\right)-\delta I\hfill \\ & \le & b-\mu N\left(t\right).\hfill \end{array}$$

(5)

Then we can get

$${\lambda }^{\alpha }\mathcal{L}\left\{N\left(t\right)\right\}-{\lambda }^{\alpha -1}N\left(0\right)\le {\displaystyle \frac{b}{\lambda }}-\mu \mathcal{L}\left\{N\left(t\right)\right\}$$

(6)

from (1) and other properties of Laplace transform. If we take $N\left(0\right)\le b/\mu$ and consider the inequality (6), we get

$$\left({\lambda }^{\alpha }+\mu \right)\mathcal{L}\left\{N\left(t\right)\right\}\le {\displaystyle \frac{b}{\lambda }}+{\lambda }^{\alpha -1}N\left(0\right)$$

and

$$\begin{array}{ccc}\hfill \mathcal{L}\left\{N\left(t\right)\right\}& \le & {\displaystyle \frac{b}{\lambda \left({\lambda }^{\alpha }+\mu \right)}}+\left({\displaystyle \frac{{\lambda }^{\alpha -1}}{{\lambda }^{\alpha }+\mu }}\right){\displaystyle \frac{b}{\mu }}\hfill \\ & \le & {\displaystyle \frac{b}{\lambda \mu }}+{\displaystyle \frac{b}{\mu }}\mathcal{L}\left\{t^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\mu t^{\alpha }\right)\right\}\hfill \end{array}$$

(7)

Applying the inverse Laplace transform to (7) and considering the boundedness of Mittag–Leffler function $E_{\alpha ,\alpha }\left(-\mu t^{\alpha }\right)$ for all $t>0$, we can write

$$\begin{array}{ccc}\hfill N\left(t\right)& \le & {\displaystyle \frac{b}{\mu }}\left(1+t^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\mu t^{\alpha }\right)\right)\hfill \\ & \le & {\displaystyle \frac{b}{\mu }}\left(1+t^{\alpha -1}C\right),\hfill \end{array}$$

(8)

where C is a constant such that $E_{\alpha ,\alpha }\left(-\mu t^{\alpha }\right)\le C$ hold for all $t>0.$

So, if $\alpha \in \left(0,1\right),$ we get $N\left(t\right)\le b/\mu$ for all $t>0$ as $t\to \infty$ from (8). Also, if $\alpha =1$ then solution of (5) comes with the solution of differential equation $N^{\prime }\left(t\right)+\mu N\left(t\right)=b.$ If the integrating factor method is processed, it is obtained that $N\left(t\right)=b/\mu +c{e}^{-\mu t}$ and so

$$N\left(t\right)=N\left(0\right)e^{-\mu t}+{\displaystyle \frac{b}{\mu }}\left(1-e^{-\mu t}\right)$$

(9)

is obtained for the initial condition $t=0.$ Standard Comparison Theorem [14] says that the right side of (9) is the maximal solution of Equation (5). Thus it is reached to the inequality

$$N\left(t\right)\le N\left(0\right)e^{-\mu t}+{\displaystyle \frac{b}{\mu }}\left(1-e^{-\mu t}\right)$$

for all $t\ge 0.$

Hence, when $N\left(0\right)\le b/\mu ,$ we obtain $N\left(t\right)\le b/\mu$ for all $t>0$ and $\alpha \in \left(0,1\right]$. This means that $\Omega$ is positively invariant for the system (3 ). □

Because the populations $V\left(t\right)$ and $R\left(t\right)$ do not feature in remainder equations of system (3), it will be enough to consider with the reduced system (10)

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }S\left(t\right)& =& b-\beta S\left(t\right)I\left(t\right)-\left(q+\mu \right)S\left(t\right),\hfill \\ \hfill {}^{C}D^{\alpha }I\left(t\right)& =& \beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I\left(t\right).\hfill \end{array}$$

(10)

3.0.2. Existence of Solutions of the System

In this section, we focus on the existence of solution of the problem

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }x\left(t\right)& =& h\left(x^t\right),t\ge 0\hfill \\ \hfill x_0\left(t\right)& =& g\left(t\right),-\theta \le t<0.\hfill \end{array}$$

(11)

Where $g=\left(g_1\left(t\right),g_2\left(t\right)\right)$ represents the initial functions of system (10) and $g\in C\left(\left[-\theta ,0\right),{\left[0,{\displaystyle \frac{b}{\mu }}\right]}^2\right).$ Also let us take $h:\Psi \subset C\left(\left[-\theta ,0\right),{\left[0,{\displaystyle \frac{b}{\mu }}\right]}^2\right)\to {\left[0,{\displaystyle \frac{b}{\mu }}\right]}^2,$ $x^t\left(s\right)=x\left(t+s\right)$ and $x\in \left\{\left(x_1,x_2\right):\left(x_1,·,x_2,·\right)\in \Omega \mathrm{defined}\mathrm{by}\right\}\subset C\left(\left[-\theta ,\infty \right),{\left[0,{\displaystyle \frac{b}{\mu }}\right]}^2\right)$. If we choose the function h as $h=\left(h_1,h_2\right)$ such that

$$\begin{array}{ccc}\hfill h_1\left(x\right)& =& b-\beta x_1\left(0\right)x_2\left(0\right)-\left(q+\mu \right)x_1\left(0\right)\hfill \\ \hfill h_2\left(x\right)& =& \beta x_1\left(0\right)x_2\left(0\right)+\beta q{x}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)x_2\left(-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)x_1\left(0\right)\hfill \end{array}$$

(12)

and $x=\left(x_1,x_2\right)=\left(S,I\right)$ then we can say that the finding the solution of the problem (11) is equivalent to solving the system (10) or equivalently the following problem.

$$\begin{array}{cc}\begin{array}{c}{}^{C}D^{\alpha }S\left(t\right)=b-\beta S\left(t\right)I\left(t\right)-\left(q+\mu \right)S\left(t\right)\\ {}^{C}D^{\alpha }I\left(t\right)=\beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I\left(t\right)\end{array},& t\in \left[0,\infty \right)\\ \begin{array}{c}S\left(t\right)=g_1\left(t\right)\\ I\left(t\right)=g_2\left(t\right)\end{array},& t\in \left[-\theta ,0\right)\end{array}$$

(13)

Also we can say that the Equation (11) has a unique solution if h is Lipschitz continuous in every compact subset $M\subset \Psi$. Indeed this result depends on Schauder fixed point theorem, [15].

In this study, the fact that the set $C=C\left(\left[-\theta ,0\right),{\mathbb{R}}^2\right)$ is a Banach space with the norm

$${∥x∥}_C=sup\left\{\left|x_1\left(t\right)\right|+\left|x_2\left(t\right)\right|:-\theta \le t\le 0\right\}$$

is also taken into consideration.

Theorem 3.

There exists a unique solution of the system (13), or equivalently, of the Equation (11) with $h:\Psi \to {\left[0,{\displaystyle \frac{b}{\mu }}\right]}^2$, defined by (12).

Proof. 

The proof depends on the result in [15]. So it is sufficient to show that h is Lipschitz continuous in every compact subset $M\subset \Psi .$ Let $x=\left(x_1,x_2\right)\in M$, and $y=\left(y_1,y_2\right)\in M$. Then considering $\left|x_i\left(t\right)\right|\le {\displaystyle \frac{b}{\mu }}$ and $\left|y_i\left(t\right)\right|\le {\displaystyle \frac{b}{\mu }}$ for $-\tau \le t\le 0,$ $i=1,2,$ we can write from the description of h

$$\begin{array}{ccc}& & ∥h\left(x\right)-h\left(y\right)∥\hfill \\ & \le & \left|h_1\left(x\right)-h_1\left(y\right)\right|+\left|h_2\left(x\right)-h_2\left(y\right)\right|\hfill \\ & =& \left|\beta \left[y_1\left(0\right)y_2\left(0\right)-x_1\left(0\right)x_2\left(0\right)\right]+\left(q+\mu \right)\left[y_1\left(0\right)-x_1\left(0\right)\right]\right|\hfill \\ & & +\left|\beta \left[x_1\left(0\right)x_2\left(0\right)-y_1\left(0\right)y_2\left(0\right)\right]+\left(\gamma +\delta +\mu \right)\left[y_1\left(0\right)-x_1\left(0\right)\right]\right|\hfill \\ & & +\left|\beta q{x}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)x_2\left(-\theta \right)e^{-\mu \theta }d\theta -\beta q{y}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)y_2\left(-\theta \right)e^{-\mu \theta }d\theta \right|\hfill \\ & =& \left|\beta \left[y_1\left(0\right)y_2\left(0\right)-y_1\left(0\right)x_2\left(0\right)+y_1\left(0\right)x_2\left(0\right)-x_1\left(0\right)x_2\left(0\right)\right]+\left(q+\mu \right)\left[y_1\left(0\right)-x_1\left(0\right)\right]\right|\hfill \\ & & +\left|\beta \left[x_1\left(0\right)x_2\left(0\right)-x_1\left(0\right)y_2\left(0\right)+x_1\left(0\right)y_2\left(0\right)-y_1\left(0\right)y_2\left(0\right)\right]+\left(\gamma +\delta +\mu \right)\left[y_1\left(0\right)-x_1\left(0\right)\right]\right|\hfill \\ & & +\left|\beta q{x}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)x_2\left(-\theta \right)e^{-\mu \theta }d\theta -\beta q{y}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)x_2\left(-\theta \right)e^{-\mu \theta }d\theta \right.\hfill \\ & & \left.+\beta q{y}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)x_2\left(-\theta \right)e^{-\mu \theta }d\theta -\beta q{y}_1\left(0\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)y_2\left(-\theta \right)e^{-\mu \theta }d\theta \right|\hfill \\ & \le & \beta \left(\left|y_1\left(0\right)\right|\left|x_2\left(0\right)-y_2\left(0\right)\right|+\left|x_2\left(0\right)\right|\left|x_1\left(0\right)-y_1\left(0\right)\right|\right)\hfill \\ & & +\beta \left(\left|x_1\left(0\right)\right|\left|x_2\left(0\right)-y_2\left(0\right)\right|+\left|y_2\left(0\right)\right|\left|x_1\left(0\right)-y_1\left(0\right)\right|\right)\hfill \\ & & +\left(q+2\mu +\gamma +\delta \right)\left|y_1\left(0\right)-x_1\left(0\right)\right|\hfill \\ & & +\beta q\left|x_1\left(0\right)-y_1\left(0\right)\right|\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)\left|x_2\left(-\theta \right)\right|e^{-\mu \theta }d\theta \hfill \\ & & +\beta q\left|y_1\left(0\right)\right|\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)\left|x_2\left(-\theta \right)-y_2\left(-\theta \right)\right|e^{-\mu \theta }d\theta \hfill \end{array}$$

$$\begin{array}{ccc}& \le & 2\beta {\displaystyle \frac{b}{\mu }}\left(\left|x_1\left(0\right)-y_1\left(0\right)\right|+\left|x_2\left(0\right)-y_2\left(0\right)\right|\right)\hfill \\ & & +\left(q+2\mu +\gamma +\delta \right)\left|y_1\left(0\right)-x_1\left(0\right)\right|\hfill \\ & & +\beta q{\displaystyle \frac{b}{\mu }}F\left(\left|x_1\left(0\right)-y_1\left(0\right)\right|+\left|x_2\left(-\theta \right)-y_2\left(-\theta \right)\right|\right)\hfill \\ & \le & \left(2\beta {\displaystyle \frac{b}{\mu }}+q+2\mu +\gamma +\delta +\beta q{\displaystyle \frac{b}{\mu }}F\right){∥x-y∥}_C.\hfill \end{array}$$

and so we conclude

$$∥h\left(x\right)-h\left(y\right)∥\le \left({\displaystyle \frac{\beta b}{\mu }}\left(2+qF\right)+q+2\mu +\gamma +\delta \right){∥x-y∥}_C.$$

Thus if we take

$$K\ge {\displaystyle \frac{\beta b}{\mu }}\left(2+qF\right)+q+2\mu +\gamma +\delta$$

then the inequality

$$∥h\left(x\right)-h\left(y\right)∥\le K{∥x-y∥}_C$$

holds in every compact subset $M\subset \Psi .$ As a result, it is concluded that there is only one solution of system (13) since h satisfies the Lipschitz inequality in every compact M subset of $\Psi$. □

3.0.3. Disease-Free Equilibria, Basic Reproduction Number

Since the equilibria of the proposed model are the solutions of the system (10) then disease-free equilibrium point which we will show with $P_{DF}=\left(S_0,I_0\right)$ provides the equations in the system.

Therefore, system (10) always has a disease-free equilibria

$$P_{DF}=\left(S_0,I_0\right)=\left({\displaystyle \frac{b}{q+\mu }},0\right).$$

(14)

Now we will determine the threshold value ${\mathcal{R}}_0$ referred as basic reproduction number is the number of secondary infections caused by one infectious individual. This parameter allows us to have an idea about the dynamics of the epidemic and to make predictions. For this reason it is a very important value in epidemiology. Using the next generation matrix method, the value of ${\mathcal{R}}_0$ related to the model (10) is calculated by following terminology below, [16,17].

First of all, let us specify that for the value of the integral $\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)e^{-\mu \theta }d\theta$ briefly will be used notation F, throughout the remainder of the text.

Let $X={(I,S)}^T$. Then we can write the system (10) in the form

$$\left[\begin{array}{c}{}^{C}D^{\alpha }I\left(t\right)\\ {}^{C}D^{\alpha }S\left(t\right)\end{array}\right]=\underset{\mathcal{M}\left(X\right)}{\underbrace{\left[\begin{array}{c}\beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta \\ 0\end{array}\right]}}-\underset{\mathcal{N}\left(X\right)}{\underbrace{\left[\begin{array}{c}\left(\gamma +\delta +\mu \right)I\left(t\right)\\ \beta S\left(t\right)I\left(t\right)+\left(q+\mu \right)S\left(t\right)-b\end{array}\right]}},$$

that is

$${}^{C}D^{\alpha }X\left(t\right)=\mathcal{M}\left(X\right)-\mathcal{N}\left(X\right).$$

The values at $P_{DF}$ of the derivatives of $\mathcal{M}\left(X\right)$ and $\mathcal{N}\left(X\right)$ with respect to $I,$ $S,$ respectively, come in sight with the following Jacobian matrices:

$$d\mathcal{M}\left(P_{DF}\right)=\left[\begin{array}{cc}\beta S_0\left(1+qF\right)& \beta I_0\left(1+qF\right)\\ 0& 0\end{array}\right],$$

$$d\mathcal{N}\left(P_{DF}\right)=\left[\begin{array}{cc}\gamma +\delta +\mu & 0\\ \beta S_0& \beta I_0+q+\mu \end{array}\right].$$

Thus we can obtain

$$M={\mathcal{M}}_{1x1}=\left[\beta S_0\left(1+qF\right)\right],$$

$$N={\mathcal{N}}_{1x1}=\left[\gamma +\delta +\mu \right]$$

and

$$M{N}^{-1}=\left[{\displaystyle \frac{\beta S_0\left(1+qF\right)}{\gamma +\delta +\mu }}\right].$$

So by using the characteristic polynomial of $M{N}^{-1}$ we get the spectral radius as

$$\rho \left(M{N}^{-1}\right)={\displaystyle \frac{\beta S_0\left(1+qF\right)}{\gamma +\delta +\mu }}.$$

Taking into account that $S_0$, the basic reproduction number of the system (10) is calculated in the form of

$${\mathcal{R}}_0=\rho \left(M{N}^{-1}\right)={\displaystyle \frac{b\beta \left(1+qF\right)}{\left(q+\mu \right)\left(\gamma +\delta +\mu \right)}}.$$

In the following part, the results about stability dynamics belonging to the epidemic model (10) have been discussed under the titles "Disease-free case" and "Endemic case", respectively.

3.0.4. Disease-Free Case

Under this title, we consider local and global stabilities of $P_{DF}=\left(S_0,I_0\right)$ for the system (10).

Theorem 4.

$P_{DF}$ is locally asymptotically stable for ${\mathcal{R}}_0<1.$

Proof. 

The Jacobian matrix at $P_{DF}=\left(S_0,I_0\right)$ of system (10) is

$$J\left(P_{DF}\right)=\left[\begin{array}{ccc}-\beta I_0-q-\mu & & -\beta S_0\\ & & \\ \beta I_0\left(1+qF\right)& & \beta S_0\left(1+qF\right)-\left(\gamma +\delta +\mu \right)\end{array}\right].$$

For the point $\left(S_0,I_0\right)=\left({\displaystyle \frac{b}{q+\mu }},0\right),$ the characteristic equation for this matrix is

$$\left(\lambda +q+\mu \right)\left(\lambda -{\displaystyle \frac{b\beta \left(1+qF\right)-\left(q+\mu \right)\left(\gamma +\delta +\mu \right)}{q+\mu }}\right)=0.$$

(15)

Using the notations ${\lambda }_1$ and ${\lambda }_2$ for the roots of Eq. ( 15) then

$${\lambda }_1=-\left(q+\mu \right)$$

and

$$\begin{array}{ccc}\hfill {\lambda }_2& =& {\displaystyle \frac{\left(q+\mu \right)\left(\gamma +\delta +\mu \right)\left(\frac{b\beta \left(1+qF\right)}{\left(q+\mu \right)\left(\gamma +\delta +\mu \right)}-1\right)}{q+\mu }}\hfill \\ & =& \left(\gamma +\delta +\mu \right)\left({\mathcal{R}}_0-1\right).\hfill \end{array}$$

When ${\mathcal{R}}_0<1,$ two roots of Equation (15) are negative. If ${\mathcal{R}}_0=1,$ then one of roots of Eq.(15) is zero. On the other hand, for ${\mathcal{R}}_0>1$ one of roots of Equation (15) has positive real parts. So $P_{DF}$ is locally asymptotically stable for ${\mathcal{R}}_0<1;$ is stable for ${\mathcal{R}}_0=1,$ and is unstable for ${\mathcal{R}}_0>1.$ □

Theorem 5.

$P_{DF}$ is uniformly asymptotically stable for ${\mathcal{R}}_0<1.$

Proof. 

Consider the nonnegative function L defined as

$$L\left(t,S,I\right)=I\left(t\right).$$

Computing the time derivative of $L,$ we obtain

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }L\left(t,S,I\right)& =& \beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I\left(t\right)\hfill \\ & =& \left(\beta S\left(t\right)+\beta q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)\right)I\left(t\right)\hfill \\ & \le & \left({\displaystyle \frac{b\beta }{q+\mu }}+{\displaystyle \frac{b\beta qF}{q+\mu }}-\left(\gamma +\delta +\mu \right)\right)I\left(t\right)\hfill \\ & =& \left(\gamma +\delta +\mu \right)\left({\displaystyle \frac{b\beta \left(1+qF\right)}{\left(q+\mu \right)\left(\gamma +\delta +\mu \right)}}-1\right)I\left(t\right)\hfill \\ & =& \left(\gamma +\delta +\mu \right)\left({\mathcal{R}}_0-1\right)I\left(t\right).\hfill \end{array}$$

When ${\mathcal{R}}_0<1,$ we say that ${}^{C}D^{\alpha }L\left(t,S,I\right)\le 0$ and ${}^{C}D^{\alpha }L\left(t,S,I\right)=0$ for $S=S_0,$ $I=I_0.$ By Theorem 1, the disease free equilibrium $P_{DF}$ is uniformly asymptotically stable in the interior of $\Omega$, when ${\mathcal{R}}_0<1$. □

3.0.5. Endemic Case

In this part, firstly the existence and uniqueness of the endemic equilibrium point of the model has been shown. Then, in addition to the course of the disease when ${\mathcal{R}}_0>1$ the local and global stability of this equilibrium point has been investigated.

3.0.6. Existence of the Endemic Equilibria

Since ${}^{C}D^{\alpha }{I}^{*}{=}^C{D}^{\alpha }{S}^{*}=0,$ the endemic equilibria denoted by $P_E=\left(S^{*},I^{*}\right)$ satisfies the algebraic equations with $I^{*}\ne 0$

$$\begin{array}{ccc}\hfill 0& =& b-\beta S^{*}{I}^{*}-\left(q+\mu \right)S^{*},\hfill \\ \hfill 0& =& \beta S^{*}{I}^{*}+\beta q{S}^{*}{I}^{*}\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I^{*}.\hfill \end{array}$$

(16)

From the second equation of (16) we write

$$0=I^{*}\left(\beta S^{*}+\beta qF{S}^{*}-\left(\gamma +\delta +\mu \right)\right).$$

Since $I^{*}\ne 0$, it must be $\beta S^{*}+\beta qF{S}^{*}-\left(\gamma +\delta +\mu \right)=0.$ Then the endemic equilibrium point comes with

$$\begin{array}{ccc}\hfill P_E& =& \left(S^{*},I^{*}\right)\hfill \\ & =& \left({\displaystyle \frac{\gamma +\delta +\mu }{\beta \left(1+qF\right)}},{\displaystyle \frac{b-\left(q+\mu \right)S^{*}}{\beta S^{*}}}\right).\hfill \end{array}$$

By the way, we say that $P_E$ can be expressed as

$$P_E=\left({\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}},{\displaystyle \frac{\left(q+\mu \right)\left({\mathcal{R}}_0-1\right)}{\beta }}\right)$$

in terms of ${\mathcal{R}}_0.$ Also it is seen clearly that the system (10) has a unique endemic equilibria iff ${\mathcal{R}}_0>1.$

3.0.7. Course of the Disease for ${\mathcal{R}}_0>1$

Now, we will focus on how the disease will progress (whether it will disappear or not) in the population when ${\mathcal{R}}_0>1$. For this, we will assume the infectious population at the beginning exists, and we will try to see course of the infectious population as time progresses for ${\mathcal{R}}_0>1$.

To do this, let us suppose $I\left(0\right)>0$ and $\underset{t\to \infty }{lim}I\left(t\right)=0.$ Then, for any given $ϵ>0,$ there exists a $T_1>0$ such that $I\left(t\right)<ϵ$ holds for $t>T_1.$ Specially, we can choose an $ϵ$ such that

$$0<ϵ<b\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_0}}\right),$$

(17)

and

$$I\left(t\right)\le {\displaystyle \frac{ϵ\mu }{b\beta }}$$

holds for all $t>T_1$.

Then

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }S\left(t\right)& =& b-\beta S\left(t\right)I\left(t\right)-\left(q+\mu \right)S\left(t\right)\hfill \\ & \ge & b-{\displaystyle \frac{\beta b}{\mu }}I\left(t\right)-\left(q+\mu \right)S\left(t\right)\hfill \\ & \ge & b-{\displaystyle \frac{\beta bϵ\mu }{\mu b\beta }}-\left(q+\mu \right)S\left(t\right)\hfill \\ & =& b-\left(q+\mu \right)S\left(t\right)-ϵ.\hfill \end{array}$$

Therefore we say that

$${}^{C}D^{\alpha }S\left(t\right)\ge b\left(1-{\displaystyle \frac{q+\mu }{b}}\right)S\left(t\right)-ϵ$$

for $t>T_1.$ This requires

$$\underset{t\to \infty }{liminf}S\left(t\right)\ge {\displaystyle \frac{b}{q+\mu }}\left(1-{\displaystyle \frac{ϵ}{b}}\right).$$

If we take into account the basic properties of limit inferior and remember the choosing (17) then we say that there exists a $T_2>0$ such that

$$S\left(t\right)\ge {\displaystyle \frac{b}{q+\mu }}\left(1-{\displaystyle \frac{b\left(1-\frac{1}{{\mathcal{R}}_0}\right)}{b}}\right)={\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}}$$

for all $t>T_2$. Thus

$$S\left(t\right)\ge {\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}}\mathrm{and}S(t-\theta )\ge {\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}}$$

hold for all $t>T_2+\theta$. If T is chosen as

$$T=max\left\{T_1,T_2+\theta \right\}$$

then

$$\begin{array}{ccc}\hfill {}^{C}D^{\alpha }I\left(t\right)& =& I\left(t\right)\left(\beta S\left(t\right)+\beta q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)\right)\hfill \\ & >& I\left(t\right)\left(\beta {\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}}+\beta q{\displaystyle \frac{b}{\left(q+\mu \right){\mathcal{R}}_0}}F-\left(\gamma +\delta +\mu \right)\right)\hfill \\ & =& I\left(t\right)\left(\gamma +\delta +\mu \right)\left({\displaystyle \frac{b\beta \left(1+qF\right)}{\left(q+\mu \right)\left(\gamma +\delta +\mu \right){\mathcal{R}}_0}}-1\right)\hfill \\ & =& 0\hfill \end{array}$$

for $t>T$. This result contradicts with $\underset{t\to \infty }{lim}I\left(t\right)=0$. Moreover, the fact ${}^{C}D^{\alpha }I\left(t\right)>0$ for $0<\alpha <1$ means that $I\left(t\right)$ is increasing, [18]. But, when ${\mathcal{R}}_0>1$ it is not possible that $\underset{t\to \infty }{lim}I\left(t\right)=0$ with the initial assumption $I\left(0\right)>0$. All these imply that the disease will not disappear in the population when ${\mathcal{R}}_0>1.$ This case is what should be happen for a consistent model and a meaningful threshold ${\mathcal{R}}_0.$

Now let’s analyze the behavior of function I which do not converge to zero when ${\mathcal{R}}_0>1$ and of the function S.

3.0.8. Local and Global Asymptotic Stability of the Endemic Equilibrium Point

Theorem 6.

$P_E$ is locally asymptotically stable in Ω for ${\mathcal{R}}_0>1$ and $\alpha \in \left(0,1\right].$

Proof. 

The Jacobian matrix of system (10) at $P_E=\left(S^{*},I^{*}\right)$ is

$$J\left(P_E\right)=\left[\begin{array}{ccc}-\beta I^{*}-q-\mu & & -\beta S^{*}\\ & & \\ \beta I^{*}\left(1+qF\right)& & \beta S^{*}\left(1+qF\right)-\left(\gamma +\delta +\mu \right)\end{array}\right].$$

Herefrom, the characteristic equation of the matrix $J\left(P_E\right)$ forms with the determinant

$$det\left(\lambda I_2-J\left(P_E\right)\right)=\left|\begin{array}{ccc}-\left(q+\mu \right){\mathcal{R}}_0-\lambda & & -{\displaystyle \frac{b\beta }{\left(q+\mu \right){\mathcal{R}}_0}}\\ & & \\ \left(q+\mu \right)\left(1+qF\right)\left({\mathcal{R}}_0-1\right)& & -\lambda \end{array}\right|=0.$$

After the rearrangements, this equation is written as

$${\lambda }^2+\left(q+\mu \right){\mathcal{R}}_0\lambda +b\beta \left(1+qF\right)\left({\mathcal{R}}_0-1\right)=0.$$

(18)

Let two roots of Equation (18) be ${\lambda }_1$ and ${\lambda }_2.$ Then

$${\lambda }_1+{\lambda }_2=-\left(q+\mu \right){\mathcal{R}}_0$$

and

$${\lambda }_1{\lambda }_2=b\beta \left(1+qF\right)\left({\mathcal{R}}_0-1\right).$$

So we say that ${\lambda }_1{\lambda }_2>0$ and ${\lambda }_1+{\lambda }_2<0$ for ${\mathcal{R}}_0>1$. In that case, two roots of Eq. (18) are negative. Since all eigenvalues of the Jacobian matrix of system (10) at $P_E$ have negative real parts, $P_E=\left(S^{*},I^{*}\right)$ is locally asymptotically stable. □

Now, we will prove that $P_E$ is globally stable for special case $\alpha =1$ by using Dulac criterion and the Poincaré-Bendixson theorem.

Theorem 7.

$P_E$ is globally asymptotically stable for ${\mathcal{R}}_0>1$.

Proof. 

We will show that (10) does not have any periodic solutions in the positive quadrant of the $SI$ plane. Let us establish a continuous function $\Phi \left(S,I\right)$ defined by

$$\Phi \left(S,I\right)={\displaystyle \frac{1}{\beta SI}},S>0,I>0$$

for system (10). Setting G and H as

$$\begin{array}{ccc}\hfill G\left(S,I\right)& =& b-\beta S\left(t\right)I\left(t\right)-\left(q+\mu \right)S\left(t\right),\hfill \\ \hfill H\left(S,I\right)& =& \beta S\left(t\right)I\left(t\right)+\beta qI\left(t\right)\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta -\left(\gamma +\delta +\mu \right)I\left(t\right),\hfill \end{array}$$

we obtain

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \left(\Phi G\right)}{\partial S}}+{\displaystyle \frac{\partial \left(\Phi H\right)}{\partial I}}& =& {\displaystyle \frac{\partial }{\partial S}}\left\{{\displaystyle \frac{b}{\beta SI}}-1-{\displaystyle \frac{q+\mu }{\beta I}}\right\}\hfill \\ & +& {\displaystyle \frac{\partial }{\partial I}}\left\{1+{\displaystyle \frac{q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)S\left(t-\theta \right)e^{-\mu \theta }d\theta }{S}}-1-{\displaystyle \frac{\gamma +\delta +\mu }{\beta S}}\right\}\hfill \\ & =& -{\displaystyle \frac{b\beta I}{{\left(\beta SI\right)}^2}}\hfill \end{array}$$

for all $S>0,$ $I>0$. We see the last expression is clearly negative. Since it has the same sign almost everywhere in the positive quadrant of the $SI$ plane for appropriate Dulac function $\Phi$. According to the Bendixson-Dulac theorem [19], we can see that system (10) has no periodic orbits in the interior of the first quadrant. Thus, all solutions of system (10) tend to one of the equilibria. Here, this point is the endemic equilibrium point that only exists when ${\mathcal{R}}_0>1$. Hence $P_E$ is globally asymptotically stable in the interior of the first quadrant.

Consequently, while ${\mathcal{R}}_0>1,$ system (10) has a unique endemic equilibrium $P_E=\left(S^{*},I^{*}\right)$, which is globally asymptotically stable. □

4. Sensitivity Analysis

It is difficult to completely eradicate an epidemic in a population in a short period of time. Considering that lots of negative situations brought about by the disease, seeking and attempts to reduce the spread of the disease have great importance. In this context, with various control measures to be implemented; lowering the ${\mathcal{R}}_0$ value is one of the important goals. Therefore, it has a great importance to examine the effect of parameters on the change of ${\mathcal{R}}_0$ and to apply control measures in this direction.

Now we will focus on the sensitivity analysis of ${\mathcal{R}}_0$. Sensitivity analysis clarifies us about what direction effective each parameter is to disease transmission. To observe whether the parameters that affect the basic reproduction number have a positive or negative effect, we will explore the normalized forward sensitivity index of ${\mathcal{R}}_0.$ The normalized forward sensitivity index of the variable ${\mathcal{R}}_0$ with respect to the parameter $\varsigma$ is defined as

$${\chi }_{\varsigma }^{{\mathcal{R}}_0}={\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \varsigma }}\times {\displaystyle \frac{\varsigma }{{\mathcal{R}}_0}},$$

by using partial derivatives. Where $\varsigma$ represents the basic parameters constituting ${\mathcal{R}}_0$.

Then

$${\chi }_{\beta }^{{\mathcal{R}}_0}={\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \beta }}\times {\displaystyle \frac{\beta }{{\mathcal{R}}_0}}=1>0,$$

$$\begin{array}{ccc}\hfill {\chi }_q^{{\mathcal{R}}_0}& =& {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial q}}\times {\displaystyle \frac{q}{{\mathcal{R}}_0}}\hfill \\ & =& -{\displaystyle \frac{\mu q\left(1+\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)e^{-\mu \theta }d\theta \right)}{\left(q+\mu \right)\left(1+q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)e^{-\mu \theta }d\theta \right)}}<0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\chi }_{\gamma }^{{\mathcal{R}}_0}& =& {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \gamma }}\times {\displaystyle \frac{\gamma }{{\mathcal{R}}_0}}\hfill \\ & =& -{\displaystyle \frac{\gamma }{\gamma +\delta +\mu }}<0\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\chi }_{\delta }^{{\mathcal{R}}_0}& =& {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \delta }}\times {\displaystyle \frac{\delta }{{\mathcal{R}}_0}}\hfill \\ & =& -{\displaystyle \frac{\delta }{\gamma +\delta +\mu }}<0.\hfill \end{array}$$

With increasing of values of the parameters that have positive indices $\left(\beta \right)$, ${\mathcal{R}}_0$ increases and so the spread of the disease progresses in the population. On the other hand, the parameters in which its sensitivity indices are negative $\left(q,\gamma \mathrm{and}\delta \right)$ cause decreasing of ${\mathcal{R}}_0$. That is, the average number of secondary infections cases decreases while these parameters increase and so the spread of the disease starts to decrease.

5. Conclusion

Stability is one of the substantial problems in the design and analysis of control systems. Since an unstable control system is considered useless and potentially dangerous, the most basic matter related to the examined control systems is whether it is stable or not.

Nowadays, with the advancement of science, the desires and efforts of individuals have been increased to solve more complex problems. The dynamics of these complex problems often requires nonlinear equations. Furthermore, although there are many methods to investigate the stability of linear systems, examining the stability of nonlinear systems is an important problem. Because nonlinear systems can exhibit different types of behavior which are not visible in linear systems.

Additionally, the existence of delay in a system may lead to being much more complex of analysis and control of the system. Besides, it may cause instability and poor performance in the system. Therefore, the stability analysis of the systems with delay (especially distributed delay) has crucial important theoretically. However, let us immediately state that determining the stability of time-delayed systems is not as easy as non-delayed systems.

In mathematical epidemic models, the course of the disease in population depends on whether ${\mathcal{R}}_0$ is greater than 1 or not. If ${\mathcal{R}}_0>1,$ spread of the disease increasingly goes on. If ${\mathcal{R}}_0<1,$ there is a decrease in outbreak speed and disease dies out gradually. Also if ${\mathcal{R}}_0=1$, the outbreak speed is stable. Therefore, ${\mathcal{R}}_0$ determines how much the outbreak will be severe in a population, and it can tell us whether the population is at risk in the face of disease. Moreover, it is very important to calculate the value ${\mathcal{R}}_0$ because the size of ${\mathcal{R}}_0$ allows to determine the amount of required effort to prevent an outbreak or to eliminate an infection from a population.

On the other hand, vaccination is one of the most important mechanisms in the struggle against epidemics. Especially in cases where immunity obtained by vaccination may not be permanent, it should also be taken into account that the duration of immunization of vaccinated individuals may not be the same. Even if they were vaccinated at the same time, while immunity of some of the individuals may continue, the others may lose their immunity and become susceptible to the epidemic.

In this study, prepared by taking all these facts into consideration, a mathematical epidemic model reflecting that the protection period provided by the vaccine effect which may vary from person to person, is presented. This new $SVIR$ fractional epidemic model is formed by aid of a system of distributed delay nonlinear fractional integro-differential equations.

Firstly, equilibrium points of formed system are found and the basic reproduction number of the model are determined as

$${\mathcal{R}}_0={\displaystyle \frac{b\beta \left(1+q\underset{0}{\overset{\infty }{\int }}f\left(\theta \right)e^{-\mu \theta }d\theta \right)}{\left(q+\mu \right)\left(\gamma +\delta +\mu \right)}}.$$

Then it is proved that system (10) has a unique disease-free equilibrium point $P_{DF}$ and a unique endemic equilibrium $P_E$, which are locally asymptotically stable when ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_0>1,$ respectively. Moreover, it has been shown that in a population with infectious individuals, the disease will never disappear when ${\mathcal{R}}_0>1$. Thus, it has been reinforced that the introduced model is consistent and ${\mathcal{R}}_0$ is meaningful.

Also, the global behavior of the system (10) is examined. Firstly, it is proved that all solutions of the system tend to $P_{DF}$ by constituting the appropriate Lyapunov function for ${\mathcal{R}}_0<1.$ Then, when ${\mathcal{R}}_0>1$, the proof of global stability of $P_E$ has been presented by means of the Dulac’s criterion and the Poincaré-Bendixson theorem, which are commonly used in two-dimensional population models.