Fractional order SEIR model with generalized incidence rate

1.   Introduction

Mathematical epidemiology is concern with the dynamical evolution of diseases and its control within a population. It is an important field and attracted much interest in the recent years. Mathematical models of various communicable diseases have been used as a powerful tool to explore the realistic aspects of disease spreading. In literature, numerous mathematical models have been developed in order to analyze the spread and possible control strategies of various diseases such as ^[1,2,3].

In epidemic models the incidence rate function plays an essential role and ensure that the model provides a reasonable qualitative analysis of the disease dynamics. In the literature of infectious disease models, a verity of incidence rate functions have been implemented. The bilinear incidence rate of the type $\beta S I$ has been frequently adopted, where $\beta$ indicates the per unit contact rate, $S$ and $I$ are respectively susceptible and infected individuals ^[4,5]. In ^[6] the authors introduced a simple epidemic model using bilinear incidence rate under the assumptions that there are no births or deaths from the infection. Later on in 1973, Capasso et al., ^[7] introduced the nonlinear incidence rate $\frac{\beta SI}{1+\alpha I }$ to formulate the mathematical model of cholera epidemic in Bari. The global analysis of the SIR and SIRS models with nonlinear incidence rate have been investigated in ^[7,8]. In ^[9], the authors investigated the dynamics of SIR epidemic model with the incidence function of the form $f(S, I)$. Recently, Gao et al. ^[10] formulated a more general epidemic model using $f(S, I)$ as transmission function and explored the global dynamics of the model. Most of these models are formulated using ordinary integer-order differential equations (IDEs). In epidemiology, the previous experiences and history of an epidemic have an essential role to study its dynamics in a more realistic way. These classical models have some serious drawbacks such as they are local in nature and can not explore the dynamics of phenomena in between two integer values. Further, since the classical models do not possess the memory effects, therefore, these models can not replicate the dynamics of many real world phenomena including infectious diseases.

Fractional calculus (FC) is the generalization of classical integer-order calculus. Mathematical models with fractional-order (FO) derivative can be used to model universal phenomena with greater degree of accuracy and its applications can be found in various fields such as engineering, economics, control theory, finance and in epidemiology ^{[11,12,13,14,15,16,17]}. The increasing interest of using FDEs in modeling of real world complex problems is due to its various properties which are not found in IDEs. In order to overcome the aforesaid limitations of integer-order derivative, different FO operators have been introduced in existing literature ^[18,19,20]. A number of epidemic models using FO derivative with singular and non-singular Kernel have been proposed in the literature. Most of these models are based on either bilinear or non-linear incidence rate. Such as Saeedian et al. ^[21] proposed a simple SIR epidemic model with bilinear incidence function and explored the importance of memory effects and previous history on the disease dynamics. Further, they concluded that the precise information about the past events plays a key role in the disease eradication. Mouaouine et al. ^[22] developed a non-integer order SIR epidemic model incorporating non-linear incidence rate function.

Motivated by the previous literature in the present paper, we develop a fractional SEIR epidemic model with generalized incidence rate function of the form $f(I) S$. Further, we provide a detail stability analysis of both disease free and endemic equilibriums of the model and numerical simulations for various values of fractional order $\alpha$. The remaining sections of the manuscript are organized as: The basic definitions and relevant results are provided in section 2. Model description and its basic properties are presented in section 3. The stability results of the model are discussed in section 4. The numerical simulations and concluding remarks are given in sections 5 and 6 respectively.

2.   Preliminaries

First we recall the basic definitions and some relevant results regarding the Caputo fractional derivative ^[18,23].

Definition 2.1. For a given function $h\in \mathbb{C}^n$ and $t, \vartheta\in \mathbb{R}$, then the FO derivative having order $\vartheta$ in Caputo sense is given by

where $n-1 < \vartheta < n\in \mathbb{N}$.

Definition 2.2. For a function $h:\mathbb{R}^+\rightarrow \mathbb{R}$, the fractional integral having order $\vartheta$ is given by

Definition 2.3. ^[24] For a given dynamical system with Caputo fractional operator given by

