*Article* **A Comparative Numerical Study and Stability Analysis for a Fractional-Order SIR Model of Childhood Diseases**

**Mohamed M. Mousa 1,2 and Fahad Alsharari 1,3,\***


**Abstract:** The objective of this work is to examine the dynamics of a fractional-order susceptibleinfectious-recovered (SIR) model that simulate epidemiological diseases such as childhood diseases. An effective numerical scheme based on Grünwald–Letnikov fractional derivative is suggested to solve the considered model. A stability analysis is performed to qualitatively examine the dynamics of the SIR model. The reliability and robustness of the proposed scheme is demonstrated by comparing obtained results with results obtained from a fourth order Runge–Kutta built-in Maple syntax when considering derivatives of integer order. Graphical illustrations of the numerical results are given. The inaccuracy of some results presented in two studies exist in the literature have been clearly explained. Generalizing of the cases examined in another study, by considering a model with fraction-order derivatives, is another objective of this work as well.

**Keywords:** SIR model; fractional derivatives; Grünwald–Letnikov method; stability analysis

#### **1. Introduction**

The investigation of the prevalence of the contagious diseases with the aid of mathematical modeling has turn into an essential means to realize epidemiological patterns diseases. A lot of procedures can be studied in the formulation of the mathematical model to reflect the characteristics in spread types of a disease. Generally, there are two main types:


Childhood infections are the highly popular form of infections contagious. These infections diseases such as chicken pox, mumps, measles, etc., to which infants start their life susceptible, and mostly develop within next 5 years. Because children are very close together with their peers, such diseases can prevalence in a fast way. Because the vaccination is believed to be the best efficient approach counter to such diseases, the progress of a framework that would develop the optimum vaccine treatment stage needed to avoid the prevalence of childhood diseases is necessary. The (susceptible–infected– recovered) SIR models are traditional models that have been used to illustrate many epidemiological diseases [1–3,10–14].

The considered version of SIR model is the model given in [1–3]. This model that applied and used to achieve a well understanding of exactly how the childhood disease spreads within populations of various individuals with time, as well as the risk of surges in the susceptible individuals. The SIR is a dynamical system that is consists of three

**Citation:** Mousa, M.M.; Alsharari, F. A Comparative Numerical Study and Stability Analysis for a Fractional-Order SIR Model of Childhood Diseases. *Mathematics* **2021**, *9*, 2847. https://doi.org/ 10.3390/math9222847

Academic Editors: Mihaela Neamt,u, Eva Kaslik and Anca Rădulescu

Received: 30 September 2021 Accepted: 2 November 2021 Published: 10 November 2021

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

coupled ordinary differential equations (ODEs) that refer to the time evolution of the next three individuals:


The normalized dynamical system that describes the SIR model which considered in [3] is given below:

$$\begin{array}{l} D\_t S(t) = (1 - p)\omega - \omega \ S(t) - \beta \ S(t) \ I(t), \\ D\_t I(t) = -(\sigma + \omega) \ I(t) + \beta \ S(t) \ I(t), \\ D\_t \mathcal{R}(t) = p\,\,\omega - \omega \ \mathcal{R}(t) + \sigma \ I(t), \end{array} \tag{1}$$

