Analytical and Numerical Boundedness of a Model with Memory Effects for the Spreading of Infectious Diseases

Abstract

In this study, an integer-order rabies model is converted into the fractional-order epidemic model. To this end, the Caputo fractional-order derivatives are plugged in place of the classical derivatives. The positivity and boundedness of the fractional-order mathematical model is investigated by applying Laplace transformation and its inversion. To study the qualitative behavior of the non-integer rabies model, two steady states and the basic reproductive number of the underlying model are worked out. The local and global stability is investigated at both the steady states of the fractional-order epidemic model. After analytic treatment, a structure-preserving numerical template is constructed to numerically solve the fractional-order epidemic model. Moreover, the positivity, boundedness and symmetry of the numerical scheme are examined. Lastly, numerical experiment and simulations are accomplished to substantiate the significant traits of the projected numerical design. Consequences of the study are highlighted in the closing section.

1. Introduction

Rabies is a deadly infectious disease caused by a virus, known as rabies. This virus belongs to the genus Lyssavirus, family Rhabdoviridae, and the order Mononegavirales. The rabies virus mainly affects the brains of the mammals [1,2]. The virus has become a serious risk for more than 150 countries of the globe. It is estimated that 55,000 casualties occur every year due to this deadly virus. The main cause of the rabies virus is the dogs’ bite [2,3]. An individual who becomes infected with the rabies virus and is left untreated dies in a short period of time after the appearance of symptoms [4]. The mathematical models for the rabies disease dynamics contain integer-order differential equations. Dogs, raccoons, bats and foxes are the animals that carry the rabies virus in their saliva and they transmit it to humans via bites [5]. The role of bites by wild animals and human causalities has been studied in various countries of the world [6,7]. Dogs, being the main source of the rabies virus in Africa and Asia, are responsible for its transmission around the globe.

The disease is transmitted to domestic animals and people through contact with tainted saliva. In previous years, bat rabies has arisen as a general medical issue in the United states and Europe [8]. In 2003, a large number individuals in South America passed away from rabies following the scratches of untamed animals, especially bats, then from dogs [9]. Rabies is a global issue, although numerous nations today are considered rabies-free (with no human cases) because of mass inoculation programs. Greece had been rabies from 1987 onward, until it reappeared in 2012 [10]. Rabies is a primary health problem in many populations dense with dogs, specifically in regions where fewer or no preventive measures are adopted (vaccination and remedy) for puppies and human beings. Remedy after exposure to the rabies virus is referred to as post-exposure prophylaxis (PEP), and vaccination earlier than exposure to the infection is known as pre-exposure prophylaxis. Many countries have developed the key concept of introducing the adjoint function to a differential equation by introducing an objective function [11].

The fractional-order derivatives have some important physical features, unlike the integer-order derivatives. One of the main features is the memory effect, which is an important part of biological models. This memory effect contains the history of the past that makes the immune system ready for invaders, such as bacteria, viruses, fungi and much more. Mathematically, a fractional-order derivative has the ability to capture a long range of dynamics of the physical phenomenon by adjusting a suitable value to the fractional-order parameter. The rate of disease dynamics is not same throughout the globe. These rates vary region-wise, depending upon the different conditions of environment, health, food, hygienic and so on. So, the use of fractional-order derivatives becomes inevitable to cover the infinitely many rates of infection according to the real situation. Further, these types of derivatives increase the stability of solutions [12,13]. Shi et al. (2022) studied a stochastic SEIRS rabies model with two patches of dog population. They investigated the existence and uniqueness of the global positive solutions. They also described the criteria for the ergodic stationary distribution [14]. Similarly, Ewald et al. (2020) reviewed various plans to alleviate the pathogens that will lead to the optimal control of the infection. They discussed the current trends in mathematical modeling and highlighted the possible features of the modeling for the human–pathogen interaction in the future [15].

It is not a simple and easy task to find the exact solutions of the nonlinear dynamical systems [16]. In some cases, it becomes more complex to find an exact solution due to the involvement of a number of parameters in the model. Additionally, in some models, the exact solution becomes almost impossible. So, the need for a numerical solution arises. In the literature, many numerical schemes exist for the classical or integer-order models [17,18,19,20,21,22]. On the other hand, a few schemes can be seen for fractional-order epidemic models. These schemes also have some flaws and deficiencies, such as divergence from the exact equilibrium points, negativity of the state variables and lack of stability. By considering these facts, an efficient numerical method is formulated to obtain the reliable numerical solutions. It is determined analytically and numerically that the new scheme converges toward the exact equilibrium state and provides positive solutions.

This paper is organized as follows: The model formulation, mathematical assumptions, mathematical flowchart and model equations are included in Section 2. The model analysis, basic reproduction number, and numerical simulations of the equilibria are mostly detailed in Section 3. We discuss the parameter values that lead to the numeric value of initial sequence ratio in Section 4. The outcomes and effectiveness of the study are discussed in Section 5.

2. Preliminaries

The purpose of the present section is to introduce crucial concepts and present the mathematical model under consideration in this work. To start with, we provide firstly the definition of Caputo fractional derivatives and their properties.

Definition 1

(Caputo [23]). Let m be a positive integer; suppose that $m-1<\lambda \le m$ and let $f:[0,\infty )\to \mathbb{R}$. The Caputo fractional derivative of order λ of a function $\mathit{f}$ is defined as

(1)$\begin{array}{c}\hfill {}_0^c{D}_t^{\lambda }\mathit{f}\left(t\right)=\frac{1}{\Gamma (m-\lambda )}{\int }_{t_0}^t{(t-s)}^{m-\lambda -1}{\mathit{f}}^m\left(s\right)ds,\forall t>0.\end{array}$

Following the notation and conventions in the previous definition, it is well-known that the Laplace transform operator for the Caputo fractional derivative of the function f is given by the following expressions, for each $s\in \mathbb{C}$ as described in [24,25]

