Investigating a Fractal–Fractional Mathematical Model of the Third Wave of COVID-19 with Vaccination in Saudi Arabia

: The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is responsible for coronavirus disease-19 (COVID-19). This virus has caused a global pandemic, marked by several mutations leading to multiple waves of infection. This paper proposes a comprehensive and integrative mathematical approach to the third wave of COVID-19 (Omicron) in the Kingdom of Saudi Arabia (KSA) for the period between 16 December 2022 and 8 February 2023. It may help to implement a better response in the next waves. For this purpose, in this article, we generate a new mathematical transmission model for coronavirus, particularly during the third wave in the KSA caused by the Omicron variant, factoring in the impact of vaccination. We developed this model using a fractal-fractional derivative approach. It categorizes the total population into six segments: susceptible, vaccinated, exposed, asymptomatic infected, symptomatic infected, and recovered individuals. The conventional least-squares method is used for estimating the model parameters. The Perov fixed point theorem is utilized to demonstrate the solution’s uniqueness and existence. Moreover, we investigate the Ulam–Hyers stability of this fractal–fractional model. Our numerical approach involves a two-step Newton polynomial approximation. We present simulation results that vary according to the fractional orders ( γ ) and fractal dimensions ( θ ), providing detailed analysis and discussion. Our graphical analysis shows that the fractal-fractional derivative model offers more biologically realistic results than traditional integer-order and other fractional models.