where *D* = *d*/*dt*, *β* and *σ* are the natural birth rate and the natural death rate, respectively; however, the childhood disease mortality rate is considered extremely low. The factor *p* is a fraction of population that have been vaccinated at the birth every year (0 ≤ *p* < 1) with taking into consideration that the remainder of the citizens is susceptible. The factor *ω* describes the rate of infection of the susceptible citizens due to a contact with infected peoples. Here, the total population is normalized. Recently, a big attention has been drawn to examine mathematical models that defined by fractional differential systems specially in the field of epidemiological diseases [19–21]. The greatest important property of such models is their memory impact [22], which is not occurring in the differential systems of integer derivatives. This memory impact is quite described by the mobility of the differentiation order for the fractional derivatives and realized as an inherited property on viruses' genomes and strains (as mentioned in [21]), that is helpful for modeling epidemiological diseases. It is known that the SIR model is highly impacted by the initial data and the integer order SIR model may not perfectly explain the spread of the diseases because of the local nature of the integer order derivatives. The fractional-order derivative is non-local in its nature and depends on the initial data. Thus, for well understanding of the epidemiological disease's models, it is helpful to replace the integer order SIR model by a fractional-order one. In the present study, we will replace time derivatives in the system (1) with fractional time derivatives of order *γ*. In the last few years, many researchers simulated and analyzed various epidemiological diseases using fractionalorder dynamical models [23–27]. In this work we generally analyze and simulate the dynamics of the fractional-order SIR model that has been studied in [1,2],

$$\begin{array}{ll} \stackrel{\gets}{D}\_{t\_0}^{\gamma} D\_t^{\gamma} S(t) = (1 - p)\omega - \omega \, S(t) - \beta \, S(t) \, I(t),\\ \stackrel{\leftarrow}{D}\_{t\_0}^{\gamma} D\_t^{\gamma} I(t) = -(\sigma + \omega) \, I(t) + \beta \, S(t) \, I(t),\\ \stackrel{\leftarrow}{D}\_{t\_0}^{\gamma} D\_t^{\gamma} R(t) = p\,\omega - \omega \, R(t) + \sigma \, I(t), \end{array} \tag{2}$$

where *<sup>C</sup> t*0 *Dγ <sup>t</sup>* is the Caputo fractional-order derivative operator of order *γ* ∈ (0, 1] which can be defined as,

$${}^{C}\_{t\_{0}}D\_{t}^{\gamma}F(t) = (1/\Gamma(n-\gamma))\left(\int\_{t\_{0}}^{t} (t-\tau)^{n-\gamma-1} \frac{d^{n}F(\tau)}{d\tau^{n}} \,d\tau\right),\tag{3}$$

when it effects on a continuous function *F* on an interval [0, *T*], where *t*<sup>0</sup> is the start time, *n* = *γ*! + 1 and *γ*! is the integer part of the fraction *γ*. The Caputo fractional-order derivative has advantages for solving initial value problems only when it applied with

analytical or semi-analytical approaches. The Caputo derivative can be approximated by the Grünwald–Letnikov (GL) method using finite differences of the fractional order, similar to the Euler method, to handle numerical solutions of initial value problems. The GL method is proceeding iteratively but the sum in the scheme becomes longer and longer, which reflects the memory effect.

This paper is structured as follows. An explanation of the Grünwald–Letnikov fractional derivative representation is given in Section 2. In Section 3, the GL discretized scheme for the considered fractional-order SIR model is developed. In Section 4, the qualitative stability analysis of the fractional-order SIR model is investigated, and the free equilibrium fixed point of the disease is examined. Section 5 contains the numerical simulation of the present results including comparisons with results displayed in [1–3] and with results of the Maple built-in scheme of the fourth order Runge–Kutta (RK4) method along with graphical illustrations. A sensitivity analysis, which demonstrates how the model dynamics differently perform as the values of the model parameters changed, is including in Section 5 as well. Finally, in Section 6, conclusions and observations are provided.

#### **2. Explanation of the Grünwald–Letnikov Fractional Derivatives**

Here, some fundamental concepts of fractional Grünwald–Letnikov derivatives are presented. The form of GL fractional derivative with order *γ* ∈ (0, 1] is expressed as:

$$\,\_{t\_0}D\_t^\gamma F(t) = \lim\_{h \to 0} (1/h^\gamma) \sum\_{j=0}^{(t-t\_0)/h} (-1)^j \binom{\gamma}{j} F(t-jh),\tag{4}$$

where *<sup>h</sup>* is the calculations time step and *<sup>γ</sup> j* = <sup>Γ</sup>(*γ*+1) <sup>Γ</sup>(*j*+1)Γ(*γ*−*j*+1) is binomial coefficients.

The explicit numerical approximation of the fractional derivative of order *γ*-th at the points *kh*, (*k* = 1, 2, . . . ) is defined by the next formula [28–30]:

$$(k - (Z/h))\_{t\_0} D\_{t\_k}^{\gamma} F(t) \approx \frac{1}{h^{\gamma}} \sum\_{j=0}^{k} c\_j^{(\gamma)} F\left(t\_{k-j}\right),\tag{5}$$

where *Z* is noted as memory length, *tk* = *kh* and *c* (*γ*) *<sup>j</sup>* , (*j* = 0, 1, ...) are the binomial coefficients that easily computed from the following equations [28]:

$$c\_0^{(\gamma)} = 1, \quad c\_j^{(\gamma)} = \left(1 - \frac{1+\gamma}{j}\right) c\_{j-1\prime}^{(\gamma)} \quad (j = 1, 2, \dots). \tag{6}$$

Therefore, a general numerical approximation of the fractional-order differential equation:

$$\_{t\_0}D\_t^\gamma F(t) = G(F(t), t), \tag{7}$$

where *G* is a nonlinear function, that can be written as,

$$F(t\_k) = G(F(t\_{k-1}), t\_{k-1})h^\gamma - \sum\_{j=\upsilon}^k c\_j^{(\gamma)} F\left(t\_{k-j}\right),\tag{8}$$

With the aid of the principle of the short memory, the summation index *v* in Equation (8) is considered as *v* = 1 for all *k* < *Z*/*h* and *v* = *k* − (*Z*/*h*) for *k* > *Z*/*h*. However, we set *v* = 1 for all *k* when the principle of the short memory is not used. Clearly, the consequence for this approximation makes us pay a penalty in the shape of a few inaccuracies. If *F* (*t*) ≤ *M*, we may easily estimate the memory length *Z*, providing the needed accuracy *ξ* as [31]:

$$Z \ge \left[ M / (\xi \cdot |\Gamma(1 - \gamma)|) \right]^{1/\gamma} \tag{9}$$

An investigation of the short memory impact and the differences between long and short memory were examined and illustrated in [30]. Generally, we can easily apply the GL fractional derivative approximation to systems of fractional-order in the following form,

$$\,\_{t\_0}D\_t^\gamma F\_i(t) = G\_i(F\_1(t), F\_2(t), \dots, F\_N(t), t), \; i = 1, 2, \dots, N. \tag{10}$$

#### **3. Grünwald–Letnikov Discretization of Fractional-Order SIR Model**

For numerical simulation of the fractional-order SIR model described in the system (2), we will implement the introduced GL discretization technique to solve it. This approach is established on the fact that for many classes of continuous functions, the fractional-order derivatives of Grünwald–Letnikov and Caputo are equivalent [31]. Using the fractional GL numerical approximation in the Equation (8), the numerical solution of the system (2), when *<sup>C</sup> t*0 *Dγ t* ≡ *t*<sup>0</sup> *Dγ <sup>t</sup>* , can be obtained using the following iterative formulas,

$$\begin{split} S(t\_k) &= [(1-p)\,\omega - \omega\,\,S(t\_{k-1}) - \beta\,\,S(t\_{k-1})\,\,I(t\_{k-1})]h^\gamma - \sum\_{j=\upsilon}^{k} c\_j^{(\gamma)}\,\,S(t\_{k-j}),\\ I(t\_k) &= [-(\sigma+\omega)\,\,I(t\_{k-1}) + \beta\,\,S(t\_k)\,\,I(t\_{k-1})]h^\gamma - \sum\_{j=\upsilon}^{k} c\_j^{(\gamma)}\,\,I(t\_{k-j}),\\ R(t\_k) &= [p\omega - \omega\,\,R(t\_{k-1}) + \sigma\,\,I(t\_k)]h^\gamma - \sum\_{j=\upsilon}^{k} c\_j^{(\gamma)}\,\,R\left(t\_{k-j}\right),\end{split} \tag{11}$$

where, *k* = 1, 2, 3, ... , *N*, for *N* = *T*/*h* and *T* is the simulation time. The binomial coefficients *c* (*γ*) *<sup>j</sup>* are calculated using the formula in Equation (6). Along this paper we set the start data (initial conditions) as follows: *S*(*t*0) = *k*1, *I*(*t*0) = *k*2, *R*(*t*0) = *k*<sup>3</sup> and the start time *t*<sup>0</sup> = 0.

#### **4. Stability Analysis of the Fractional-Order SIR Model**

In this section, we examine the stability of the considered SIR model. We proceed to find the fixed points (steady solutions) that occur when *<sup>t</sup>*<sup>0</sup> *Dγ <sup>t</sup> S*(*t*) = *<sup>t</sup>*<sup>0</sup> *Dγ <sup>t</sup> I*(*t*) = *<sup>t</sup>*<sup>0</sup> *Dγ <sup>t</sup> R*(*t*) = 0, and evaluate their stability. The fixed point that occurs at a long-term behavior of the fractional system (2) are *Pf* = *S*∗ *<sup>f</sup>* = 1 − *p*, *I*<sup>∗</sup> *<sup>f</sup>* = 0, *R*<sup>∗</sup> *<sup>f</sup>* = *p* . This fixed point is corresponding to the disease-free equilibrium which is described by the nonexistence of the infected nodes. The Jacobian matrix for the SIR system (2) is defined as:

$$J = \begin{bmatrix} -\beta I - \omega & -\beta S & 0\\ \beta I & \beta S - \omega - \sigma & 0\\ 0 & \sigma & -\omega \end{bmatrix} \tag{12}$$

At disease free equilibrium fixed point, the eigenvalues relating to the matrix *J* are,

$$
\lambda\_1 = \beta(1-p) - \omega - \sigma,\\
\lambda\_2 = -\omega,\\
\lambda\_3 = -\omega. \tag{13}
$$

This disease-free equilibrium fixed point is asymptotically stable if all the eigenvalues in Equation (13) satisfy the following condition [32]:

$$|\arg(\lambda\_{1,2,3})| > \frac{\pi}{2}\gamma.\tag{14}$$

From Equations (13) and (14), the condition that guarantees the stability of SIR dynamical system at the disease-free equilibrium point is,

$$
\beta(1-p) < \omega + \sigma, \forall \, \gamma \in (0,1]. \tag{15}
$$

Therefore, the vaccination reproduction number *Vr* can be defined as,

$$V\_r = \frac{\beta (1 - p)}{\omega + \sigma}.\tag{16}$$

If *Vr* is higher than a limit value, an infectious disease can propagate in a susceptible group. Moreover, the condition in Equation (15) shows that a critical vaccination fraction *pc* can defined as,

$$p\_{\mathcal{C}} = 1 - \frac{\omega + \sigma}{\beta}.\tag{17}$$

It is worth to note that above this critical fraction, the disease-free equilibrium is stable i.e., *p* > *pc*. Therefore, the vaccination fraction must be large enough in order to effectively avoid surge period of the disease.

The other fixed point of the dynamical system (2) is the point corresponding to the endemic equilibrium point, which is described by the existence of the infected nodes. This point is obtained as *Pe* = *S*∗ *<sup>e</sup>* <sup>=</sup> <sup>1</sup>−*<sup>p</sup> Vr* , *I*<sup>∗</sup> *<sup>e</sup>* = *<sup>ω</sup> <sup>β</sup>* (*Vr* − 1), *R*<sup>∗</sup> *<sup>e</sup>* = <sup>1</sup> *βS*∗ *e <sup>ω</sup><sup>p</sup>* <sup>+</sup> *<sup>σ</sup>* <sup>−</sup> *<sup>σ</sup> <sup>β</sup>* (*ω* + *σ*) . It is clear that the endemic equilibrium will only occur when *Vr* > 1. At this point, the eigenvalues relating to the matrix *J* are obtained as,

$$
\lambda\_1 = -\omega,\ \lambda\_{2,3} = -\frac{\omega V\_r}{2} \left( 1 \pm \sqrt{1 - 4\frac{(\omega + \sigma)}{\omega V\_r}} \right). \tag{18}
$$

From Equations (14) and (18), the endemic equilibrium point *Pe* is asymptotically stable if the following condition is fulfilled,

$$1 < V\_r \le \frac{4}{\omega} (\omega + \sigma), \; \forall \; \gamma \in (0, 1]. \tag{19}$$

#### **5. Numerical Results and Discussion**

In the current section, we present numerical results with many visualizations to simulate the dynamics of the present fractional-order SIR model subjected to several fractional order *γ*, initial conditions and parameters. Along this section, we set the time step *h* = 0.01 for the fractional GL iterative scheme and for the Maple RK4 built-in scheme. Firstly, we consider the numerical values of the case considered in [1,2] to show the inaccuracy of the results displayed in these papers. As mentioned in [1,2], we assume:

$$k\_1 = 1, \; k\_2 = 0.5, \; k\_3 = 0, \; p = 0.9, \; \not p = 0.8, \; \omega = 0.4, \; \sigma = 0.03. \tag{20}$$

According to the parameters in Equation (20), the critical vaccination fraction *pc* = 0.4625. Therefore, this case falls under the disease-free equilibrium case at which the steady state solution asymptotically approaches to the following fixed point,

$$P\_f = \left( S\_f^\* = 0.1, \ I\_f^\* = 0, \ R\_f^\* = 0.9 \right) \tag{21}$$

Looking at the results that were presented in [1,2], we found that the graphical representations of the solutions are inaccurate and do not agree with the corresponding steady state solution that shown in Equation (21). In this regard, we compared our numerical results with the results presented in [1,2] for the susceptible group *S*(*t*) only at several values of the fractional order *γ*. This comparison is displayed in Figure 1. For the validation purpose of the present numerical results, we plot the solution of Maple RK4 built-in scheme along with our solution using the GL scheme at *γ* = 1.

**Figure 1.** Comparison of *S*(*t*) between (**a**) present results, (**b**) results of [1] and (**c**) results of [2].

From Figure 1, it is shown that the graphical representation of *S*(*t*) in [1,2] is incorrect even in the transient region of the solution. From Figure 1c that was presented in [2], it can be shown that *S*(0) = 0 although the authors considered *S*(0) = 1 in their study. Moreover, the graphical representations of *I*(*t*) and *R*(*t*) that displayed in [1,2] are inaccurate as well because they do not agree with the corresponding fixed point in Equation (21). For *γ* = 1, the excellent agreement between the approximate solution using RK4 scheme and the GL scheme validates the approximate solution of GL scheme. Figure 1a illustrates that our approximate solution agrees with the steady state fixed point shown in Equation (21). This note is considered as another validation of the present results.

The rest of this section is related to the generalization of the cases presented in [3] to examine the analyze of those cases subject to various fractional orders, mainly *γ* = 1, *γ* = 0.75 and *γ* = 0.5 for the same initial conditions and parameters. Because of the considered SIR model encounters long interepidemic periods, we set the simulation time *T* = 250. For the purpose of the comparison with the results of [3], two versions of the solutions graphs were displayed. One at a small-time behavior *T* = 10, to show how the solutions at the transient region behave; and the other at a long-time behavior *T* = 250, to show how the solutions behave at long interepidemic period (steady state solution). Here, we set *β* = 0.8, *ω* = 0.4 and *σ* = 0.03 for all the considered cases. According to these parameters, the critical value of the vaccination portion is obtained as *pc* = 0.4625. The results of the RK4 method at *γ* = 1 are plotted by circle points (◦).

• **Case 1**

We consider the initial conditions and case parameters as *k*<sup>1</sup> = 1, *k*<sup>2</sup> = 0, *k*<sup>3</sup> = 0 and *p* = 0.9. Here, we have a stable disease-free equilibrium case with a vaccination reproduction number *Vr* = 0.186047.

This case demonstrates the effect of high-level vaccination exposure on the initial individuals with no infective individuals. Figure 2 displays a comparison between present numerical results for *S*, *I* and *R* versus time for various fractional orders, *γ* = 1, *γ* = 0.75 and *γ* = 0.5, against the results obtained in [3] when *γ* = 1. It can be shown that there is an excellent agreement between the results of GL scheme, RK4 scheme and the obtained in at *γ* = 1. Furthermore, the graphical representation of the numerical results demonstrates that the solutions are monotone dependent on the order of the fractional derivative *γ*. The susceptible, infective, and removed individuals decrease as the fractional order *γ* decreases in both transient and steady state stage. Generally, as *γ* diminishes, the fixed equilibrium solution reduces as well. For this case, the disease has been successfully eradicated because of the effect of high-level vaccination fraction *p*. It is worth to notice that all individuals normally remain disease-free for all times.

**Figure 2.** Behavior of various individual fractions (*S* in black, *I* in red and *R* in blue) versus time for Case 1 at different fractional order values (**a**) till *T* = 10, (**b**) results of [3] and (**c**) till *T* = 250.

#### • **Case 2**

In this case we make a slight change in the initial condition to study the influence of high-level vaccination exposure on the initial low-level infective individuals. Here we consider *k*<sup>1</sup> = 0.8, *k*<sup>2</sup> = 0.2, *k*<sup>3</sup> = 0 and *p* = 0.9. This case is also a stable disease-free equilibrium one, with a vaccination reproduction number *Vr* = 0.186047, in which the disease will be eradicated.

Figure 3 displays a comparison between present numerical results for susceptible, infective and removed individuals versus time for various fractional order *γ*, along with the results obtained in [3] when *γ* = 1. When considering an integer order, an excellent agreement between the results of GL scheme, RK4 scheme and the results obtained in [3] is realized. As shown from Figure 3, the susceptible and infective individuals reduce with the passage of time whereas the removed individual raises due to the appearance of recovered and vaccinated groups with long-lasting resistance against disease and hence the disease will be eradicated. Here, the disease-free equilibrium is reachable as soon as the vaccination amount level is more than the critical vaccination value *pc* (i.e., *p* > *pc*) as occurred in Case 1. The impact of the fractional order *γ* on the dynamic of the solutions is the same as its effect in Case 1.

**Figure 3.** Behavior of various individual fractions (*S* in black, *I* in red and *R* in blue) versus time for Case 2 at different fractional order values (**a**) till *T* = 10, (**b**) results of [3] and (**c**) till *T* = 250.

#### • **Case 3**

Here we consider an endemic-equilibrium case by setting *k*<sup>1</sup> = 0.8, *k*<sup>2</sup> = 0.2, *k*<sup>3</sup> = 0 and *p* = 0.3. This case explains the impact of low-level vaccination exposure on the initial individuals with low-level infective population. Figure 4 shows a graphical representation of numerical results for the susceptible, infective, and removed population versus time for various fractional order *γ*, alongside the results obtained in [3]. As expected, when *γ* = 1, an excellent agreement between the obtained results and the results displayed in [3] is achieved.

From Figure 4, it can be notice that the number of susceptible individuals reduces, while the removed individuals increase with a small amount as time expands. Nevertheless, it is notable that the infective individuals will never vanish as time expands and hence the endemic case continues with no disease eradication. This proves that a disease-free equilibrium could not be accomplished when the vaccination amount is less than the critical vaccination value (i.e., *p* < *pc*). In this case, the endemic equilibrium stays stable because a vaccination reproduction number *Vr* = 1.302326 satisfies the stability condition in Equation (19). In such an endemic-equilibrium case, the fractional order *γ* plays an influential role. The fractional order affects the susceptible, infective, and removed individuals by a different manner from disease-free equilibrium cases. When *γ* < 1, we found that as the fractional order *γ* decreases, the susceptible and infective groups rapidly decrease in the transient stage and then begin to increase once again in the steady state stage. In addition, for *γ* < 1, the susceptible individuals will be more than the amount for *γ* = 1 in the transient stage. While the infective individuals remain less than the amount for *γ* = 1 in both transient and steady state stages.

**Figure 4.** Behavior of various individual fractions (*S* in black, *I* in red and *R* in blue) versus time for Case 3 at different fractional order values (**a**) till *T* = 10, (**b**) results of [3] and (**c**) till *T* = 250.

#### • **Case 4**

Here we consider another endemic-equilibrium case by setting *k*<sup>1</sup> = 0.8, *k*<sup>2</sup> = 0.2, *k*<sup>3</sup> = 0 and *p* = 0. This case describes the influence of low-level vaccination exposure on the initial individuals with no vaccinated people. Figure 5 shows the results visualization for the susceptible, infective, and removed population with time for various values of *γ*, beside the results obtained in [3] for *γ* = 1.

As expected, an excellent agreement between our results and the results obtained in [3] is realized. Similar to Case 3, we have a stable endemic-equilibrium situation with a vaccination reproduction number *Vr* = 1.860465. Here, the number of susceptible individuals decreases while the infective individuals increase by the time for the steady state stage. The only role of the removed individuals is the very little fraction of the recovered population with long-lasting immunity. Here, the disease is quickly transported to the bulk of people. The impact of the fractional order *γ* on the dynamic of the solutions is similar in its impact in Case 3.

A sensitivity analysis of the considered cases is displayed in Table 1 which illustrates how the model dynamics differently behave as per parameters values changed. From the data presented in the table, it is clear that the numerical results are in excellent agreement with the results of the system fixed points in both disease-free equilibrium and endemic equilibrium case. For the disease-free equilibrium cases, the fraction of the infected population is vanishing for all values of the fractional order *γ* by the time and hence a disease has been eradicated. While for the endemic equilibrium cases, the fraction of the infected population decreases with the decreasing of the fractional order *γ* by the time, but the disease is not eradicated.

**Figure 5.** Behavior of various individual fractions (*S* in black, *I* in red and *R* in blue) versus time for Case 4 at different fractional order values (**a**) till *T* = 10, (**b**) results of [3] and (**c**) till *T* = 250.


**Table 1.** Impact of various parameter values on the population fractions.

#### **6. Conclusions**

This paper deals with a new application of a Grünwald–Letnikov method for the achievement of a numerical simulation of a fractional-order SIR epidemiological diseases model. The significant advantage of the considered technique is that it can be used without system linearization or any other restrictions for any fractional-order dynamical system. A stability analysis is performed to qualitatively investigate the dynamics of the model in both disease-free and endemic equilibrium case. It is explained that the graphical representation presented in [1,2] are inaccurate. One of the issues that validates the present results is the excellent agreement between them with the results presented in [3] and RK4 results for integer order *γ* = 1. The excellent agreement between the asymptotic behavior of the solutions with the numerical results is another evident on the validation and the reliability of the considered numerical technique in studying such dynamical models in epidemiology. The proposed technique is reliable and accurate in simulating the behavior of the solutions for a long-time interval. The disease-free equilibrium is stable if the vaccination amount above a critical value *pc*. The dynamics of the model is strongly dependent on the value of the fractional order *γ*. For the endemic equilibrium cases, the fraction of the infected population decreases by the reduction of the fractional order *γ* as time passes.

**Author Contributions:** Both authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Majmaah University, grant number R-2021-260.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the Deanship of Scientific Research at Majmaah University for supporting this work under Project Number No.: R-2021-260.

**Conflicts of Interest:** The authors declare that there are no conflicts of interest.

#### **References**