(2)$\begin{array}{cc}\hfill L\left\{{}_0^c{D}_t^{\lambda }f\left(t\right)\right\}\left(s\right)& =s^{\lambda }F\left(s\right)-\sum _{k=0}^{n-1}{s}^{k-n-1}{f}^{\left(k\right)\left(0\right)}\hfill \\ \hfill & =\frac{s^{n}F\left(s\right)-s^{n-1}f\left(0\right)-···-f^{n-1}\left(0\right)}{s^{n-\lambda }}.\hfill \end{array}$

Definition 2.

The classical Mittag–Leffler function with a single parameter $\alpha >0$ is the function $E_{\alpha }:\mathbb{C}\to \mathbb{C}$ given in closed form by

(3)$E_{\alpha }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+1)},\forall z\in \mathbb{C}.$

Meanwhile, the Mittag–Leffler function with two parameters $\alpha >0$ and $\beta >0$ is provided by the formula (see Section 2 of [26])

(4)$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},\forall z\in \mathbb{C}.$

Definition 3.

A value $v^{*}$ is a steady state of the system ${}_0^C{D}_t^{\lambda }v\left(t\right)=\mathit{f}(t,v\left(t\right))$, $v(t_0>0)$ if $\mathit{f}(t,v^{*})=0$, for all $t\ge 0$ (see [27,28]).

In our compartmental deterministic SIR model for rabies in human, the human population is divided into susceptible (S), infectious (I), and recovered (R) classes. These classes or compartments are functions that depend only on time. Moreover, this model may be represented as the system of ordinary differential equations [29]:

(5)$\begin{array}{ccc}& & \frac{dS}{dt}=\beta -\alpha (1-\theta )SI-\rho S-m_{1}S,\hfill \end{array}$

(6)$\begin{array}{ccc}& & \frac{dI}{dt}=\alpha (1-\theta )SI-m_{2}I-(\rho +\delta )I,\hfill \end{array}$

(7)$\begin{array}{ccc}\hfill & & \frac{dR}{dt}=m_{1}S+m_{2}I-\rho R.\hfill \end{array}$

The parameters involved in systems (5)–(7) possess concrete physical connotations. For example, $\beta$ is the human recruitment rate, $\delta$ is the death due to rabies, $\alpha$ is the constant rate, $\rho$ is the natural death rate, $\theta$ is the rate of rabies awareness, $m_1$ is the human vaccination rate, and $m_2$ is the rate of rabies treatment. Throughout, we will assume that all of them are positive constants.

For purposes of this work, the classic rabies model (5)–(7) is transformed in the scenario of fractional calculus as

(8)$\begin{array}{ccc}& & {}_0^C{D}_t^{\lambda }S\left(t\right)=\beta -\alpha (1-\theta )SI-\rho S-m_{1}S,\hfill \end{array}$

(9)$\begin{array}{ccc}& & {}_0^C{D}_t^{\lambda }I\left(t\right)=\alpha (1-\theta )SI-m_{2}I-(\rho +\delta )I,\hfill \end{array}$

(10)$\begin{array}{ccc}& & {}_0^C{D}_t^{\lambda }R\left(t\right)=m_{1}S+m_{2}I-\rho R,\hfill \end{array}$

and we will consider initial conditions of the form $S\left(0\right)=S_0\ge 0$, $I\left(0\right)=I_0\ge 0$ and $R\left(0\right)=R_0\ge 0$. We establish now some important properties on the solutions of this Caputo fractional model.

Theorem 1.

For any initial values $(S\left(0\right),I\left(0\right),R\left(0\right))\in {\mathbb{R}}_{+}^3$, the solution $\left(S\right(t),I(t),R(t\left)\right)$ for the system (8)–(10) is positive-invariant in ${\mathbf{R}}_{+}^3$.

Theorem 2.

If $S\left(0\right)\ge 0$, $I\left(0\right)\ge 0$ and $R\left(0\right)\ge 0$, then the solution of systems (8)–(10) is uniformly bounded at all times.

The proofs of Theorems 1 and 2 are given in Appendix A.

3. Analytical Results

The present section is devoted to establishing some properties on the solutions of epidemic model (8)–(10). To start with, it is important to notice that this epidemic model possesses two types of steady-state solutions, namely, a virus-free steady state (VFSS) and a virus-existing steady state (VESS).

To calculate the steady states of model (8)–(10), we assume that the solutions satisfy ${}^{c_0}{D}_t^{\lambda }S\left(t\right){=}^{c_0}{D}_t^{\lambda }I\left(t\right){=}^{c_0}{D}_t^{\lambda }R\left(t\right)=0$, for each $t\ge 0$. After some simplifications, it is possible to obtain the VFSS and the VESS. They are respectively represented by $(S^0,I^0,R^0)$ and $(S^{*},I^{*},R^{*})$, and they are provided by the following formulas:

(11)$\begin{array}{cc}\hfill (S^0,I^0,R^0)& =\left(\frac{\beta }{\rho +m_1},0,\frac{m_{1}S}{\rho }\right),\hfill \end{array}$

(12)$\begin{array}{cc}\hfill (S^{*},I^{*},R^{*})& =\left(\frac{m_2+\rho +\delta }{\alpha (1-\theta )},\frac{\beta S^{*}(\rho +m_1)}{\alpha (1-\theta )S^{*}},\frac{m_1{S}^{*}+m_2{I}^{*}}{\rho }\right).\hfill \end{array}$

Notice that when an infected individual appears among the individuals, then epidemics will exist in the population. Mathematically, this situation implies that the derivative of the infected population with respect to time is positive at the time that the infected individual is present. As a consequence, $I\left[(\alpha (1-\theta )S-m_2-(\rho +\delta ))\right]>0$, which implies that $\alpha (1-\theta )S-m_2-(\rho +\delta )>0$ in view that $T>0$. Equivalently,

$\begin{array}{c}\hfill \frac{\alpha (1-\theta )S}{m_2+\rho +\delta }>1.\end{array}$

