A Mathematical Model of Vaccinations Using New Fractional Order Derivative

Purpose: This paper studies a simple SVIR (susceptible, vaccinated, infected, recovered) type of model to investigate the coronavirus’s dynamics in Saudi Arabia with the recent cases of the coronavirus. Our purpose is to investigate coronavirus cases in Saudi Arabia and to predict the early eliminations as well as future case predictions. The impact of vaccinations on COVID-19 is also analyzed. Methods: We consider the recently introduced fractional derivative known as the generalized Hattaf fractional derivative to extend our COVID-19 model. To obtain the fitted and estimated values of the parameters, we consider the nonlinear least square fitting method. We present the numerical scheme using the newly introduced fractional operator for the graphical solution of the generalized fractional differential equation in the sense of the Hattaf fractional derivative. Mathematical as well as numerical aspects of the model are investigated. Results: The local stability of the model at disease-free equilibrium is shown. Further, we consider real cases from Saudi Arabia since 1 May–4 August 2022, to parameterize the model and obtain the basic reproduction number R0v≈2.92. Further, we find the equilibrium point of the endemic state and observe the possibility of the backward bifurcation for the model and present their results. We present the global stability of the model at the endemic case, which we found to be globally asymptotically stable when R0v>1. Conclusion: The simulation results using the recently introduced scheme are obtained and discussed in detail. We present graphical results with different fractional orders and found that when the order is decreased, the number of cases decreases. The sensitive parameters indicate that future infected cases decrease faster if face masks, social distancing, vaccination, etc., are effective.


Introduction
Mathematical models are recognized as crucial in epidemiology for understanding the dynamics of diseases and making predictions about their long-term behavior. With the passage of time and the emergence of new infectious diseases in humans populations, mathematical models have been used to determine the peak infection curve, the days to eradication, and the number of possible future cases. In the case of the coronavirus disease, the early prediction of the peak of the infection, the basic reproduction number, and the individuals and protect them from future illness. It is safe and most people can use it without any fear [34][35][36][37].
This paper investigates the dynamics of the coronavirus infection in Saudi Arabia using recently reported cases through the fractional-order vaccination model. We use the recent cases of the fourth wave in Saudi Arabia and implement a mathematical model first in the integer order and then extend it to the generalized order model. The model is directly fitted to the cases in which the vaccine is present and studies the equilibrium points analysis. We observe the possibility of a backward bifurcation phenomenon, where the disease-free equilibrium coexists with the endemic state and hence the global asymptotical stability of the disease-free equilibrium does not exist. We divide the work section-wise: Section 2 gives details of the newly fractional derivative considered by Hattaf [38]. Further, it discusses the formulation of the problem and further extends the model into the fractional order system. Furthermore, we study the equilibrium points, the basic reproduction, backward bifurcation, and the local asymptotical stability of the disease-free case and explain the algorithm for the numerical simulation of the fractional model. In Section 3, we discuss the graphical results and present results regarding disease controls. The results are summarized briefly in Section 4.
Below, for Definition (1), we can write the corresponding fractional integral as Definition 2 ([38]). For the newly fractional derivative D q,q 1 l 1 ,φ , the corresponding fractional integral can be expressed as where RL I q 1 l 1 ,φ of order q 1 denotes the standard weighted Riemann-Liouville fractional integral and is defined by Theorem 1 ([39]). Suppose y = 0 is an equilibrium point of and V(y) is a continuously differentiable function in a neighborhood U ∈ R n of the origin holds the conditions below: (i) V(0) = 0 and V(y) > 0 for all y ∈ U\{0}; (ii) D q,q 1 0,φ V(y) ≤ 0 for all y ∈ U\{0}.