Introduction
The onset of the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in late 2019 marked the beginning of a global pandemic that profoundly impacted all aspects of human society.Coronavirus disease-19 (COVID-19) has caused widespread illness, loss of life, economic disruptions, and a significant strain on healthcare systems worldwide.Typical symptoms include fever, cough, fatigue, and difficulty breathing, with the disease being more severe for the elderly and those with underlying chronic conditions.
Mathematical modeling has been instrumental in understanding and combating the pandemic.Among the various models, fractal-fractional modeling stands out, shedding light on the intricate transmission patterns of COVID-19.Rooted in fractal geometry and fractional calculus, this approach captures the self-similar patterns in the virus's spread.It provides a mathematical framework mirroring the nonlinearity and variability observed in actual epidemic data [1][2][3].Multiple mathematical frameworks have been developed and scrutinized to understand COVID-19 transmission globally.For instance, Ref. [4] proposed a model incorporating lock-down measures for COVID-19.Additionally, the impact of unidentified cases was examined in [5] through a distinct mathematical framework.Specific preventive strategies' influence on the trajectory and management of COVID-19 in Pakistan was detailed in [6].Furthermore, Ref. [7] considered the environmental dissemination of the virus in a transmission model focused on a Saudi Arabian case study.Notably, the majority of COVID-19 models rely on integer-order derivatives, which may sometimes inadequately represent certain realistic features.
To overcome these challenges, fractional-order derivatives serve as a valuable tool for understanding disease dynamics, providing additional crucial insights.Models using fractional-order mathematics inherently possess memory properties, offering a more suitable approach to characterize epidemic models.The application of fractional-order derivatives in disease dynamics has been explored in works such as [8,9] and related studies.
In recent times, a surge in the development of fractional-order mathematical models aimed at unraveling the dynamic behaviors of the emerging COVID-19 pandemic has been observed.For instance, Ref. [10] presented a COVID-19 model utilizing Caputo fractional operators, focusing on pandemic data from China.Additionally, a model centered on the fractional Atangana-Baleanu operator was introduced to depict the spread of COVID-19.The analysis prompted a discussion on the represented transmission dynamics and their implications for societal infection patterns [11].
In [12], a model based on the SEIR compartment, utilizing non-singular derivatives, examined the progression of COVID-19 in India.The q-Homotopy analysis transform method was employed for solution derivation.Additionally, Ref. [13] introduced dynamic and numerical approximations for an arbitrary-order COVID-19 system, considering three specific derivative operators.Furthermore, Ref. [14] advanced a fractional model within the Atangana-Baleanu-Caputo context, emphasizing the advantages of non-integer orders over integer orders.
This paper aims to explore the transmission dynamics of COVID-19 in Saudi Arabia by employing a compartmental model grounded in Caputo fractional-order derivatives.The progression of infectious individuals is categorized into reported cases (identified through COVID-19 testing) and unreported cases (those who have not undergone testing).As a result, the initial SAIR (Susceptible-Exposed-Infected-Recovered) model with vaccination evolves into SVAI u I r R (Susceptible-Vaccinated-Exposed-Asymptomatic infected-Symptomatic infected-Recovered) (refer to Figure 1 for the virus's dynamics in each compartment).Utilizing actual data from Saudi Arabia, we estimate the parameter values for the suggested mathematical model of COVID-19.The paper's structure is outlined as follows: Section 2 presents preliminary findings utilized throughout the paper.Section 3 details the Caputo fractional-order model and the estimation methodology.Section 4 delves into fundamental model analysis, including solution existence and uniqueness, the basic reproduction number indicating disease transmissibility, and the stability of model equilibria.Section 5 introduces the numerical scheme for the FF model.Numerical simulations and model discussions are covered in Section 6, and Section 7 concludes with final observations.

Preliminary Results
This section gives some basic definitions and fundamental concepts related to fractional and fractal-fractional calculus and generalized Banach spaces; for more details, we suggest Refs.[1,[15][16][17].
Definition 1. Recall that the Riemann-Louiville fractional integral of order γ of a function f : R + −→ R is defined by where γ ∈ R + is the order of integration, and Γ(γ) = +∞ 0 e −t t γ−1 dt is the Gamma function.

Definition 2. The Caputo fractional derivatives of order
Here, n = [γ] + 1, where [γ] denotes the integer part of γ from the set of positive real numbers R + .The Caputo derivative and the Riemann-Liouville integral adhere to the subsequent properties: For C 0 D γ t (C), where C is an element of the real numbers R, the result is zero;

Definition 3 ([15]
).Consider a function ϱ that is differentiable within the open interval (a, b).Assuming ϱ possesses fractal differentiability of order θ over this interval, we can then define its fractal-fractional derivative (abbreviated as FF derivative) of order γ following the Caputo approach, characterized by a power law.
where dϱ(t) Lemma 1. Equation (3) can be expressed in the following manner: The Perov fixed point theorem is used to demonstrate the existence and the uniqueness.To proceed, we shall explore several key topics from generalized Banach spaces that are linked to this theorem.Definition 4. Let E be a vector space over K = R or C. A generalized norm on E is a map: Banach space (in short, GBS) if the vector-valued metric space generated by its vector-valued metric where O n is denoted for the zero n × n matrix.Definition 6.Given a generalized metric space (E , δ G ) and an operator N mapping E onto itself, the operator N is defined as a Υ-contraction.This is characterized by a matrix Υ in the set M n×n (R + ), which converges to the zero matrix O n .It holds that for any elements ϱ, v in E , the following inequality is satisfied: The Perov fixed point theorem, which is an extension of the Banach contraction principle, follows.
Theorem 1 ([17]).Let E be a complete generalized metric space and let N : E −→ E be an Υ-contraction operator.Then, N has a unique fixed point in E .

Mathematical Model Description
To formulate a mathematical model for the spread of coronavirus, considering the impact of vaccination, we divide the total population N(t) into six subclasses: susceptible S(t), vaccinated V (t), exposed A(t), asymptomatic infected I u (t), symptomatic infected I r (t), and recovered R(t).The susceptible class is generated by newborn babies at a rate of Λ.The number of susceptible individuals is decreased on the one hand by transmission to the exposed class at the rate , where ς is the effective contact rate with infected people and ϵ is the coefficient of symptomatic individuals, and on the other hand by vaccination at the rate ω V and natural mortality m, while it increases through the migration of both vaccinated and recovered people who lose immunity at the rates ω S and ω 2 , respectively.The group of exposed people grows by acquiring susceptible and vaccinated humans who have been in contact with infected individuals, while it decays by losing exposed people to the symptomatic and asymptomatic infected subclasses at the rates of ρσ and ρ(1 − σ), respectively.People displaying clear signs of illness are transmitted to the recovery group at a rate v u , while they die naturally or due to the coronavirus illness at the rates m and d 1 , respectively.Conversely, asymptomatic people become recovered at the rate v r and die due to the disease at the rate d 2 .Based on the above discussions and diagram (Figure 1): Thus, the mathematical model is represented by the following system of differential equations: Furthermore, the initial conditions considered are as follows:

Parameter Estimation Procedure
The method for estimating the model parameters is detailed in this section of the article.The conventional least-squares method is used for this objective.This method has been used in several recent publications and is employed to determine the optimal fit for a collection of data points by reducing the total of the offsets or residuals of points from the contoured curve [6,9,18].Sources cited in [19] were used to determine the birth rate Λ and the natural mortality rate m. Figure 2 illustrates the curve of the model fit to the actual data points, while Table 1 enumerates the estimated and fitted parameters.The most crucial threshold parameter, the basic reproduction number, has an estimated numerical value of 1.20.

Mathematical Model Description by Fractal-Fractional Derivative
The following initial conditions apply: The system ( 6) and ( 7) can be rewritten in the following form: Here, the vector Y (t) = (S(t), V (t), A(t), I u (t), I r (t), R(t)) T , and the operator F is defined as follows: The initial condition is

Existence and Uniqueness Results
Next, we convert the system (6) into the following integral equation by using initial conditions (10) and Lemma 1: The existence and uniqueness of the solution of systems ( 6) and ( 7) are discussed in the generalized Banach space Assume the existence of a vector, each of whose components is positive, that satisfies In addition, if the matrix , and given that H(γ, θ) represents the beta function for γ and θ, the system outlined in (6) and (7) possesses a unique solution within the specified space ∏ 6 1 C([0, T]), and the solution holds for all t > 0 in B(Y 0 , Υ).Here, B(Y 0 , Υ) is the closed ball in the generalized Banach space Proof.It is clear that Y is a solution to the systems ( 6) and (7) if and only if it satisfies the following equation: where M is the operator defined by M : We begin by proving that the operator defined above satisfies the following inequality: To this end, let ℜ ∈ B(Y 0 , Υ), and we have where Then, the operator M maps B(Y 0 , Υ) onto itself.Further, let Y = (S, V, A, I u , I r , R), Ȳ = ( S, V, Ā, Īu , Īr , R) ∈ B(Y 0 , Υ), and when t falls within the interval [0, T], we observe By determining the supremum across t, we derive In the same manner, we have By the linearity of F 4 , F 5 , F 6 we can directly find By rewriting Equations ( 15)-( 20) in matrix form, we find , where ℧ is defined in (2).Then, F is a G-Lipschitz operator.By using this fact, we conclude that for any Y, Ȳ ∈ B(Y 0 , Υ): Since the matrix hence, by applying the Perov fixed point theorem 1, the systems ( 6) and ( 7) have a unique solution.

Stability Analysis of the Fractal-Fractional COVID-19 Model in the Ulam-Hyers Context
In this subsection, we aim to establish the Ulam-Hyers stability for our suggested model, drawing on definitions from [20].Definition 7. Let (X, d G ) be a generalized metric space and F : X → X be an operator.Then, the fixed point equation is said to be generalized Ulam-Hyers stable (GUHS) if there exists an increasing function ψ : R m + → R m + , continuous in 0 R m with ψ(0) = 0, such that, for any ε := (ε 1 , . . . ,ε m ) with ε i > 0 for i ∈ {1, . . ., m} and any solution Y * ∈ X of the inequalities Let the following hold: Lemma 2. The solution of the perturbed model fulfills the relation given below: Proof.The resolution of ( 23) is expressed as Then, we have sup Repeating the same procedure as the rest of the equations of system (23), we have sup sup sup sup Hence, the proof is completed.
Proof.Assuming X = (S, V, A, I u , I r , R) represents a possible resolution of the inequality (24), and X * = (S * , V * , A * , I * u , I * r , R * ) is the distinct solution of (6), then In the same manner, we find and Considering the convergence of the matrix ℧ to zero, it then follows that Therefore, (6) exhibits GUHS.

Numerical Scheme for the FF Model
In this section, we employ a two-step Newton polynomial approximation [21] to establish the numerical schemes for our fractal-fractional model describing coronavirus dynamics.We reorder the fractal-fractional problem as follows: . At the point t n+1 , we have the following: Thus, we can write Using the two steps of the Newton polynomial, we obtain After some calculation, we obtain the numerical scheme for Model (6) of the coronavirus epidemic as follows: , .

Simulations and Discussion
In this section, we employ the numerical scheme (28) to explore simulations of system (1).Given the parameter values listed in Table 1 and the initial conditions S(0) = 34,813,051, V(0) = 2000, A(0) = 1500, I u = 204, I r = 80, and R(0) = 0, we discuss the influence of fractal and fractional orders, as well as various values of the contact and vaccination parameters, on the dynamics of COVID-19 disease transmission.
Figures 3-5 illustrate the numerical results depicting the impact of fractal and fractional orders on the dynamics of COVID-19 disease transmission.It is noted that the populations in the susceptible subclass decrease, while they increase in the remaining five groups over time.
In Figure 3, we fix the fractional order γ = 1 and observe the behavior of the six groups of the model under the impact of different values of the fractal dimension θ = 1, 0.97, 0.95, 0.93.When the fractal dimension is smaller, both the decrease and increase rates in the susceptible group and the other five groups are slower when compared to higher values of θ. Figure 4 displays the numerical graphics of the six groups of the model when fixing the fractal dimension θ = 1 and varying the fractional order γ = 1, 0.97, 0.95, 0.93; we see that for smaller values of the fractional order γ, the rates of decline in the susceptible group and the growth in the vaccinated group are slower compared with large values of gamma, while smaller values of γ lead to rapid growth in the subclasses of exposed, asymptomatic infected, symptomatic infected, and recovered.
We observe the same behavior in Figure 5, when we vary both the fractal dimension and fractional order together.However, the deceleration of decay and growth rates in the six subclasses is more pronounced compared to the previous case.
This observation demonstrates the feasibility of using fractal-fractional derivatives to enhance the comprehension of the dynamics of a disease.In this part of the simulation, we explore effective strategies to control the spread of COVID-19.In this regard, we set the values of the fractal dimension and the fractional order to θ = 0.93 and γ = 0.92, respectively, and observe the influence of varying values of ς, ω V , ω S , and k on the infected populations.

A(t)
First, we examine the effect of lock-down on the transmission of the coronavirus by changing the contact rate ς.As shown in Figure 6, we note that restricting interactions between susceptible individuals and both asymptomatic and symptomatic infected individuals can decrease the number of future infection cases.Next, to demonstrate the influence of vaccination, we present the numerical results of the proposed model in Figure 7, both with and without the vaccination parameters.It is observed that there is a clear decrease in the number of infections when vaccine parameters are included.Figures 8-10 depict the influence of the k, ω V and ω S parameter values individually, in order to investigate the optimal approach for implementing vaccination within the community.
In Figure 7, we observe the impact of an increasing vaccination rate (ω V ) on the size of the infected groups.It becomes evident that as the vaccination rate rises, the number of infections decreases over time.
Figure 9 displays the evolution of the number of infections under improved vaccine effectiveness.It must be noted first that if k is close to 0, the vaccine is more effective.It is observed that as the vaccine's effectiveness increases, the number infected individuals decreases.
Figure 10 depicts the impact of the vaccine's waning immunity rate, and we see that the number of infections decreases significantly when we take smaller values of ω S .
Then, we can conclude that by increasing both the vaccination rate (ω V ) and vaccine efficacy (k) while decreasing the waning immunity rate (ω S ), we can achieve a rapid reduction in the number of infected cases within the population.
To validate this conclusion, Figure 11 presents numerical results obtained by increasing the vaccination rate from 0.0428 to 0.0488 and decreasing the waning immunity rate from 0.01 to 0.009.It is evident that the reduction in the number of infected cases occurs more rapidly when compared to the results shown in Figure 7.

Conclusions
This research presents a developed model of COVID-19 virus transmission.We used fractalfractional derivatives to create a suitable model of the third wave of COVID-19 in the KSA.The qualitative examination of the suggested model began with showing the existence and uniqueness of solutions as provided by the Perov fixed point theorem.Furthermore, we demonstrated the stability of the suggested system's solutions using Ulam-Hyers principles.In addition, the numerical approach for the model was provided using the two-step Newton polynomial approximation.Then, by adjusting the fractal dimension values, fractional order, contact rate, vaccination rate, and loss of immunity rate, we graphed the convergence and stability of the system's six solutions.Using our results, better responses will be implemented in the following waves, which will help eliminate the disease forever from society.
Our study did, however, have some limitations.First, the hospitalization class was not taken into account in this study.Additionally, the vaccination process typically consists of three stages: the first dose, the second dose, and the booster dose.Therefore, it is necessary to divide the vaccinated segment into three subclasses.Furthermore, while our study provides insights into the dynamics of infection rates in relation to vaccination strategies, we acknowledge its limitations.Primarily, it does not encompass the broader implications of vaccination, such as the impact on mortality rates.Further research is explore this critical aspect, especially in light of emerging data that suggest complex outcomes of vaccination programs in terms of public health.Future studies should aim to integrate comprehensive mortality rate analyses, comparing vaccinated and unvaccinated populations, to provide a more holistic understanding of vaccine efficacy and public health implications.

Figure 1 .
Figure 1.Flow diagram showing the transmission between different compartments.

Figure 3 .
Figure 3. Numerical results for model (6) under various values of the fractal dimension θ.

Figure 8 .
Figure 8.The impact of the vaccination rate on the components of the infected classes.

Figure 9 .
Figure 9.The impact of the vaccination efficacy rate on the components of the infected classes.

Figure 10 .
Figure 10.The impact of the loss of immunity rate on the components of the infected classes.

Table 1 .
Definitions of parameters used in model (4) and their values in numerical simulations.