At the outset of epidemic, we have that S is approximately equal to $\frac{\beta }{\rho +m_1}$. All this together yields that the basic reproductive number is given by

(13)$R_0=\frac{\beta \alpha (1-\theta )}{(m_2+\rho +\delta ).(\rho +m_1)}$

It is well-known that the basic reproductive number determines whether the disease will die out or persist in the population. Indeed, recall that when $R_0<1$, then the disease will die out (no epidemic). Meanwhile, if $R_0>1$, then the disease will become an epidemic. Moreover, the value of $R_0$ indicates that how contagious the disease is.

Next, we will investigate the stability of the steady states associated with the fractional model (8)–(10). To this end, two theorems are established.

Theorem 3.

The VFSS is locally asymptotically stable if $R_0<1$.

The proof of Theorem 3 is given in Appendix A.

Theorem 4.

If $R_0>1$, then the VESS is locally asymptotically stable.

Proof. 

We proceed in the similar fashion as in the proof of the previous theorem. We calculate the matrix $J_{S,I}$ as in the previous result, and evaluate it now at the VESS. Then the eigenvalues of the matrix $J_{S^{*},I^{*}}$ satisfy the quadratic equation

(14)${\lambda }^2+(\rho +m_1+\beta (\rho +m_1))\lambda +\beta (\rho +m_1)(m_2+\rho +\delta )=0.$

Under the hypothesis, notice that the independent coefficient as well as the coefficient of the linear term in this quadratic expression are positive. By the Routh–Hurwitz criterion, we conclude that the system is stable at the VESS.    □

In order to prove the global stability of fractional system (8)–(10), we present an important lemma, which is as follows.

Lemma 1

([30]). Let $W\left(t\right)\in {\mathbb{R}}^{+}$ be a continuous function. Then, for any time $t\ge t_0$:

$\begin{array}{c}\hfill \begin{array}{cc}& {}_0^c{D}_t^{\theta }\left[W\left(t\right)-W^{*}-W^{*}log\frac{W\left(t\right)}{W^{*}}\right]\le \left(1-\frac{W^{*}}{W\left(t\right)}\right){}_0^c{D}_t^{\theta }W\left(t\right).\hfill \\ & W^{*}\in {\mathbb{R}}^{+},forall\theta \in (0,1).\hfill & \hfill \end{array}\end{array}$

Theorem 5.

If $R_0<1$, then the VFSS is globally asymptotically stable.

Proof. 

We consider the Volterra-type Lyapunov function as

(15)$\begin{array}{ccc}\hfill W& =& \left(S-S^0-S^{0}log\frac{S}{S^0}\right)+I+\left(R-R^0-R^{0}log\frac{R}{R^0}\right)\hfill \end{array}$

After applying the Lemma 1, we have

$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\lambda }W& \le & \left(1-\frac{S^0}{S}\right){}_0^C{D}_t^{\lambda }S+{}_0^C{D}_t^{\lambda }I+\left(1-\frac{R^0}{R}\right){}_0^C{D}_t^{\lambda }R,\hfill \\ \hfill {}_0^C{D}_t^{\lambda }W& \le & \left(1-\frac{S^0}{S}\right)(\beta -\alpha (1-\theta )SI-\rho S-m_{1}S)+\alpha (1-\theta )SI\hfill \\ & -& m_{2}I-(\rho +\delta )I+\left(1-\frac{R^0}{R}\right)(m_{1}S+m_{2}I-\rho R),\hfill \end{array}$

After some computations, we obtain

(16)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_0^C{D}_t^{\lambda }W\le & \frac{-\beta {(S-S^0)}^2}{S{S}^0}-\alpha (1-\theta )(S-S^0)(I-I^0)-(m_2+\rho +\delta )I\hfill \\ & \left(1-\frac{\alpha (1-\theta )S}{(m_2+\rho +\delta )}\right)-\frac{m_{1}S}{R{R}^0}{(R-R^0)}^2-\frac{m_{2}I}{R{R}^0}{(R-R^0)}^2.\hfill \end{array}\end{array}$

From Equation (16), ${}_0^C{D}_t^{\lambda }W\le 0$ for $R_0<1$, and ${}_0^C{D}_t^{\lambda }W=0$ only if $S=S^0$, $I=0,R=R^0$. Therefore, it can be concluded that system (8)–(10) is globally asymptotically stable.    □

Theorem 6.

If $R_0>1$, then the VESS is globally asymptotically stable.

Proof. 

In the same case as the disease-free equilibrium point, we assume to consider the Volterra-type Lyapunov function defined as

(17)$\begin{array}{c}\hfill \begin{array}{cc}\hfill W& =K_1\left(S-S^{*}-S^{*}ln\frac{S}{S^{*}}\right)+K_2\left(I-I^{*}-I^{*}ln\frac{I}{I^{*}}\right)\hfill \\ & +K_3\left(R-R^{*}-R^{*}ln\frac{R}{R^{*}}\right),\hfill \end{array}\end{array}$

where, $K_i,i=1,2,3$ are all positive constants that can be chosen later. Then, by substituting Equation (17) into main model (8)–(10) and using Lemma 1, we reach the following:

$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\lambda }W& \le & K_1\left(1-\frac{S^{*}}{S}\right){}_0^C{D}_t^{\lambda }S+K_2\left(1-\frac{I^{*}}{I}\right){}_0^C{D}_t^{\lambda }I+K_3\left(1-\frac{R^{*}}{R}\right){}_0^C{D}_t^{\lambda }R,\hfill \\ \hfill {}_0^C{D}_t^{\lambda }W& \le & K_1\left(\frac{S-S^{*}}{S}\right)(\beta -\alpha (1-\theta )SI-\rho S-m_{1}S)+K_2\left(\frac{I-I^{*}}{I}\right)\hfill \\ & & [\alpha (1-\theta )SI-m_{2}I-(\mu +\delta )I]+K_3\left(\frac{R-R^{*}}{R}\right)(m_{1}S+m_{2}I-\rho R),\hfill \end{array}$