the constant $x^*$ is an equilibrium point if and only if $f(t, x^*) = 0$.

In order to implement Lyapunov stability method for a system involving Caputo derivative, we re-call the relevant result from ^[24,25].

Theorem 2.4. For an equilibrium point given by $x^*$ for the system in Caputo sense (2.1) and $\Omega\in \mathbb{R}^n$ be the domain such that $x^*\in \Omega$ and let $\mathcal{G}:[0, \infty)\times\Omega \rightarrow \mathbb{R}$, be a continuously differentiable function and if

and

$\forall$ $\vartheta\in(0, 1)$ and $x\in\Omega$. Where $V_1(x)$, $V_2(x)$ and $V_3(x)$ are continuously positive definite functions over $\Omega$, then the point $x$ of (2.1) is stable uniformly asymptotically.

3.   Mathematical formulation of the model

To construct the propose model we divide the total population into four subclasses i.e., suspectable $S(t)$, exposed $E(t)$, infected $I(t)$ and those who recovered are denoted by $R(t)$. The suspectable population is recruited at the rate $\Lambda$. The natural death rate in all classes is denoted by $d$. The function $Sf(I)$ denotes the generalized incidence rate function describing the rate at which the suspectable become infectious and join the exposed class. The exposed class become infected and enter to infected class at the rate $\tau$. The parameter $\mu_2$ is the recovery rate of infected class whereas $\delta$ is the disease related death rate of infected class. The waning of immunity of the recovered class is denoted by the parameter $\theta$ and join the suspectable class again. The dynamics in each class is governed by a differential equation. Hence, the proposed fractional SEIR with generalized incidence rate is given by the below nonlinear system of FDEs.

In (3.1), ${}^CD_t^\alpha$ denotes the Caputo derivative having order $\alpha\in(0, 1]$ in order to describe the memory effects in the proposed epidemic model.

To make our study more effective we have assumed that $f(I)$ is to be non-negative and continuously differentiable in the interior of $\mathbb{R}_+$ and further, the same hypotheses is taken in account as mentioned in ^[26] i.e., $f(I)$ is a real locally Lipschitz function on the interval $[0, \infty)$ such that

ⅰ. $f(0) = 0$ and $f(I) > 0$ for $I > 0$;

ⅱ. $f(I)/I$ is continuous and monotonously non-increasing for $I > 0$ and $\lim_{I\rightarrow 0^+}\frac{f(I)}{I}$ exists, denoted by $\beta$ with $\beta > 0$.

It is obvious from condition (ⅱ)

Then, $\int_{0+}^1 1/f(u)du = 1$ and thus the assumption (ⅱ) in ^[26] is redundant. The incidence rate $f(I)$ is dependent on the concentration of infection. Furthermore, it most famous generalized forms satisfying the above hypothesis are given in the Table 1.

3.1.   Existence and positivity of the solution

To present the non-negativity of the system solution, let

To proceeds further, we recall the generalized mean values theorem ^[27].

Lemma 3.1. Let suppose that $h(y) \in C[a_1, a_2]$ and $^CD^\alpha_t h(y) \in (a_1, a_2]$, then

where $a_1\leq\chi\leq t,$ $\forall$ $t\in (a_1, a_2]$.

Corollary 3.2. Suppose that $h(y) \in C[a_1, a_2]$ and $^CD^\alpha_t h(y) \in (a_1, a_2]$, where $\alpha \in (0, 1]$. Then if

Theorem 3.3. A unique solution $y(t)$ of (3.1) exists with a view to remain in $\mathbb{R}^4_+$. Furthermore, the solution is positive.

Proof. The exitance of the Caputo fractional SEIR model can be shown with the help of theorem 3.1 from ^[28], while the uniqueness of the solution can be easily obtained by making use of the remark 3.2 in ^[28] for all positive values of $t$. In order to explore the solution positivity, it is necessary to show that on each hyperplane bounding the non-negative orthant, the vector field points to $\mathbb{R}^4_+$. Utilizing the aforesaid conditions on incidence function $f(I)$, we deduced form the system (3.1)