Model Formulation
We consider an SVIR model and denote its total population by N(t). The model consists of four components: the healthy individuals that have the ability to become infected after close contact with infected COVID-19 people is shown by S(t); individuals that are vaccinated are given by V(t); individuals that are infected are given by I(t); and those recovered from infection of COVID-19 or vaccination are given by R(t). We write N(t) = S(t) + V(t) + I(t) + R(t). The population of healthy individuals is obtained through the birth rate Λ, while the natural mortality rate in each compartment is given by µ. Healthy individuals become infected when they have close contact with infected people, and hence the route of the transmission is βSI/N, while vaccinated individuals after close contact with infected people are shown through the route β 1 V I/N. The portion of healthy individuals to be vaccinated is shown by ω. The vaccinated and the infected individuals are recovered, respectively, by the rate γ 1 and γ. The disease mortality of the COVID-19 infected people in the infected compartment is given by d 1 . With these assumptions, the COVID-19 model with vaccination is given by the following nonlinear differential equations: with the non-negative initial conditions We consider the following biologically feasible region for the model (5), Γ = (S, V, I, R) ∈ R 4 + : S, V, I, R ≥ 0, and N ≤ Λ/µ , which is positively invariant for any trajectory of the system for an initial condition, which will remain in Γ for every time t ≥ 0. Therefore, the region is positively invariant, and its dynamical results can be studied within Γ. It can be observed from model (5) that the equation R can be eliminated without any loss of generality, as it does not appear in the rest of the equation. The results of R can be easily obtained using the relation R = N − S − V − I. Using this fact, in the following, we focus our study to analyze the fractional model without the last equation. There are no transmission rates from equations R to the rest of the equations, so one can ignore and reduce it, while the results of recovery cases can be obtained using the equation

A Fractional Model
We apply the recently introduced fractional derivative by Hattaf given in Definition 1 to our model (5) and obtain the following generalized fractional order model: and the related initial conditions

Analysis of the Model
This section considers the mathematical results involved in the fractional order system (6). In a dynamical system, first, we obtain the possible equilibrium points of the disease model (6). In general, the models often formulated for the disease-related human population consist of two equilibrium points, the infection-free and the infected. The infection-free equilibrium can be denoted by P 0 of the model (6), which one can obtain as follows: , 0 .
Another important concept in disease epidemiology is the computation of the basic reproduction number. The basic reproduction tells us about the disease's progress, and whether it can be controlled or spread within the population. For our SVIR-type model (6), we can determine this number using the last equation of the system (6) within the disease-free case P 0 and obtain the following result: The basic reproduction, or in this case the vaccine reproduction number R v 0 , consists of two parts: the first part R 1 is associated with vaccine cases, while the other one R 2 is related to cases without vaccination. It is obvious that the vaccine reduces the basic reproduction number, as vaccines for any disease in the literature prove that vaccines are the best control of disease. We obtain the basic reproduction number with no vaccination by putting ω = 0, and obtain the following:

Endemic Equilibria
Here, we shall investigate the endemic equilibrium of the vaccine model (6) given by P 1 P 1 = (S, V, I) = (S * , V * , I * ) (8) and can determined by equating H D 0,φ I(t) = 0, and obtain the following, Inserting (9) into the following expression, and obtain We obtain the following, where a 0 = ββ 1 , Here, it can be seen that a 0 > 0 and a 2 can be positive if The result for the endemic equilibria can be summarized in the following form: The coronavirus model (6) has the following: There exists a unique endemic equilibrium if a 1 < 0 and a 2 = 0 → R v 0 = 1, 3. We can have two endemic equilibria if a 2 > 0 → R v 0 < 1, a 1 < 0 and its related discriminant is positive 4. Above the other cases, there is no possible equilibria.
It is clear from the first part (i) of the Theorem 2 that we have a unique positive endemic equilibrium whenever R v 0 > 1. Further, the third part of Theorem 2 tells us about the occurrence of the phenomenon of backward bifurcation in the COVID-19 infection model (6). This means that disease-free equilibrium coexists with endemic equilibrium, and the model will not be globally asymptotically stable. In such a case, the disease will persist in the population for a long time and need vaccination and other control measures for its elimination and control. To achieve the mathematical expression and its graphical result, we set the discriminant a 2 1 − 4a 0 a 2 = 0 and then solve further for the critical values of R v 0 denoted by R c , which is given by Therefore, the backward bifurcation may occur for the values of R v 0 such that R c < R v 0 < 1. Consider the listed value Λ = 1273.94, β = 0.6, β 1 = 0.2, γ = 0.05, γ 1 = 0.04, µ = 1/(74.87 × 365), d 1 = 0.024 and ω = 0.15. The related bifurcation plot is given in Figure 1. In Figure 1, one can see that β is the bifurcation parameter that can cause the backward bifurcation. In such cases, the model may or may not be globally asymptotically stable at the disease-free equilibrium.