After some computations, we obtain

(18)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_0^C{D}_t^{\lambda }W& \le -K_1\frac{-\beta {(S-S^{*})}^2}{S{S}^{*}}-K_1\alpha (1-\theta )(S-S^{*})(I-I^{*})-K_3\frac{m_{1}S{(R-R^{*})}^2}{R{R}^{*}}\hfill \\ & -K_3\frac{m_{2}I{(R-R^{*})}^2}{R{R}^{*}}.\hfill \end{array}\end{array}$

Finally, by setting $K_1=K_2=K_3=1$, we have

$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\lambda }W& \le & -\frac{\beta {(S-S^{*})}^2}{S{S}^{*}}-\alpha (1-\theta )(S-S^{*})(I-I^{*})-\frac{m_{1}S{(R-R^{*})}^2}{R{R}^{*}}-\frac{m_{2}I{(R-R^{*})}^2}{R{R}^{*}}\le 0\hfill \end{array}$

where ${}_0^C{D}_t^{\lambda }W\le 0$ for $R_0>1$, and ${}_0^C{D}_t^{\lambda }W=0$ only if $S=S^{*}$, $I=I^{*},R=R^{*}$. Therefore, it is concluded that system (8)–(10) is globally asymptotically stable.    □

4. Numerical Results

This section is devoted to present a reliable structure-preserving numerical method for the solution of the continuous model (8)–(10). To that end, let $t>0$ and fix the time interval $(0,t)$. Let us divide this temporal interval into $n\in \mathbb{N}$ identical sub-intervals with grid points given by

(19)$0=t_0<t_1<···<t_{n+1}=t=(n+1)h.$

Here, $h=t_{j+1}-t_j$, for each $j=1,2,\dots ,n$. For the sake of convenience, let us convey that $S_j=S\left(t_j\right)$, $I_j=I\left(t_j\right)$ and $R_j=R\left(t_j\right)$, for each $j=0,1,\dots ,n$.

Definition 4.

Let Z be any of the functions S, I or R, and assume that $j=0,1,\dots ,n-1$. Then the nonstandard Grünwald–Letnikov approximation for the Caputo fractional derivative of Z of order λ at the time $t_{n+1}$ is given by

(20) ${}_0^C{\Delta }_t^{\lambda }Z\left(t_{n+1}\right)=\frac{1}{\psi {\left(h\right)}^{\lambda }}\left(Z\left(t_{n+1}\right)-\sum _{u=1}^{n+1}{c}_u^{\lambda }Z\left(t_{n+1-u}\right)-{\gamma }_{n+1}^{\lambda }Z\left(t_0\right)\right).$

Here, $\psi :\mathbb{R}\to \mathbb{R}$ is a suitable function, $c_1^{\lambda }=\lambda$ and

(21) $c_v^{\lambda }=\left(1-\frac{\lambda +1}{v}\right)c_{v-1}^{\lambda },\forall v=2,3,\dots$

If $0<\lambda <1$, then it is easy to check that the coefficients $c_u^{\lambda }$ and ${\gamma }_u^{\lambda }$ satisfy the following chains of inequalities for each $u\ge 1$:

• $0<c_{u+1}^{\lambda }<c_u^{\lambda }<\dots <c_1^{\lambda }=\lambda <1$.

• $0<{\gamma }_{u+1}^{\lambda }<{\gamma }_u^{\lambda }<\dots <{\gamma }_1^{\lambda }=\frac{1}{\Gamma (1-\lambda )}$.

Moreover, if $\lambda \in (0,1)$, then ${lim}_{u\to \infty }{c}_u^{\lambda }={lim}_{u\to \infty }{\gamma }_u^{\lambda }=0$ and ${\sum }_{u=1}^{\infty }{c}_u^{\lambda }=1$.

Using these conventions, the Grünwald–Letnikov nonstandard finite-difference (GL-NSFD) method employed to approximate the solutions of the epidemic model (8)–(10) is given by the system of discrete equations, for each $n=0,1,2,\dots$:

(22)$\begin{array}{cc}\hfill {}_0^C{\Delta }_t^{\lambda }{S}_{n+1}& =\beta -\alpha (1-\theta )S_{n+1}{I}_n-\rho S_{n+1}-m_1{S}_{n+1},\hfill \end{array}$

(23)$\begin{array}{cc}\hfill {}_0^C{\Delta }_t^{\lambda }{I}_{n+1}& =\alpha (1-\theta )S_{n+1}{I}_n-m_2{I}_{n+1}-(\rho +\delta )I_{n+1},\hfill \end{array}$

(24)$\begin{array}{cc}\hfill {}_0^C{\Delta }_t^{\lambda }{R}_{n+1}& =m_1{S}_{n+1}+m_2{I}_{n+1}-\rho R_{n+1}.\hfill \end{array}$

By using Formula (21) in the system of Equations (22)–(24), we obtain

(25)$\begin{array}{c}\hfill \frac{1}{\psi {\left(h\right)}^{\lambda }}(S_{n+1}-\sum _{u=1}^{n+1}{c}_{\mu }^{\lambda }{S}_{n+1-u}-{\gamma }_{n+1}^{\lambda }{S}_0)=\beta -\alpha (1-\theta )S_{n+1}{I}_n-\rho S_{n+1}-m_1{S}_{n+1},\end{array}$

(26)$\begin{array}{c}\hfill \frac{I}{\psi {\left(h\right)}^{\lambda }}(I_{n+1}-\sum _{u=1}^{n+1}{c}_{\mu }^{\lambda }{I}_{n+1-u}-{\gamma }_{n+1}^{\lambda }{I}_0)=\alpha (1-\theta )S_{n+1}{I}_n-m_2{I}_{n+1}-(\rho +\delta )I_{n+1},\end{array}$