Hence, using the above corollary, we got the target that is, the solution will stay in $\mathbb{R}^4_+$ and hence, the biologically feasible region is constructed as:

Next we explore the equilibria and basic threshold quantity ${\mathcal{R}}_0$ of the model in the following subsection.

3.2.   Equilibria and basic reproduction number

To evaluate the equilibria of the proposed model (3.1) we need to solve the following linearized system:

Thus, we have

Theorem 3.4. The fractional SEIR model (3.1) have at the most two equilibria which are:

The disease free equilibrium (DEF) given by

The endemic equilibrium (EE) denoted by $E_1 = \Big(S^{*}, E^{*}, I^{*}, R^{*}\Big),$

where,

where $I^{*}$ is a positive zero of the function $H$ defined below

Further,

Hence, $H$ has a positive zero (and only one in this case) if and only if $\frac{d(d+\tau)(d+\mu_1+\mu_2+\delta)}{\beta \tau}-\Lambda < 0$ or equivalently,

The expression for the most biologically important threshold parameter known as the basic reproduction number ${\cal{R}}_0$ obtained by the next generation approach which is given as below:

Where, $\frac{\partial f(I_0)}{\partial I}$ is the partial derivative of incidence function $f$ at DFE. Clearly, the EE $E_{1}$, exist if ${\mathcal{R}}_0 > 1$.

4.   Stability of DEF

The Jacobian matrix $J_{E_{0}}$ of the fractional model SEIR (3.1), evaluated around the DEF, $E_{0}$ is as follow:

Theorem 4.1. For any two positive integers $r_1$ and $r_2$ with $gcd(r_1, r_2) = 1$. Let $\alpha = (\frac{r_1}{r_2})$ and define $M = r_2$, then the model DEF denoted by $E_{0}$ is stable locally asymptotically provided that $|arg(\lambda)| > \frac{\pi}{2M}$, where $\lambda$ denotes the possible roots of the characteristic equation (4.2) of the matrix $J_{E_0}$.

Proof. By expansion of (4.2), we get below equation in term of $\lambda$.

where the coefficients are given below:

The arguments of the roots of $\lambda^{p_1}+d_1 = 0$ are as follow:

In similar pattern, it can be shown that argument of the roots of $\lambda^{p_1}+d_2 = 0$ are also greater than $\frac{\pi}{2M}$. Further, if ${\cal{R}}_{0} < 1,$ then the desired condition ($|arg(\lambda)| > \frac{\pi}{2M}$) is satisfied for all roots of polynomial (4.3). For ${\cal{R}}_{0} > 1,$ then with the help of Descartes rule of signs, there exits exactly one root of characteristic equation with $|arg(\lambda)| < \frac{\pi}{2M}$. Thus the DEF is stable locally asymptotically if ${\cal{R}}_{0} < 1$ and unstable otherwise.

Theorem 4.2. The DEF, $E_0$, of the fractional SEIR model (3.1) is globally asymptotically stable (GAS) within the region $\Phi$ if $\mathcal{R}_0 < \mathit{1}$.

Proof. To present the proof, we consider the following appropriate Lyapunov function:

where $\mathcal{A}_j$, for j = 1, 2, which are positive constant to be decide later. Evaluating the time Caputo fractional derivative of $\mathcal{L}(t)$ we obtain

Utilizing (3.1), we obtain

Let the constants be $\mathcal{A}_1 = \tau$ and $\mathcal{A}_2 = (d+\tau)$, then simplifying, we have,

It is clear that when ${\cal{R}}_0 < 1$ then ${}^CD_t^\alpha \mathcal{L}(t)$ is $-ve$, therefore, we conclude that the disease free case $E_{0}$ is GAS in the region $\Phi$.

Next, we present the global stability of the fractional order $SEIR$ model (3.1) at the endemic case.

4.1.   Global stability of endemic equilibrium

Theorem 4.3. If ${\cal{R}}_0 > 1$, then the EE, $E_1$ of the system (3.1) is GAS.