Stability Analysis
The stability analysis of Model (5) can be studied for the disease-free case. We present these results in the following theorem: Theorem 3. The fractional-order SVIR model is locally asymptotically stable, provided that R v 0 < 1.

Proof.
We have the Jacobian matrix of the system (6) evaluated for the disease-free case P 0 , and it is given by The Jacobian matrix J(P 0 ) can be expanded, and we can obtain the eigenvalues as follows: . The first two eigenvalues are obviously negative, while the third eigenvalue can be negative if R v 0 < 1. Therefore, all the roots of the Jacobian matrix contain negative real parts, so the fractional-order system (6) at the disease-free equilibrium point P 0 is locally asymptotically stable provided that

Global Stability
To show the global stability of the model at P 0 when R v 0 ≤ 1, we construct the Lyapunove function given by Therefore, the above result can be stated as follows:

Global Stability at Endemic State
We need the following results in the proof of the following theorem, with the assumption βSI ≤ N and β 1 V I ≤ N: Proof. We define the Lyapunove function given by where Ψ(y) = y − 1 − ln y, for y > 0. It is clear that Ψ(y) attains its global minimum at y = 1 and Ψ(1) = 0. Therefore, Ψ(y) ≥ 0 for every y > 0. Thus, L(S, V, I) ≥ 0 with L(S * , V * , I * ) = 0. Applying the Corollary 2 given in [40], we have Using the equation from system (6) into (15), and calculating the terms of Equation (15), we have Using Equations (16)- (18) in Equation (15), we have Therefore, H D q,q 1 0,φ L(t) ≤ 0 when R v 0 > 1. It follows from Theorem 1 that the endemic equilibrium P 1 of the fractional-order model (6) when R v 0 > 1 is globally asymptotically stable.

Numerical Scheme and Its Results
This section discusses the numerical scheme for the new fractional generalized derivative and obtains numerical simulations. We follow the same stepping as mentioned in [41]. The generalized fractional derivative is given by D q,q 1 0,φ y(t) = g(t, y(t)).
Equation (20) is converted into the following form: Considering that t n = n∆t, with n ∈ N, we have which leads to the following, where l(ξ, y(ξ)) = φ(ξ)g(ξ, y(ξ)). One can approximate the function l in the interval [t k , t k+1 ] as it is given in [42]. The Lagrange polynomial interpolation passing these points (t k−1 , l(t k−1 , y k−1 )) and (t k , l(t k , y k )) is as follows: Thus, The integrals inside Equation (28) can be determined as follows: where Finally, we achieve the required scheme as follows: n,k,q 1 +φ(t k−1 )g(t k−1 , y k−1 )A 2 n,k,q 1 (28)