(27)$\begin{array}{c}\hfill \frac{1}{\psi {\left(h\right)}^{\lambda }}(R_{n+1}-\sum _{u=1}^{n+1}{c}_{\mu }^{\lambda }{R}_{n+1-u}-{\gamma }_{n+1}^{\lambda }{R}_0)=m_1{S}_{n+1}+m_2{I}_{n+1}-\rho R_{n+1}.\end{array}$

After some algebraic steps, it is possible to check that this system is equivalent to

(28)$\begin{array}{cc}\hfill S_{n+1}& =\frac{\psi {\left(h\right)}^{\lambda }\beta +{\sum }_{u=1}^{n+1}{c}_u^{\lambda }{S}_{n+1-u}+{\gamma }_{n+1}^{\lambda }{S}_0}{(1+\alpha (1-\theta )I_n+\rho +m_1)\psi {\left(h\right)}^{\lambda }},\hfill \end{array}$

(29)$\begin{array}{cc}\hfill I_{n+1}& =\frac{\psi {\left(h\right)}^{\lambda }\alpha (1-\theta )S_{n+1}{I}_n+{\sum }_{u=1}^{n+1}{c}_u^{\lambda }{I}_{n+1-u}-{\gamma }_{n+1}^{\lambda }{I}_0}{1+(m_2+\left({\rho }^{\lambda +\delta }\right))\psi {\left(h\right)}^{\lambda }},\hfill \end{array}$

(30)$\begin{array}{cc}\hfill R_{n+1}& =\frac{\psi {\left(h\right)}^{\lambda }{m}_1{S}_{n+1}+\psi \left(h\right)m_2{I}_{n+1}+{\sum }_{u=1}^{n+1}{c}_u^{\lambda }{R}_{n+1-u}+{\gamma }_{n+1}{R}_0}{1+\psi {\left(h\right)}^{\lambda }\rho }.\hfill \end{array}$

The next results summarize the capability of the numerical method to preserve the positivity and the boundedness of approximations.

Theorem 7.

Suppose that $S_0>0$,$I_0>0$,$R_0>0$ and $\psi {\left(h\right)}^{\lambda }>0$. Then $S_n>0$,$I_n>0$ and $R_n>0$, for each $n=0,1\dots$

Proof. 

The conclusion is true for $n=0$ by hypothesis. So, assume that the conclusion of this result is satisfied. Using Formulas (28)–(30), the induction hypothesis and the fact that all the coefficients are positive, it obviously follows that $S_{n+1}>0$, $I_{n+1}>0$ and $R_{n+1}>0$. The conclusion of this result follows now by mathematical induction.    □

Theorem 8.

Suppose that $S_0$, $I_0$ and $R_0$ are positive, and let $\psi {\left(h\right)}^{\lambda }>0$. If $\alpha =0$ or $\theta =0$, then there are constants $N_{n,\lambda }>0$ for each $n=0,1,2,\dots$, such that

(31) $max\{S_n,I_n,R_n\}\le N_{n+1,\lambda },\forall n=0,1,2,\dots$

Proof. 

For illustration purposes, we proceed recursively, though the proof could be carried out also using mathematical induction. Beforehand, add all three Equations (28)–(30) to obtain the following identity:

(32)$\begin{array}{cc}\hfill & (1+\rho )\psi {\left(h\right)}^{\lambda }{S}_{n+1}+(1+(\rho +\delta )\psi {\left(h\right)}^{\lambda })I_{n+1}+(1+\psi {\left(h\right)}^{\lambda }\rho )R_{n+1}\hfill \\ \hfill & =\sum _{u=1}^{n+1}{c}_u^{\lambda }(S_{n+1-u}+I_{n+1-u}+R_{n+1-u})+{\gamma }_{n+1}^{\lambda }(S_0+I_0+R_0)+\beta \psi {\left(h\right)}^{\lambda }.\hfill \end{array}$

Notice that $N_{0,\lambda }=max\{S_0,I_0,R_0\}$ obviously satisfies the conclusion of the theorem. Letting now $n=0$, it is easy to check that

(33)$(1+\rho )\psi {\left(h\right)}^{\lambda }{S}_1+(1+(\rho +\delta )\psi {\left(h\right)}^{\lambda })I_1+(1+\psi {\left(h\right)}^{\lambda }\rho )R_1<N_{1,\lambda }^{*},$

where

(34)$N_{1,\lambda }^{*}=(S_0+I_0+R_0)\left(\lambda +\frac{1}{\Gamma (1-\lambda )}\right)+\beta \psi {\left(h\right)}^{\lambda }.$

It is obvious then that $S_1\le N_{1,\lambda }$, $I_1\le N_{1,\lambda }$ and $R_1\le N_{1,\lambda }$, where

(35)$N_{1,\lambda }=N_{1,\lambda }^{*}max\left\{\frac{1}{(1+\rho )\psi {\left(h\right)}^{\lambda }},\frac{1}{1+(\rho +\delta )\psi {\left(h\right)}^{\lambda }},\frac{1}{1+\psi {\left(h\right)}^{\lambda }\rho }\right\}.$

The conclusion can be reached then recursively.    □

The pseudocode of proposed numerical technique (28)–(30) is given in Algorithm 1.

Algorithm 1: Pseudocode of proposed scheme (28)–(30).

procedureFunction$(S,I,R)$     Get$\left(\lambda \right)$     INIT → $\beta ,\alpha ,\theta ,\delta ,\rho ,m_1,m_2$     SET → t     Initial values $\leftarrow S,I,R$     for each n until length(t)+1 do           ${\gamma }_n\leftarrow h^{\lambda ,\Gamma }$     end for     repeat     $n\leftarrow$ length(t)-1     for $SIR$ do           for each u until n do                $c_1\to \lambda$                if $u>1$ then                    $c\to {\lambda }^u$                end if                $\sum c^u\times {(S,I,R)}^{n+1-u}$                Update $S,I,R$           end for           ${[S,I,R]}^{n+1}=$ right hand side of (28)–(30)     end for end procedure