Proof. Before to start the proof, at the steady state $E_1$ of the model (3.1) we derive the following relations

Now, we define the following Lyapunov function

The differentiation of $\mathcal{V}(t)$ along with the solution of (3.1) is

By direct calculations, we have that:

After some arrangement we have

It follows from the property arithmetic mean and we have

and if

then, by Lyapunov stability theorem, it ensures that the model is GAS at $E_1$ when ${\cal{R}}_0 > 1$.

5.   Numerical results

To obtain the numerical solution of the fractional SEIR model (3.1), we take the general incidence function particularly as bilinear $f(I) = \beta I$ and saturated i.e., $f(I) = \frac{\beta I}{1+\alpha_1 I}$. The numerical values of the parameters are $\Lambda = 0.8$ $\beta = 0.00004$, $\alpha_1 = 0.1$, $d = 0.001$, $\mu_1 = 0.02$, $\theta = 0.02$, $\tau = 0.02$, $\mu_2 = 0.02$, and $\delta = 0.004$. For the numerical solution of fractional order model (3.1), we use the predictor corrector method. The time level considered in the numerical simulations of the model is 400 days. We present the obtained results in the form of graphics as well as Tables. For the Tables 2 and 3, we use the step size $h = 0.1$ while for the graphical results we considered $h = 0.05$. The graphical results of bilinear case are depicted in 1 and 2. We considered $\beta = 0.004$ on which the reproduction number $\mathcal{R}_0 = 67.7249 > 1$ which shows that the population is highly endemic and present the graphical results in Figure 1. Further, we obtain Figure 2, in which we considered $\beta = 0.00004$ on which $\mathcal{R}_0 = 0.67721 < 1$. In Figure 1 and 2 we used the order of the fractional parameter $\alpha = 1, 0.95, 0.9, 0.85$. In Figure 1, the subetaaphs show, the population of susceptible, exposed, infected and recovered individuals. By decreasing the value of $\alpha$, the population of exposed, infected and recovered individuals decreases efficiently while the population of susceptible increases. For such a high endemic case the population of infected compartments and their decrease at fractional order parameter is very important. To make such model with the realistic data could be useful for the data fitting. The sub-graphs in Figure 2 shows that decreasing the fractional order parameter $\alpha = 1, 0.95, 0.9, 0.85$, we can see that the population of infected compartments are increases while the population of susceptible and recovered individuals decreases. The graphical interpretations of the model (3.1) with nonlinear saturated incidence rate are presented in Figures 3 and 4 for $\mathcal{R}_0 > 1$ and $\mathcal{R}_0 < 1$ respectively. Furthermore, for the saturated incidence rate we considered the fractional order parameter $\alpha = 1, 0.95, 0.9, 0.85$ and give the Tables 2 and 3, where each variable for the step-size $h = 0.1$ their values are presented.

6.   Conclusion

The present work investigate the dynamics of a fractional SEIR epidemic model with generalized incidence rate. We presented the fractional SEIR model in Caputo sense and explored its basic mathematical results. The stability results of the disease free and endemic equilibriums are investigated. The fractional model at the disease free case is locally as well as globally asymptotically stable when $\mathcal{R}_0 < 1$. Further, we presented the global stability of the fractional order model at the endemic case when $\mathcal{R}_0 > 1$ by using the extended Lyapunov function theory. Then, the numerical results via Tables and graphically are obtained for the fractional order model when the basic reproduction number less or greater than unity. To obtained the model simulations we particularly considered the most widely used bilinear and situated incidence rates. The graphical result for the high endemic value of the fractional order parameter $\alpha$, suggest that the at high endemic case the infected compartments are decreases well by decreasing the fractional order parameter $\alpha$. From the graphical results we concluded that the saturated incidence rate is more appropriate and biologically feasible than the bilinear case. Further, we concluded that the fractional order model is the generalization of integer order model and it gives useful information at each instant of time of interest. In future we will study the dynamics of the present model using Atangana-Baleanu-Caputo and Caputo-Fabrizio fractional derivatives.

Conflict of interest

All authors declare no conflict of interest.