Sensitivity Analysis
Sensitivity analysis is critical for determining how best to minimize coronavirus mortality and morbidity, as well as the relative relevance of the many factors responsible for its transmission and prevalence. In this subsection, we will find the model parameters that have a large influence on R v 0 . The following formula should be used to determine the sensitivity analysis of the parameters involved in the basic reproduction number R v 0 [43].
Definition 3. The normalized forward sensitivity index of a variable, w, for which differentiability depends on a parameter, q, is defined as Using the formula mentioned in the above definitions, we calculate the analytical for each of the different parameters β = 0.28, β 1 = 0.2, γ = 0.05, γ 1 = 0.04, µ = 1/(74.87 × 365), d 1 = 0.024, and ω = 0.15. We now calculate the sensitivity index of R v 0 with respect to the parameter β as In a similar way, we can calculate the rest of the indices as shown in Table 1. It can be observedfrom the values given in Table 1 that the parameter β 1 , followed by γ, d 1 , β, and so on, can increase or decrease the basic reproduction number. Upon decreasing the contact among the susceptible and infected, and increasing the recovery rate by vaccinating the individuals, the number of infected individuals shall decrease.

Results and Discussion
We consider the following numerical values and the initial conditions in our numerical simulation of the model (6) Figure 3 shows the behavior of the model variables for various values of q and keeping q 1 fixed. Figure 4 is given to show the behavior of the model when q is fixed and q 1 varies for various values. Varying both the values of q and q 1 , we have plotted the results graphically in Figure 5. The contact rates β and β 1 have a great impact on the disease spread and control, which is shown graphically in Figure 6. When the value of β (the contact among the healthy and the infected compartments, such as social distancing, avoiding gatherings, using face masks, etc.) and β 1 (the contact among healthy and vaccinated people) decrease, the number of infected people decreases. We give the results of the model numerically in Figures 7 and 8 when φ(t) = (1 + q) q 1 φ(t) = (1 + exp(−t)) q 1 , respectively, for many values of the fractional-order parameters q and q 1 . One of the advantages of this new fractional derivative is the use of the function φ(t), where one can fit the data well using an appropriate value for φ(t). We also compare the present scheme given in [41] with [42]. We provide such a comparison in Figures 9 and 10. In Figure 9, we fix φ(t) = 1, q = q 1 = 1 and show the comparison of the present method with the Atangana-Taufik scheme shown in [42]. Similarly, when q = 0.9, 0.7, we show the comparison of the results in Figure 10. It should be noted that the present scheme generalizes the Atangana-Taufik method, so we put φ(t) = 1. The comparison of the basic reproduction number with and without vaccination is shown in Figure 11. It is clear from the comparison results that the present scheme is matched perfectly with the scheme given in [42].          (d) Figure 10. Comparison of the present method with Atangana-Taufik method, for q = 0.9, 0.7, q 1 = 1 and φ(t) = 1. Subfigure (a) shows the comparison of the schemes for healthy population when φ(t) = 1 and q = 0.9, 0.7, q 1 = 1. Subfigure (b) shows the comparison of the schemes for vaccinated population when φ(t) = 1 and q = 0.9, 0.7, q 1 = 1. Subfigure (c) shows the comparison of the schemes for infected population when φ(t) = 1 and q = 0.9, 0.7, q 1 = 1. Subfigure (d) shows the comparison of the schemes for recovered population when φ(t) = 1 and q = 0.9, 0.7, q 1 = 1.

Conclusions
We investigated the dynamics of coronavirus infection cases from Saudia Arabia using the fractional model. The new fractional derivative considered was recently reported in the literature and is known as the generalized Hattaf fractional derivative. We presented the background results for the new fractional derivative and then considered the model using the Hattaf derivative. The equilibrium points are obtained, and their stability is discussed. It can be observed that the fractional model is locally asymptotically stable when R 0 < 1, and it is unstable otherwise. Further, we obtain the global stability of the model when the basic reproduction number R 0 > 1. Considering the reported cases of the COVID-19 in Saudi Arabia, the model is fitted to the data, and their results are obtained for both integer and non-integer cases. For the behaviors of the fractional-order parameters q and q 1 , we have presented some numerical results that demonstrate the effectiveness of the new fractional derivative. The basic reproduction number without vaccination is R 0 = 3.78, while with vaccination, it is R v 0 = 2.92. One of this generalized fractional derivatives is the use of the new function that results when dealing with data of a different nature, which are difficult to fit using other fractional operators. The parameters' values that decrease the future cases have been shown graphically.