Finally, we present some numerical simulations which assess the capability of the numerical scheme to preserve the positivity and the boundedness. To that end, we employ the parameters summarized in Table 1. Throughout, the initial conditions are given as $S\left(0\right)=0.09$, $I\left(0\right)=0.02$, $R\left(0\right)=0.02$. The results of our numerical simulations for both the VFSS and the VESS are provided in Figure 1 for various values of $\lambda \in (0,1)$. The results show that the scheme is capable of preserving the positivity and the boundedness of the approximations. Moreover, they confirm the capability of the numerical model to preserve the stability of the steady-state solutions of the epidemiological system.

5. Conclusions

In this article, a classical rabies model is transformed into the fractional-order rabies model. As the fractional differential equation represents the class of differential equations by letting the different values of fractional-order parameter. So, the disease dynamics be understood deeply and concisely. The positivity and the boundedness properties of the fractional rabies model are explored with the help of the Laplace theory. Two steady states of the underlying model, i.e., the virus-free and virus-existing steady states, are worked out. The basic reproductive number of the fractional system is depicted. It is investigated that the non-integer order system is locally and globally stable at the virus-free steady state when $R_0<1$, and for $R_0>1$ the system is locally and globally stable at the virus-existence steady state. For the numerical analysis, the robust numerical template is constructed. The significant features of this numerical design are the positivity, boundedness and convergence toward accurate steady states. These traits of the numerical design are identified by establishing some standard results. Moreover, simulations are presented to validate all of the key features of the novel numerical design.

Author Contributions

Conceptualization, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; methodology, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; software, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; validation, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; formal analysis, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; investigation, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; resources, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; data curation, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; writing—original draft preparation, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; writing—review and editing, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; visualization, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; supervision, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; project administration, Z.I., J.E.M.-D., N.A., A.J., M.R. and A.R.; funding acquisition, J.E.M.-D. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data will be available from the corresponding author (J.E.M.-D. upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figure and Table

Figure 1. Numerical simulations of the solutions for the continuous epidemic model investigated in this work, for the VFSS (left column) and the VESS (right column). The graphs present the dynamics for [Forumla omitted. See PDF.] (top row), [Forumla omitted. See PDF.] (middle row) and [Forumla omitted. See PDF.] (bottom row). Various values of [Forumla omitted. See PDF.] were employed in each case (see the legends). (a) [Forumla omitted. See PDF.] for VFSS, (b) [Forumla omitted. See PDF.] for VESS, (c) [Forumla omitted. See PDF.] for VFSS, (d) [Forumla omitted. See PDF.] for VESS, (e) [Forumla omitted. See PDF.] for VFSS and (f) [Forumla omitted. See PDF.] for VESS.

Table 1

Parameter values used to simulate the dynamics of model (8)–(10).

Parameters	Values (VFSS)	Values (VESS)
$\beta$	$0.5$	$0.5$
$\delta$	$0.01$	$0.01$
$\rho$	$0.25$	$0.25$
$m_1$	$0.25$	$0.25$
$m_2$	$0.001$	$0.001$
$\theta$	$0.0001$	$0.0001$
$\alpha$	$0.2$	$0.3$

Appendix A

Theorem A1.

For any initial values $(S\left(0\right),I\left(0\right),R\left(0\right))\in {\mathbb{R}}_{+}^3$, the solution $\left(S\right(t),I(t),R(t\left)\right)$ for the system (8)–(10) is positive-invariant in ${\mathbf{R}}_{+}^3$.

Proof. 

Let us define the norm function $f_{\infty }={sup}_{t\in D_f}\left|f\left(t\right)\right|$, and consider Equation (8). Bounding from below, it is easy to check that the following inequalities are satisfied: (A1)$$\begin{array}{cc}\hfill {}_0^C{D}_t^{\lambda }S\left(t\right)& =\beta -\alpha (1-\theta )SI-\rho S-m_{1}S\ge -\left(\alpha \right|I|+\rho +m_1)S\hfill \\ \hfill & \ge -(\alpha \underset{t\in D_f}{sup}\left|I\right|+\rho +m_1)S\ge -(\alpha I_{\infty }+\rho +m_1)S\ge -K_{1}S,\forall t\ge 0.\hfill \end{array}$$ Here, $K_1=\alpha I_{\infty }+\rho +m_1$. As a consequence, ${}_0^C{D}_t^{\lambda }S\left(t\right)+K_{1}S\ge 0$, for each $t\ge 0$. Applying the Laplace operator on both sides of this inequality, using properties for this transform and rearranging terms to solve for $L\left\{S\right(t\left)\right\}$, we readily obtain (A2)$$L\left\{S\left(t\right)\right\}\ge \frac{S\left(0\right)S^{\lambda -1}}{S^{\lambda }+K_1}.$$ After applying the inverse Laplace operator on (A2), we reach that $S\left(t\right)\ge 0$. Similarly, it is possible to prove that $I\left(t\right)\ge 0$ and $R\left(t\right)\ge 0$, for each $t\ge 0$. □

Theorem A2.

If $S\left(0\right)\ge 0$, $I\left(0\right)\ge 0$ and $R\left(0\right)\ge 0$, then the solution of the system (8)–(10) is uniformly bounded at all times.

Proof. 

Consider the system of equations for the fractional-order rabies model, and define $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$, for each $t\ge 0$. Adding the three fractional differential equations of the mathematical model and rearranging terms, we obtain that (A3)$${}_0^C{D}_t^{\lambda }N\left(t\right)+\rho S\left(t\right)+(\rho +\delta )I\left(t\right)+\left(\rho \right)R\left(t\right)=\beta .$$ Applying Laplace transform on both sides of this equation, using the formula for the Laplace transform of Caputo derivatives and rearranging terms yields (A4)$$L\left\{N\left(t\right)\right\}+\frac{\delta }{S^{\lambda }+\rho }L\left\{I\left(t\right)\right\}=\frac{\beta S^{-1}}{S^{\lambda }+\rho }+\frac{S^{\lambda -1}}{S^{\lambda }+\rho }N\left(0\right).$$ Let $M=max\{N_0,\frac{\beta }{\rho }\}$. Applying now the inverse Laplace transform, we readily obtain (A5)$$\begin{array}{cc}\hfill N\left(t\right)+I\left(t\right)\delta t^{\lambda -1}{E}_{\lambda ,\lambda }(-\rho t^{\lambda })& =N_0{E}_{\lambda ,1}(-\rho t^{\lambda })+\frac{\beta }{\rho }·\rho t^{\lambda }{E}_{\lambda ,\lambda +1}+(-\rho t^{\lambda })\hfill \\ \hfill & \le M\{E_{\lambda ,1}(-\rho t^{\lambda })+\rho t^{\lambda }{E}_{\lambda ,\lambda +1}(-\rho t^{\lambda })\}.\hfill \end{array}$$ As a consequence, we obtain that $N\left(t\right)\le M$. This implies that $N\left(t\right)$ is bounded by M, which means that the solution is bounded in ${\mathbb{R}}_{+}^3$, for all $t\in [0,\infty )$. □

Theorem A3.

The VFSS is locally asymptotically stable if $R_0<1$.

Proof. 

Beforehand, observe that the function R does not appear in the expressions of the first two equations of system (8)–(10). So, for the analysis, we may consider only the first two equations of the epidemic model. Let us define then f and g as the following two functions of S and I: (A6)$$\begin{array}{cc}\hfill f(S,I)& =\beta -\alpha (1-\theta )SI-\rho S-m_{1}S,\hfill \end{array}$$ (A7)$$\begin{array}{cc}\hfill g(S,I)& =\alpha (1-\theta )SI-m_{2}I-(\rho +\delta )I.\hfill \end{array}$$ Notice that the Jacobian matrix of this system of equations is given by (A8)$$J_{S,I}=\left(\begin{array}{cc}{f}_S& f_I\\ g_S& g_I\end{array}\right)=\left(\begin{array}{cc}-\alpha (1-\theta )I-(\rho +m_1)& -\alpha (1-\theta )S\\ \alpha (1-\theta )I& \alpha (1-\theta )S-(m_2+(\rho +\delta ))\end{array}\right).$$ Evaluating now the Jacobian matrix on the VFSS, we readily obtain (A9)$$J_{S^0,I^0}=\left(\begin{array}{cc}-(\rho +m_1)& -\alpha (1-\theta )\left(\frac{\beta }{\rho +m_1}\right)\\ 0& \alpha (1-\theta )\left(\frac{\beta }{\rho +m_1}\right)-(m_2+(\rho +\delta ))\end{array}\right).$$ After some algebraic calculations, we may readily check that the eigenvalues of the Jacobian matrix are the numbers (A10)$$\begin{array}{cc}\hfill {\lambda }_1& =-(\rho +m_1),\hfill \end{array}$$ (A11)$$\begin{array}{cc}\hfill {\lambda }_2& =\alpha (1-\theta )\frac{\beta }{\rho +m_1}-(m_2+(\rho +\delta )).\hfill \end{array}$$ The number ${\lambda }_1$ is obviously negative. Meanwhile, ${\lambda }_2$ is likewise negative in view that $R_0<1$. We conclude then that the VFSS is locally asymptotically stable, as desired. □

References

1. Hayman, D.T.S.; Johnson, N.; Horton, D.L.; Hedge, J.; Wakeley, P.R.; Banyard, A.C.; Zhang, S.; Alhassan, A.; Fooks, A.R. Evolutionary history of rabies in Ghana. PLoS Neglected Trop. Dis.; 2011; 5, e1001. [DOI: https://dx.doi.org/10.1371/journal.pntd.0001001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21483707]

2. Kiddy, J.K.; Asamoah, K.; Oduro, F.T.; Bonyah, E.; Seidu, B. Modelling of Rabies Transmission Dynamics Using Optimal Control Analysis. J. Appl. Math.; 2017; 2017, 2451237.

3. Ega, T.T.; Luboobi, L.S.; Kuznetsov, D. Modeling the Dynamics of Rabies Transmission with Vaccination and Stability Analysis. Appl. Comput. Math.; 2015; 4, pp. 409-419.

4. Demirci, E. A New Mathematical Approach for Rabies Endemy. Appl. Math. Sci.; 2014; 8, pp. 59-67. [DOI: https://dx.doi.org/10.12988/ams.2014.39525]

5. Balcha, C.; Abdela, N. Review of Rabies Preventions and Control. Int. J. Public Health Sci.; 2017; 6, pp. 343-350. [DOI: https://dx.doi.org/10.11591/ijphs.v6i4.10781]

6. Keeling, M.J.; Rohani, P. Modeling Infectious Diseases in Humans and Animals; Princeton University Press: Princeton, NJ, USA, 2008.

7. Bartoszynki, R. On the Risk of Rabies; Mathematical Institute, Polish Academy of Sciences: Warsaw, Poland, 1975; pp. 355-377.

8. Zarin, R.; Ahmed, I.; Kumam, P.; Zeb, A.; Din, A. Fractional modeling and optimal control analysis of rabies virus under the convex incidence rate. Results Phys.; 2021; 28, 104665. [DOI: https://dx.doi.org/10.1016/j.rinp.2021.104665]

9. Wiraningsih, E.D.; Widodo,; Aryati, L.; Toaha, S.; Lenhart, S. Optimal Control for SEIR Rabies Model between Dogs and Human with Vaccination Effect in dogs. Proceedings of the 6th IMT-GT Conference on Mathematics, Statistics and Its Applications (ICMSA2010), Universiti Tunku Abdul Rahman; Kuala Lumpur, Malaysia, 3–4 November 2018.

10. Tsiodras, S.; Dougas, G.; Baka, A.; Billinis, C.; Doudounakis, S.; Balaska, A.; Georgakopoulou, T.; Rigakos, G.; Kontos, V.; Tasioudi, K.E. et al. Re-emergence of animal rabies in northern Greece and subsequent human exposure. Eurosurveillance; 2013; 18, 20474.

11. Warrel, M.J. Current rabies vaccines and prophylaxis schedules: Preventing rabies before and after exposure. Travel Med. Infect. Dis.; 2012; 10, pp. 1-15. [DOI: https://dx.doi.org/10.1016/j.tmaid.2011.12.005]

12. Arafa, A.A.M.; Rida, S.Z.; Khalil, M. A fractional-order model of HIV infection: Numerical solution and comparisons with data of patients. Int. J. Biomath.; 2014; 7, 1450036. [DOI: https://dx.doi.org/10.1142/S1793524514500363]

13. Rihan, F.A.; Lakshmanan, S.; Hashish, A.H.; Rakkiyappan, R.; Ahmed, E. Fractional-order delayed predator–prey systems with Holling type-II functional response. Nonlinear Dyn.; 2015; 80, pp. 777-789. [DOI: https://dx.doi.org/10.1007/s11071-015-1905-8]

14. Shi, Z.; Jiang, D.; Zhang, X.; Alsaedi, A. A stochastic SEIRS rabies model with population dispersal: Stationary distribution and probability density function. Appl. Math. Comput.; 2022; 427, 127189. [DOI: https://dx.doi.org/10.1016/j.amc.2022.127189]

15. Ewald, J.; Sieber, P.; Garde, R.; Lang, S.N.; Schuster, S.; Ibrahim, B. Trends in mathematical modeling of host–pathogen interactions. Cell. Mol. Life Sci.; 2020; 77, pp. 467-480. [DOI: https://dx.doi.org/10.1007/s00018-019-03382-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31776589]

16. Fernández, F.M. On some approximate methods for nonlinear models. Appl. Math. Comput.; 2009; 215, pp. 168-174. [DOI: https://dx.doi.org/10.1016/j.amc.2009.04.060]

17. Dayan, F.; Rafiq, M.; Ahmed, N.; Raza, A.; Ahmad, M.O. A dynamical study of a fuzzy epidemic model of Mosquito-Borne Disease. Comput. Biol. Med.; 2022; 148, 105673. [DOI: https://dx.doi.org/10.1016/j.compbiomed.2022.105673] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35803748]

18. Ahmed, N.; Macías-Díaz, J.E.; Raza, A.; Baleanu, D.; Rafiq, M.; Iqbal, Z.; Ahmad, M.O. Design, Analysis and Comparison of a Nonstandard Computational Method for the Solution of a General Stochastic Fractional Epidemic Model. Axioms; 2022; 11, 10. [DOI: https://dx.doi.org/10.3390/axioms11010010]

19. Macías-Díaz, J.E. An explicit dissipation-preserving method for Riesz space-fractional nonlinear wave equations in multiple dimensions. Commun. Nonlinear Sci. Numer. Simul.; 2018; 59, pp. 67-87. [DOI: https://dx.doi.org/10.1016/j.cnsns.2017.10.019]

20. Macías-Díaz, J.E. Numerical simulation of the nonlinear dynamics of harmonically driven Riesz-fractional extensions of the Fermi–Pasta–Ulam chains. Commun. Nonlinear Sci. Numer. Simul.; 2018; 55, pp. 248-264. [DOI: https://dx.doi.org/10.1016/j.cnsns.2017.07.012]

21. Macías-Díaz, J.E. Numerical study of the transmission of energy in discrete arrays of sine-Gordon equations in two space dimensions. Phys. Rev. E; 2008; 77, 016602. [DOI: https://dx.doi.org/10.1103/PhysRevE.77.016602]

22. Macías-Díaz, J.E.; Ahmed, N.; Rafiq, M. Analysis and Nonstandard Numerical Design of a Discrete Three-Dimensional Hepatitis B Epidemic Model. Mathematics; 2019; 7, 1157. [DOI: https://dx.doi.org/10.3390/math7121157]

23. Caputo, M. Linear model of dissipation whose q is almost frequency independent. II. Geophys. J. Int.; 1967; 13, pp. 529-539. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1967.tb02303.x]

24. Mittag-Leffler, M.G. Sur la nouvelle fonction E_α(x). Comptes Rendus Acad. Sci. Paris; 1903; 137, pp. 554-558.

25. Mittag-Leffler, M.G. Sur la representation analytique d’une branche uniforme d’une fonction monogene. Acta Math.; 1905; 29, pp. 101-181. [DOI: https://dx.doi.org/10.1007/BF02403200]

26. Podlubny, I. Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Mathematics in Science and Engineering; Book Section 2 Elsevier: Amsterdam, The Netherlands, 1999; Volume 198, pp. 41-119.

27. Li, Y.; Chen, Y.; Podlubny, I. Mittag–Leffler stability of fractional order nonlinear dynamic systems. Automatica; 2009; 45, pp. 1965-1969. [DOI: https://dx.doi.org/10.1016/j.automatica.2009.04.003]

28. Biala, T.A.; Khaliq, A.Q.M. A fractional-order compartmental model for the spread of the COVID-19 pandemic. Commun. Nonlinear Sci. Numer. Simul.; 2021; 98, 105764. [DOI: https://dx.doi.org/10.1016/j.cnsns.2021.105764] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33746459]

29. Jawaz, M.; Rehman, M.A.; Ahmed, N.; Baleanu, D.; Iqbal, M.S.; Rafiq, M.; Raza, A. Analysis and numerical effects of time-delayed rabies epidemic model with diffusion. Int. J. Nonlinear Sci. Numer. Simul.; 2022; [DOI: https://dx.doi.org/10.1515/ijnsns-2021-0233]

30. Leon, C.V. Volterra Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul.; 2015; 24, pp. 75-85. [DOI: https://dx.doi.org/10.1016/j.cnsns.2014.12.013]