Open AccessArticle

Analysis of the Solution of a Model of SARS-CoV-2 Variants and Its Approximation Using Two-Step Lagrange Polynomial and Euler Techniques

Muhammad Usman

Mujahid Abbas

^2,3

and

Andrew Omame

^1,4,*

Abdus Salam School of Mathematical Sciences, Government College University, Katchery Road, Lahore 54000, Pakistan

Department of Mathematics, Government College University, Katchery Road, Lahore 54000, Pakistan

Department of Medical Research, China Medical University Hospital, China Medical University, Taichung 40402, Taiwan

⁴

Department of Mathematics, Federal University of Technology, Owerri 460114, Nigeria

Author to whom correspondence should be addressed.

Axioms 2023, 12(5), 480; https://doi.org/10.3390/axioms12050480

Submission received: 7 March 2023 / Revised: 1 May 2023 / Accepted: 10 May 2023 / Published: 16 May 2023

(This article belongs to the Special Issue Dynamical Systems: Theory and Applications in Mathematical Biology)

Download

Browse Figures

Versions Notes

Abstract

In this paper, a vaccination model for SARS-CoV-2 variants is proposed and is studied using fractional differential operators involving a non-singular kernel. It is worth mentioning that variability in transmission rates occurs because of the particular population that is vaccinated, and hence, the asymptomatic infected classes are classified on the basis of their vaccination history. Using the Banach contraction principle and the Arzela–Ascoli theorem, existence and uniqueness results for the proposed model are presented. Two different numerical approaches, the fractional Euler and Lagrange polynomial methods, are employed to approximate the model’s solution. The model is then fitted to data associated with COVID-19 deaths in Pakistan between 1 January 2022 and 10 April 2022. It is concluded that our model is much aligned with the data when the order of the fractional derivative

ζ = 0.96

. The two different approaches are then compared with different step sizes. It is observed that they behave alike for small step sizes and exhibit different behaviour for larger step sizes. Based on the numerical assessment of the model presented herein, the impact of vaccination and the fractional order are highlighted. It is also noted that vaccination could remarkably decrease the spikes of different emerging variants of SARS-CoV-2 within the population.

Keywords:

SARS-CoV-2; delta variant; omicron variant; two-step Lagrange polynomial; fractional Euler; fractional derivative

MSC:

34C60; 92C42; 92D30; 92D25

1. Introduction

In recent years, different variants of “Severe Acute Respiratory Syndrome CoronaVirus 2 (SARS-CoV-2)” have appeared in different parts of the world [1]. The mutated variants have also been associated with an increase in SARS-CoV-2 cases [1]. In addition to the different preventive mechanisms against SARS-CoV-2, several vaccines have also been developed to fight the disease [2,3]. It was observed that many of the vaccines possess high effectiveness against emerging SARS-CoV-2 variants [3]. Recently, many epidemiological studies have investigated the efficacy of current SARS-CoV-2 vaccines against variants of concern (VOCs) [4,5,6,7,8,9]. For instance, Tang et al. [4] studied the effectiveness of the Pfizer and Moderna vaccines against the SARS-CoV-2 Delta variant in Qatar. They reported that those who are fully vaccinated with any of the two vaccines have at least 51.9% protection against the Delta variant. Similarly, on the basis of the studies conducted in the United Kingdom [5], Canada [6], the USA [7,8,9], and Israel [10], it was observed that the Pfizer and Moderna vaccines showed 39–75% effectiveness against the SARS-CoV-2 variants. Many countries have achieved more than 80% vaccination coverage; however, Pakistan is among the countries with the lowest record of complete vaccination coverage.

On the other hand, the role of mathematical modelling in proposing different viable solutions to governments and policymakers to fight against infectious diseases has become very significant. The classical mathematical models of the integer-order derivatives have been employed in studying infectious diseases [11,12,13,14,15,16]. For instance, Ali et al. [11] proposed an avian influenza model with asymptomatic carriers and analyzed it using a spectral method. Additionally, Bonyah et al. [12] investigated a mathematical model for dengue fever and Zika virus co-dynamics. Rihan and Alsakaji [13] designed a COVID-19 infection model with asymptomatic infected and interacting people. The authors [14] investigated a COVID-19 model with time delay. Most of these models have deficiencies and limitations, as their formulations are based on the classical integer-order derivative, which lacks the history of the disease. To overcome these problems, methodologies involving fractional differential operators with non-local and non-singular kernels are being studied.

In the study of epidemiological models, the use of fractional-order differential equations plays a vital role in obtaining results that are more realistic. Riemann–Liouville and Caputo introduced fractional operators that relied on a power-law kernel. These definitions have some limitations when applied to the study of biological processes [17,18,19]. To overcome these shortcomings, Caputo and Fabrizio (CF) [20] and Atangana and Baleanu (AB) [21] proposed operators whose kernels involve the exponential and Mittag–Leffler functions, respectively. These operators have been successfully applied in the modeling of different real-life phenomena [21]. Note that the Caputo–Fabrizio derivative has its own merits but fails to handle complex models, whereas the Atangana–Baleanu derivative works well for modelling complex real-life phenomena [22,23,24,25,26].

Although many methods have been developed to approximate the solutions of fractional differential equations [27,28,29], to the best of our knowledge, a comparative study of the behaviour of two-step Lagrange polynomial and fractional Euler schemes has not yet been considered. This paper aims to do so with the help of a well-formulated mathematical model for dual strains of SARS-CoV-2. The model is also investigated using real data and simulated to suggest various control mechanisms against SARS-CoV-2 variants. The AB derivative is applied to study the proposed model in this paper. We believe that our findings in this paper will provide better insights to policymakers and health agencies in overcoming not only SARS-CoV-2 variants but also other vaccine-preventable diseases in the future.

The main contributions of this paper are:

(i).: The consideration and analysis of a novel model for SARS-CoV-2 with two variants. The new model incorporates variability in transmission due to vaccination history;
(ii).: The study of necessary conditions for the existence and uniqueness of the solution of the model;
(iii).: The provision of proof of the Ulam–Hyers stability result;
(iv).: The evaluation of the fractional system numerically using the two-step Lagrange polynomial and fractional Euler methods;
(v).: Highlighting the impact of vaccination and the fractional-order derivative.

Definition 1.

Let

f \in H^{1} (a_{1}, a_{2})

with

a_{2} > a_{1},

and

ζ \in [0, 1]

. The Atangana–Baleanu fractional derivative of f of order ζ in the Caputo sense is defined by

_{a}^{A B C} D_{t}^{ζ} f (t) = \frac{F (ζ)}{(1 - ζ)} \int_{a_{1}}^{t} \frac{d f (ϱ)}{d ϱ} E_{ζ} [- ζ \frac{{(t - ϱ)}^{ζ}}{1 - ζ}] d ϱ,

where

F (ζ) = (1 - ζ) + \frac{ζ}{Γ (ζ)}

is a normalization function satisfying

F (0) = F (1) = 1

and

E_{ζ}

(.) is the Mittag–Leffler function given by

E_{ζ} (t) = \sum_{k = 0}^{\infty} \frac{t^{k}}{Γ (ζ k + 1)}, ζ > 0 .

Definition 2

([21]). The Atangana–Baleanu integral of f (in the Caputo sense) of order ζ is defined by

_{a_{1}}^{A B C} I_{t}^{ζ} f (t) = \frac{1 - ζ}{F (ζ)} f (t) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{a_{1}}^{t} f (ϱ) {(t - ϱ)}^{ζ - 1} d ϱ .

Definition 3

([21]). The Laplace transform of the Atangana–Baleanu fractional derivative of

f

of order ζ in the Caputo sense is given by

L \{_{a_{1}}^{A B C} D_{t}^{ζ} f (t)\} = F (ζ) \frac{s^{ζ} L {f (t)} - s^{ζ - 1} f (a_{1})}{s^{ζ} (1 - ζ) + ζ},

where

L

is the Laplace transform operator.

Definition 4

([30]). The Laplace transform of the Mittag–Leffler function

E_{ζ}

(.) is given by

L \{E_{ζ} (- θ t^{ζ})\} = \frac{s^{ζ - 1}}{s^{ζ} + θ}, w h e r e θ \in R .

2. Model Formulation

At any given time,

t,

the total human population

N

consists of the following mutually exclusive classes:

S (t)

: susceptible or infected individuals;

V (t)

: individuals vaccinated against SARS-CoV-2;

I_{1} (t)

: individuals infected with the Delta variant;

I_{2} (t)

: individuals with Omicron variants (asymptomatic);

I_{1}^{v} (t)

: individuals infected with the Delta variant;

I_{2}^{v} (t)

: individuals infected with the Omicron variant (asymptomatic) coming from the vaccinated class;

H_{1} (t)

: individuals infected with the Delta variant (symptomatic);

H_{2} (t)

: individuals with Omicron variants (symptomatic); and

R_{1} (t)

and

R_{2} (t)

: individuals who have recovered from the Delta and Omicron variants, respectively. Unvaccinated persons are recruited into the population at the rate

(1 - ξ) Ψ

, with

ξ

representing the vaccination rate. Individuals in an unvaccinated state contract the Delta SARS-CoV-2 variant at the rate

\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}

. They also acquire the Omicron variant at the rate

\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}

. Individuals in this class are vaccinated at the rate

γ

. Natural death occurs at the rate

ϑ

. When immunity is lost, individuals recovered from Delta or Omicron return to a susceptible state at the rate

η_{1}

η_{2}

. Vaccinated persons are recruited into the population at the rate

ξ Ψ

. Individuals in the vaccinated state could contract Delta variant at the rate

(1 - ψ) \frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}

, where

ψ

is the vaccine’s effectiveness against Delta variant. Individuals in this state could also contract an infection of the Omicron variant at the rate

(1 - ω) \frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}

, where

ω

is the vaccine’s effectiveness against the Omicron variant. The parameters describing the flows from one epidemiological state to the other are given in Table 1. In this model, strain 1 denotes the Delta variant, while strain 2 denotes the Omicron variant. Additionally, vaccinated susceptible individuals (also assumed to have completed two doses of the available vaccines) have a reduced rate of infection with both variants, although with higher efficacy against the Delta variant [4]. It is also assumed that vaccinated individuals who become infected with strain 2 have a reduced transmissibility rate relative to unvaccinated infected individuals. Furthermore, it is further assumed that immigrants into the population have completed their vaccination dosage.

The set of equations for the fractional order system is given by:

\begin{matrix} _{0}^{A B C} D_{t}^{ζ} S (t) & = (1 - ξ) Ψ - (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)} + \frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2})}{N} + γ + ϑ) S (t) + η_{1} R_{1} (t) + η_{2} R_{2} (t), \\ _{0}^{A B C} D_{t}^{ζ} V (t) & = ξ Ψ + γ S (t) - [(1 - ψ) (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) + (1 - ω) (\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) + ϑ] V (t), \\ _{0}^{A B C} D_{t}^{ζ} I_{1} (t) & = (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) S (t) - (α_{1} + ϱ_{1} + ϑ) I_{1} (t), \\ _{0}^{A B C} D_{t}^{ζ} I_{1}^{v} (t) & = (1 - ψ) (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) V (t) - (α_{2} + q_{1} + ϑ) I_{1}^{v} (t), \\ _{0}^{A B C} D_{t}^{ζ} H_{1} (t) & = α_{1} I_{1} (t) + α_{2} I_{1}^{v} (t) - (q_{2} + δ_{1} + ϑ) H_{1} (t), \\ _{0}^{A B C} D_{t}^{ζ} R_{1} (t) & = τ_{1} I_{1} (t) + q_{1} I_{1}^{v} (t) + q_{2} H_{1} (t) - (ϑ + η_{1}) R_{1} (t), \\ _{0}^{A B C} D_{t}^{ζ} I_{2} (t) & = (\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) S (t) - (α_{3} + ϱ_{2} + ϑ) I_{2} (t), \\ _{0}^{A B C} D_{t}^{ζ} I_{2}^{v} (t) & = (1 - ω) (\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) V (t) - (α_{4} + q_{3} + ϑ) I_{2}^{v} (t), \\ _{0}^{A B C} D_{t}^{ζ} H_{2} (t) & = α_{3} I_{2} (t) + α_{4} I_{2}^{v} (t) - (q_{4} + δ_{2} + ϑ) H_{2} (t), \\ _{0}^{A B C} D_{t}^{ζ} R_{2} (t) & = τ_{2} I_{2} (t) + q_{3} I_{2}^{v} (t) + q_{4} H_{2} (t) - (ϑ + η_{2}) R_{2} (t) . \end{matrix}

(1)

System (1) can be represented in a compact form given as:

\{\begin{matrix} _{0}^{A B C} D_{t}^{ζ} K (t) & = K (t, K (t)), \\ K (0) & = K_{0}, \end{matrix}

(2)

where the vector

K (t) = {(\begin{matrix} S (t) & V (t) & I_{1} (t) & I_{1}^{v} (t) & H_{1} (t) & R_{1} (t) & I_{2} (t) & I_{2}^{v} (t) & H_{2} (t) & R_{2} (t) \end{matrix})}^{T} \in R^{10}

, for

t \in J = [0, b]

. That is, we have the function

K : J \to R^{10}

. Additionally,

K : J \times R^{10} \to R^{10}

defines a function where

\begin{matrix} K_{1} (t, K (t)) & = (1 - ξ) Ψ - (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)} + \frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)} + γ + ϑ) S (t) \\ + η_{1} R_{1} (t) + η_{2} R_{2} (t), \\ K_{2} (t, K (t)) & = ξ Ψ + γ S (t) - [(1 - ψ) (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) + (1 - ω) (\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) + ϑ] V (t), \\ K_{3} (t, K (t)) & = (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) S (t) - (α_{1} + τ_{1} + ϑ) I_{1} (t), \\ K_{4} (t, K (t)) & = (1 - ψ) (\frac{β_{1} (I_{1} (t) + ϕ_{1 v} I_{1}^{v} (t) + θ_{1} H_{1} (t))}{N (t)}) V (t) - (α_{2} + q_{1} + ϑ) I_{1}^{v} (t), \\ K_{5} (t, K (t)) & = α_{1} I_{1} (t) + α_{2} I_{1}^{v} (t) - (q_{2} + δ_{1} + ϑ) H_{1} (t), \\ K_{6} (t, K (t)) & = τ_{1} I_{1} (t) + q_{1} I_{1}^{v} (t) + q_{2} H_{1} (t) - (ϑ + η_{1}) R_{1} (t), \\ K_{7} (t, K (t)) & = (\frac{β_{2} (I_{2} + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) S (t) - (α_{3} + τ_{2} + ϑ) I_{2} (t), \\ K_{8} (t, K (t)) & = (1 - ω) (\frac{β_{2} (I_{2} (t) + ϕ_{2 v} I_{2}^{v} (t) + θ_{2} H_{2} (t))}{N (t)}) V (t) - (α_{4} + q_{3} + ϑ) I_{2}^{v} (t), \\ K_{9} (t, K (t)) & = α_{3} I_{2} (t) + α_{4} I_{2}^{v} (t) - (q_{4} + δ_{2} + ϑ) H_{2} (t), \\ K_{10} (t, K (t)) & = τ_{2} I_{2} (t) + q_{3} I_{2}^{v} (t) + q_{4} H_{2} (t) - (ϑ + η_{2}) R_{2} (t) . \end{matrix}

(3)

System (3) can be transformed to the Volterra integral equation given by

K (t) = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t, K (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ .

(4)

Model’s Basic Properties

Theorem 1.

The closed set

\begin{matrix} D & = {(S (t), V (t), I_{1} (t), I_{1}^{v} (t), H_{1} (t), R_{1} (t), I_{2} (t), I_{2}^{v} (t), H_{2} (t), R_{2} (t)) \in R_{+}^{10} : \\ S (t) + V (t) + I_{1} (t) + I_{1}^{v} (t) + H_{1} (t) + R_{1} (t) + I_{2} (t) + I_{2}^{v} (t) + H_{2} (t) + R_{2} (t) \leq \frac{Ψ}{ϑ}} \end{matrix}

is positively invariant relative to System (1).

Proof.

If all equations in System (1) are added, then we obtain that

_{0}^{A B C} D_{t}^{ζ} N (t) = Ψ - ϑ N (t) - (δ_{1} H_{1} (t) + δ_{2} H_{2} (t)),

(5)

which could be written as

_{0}^{A B C} D_{t}^{ζ} N (t) \leq Ψ - ϑ N (t) .

(6)

Upon applying the Laplace transform of the AB derivative to (6), we obtain

F (ζ) \frac{s^{ζ} L \{N (t)\} - s^{ζ - 1} N (0)}{s^{ζ} (1 - ζ) + ζ} \leq L \{Ψ - ϑ N (t)\},

which can be written as

\begin{matrix} L \{N (t)\} & \leq \frac{F (ζ) s^{ζ} N (0) + [s^{ζ} (1 - ζ) + ζ] Ψ}{s {[F (ζ) + ϑ (1 - ζ)] s^{ζ} + ϑ ζ}} \\ = \frac{Ψ}{ϑ s} + \frac{F (ζ) N (0)}{λ} \frac{s^{ζ - 1}}{s^{ζ} + \frac{ϑ ζ}{λ}} + \frac{(1 - ζ) Ψ}{λ} \frac{s^{ζ - 1}}{s^{ζ} + \frac{ϑ ζ}{λ}} - \frac{Ψ}{ϑ} \frac{s^{ζ - 1}}{s^{ζ} + \frac{ϑ ζ}{λ}}, \end{matrix}

where

λ = [F (ζ) + ϑ (1 - ζ)] .

On taking the inverse Laplace transform of both sides of the above inequality, we have

N (t) \leq \frac{Ψ}{ϑ} + \frac{F (ζ) N (0)}{λ} E_{ζ} (- \frac{ϑ ζ}{λ} t^{ζ}) + \frac{(1 - ζ) Ψ}{λ} E_{ζ} (- \frac{ϑ ζ}{λ} t^{ζ}) - \frac{Ψ}{ϑ} E_{ζ} (- \frac{ϑ ζ}{λ} t^{ζ}) .

Due to the asymptotic nature of the Mittag–Leffler function [30], we have

N (t) \leq \frac{Ψ}{ϑ}

t \to \infty

. Thus, System (1) is positively invariant. The schematic diagram is given in Figure 1. □

3. Existence and Uniqueness of The Solution

3.1. Existence

In this section, we will study the necessary conditions for the existence of solutions of the model. Consider a Banach space

E = C [J, R^{10}]

equipped with the norm

∥ Φ ∥ =_{t \in J}^{sup} | Φ (t) |,

where,

| Φ (t) | = | Φ_{1} (t) | + | Φ_{2} (t) | + | Φ_{3} (t) | + | Φ_{4} (t) | + | Φ_{5} (t) | + | Φ_{6} (t) | + | Φ_{7} (t) | + | Φ_{8} (t) |

+ | Φ_{9} (t) | + | Φ_{10} (t) |

. Norms on

C ([J, R^{10}])

C ([J, R])

will be clear from the context.

Theorem 2

([34]). Let M be a non-empty closed, bounded, and convex subset in a Banach space

E = C ([J, R^{10}])

. If

P_{1}, P_{2} : M \to E

are two operators satisfying the following:

(i).: $P_{1} Φ_{1} + P_{2} Φ_{2} \in M$ , whenever $Φ_{1}, Φ_{2} \in M$ ;
(ii).: $P_{2}$ is a contraction;
(iii).: $P_{1}$ is compact and continuous;

then there exists

Φ \in M

such that

Φ = P_{1} Φ + P_{2} Φ

Theorem 3.

K : J \times R^{10} \to R^{10}

is continuous with

| K (t, Φ) | \leq | Ψ (t) |

and satisfies

|K (t_{1}, Φ (t_{1})) - K (t_{2}, Φ (t_{2}))| \leq L_{M} |t_{1} - t_{2}|

for all

(t, Φ) \in J \times R^{10} a n d Ψ \in C (J, R_{+})

with

∥ Ψ ∥ = {sup}_{t \in J} | Ψ (t) |

, then the proposed Model (4) has at least one solution.

Proof.

Consider

B_{η} = {Φ \in E : ∥ Φ ∥ \leq η}

, where

η \geq ∥Φ_{0}∥ + Ω ∥ Ψ ∥

Φ_{0} \in R^{10}

and

Ω = [\frac{1 - ζ}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}]

. Obviously,

B_{η}

is a closed convex and bounded subset of E. Define operators

P_{1}, P_{2} : B_{η} \to E

(P_{1} Φ) (t) = \frac{1 - ζ}{F (ζ)} K (t, Φ (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ, \forall t \in J,

(P_{2} Φ) (t) = Φ_{0} (t), \forall t \in J,

respectively. According to the given assumptions,

K : J \times R^{10} \to R^{10}

is continuous and satisfies the condition

K (t, Φ (t)) \leq | Ψ (t) |,

for each

t \in J

and

Φ \in R^{10}

. That is,

K : J \times R^{10} \to R^{10}

is point-wise bounded.

Now, for any

Φ_{1}, Φ_{2} \in B_{η},

\begin{matrix} ∥ (P_{1} Φ_{1}) (t) + (P_{2} Φ_{2}) (t) ∥ =_{t \in J}^{sup} | P_{1} Φ_{1} (t) + P_{2} Φ_{} (t) | \\ =_{t \in J}^{sup} [| Φ_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, Φ_{2} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ_{2} (ϱ)) d ϱ |] \\ \leq_{t \in J}^{sup} [| Φ_{0} (t) | + \frac{1 - ζ}{F (ζ)} | Ψ (t) | + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} | Ψ (ϱ) | d ϱ] \\ \leq ∥ Φ_{0} ∥ + \frac{1 - ζ}{F (ζ)} ∥ Ψ ∥ + \frac{ζ ∥ Ψ ∥}{F (ζ) Γ (ζ)} sup_{t \in J} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} d ϱ \\ \leq ∥ Φ_{0} ∥ + [\frac{1 - ζ}{F (ζ)}] ∥ Ψ ∥ + \frac{b^{ζ}}{F (ζ) Γ (ζ)} ∥ Ψ ∥ \\ = ∥ Φ_{0} ∥ + [\frac{1 - ζ}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}] ∥ Ψ ∥ \\ = ∥ Φ_{0} ∥ + Ω ∥ Ψ ∥ \leq η . \end{matrix}

Hence,

P_{1} Φ_{1} + P_{2} Φ_{2} \in B_{η}

Clearly,

P_{2}

is a contraction since it is a constant operator.

P_{1}

is also continuous because

K

is continuous.

Now, for any

Φ \in B_{η}

, we have

\begin{matrix} ∥ (P_{1} Φ) (t) ∥ =_{t \in J}^{sup} | P_{1} Φ (t) | \\ =_{t \in J}^{sup} [| \frac{1 - ζ}{F (ζ)} K (t, Φ (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ |] \\ \leq_{t \in J}^{sup} [\frac{1 - ζ}{F (ζ)} | Ψ (t) | + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} | Ψ (ϱ) | d ϱ] \\ \leq \frac{1 - ζ}{F (ζ)} ∥ Ψ ∥ + \frac{ζ ∥ Ψ ∥}{F (ζ) Γ (ζ)} sup_{t \in J} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} d ϱ \\ \leq [\frac{1 - ζ}{F (ζ)}] ∥ Ψ ∥ + \frac{b^{ζ}}{F (ζ) Γ (ζ)} ∥ Ψ ∥ \\ = [\frac{1 - ζ}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}] ∥ Ψ ∥ \\ = Ω ∥ Ψ ∥ \leq η . \end{matrix}

Thus,

P_{1} (B_{η}) \subset B_{η}

\bar{P_{1} (B_{η})}

is bounded and closed. In order to apply the Arzela–Ascoli theorem, it suffices to prove that

\bar{P_{1} (B_{η})}

is equicontinuous.

Now, for any

Φ \in B_{η}

, we have

\begin{matrix} | (P_{1} Φ) (t_{2}) - (P_{1} Φ) (t_{1}) | = | \frac{1 - ζ}{F (ζ)} K (t_{2}, Φ (t_{2})) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t_{2}} {(t_{2} - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ \\ - \frac{1 - ζ}{F (ζ)} K (t_{1}, Φ (t_{1})) - \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t_{1}} {(t_{1} - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ | \\ = \frac{1 - ζ}{F (ζ)} | K (t_{2}, Φ (t_{2})) - K (t_{1}, Φ (t_{1})) | + \frac{ζ}{F (ζ) Γ (ζ)} |\int_{0}^{t_{1}} [{(t_{2} - ϱ)}^{ζ - 1} - {(t_{1} - ϱ)}^{ζ - 1}] K (ϱ, Φ (ϱ)) d ϱ \\ + \int_{t_{1}}^{t_{2}} {(t_{2} - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ| \\ \leq \frac{1 - ζ}{F (ζ)} L_{M} | t_{2} - t_{1} | + \frac{ζ}{F (ζ) Γ (ζ)} | \int_{0}^{t_{1}} [{(t_{2} - ϱ)}^{ζ - 1} - {(t_{1} - ϱ)}^{ζ - 1}] K (ϱ, Φ (ϱ)) d ϱ \\ + \int_{t_{1}}^{t_{2}} {(t_{2} - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ | \\ \leq \frac{1 - ζ}{F (ζ)} L_{M} | t_{2} - t_{1} | + \frac{∥ Ψ ∥}{F (ζ) Γ (ζ)} [(t_{2}^{ζ} - t_{1}^{ζ})] . \end{matrix}

(7)

Clearly, the right-hand side of the above inequality vanishes as

t_{2} \to t_{1}

. Thus,

P_{1} B_{η}

is equicontinuous, and so is

\bar{P_{1} (B_{η})}

. Thus,

\bar{P_{1} (B_{η})}

, being closed, bounded, and equicontinuous, is compact, which implies that

P_{1}

is a compact operator. Hence, all the conditions of Theorem 2 are satisfied. Thus, there exists

Φ

E

such that

Φ = P_{1} Φ + P_{2} Φ

. That is,

Φ = Φ_{0} + \frac{1 - ζ}{F (ζ)} K (t, Φ) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ .

□

3.2. Uniqueness

Theorem 4.

K \in C ([J, R^{10}])

satisfy the Lipschitz condition

|K (t, Φ_{1} (t)) - K (t, Φ_{2} (t))| \leq L_{K} |Φ_{1} (t) - Φ_{2} (t)|,

(8)

for all

t \in

J

and each

Φ_{1}, Φ_{2} \in E

L_{K} > 0

, then an integral equation, (4), which is equivalent to System (1), has a unique solution provided that

Ω L_{K} < 1 .

Proof.

Define an operator

P : E \to E

(P Φ) (t) = Φ_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, Φ) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} K (ϱ, Φ) {(t - ϱ)}^{ζ - 1} d ϱ .

Now, for any

Φ_{1}, Φ_{2} \in E

, we obtain

\begin{matrix} ∥ (P Φ_{1}) (t) - (P Φ_{2}) (t) ∥ \leq sup_{t \in J} [| Φ_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, Φ_{1} (t)) \\ + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ_{1} (ϱ)) d ϱ \\ - (Φ_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, Φ_{2} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ_{2} (ϱ)) d ϱ) |] \\ \leq sup_{t \in J} \frac{1 - ζ}{F (ζ)} [| K (t, Φ_{1} (t)) - K (t, Φ_{2} (t)) |] + sup_{t \in J} \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} | K (ϱ, Φ_{1} (ϱ)) \\ - K (ϱ, Φ_{2} (ϱ)) | d ϱ \\ \leq & sup_{t \in J} \frac{1 - ζ}{F (ζ)} [L_{K} | Φ_{1} (t) - Φ_{2} (t) |] + sup_{t \in J} \frac{L_{K} ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} | Φ_{1} (ϱ) - Φ_{2} (ϱ) | d ϱ \\ \leq & \frac{1 - ζ}{F (ζ)} [L_{K} ∥ Φ_{1} - Φ_{2} ∥] + \frac{L_{K} ζ ∥ Φ_{1} - Φ_{2} ∥}{F (ζ) Γ (ζ)} sup_{t \in J} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} d ϱ \\ \leq & [\frac{1 - ζ}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}] L_{K} ∥ Φ_{1} - Φ_{2} ∥ \\ = & Ω L_{K} ∥ Φ_{1} (t) - Φ_{2} (t) ∥ . \end{matrix}

This implies that P is a contraction. As

P (Φ (t)) = P_{1} (Φ (t)) + P_{2} (Φ (t))

P B_{η} \subset B_{η}

. Additionally, the set

B_{η}

is closed and convex. The results follow from the Banach contraction principle. □

3.3. The Basic Reproduction Number of the Model

The model’s disease-free equilibrium (DFE) is given by

\begin{matrix} ψ_{0} & = (S^{*}, V^{*}, I_{1}^{*}, I_{1}^{v *}, H_{1}^{*}, R_{1}^{*}, I_{2}^{*}, I_{2}^{v *}, H_{2}^{*}, R_{2}^{*}) \\ = (\frac{(1 - ξ) Ψ}{γ + ϑ}, \frac{ξ Ψ (γ + ϑ) + γ Ψ (1 - ξ)}{ϑ (γ + ϑ)}, 0, 0, 0, 0, 0, 0, 0, 0) . \end{matrix}

Using the next-generation operator method [35], the transfer matrices associated with the model are

F = (\begin{matrix} \frac{β_{1} S^{*}}{N} & \frac{β_{1} ϕ_{1 v} S^{*}}{N} & \frac{β_{1} θ_{1} S^{*}}{N} & 0 & 0 & 0 \\ \frac{(1 - ψ) β_{1} V^{*}}{N} & \frac{(1 - ψ) β_{1} ϕ_{1 v} V^{*}}{N} & \frac{(1 - ψ) β_{1} θ_{1} V^{*}}{N} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{β_{2} S^{*}}{N} & \frac{β_{2} ϕ_{2 v} S^{*}}{N} & \frac{β_{2} θ_{2} S^{*}}{N} \\ 0 & 0 & 0 & \frac{(1 - ω) β_{2} V^{*}}{N} & \frac{(1 - ω) β_{2} ϕ_{2 v} V^{*}}{N} & \frac{(1 - ω) β_{2} θ_{2} V^{*}}{N} \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}),

(9)

V = (\begin{matrix} L_{1} & 0 & 0 & 0 & 0 & 0 \\ 0 & L_{2} & 0 & 0 & 0 & 0 \\ - α_{1} & - α_{2} & L_{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & L_{4} & 0 & 0 \\ 0 & 0 & 0 & 0 & L_{5} & 0 \\ 0 & 0 & 0 & - α_{3} & - α_{4} & L_{6} \end{matrix}),

(10)

where

\begin{matrix} L_{1} & = α_{1} + τ_{1} + ϑ, L_{2} = α_{2} + q_{1} + ϑ, L_{3} = q_{2} + δ_{1} + ϑ, \\ L_{4} & = α_{3} + τ_{2} + ϑ, L_{5} = α_{4} + q_{3} + ϑ, L_{6} = q_{4} + δ_{2} + ϑ . \end{matrix}

The model’s reproduction number is defined by

R_{0} = ρ (F V^{- 1}) = \max {R_{0 W}, R_{0 M}}

, where

R_{0 D}

and

R_{0 M}

are reproduction numbers for the Delta and Omicron variants, respectively, and are given by

R_{0 D} = \frac{β_{1} (L_{3} + α_{1} θ_{1}) L_{2} S^{*} + (1 - ψ) β_{1} (L_{3} ϕ_{1 v} + α_{2} θ_{1}) L_{1} V^{*}}{L_{1} L_{2} L_{3} N^{*}},

R_{0 M} = \frac{β_{2} (L_{6} + α_{3} θ_{2}) L_{5} S^{*} + (1 - ω) β_{2} (L_{6} ϕ_{2 v} + α_{4} θ_{2}) L_{4} V^{*}}{L_{4} L_{5} L_{6} N^{*}} .

3.4. Local Asymptotic Stability of the Disease-Free Equilibrium (DFE) of the Model

Theorem 5.

The DFE,

Z_{0}

, of System (1) is locally asymptotically stable (LAS) if

R_{0} < 1

and unstable as long as

R_{0} > 1

Proof.

We shall analyze the local stability of Model (1) by evaluating the Jacobian at the DFE,

Z_{0}

, given by:

(\begin{matrix} - Q_{1} & 0 & - \frac{β_{1} S^{*}}{N} & - \frac{β_{1} ϕ_{1 v} S^{*}}{N} & - \frac{β_{1} θ_{1} S^{*}}{N} & η_{1} & - \frac{β_{2} S^{*}}{N} & - \frac{β_{2} ϕ_{2 v} S^{*}}{N} & - \frac{β_{2} θ_{2} S^{*}}{N} & η_{2} \\ γ & - ϑ & - \frac{β_{1} (1 - ψ) V^{*}}{N} & - \frac{β_{1} (1 - ψ) ϕ_{1 v} V^{*}}{N} & - \frac{β_{1} θ_{1} (1 - ψ) V^{*}}{N} & 0 & - \frac{β_{2} (1 - ω) V^{*}}{N} & - \frac{β_{2} ϕ_{2 v} (1 - ω) V^{*}}{N} & - \frac{β_{2} θ_{2} (1 - ω) V^{*}}{N} & 0 \\ 0 & 0 & \frac{β_{1} S^{*}}{N} - L_{1} & \frac{β_{1} ϕ_{1 v} S^{*}}{N} & \frac{β_{1} θ_{1} S^{*}}{N} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{(1 - ψ) β_{1} V^{*}}{N} & \frac{(1 - ψ) β_{1} ϕ_{1 v} V^{*}}{N} - L_{2} & \frac{(1 - ψ) β_{1} θ_{1} V^{*}}{N} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & α_{1} & α_{2} & L_{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & τ_{1} & q_{1} & q_{2} & - Q_{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{β_{2} S^{*}}{N} - L_{4} & \frac{β_{2} ϕ_{2 v} S^{*}}{N} & \frac{β_{2} θ_{2} S^{*}}{N} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{β_{2} (1 - ω) V^{*}}{N} & \frac{β_{2} ϕ_{2 v} (1 - ω) V^{*}}{N} - L_{5} & \frac{β_{2} θ_{2} (1 - ω) V^{*}}{N} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & α_{3} & α_{4} & - L_{6} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & τ_{2} & q_{3} & q_{4} & - Q_{3} \end{matrix})

where,

Q_{1} = (γ + ϑ), Q_{2} = ϑ + η_{1}, Q_{3} = (ϑ + η_{2})

The first four eigenvalues are given by

φ_{1} = - ϑ, φ_{2} = - (ϑ + η_{1}), φ_{3} = - (ϑ + η_{2}), φ_{4} = - (ϑ + γ),

whereas the remaining six eigenvalues are the solutions of the characteristic polynomial equations and are given as follows:

φ^{3} + ξ_{11} φ^{2} + ξ_{22} φ + ξ_{33} = 0,

(11)

and

φ^{3} + γ_{11} φ^{2} + γ_{22} φ + γ_{33} = 0,

(12)

where

\begin{matrix} ξ_{11} & = (L_{1} + L_{2} + L_{3} - \frac{β_{1} [S^{*} + (1 - ψ) ϕ_{1 v} V^{*}]}{N^{*}}), \\ ξ_{22} & = (L_{1} L_{2} + L_{1} L_{3} + L_{2} L_{3} - \frac{β_{1} (L_{2} + L_{3} + α_{1} θ_{1}) S^{*}}{N^{*}} - \frac{β_{1} (1 - ψ) (α_{2} θ_{1} + L_{1} ϕ_{1 v} + L_{3} ϕ_{1 v}) V^{*}}{N^{*}}), \\ ξ_{33} & = L_{1} L_{2} L_{3} (1 - R_{0 D}), \\ γ_{11} & = (L_{4} + L_{5} + L_{6} - \frac{β_{2} [S^{*} + (1 - ω) ϕ_{2 v} V^{*}]}{N^{*}}), \\ γ_{22} & = (L_{4} L_{5} + L_{4} L_{6} + L_{5} L_{6} - \frac{β_{2} (L_{5} + L_{6} + α_{2} θ_{2}) S^{*}}{N^{*}} - \frac{β_{2} (1 - ω) (α_{4} θ_{2} + L_{4} ϕ_{2 v} + L_{6} ϕ_{2 v}) V^{*}}{N^{*}}), \\ γ_{33} & = L_{4} L_{5} L_{6} (1 - R_{0 M}) . \end{matrix}

From the Routh–Hurwitz criterion, Equations (11) and (12) have roots with negative real parts if and only if

ξ_{11} > 0

ξ_{22} > 0

ξ_{33} > 0

ξ_{11} ξ_{22} > ξ_{33}

and

γ_{11} > 0

γ_{22} > 0

γ_{33} > 0

γ_{11} γ_{22} > γ_{33}

, provided that

R_{0 D} < 1

and

R_{0 M} < 1

. Hence, the DFE,

Z_{0}

, is locally asymptotically stable if

R_{0} = \max {R_{0 D}, R_{0 M}} < 1

. □

This could also be verified by Theorem 2 in [35].

4. Ulam–Hyers Stability

The system’s Ulam–Hyers (UH) stability and generalized UH stability [36,37] are discussed in this section.

Let

E = C (J, R^{10})

be the space of continuous functions from

J

R^{10}

endowed with the norm

∥ Φ ∥ =_{t \in J}^{s u p} |Φ (t)|

, where

J = [0, b]

Definition 5.

The transformed system given below

\{\begin{matrix} _{0}^{A B C} D_{t}^{ζ} Φ (t) & = K (t, Φ (t)), \\ Φ (0) & = Φ_{0}, \end{matrix}

(13)

is UH stable if there exists

k > 0

such that for

ϵ > 0

and

\bar{Φ} \in E

with the following inequality

∥^{A B C} D^{ζ} \bar{Φ} (t) - K (t, \bar{Φ} (t) ∥ \leq ϵ, t \in J, ϵ = max {(ϵ_{i})}^{T}, i = 1, 2, \dots 10,

(14)

there exists a unique solution

Φ \in E

of the fractional system (13) such that

∥ \bar{Φ} (t) - Φ (t) ∥ \leq k ϵ, t \in J, k = max {(k_{j})}^{T}, j = 1, 2, \dots 10,

where,

\begin{matrix} \bar{Φ} (t) & = {(\begin{matrix} \bar{S} (t) & \bar{V} (t) & {\bar{I}}_{1} (t) & {\bar{I}}_{1}^{v} (t) & {\bar{H}}_{1} (t) & {\bar{R}}_{1} (t) & {\bar{I}}_{2} (t) & {\bar{I}}_{2}^{v} (t) & {\bar{H}}_{2} (t) & {\bar{R}}_{2} (t) \end{matrix})}^{T}, \\ Φ (t) & = {(\begin{matrix} S (t) & V (t) & I_{1} (t) & I_{1}^{v} (t) & H_{1} (t) & R_{1} (t) & I_{2} (t) & I_{2}^{v} (t) & H_{2} (t) & R_{2} (t) \end{matrix})}^{T}, \\ Φ (0) & = {(\begin{matrix} S (0) & V (0) & I_{1} (0) & I_{1}^{v} (0) & H_{1} (0) & R_{1} (0) & I_{2} (0) & I_{2}^{v} (0) & H_{2} (0) & R_{2} (0) \end{matrix})}^{T} . \end{matrix}

Definition 6.

The model system (13) is generalized UH stable if there exists a continuous function

ϕ : R^{+} \to R^{+}

with

ϕ (0) = 0

such that for any approximate solution

\bar{Φ} \in E

of System (14), there exists a unique solution

Φ \in E

that satisfies the following:

∥ \bar{Φ} (t) - Φ (t) ∥ \leq ϕ (ϵ), t \in J, ϕ = max {(ϕ_{j})}^{T}, j = 1, 2, \dots, 10 .

Remark 1.

A function

\bar{Φ} \in E

satisfies the inequality (14) if and only if there exists a function

h \in E

with the following properties:

(i).: $∥h (t)∥ \leq ϵ, t \in J$ .
(ii).: $^{A B C} D^{ζ} \bar{Φ} (t) = K (t, \bar{Φ} (t) + h (t)$ , $t \in J$ .

Lemma 1.

For

\bar{Φ} \in E

in the inequality (14), we have the following:

|\bar{Φ} (t) - ({\bar{Φ}}_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, \bar{Φ} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, \bar{Φ} (ϱ)) d ϱ)| \leq Ω ϵ .

(15)

Proof.

Using (ii) of Remark 1, we have

^{A B C} D^{ζ} \bar{Φ} (t) = K (t, \bar{Φ} (t)) + h (t)

t \in J

Applying the ABC integral, we obtain that

\begin{matrix} \bar{Φ} (t) & = {\bar{Φ}}_{0} + \frac{1 - ζ}{F (ζ)} K (t, \bar{Φ} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, \bar{Φ} (ϱ)) d ϱ + \frac{1 - ζ}{F (ζ)} h (t) \\ + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} h (ϱ) d ϱ . \end{matrix}

(16)

Re-arranging this, applying the norm on both sides and using (i) of Remark 1, we have

\begin{matrix} |\bar{Φ} (t) - ({\bar{Φ}}_{0} (t) + \frac{1 - ζ}{F (ζ)} K (t, \bar{Φ} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, \bar{Φ} (ϱ)) d ϱ)| \\ \leq \frac{1 - ζ}{F (ζ)} |h (t)| + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} |h (ϱ)| d ϱ \\ \leq (\frac{(1 - ζ)}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}) ϵ \leq Ω ϵ . \end{matrix}

□

Theorem 6.

Let

Φ \in E

and

K : J \times R^{10} \to R^{10}

satisfy the Lipschitz condition with the Lipschitz constant

L_{K} > 0

and

1 - Ω L_{K} > 0

, where

Ω = \frac{(1 - ζ)}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}

. Then, Model (13) is generalized UH stable.

Proof.

Suppose that

\bar{Φ} \in E

satisfies the inequality in (14) and

Φ \in E

is a unique solution of (13). Then,

\forall ϵ > 0, t \in J

together with Lemma 1, we have

\begin{matrix} ∥\bar{Φ} (t) - Φ (t)∥ =_{t \in J}^{sup} |{\bar{Φ}}_{0} + \frac{1 - ζ}{F (ζ)} K (t, \bar{Φ} (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, \bar{Φ} (ϱ)) d ϱ \\ + \frac{1 - ζ}{F (ζ)} h (t) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} h (ϱ) d ϱ \\ - (Φ_{0} + \frac{1 - ζ}{F (ζ)} K (t, Φ (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, Φ (ϱ)) d ϱ)| \\ \leq_{t \in J}^{sup} |{\bar{Φ}}_{0} - Φ_{0}| +_{t \in J}^{sup} |(\frac{1 - ζ}{F (ζ)} + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} d ϱ) h (t)| \\ +_{t \in J}^{sup} \frac{1 - ζ}{F (ζ)} ∥K (t, \bar{Φ} (t)) - K (t, Φ (t))∥ +_{t \in J}^{sup} \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} ∥K (t, \bar{Φ} (t)) - K (t, Φ (t))∥ d ϱ \\ \leq Ω ϵ + \frac{(1 - ζ) L_{K}}{F (ζ)} ∥\bar{Φ} - Φ∥ + \frac{ζ L_{K}}{F (ζ) Γ (ζ)}_{t \in J}^{sup} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} ∥\bar{Φ} - Φ∥ d ϱ \\ \leq Ω ϵ + (\frac{(1 - ζ)}{F (ζ)} + \frac{b^{ζ}}{F (ζ) Γ (ζ)}) L_{K} ∥\bar{Φ} - Φ∥ \\ = Ω ϵ + Ω L_{K} ∥\bar{Φ} (t) - Φ (t)∥ . \end{matrix}

Thus,

∥ \bar{Φ} - Φ ∥ \leq k ϵ,

(17)

where

k = \frac{Ω}{1 - Ω L_{K}}

Taking

ϕ (ϵ) = k ϵ

, we conclude that the model (13) is both UH and generalized UH stable. □

5. Numerical Scheme

5.1. Two-Step Lagrange Polynomial Method

We follow the numerical scheme similar to the one given in [38]. Applying the fundamental theorem of fractional calculus, the fractional differential equation given in (4) is transformed to an integral equation given by

K (t) = K_{0} + \frac{1 - ζ}{F (ζ)} K (t, K (t)) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t} {(t - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ,

where the vector

K (t) = {(\begin{matrix} S (t) & V (t) & I_{1} (t) & I_{1}^{v} (t) & H_{1} (t) & R_{1} (t) & I_{2} (t) & I_{2}^{v} (t) & H_{2} (t) & R_{2} (t) \end{matrix})}^{T} \in R^{10}

, for each

t \in [0, b]

At a given point,

t = t_{σ + 1} = (σ + 1) h

, where

h = t_{σ + 1} - t_{σ}

is the time step size and

σ = 0, 1, 2 \dots

. The above equation discretizes to

\begin{matrix} K (t_{σ + 1}) & = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t_{σ + 1}} {(t_{σ + 1} - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ, \\ K (t_{σ + 1}) & = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ) Γ (ζ)} \sum_{ϵ = 0}^{σ} \int_{t_{ϵ}}^{t_{ϵ + 1}} {(t_{σ + 1} - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ . \end{matrix}

(18)

A function

K (ϱ, K (ϱ))

[t_{ϵ}, t_{ϵ + 1}]

can be approximated with the help of the Lagrange two-point interpolation as follows:

\begin{matrix} K (ϱ, K (ϱ)) & ≊ \frac{ϱ - t_{ϵ - 1}}{t_{ϵ} - t_{ϵ - 1}} K (t_{ϵ}, K (t_{ϵ})) - \frac{ϱ - t_{ϵ}}{t_{ϵ} - t_{ϵ - 1}} K (t_{ϵ - 1}, K (t_{ϵ - 1})) \\ = \frac{K (t_{ϵ}, K (t_{ϵ}))}{h} (ϱ - t_{ϵ - 1}) - \frac{K (t_{ϵ - 1}, K (t_{ϵ - 1}))}{h} (ϱ - t_{ϵ}) . \end{matrix}

(19)

Substituting (19) into (18), we obtain

\begin{matrix} K (t_{σ + 1}) & = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ) Γ (ζ)} \sum_{ϵ = 0}^{σ} (\frac{K (t_{ϵ}, K (t_{ϵ}))}{h} \int_{t_{ϵ}}^{t_{ϵ + 1}} (ϱ - t_{ϵ - 1}) {(t_{σ + 1} - ϱ)}^{ζ - 1} d ϱ \\ - \frac{K (t_{ϵ - 1}, K (t_{ϵ - 1}))}{h} \int_{t_{ϵ}}^{t_{ϵ + 1}} (ϱ - t_{ϵ}) {(t_{σ + 1} - ϱ)}^{ζ - 1} d ϱ) . \end{matrix}

(20)

Upon integrating, we have

\begin{matrix} K (t_{σ + 1}) & = & K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

(21)

Adopting the numerical scheme (21) into the fractional system (4) yields the following numerical solution:

\begin{matrix} S (t_{σ + 1}) & = & S (0) + \frac{(1 - ζ)}{F (ζ)} K_{1} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{1} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{1} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} V (t_{σ + 1}) & = & V (0) + \frac{(1 - ζ)}{F (ζ)} K_{2} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{2} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{2} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} I_{1} (t_{σ + 1}) & = & I_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{3} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{3} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{3} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} I_{1}^{v} (t_{σ + 1}) & = & I_{1}^{v} (0) + \frac{(1 - ζ)}{F (ζ)} K_{4} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{4} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{3} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} H_{1} (t_{σ + 1}) & = & H_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{5} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{5} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{5} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} R_{1} (t_{σ + 1}) & = & R_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{6} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{6} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{6} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

(22)

\begin{matrix} I_{2} (t_{σ + 1}) & = & I_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{7} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{7} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{7} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} I_{2}^{v} (t_{σ + 1}) & = & I_{2}^{v} (0) + \frac{(1 - ζ)}{F (ζ)} K_{8} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{8} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{8} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} H_{2} (t_{σ + 1}) & = & H_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{9} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{9} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{9} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

\begin{matrix} R_{2} (t_{σ + 1}) & = & R_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{10} (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ)} \sum_{ϵ = 0}^{σ} \{\frac{h^{ζ} K_{10} (t_{ϵ}, K (t_{ϵ}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ} (σ - ϵ + ζ + 2) \\ - & {(σ - ϵ)}^{ζ} (σ - ϵ + 2 ζ + 2)] - \frac{h^{ζ} K_{10} (t_{ϵ - 1}, K (t_{ϵ - 1}))}{Γ (ζ + 2)} [{(σ - ϵ + 1)}^{ζ + 1} - {(σ - ϵ)}^{ζ} (σ - ϵ + ζ + 1)]\} . \end{matrix}

The error analysis, stability, and accuracy of this numerical scheme have been well studied in [38].

5.2. Fractional Euler Method

At a given time

t = t_{σ + 1} = (σ + 1) h

, where

h = t_{σ + 1} - t_{σ}

is the time step size and

σ = 0, 1, 2 \dots

, discretizing the Volterra integral (4), with the help of the Euler method [39], we move as follows:

\begin{matrix} K (t_{σ + 1}) & = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ) Γ (ζ)} \int_{0}^{t_{σ + 1}} {(t_{σ + 1} - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ \\ K (t_{σ + 1}) & = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ}{F (ζ) Γ (ζ)} \sum_{ϵ = 0}^{σ} \int_{t_{ϵ}}^{t_{ϵ + 1}} {(t_{σ + 1} - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ . \end{matrix}

(23)

With the help of the product rectangle rule [40],

\sum_{ϵ = 0}^{σ} \int_{t_{ϵ}}^{t_{ϵ + 1}} {(t_{σ + 1} - ϱ)}^{ζ - 1} K (ϱ, K (ϱ)) d ϱ = \frac{h^{ζ}}{ζ} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K (t_{ϵ}, K (t_{ϵ})) .

Thus,

K (t_{σ + 1}) = K (0) + \frac{(1 - ζ)}{F (ζ)} K (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K (t_{ϵ}, K (t_{ϵ})) .

(24)

Adopting the numerical scheme (24) into System (4) gives the following numerical solution:

\begin{matrix} S (t_{σ + 1}) & = S (0) + \frac{(1 - ζ)}{F (ζ)} K_{1} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{1} (t_{ϵ}, K (t_{ϵ})) . \\ V (t_{σ + 1}) & = V (0) + \frac{(1 - ζ)}{F (ζ)} K_{2} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{2} (t_{ϵ}, K (t_{ϵ})) . \\ I_{1} (t_{σ + 1}) & = I_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{3} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{3} (t_{ϵ}, K (t_{ϵ})) . \\ I_{1}^{v} (t_{σ + 1}) & = I_{1}^{v} (0) + \frac{(1 - ζ)}{F (ζ)} K_{4} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{4} (t_{ϵ}, K (t_{ϵ})) . \\ H_{1} (t_{σ + 1}) & = H_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{5} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{5} (t_{ϵ}, K (t_{ϵ})) . \\ R_{1} (t_{σ + 1}) & = R_{1} (0) + \frac{(1 - ζ)}{F (ζ)} K_{6} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{6} (t_{ϵ}, K (t_{ϵ})) . \\ I_{2} (t_{σ + 1}) & = I_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{7} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{7} (t_{ϵ}, K (t_{ϵ})) . \\ I_{2}^{v} (t_{σ + 1}) & = I_{2}^{v} (0) + \frac{(1 - ζ)}{F (ζ)} K_{8} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{8} (t_{ϵ}, K (t_{ϵ})) . \\ H_{2} (t_{σ + 1}) & = H_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{9} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{9} (t_{ϵ}, K (t_{ϵ})) . \\ R_{2} (t_{σ + 1}) & = R_{2} (0) + \frac{(1 - ζ)}{F (ζ)} K_{10} (t_{σ}, K (t_{σ})) + \frac{ζ h^{ζ}}{F (ζ) Γ (ζ + 1)} \sum_{ϵ = 0}^{σ} [{(σ - ϵ + 1)}^{ζ} - {(σ - ϵ)}^{ζ}] K_{10} (t_{ϵ}, K (t_{ϵ})) . \end{matrix}

(25)

The error estimates, stability analysis, and the high accuracy of this numerical scheme were investigated in [41].

6. Simulations of the SARS-CoV-2 Model (1)

6.1. Baseline Values and Initial Conditions

The life expectancy in Pakistan is 67.27 years, and the population is composed of 231,402,116 individuals [33]. We set the recruitment rate to be

\frac{231, 402, 116}{67.27 \times 365}

per day, while the natural death rate is obtained as

\frac{1}{67.27 \times 365}

per day. In Pakistan, as of 1 January 2022, the cumulative COVID-19 cases recorded stands at 1,296,527, with 28,941 cumulative deaths recorded. The initial conditions are set as

S (0) = 160, 000, 000, V (0) = 70, 000, 000, I_{1} (0) = I_{1}^{v} (0) = H_{1} (0) = 319, 895, R_{1} (0) = 0, I_{2} (0) = I_{2}^{v} (0) = H_{2} (0) = 319, 895, R_{2} (0) = 0

6.2. Model Fitting

Using the data available for the daily confirmed COVID-19 deaths in Pakistan [33] between 1 January 2022 and 10 April 2022, the fractional model is fitted to the data. During this period, both the Delta and Omicron variants were in high circulation within Pakistan. Parameter values estimated from the fitting are shown in Table 1. The fitting shown in Figure 2 reveals that our model fits very well to the data when the fractional order

ζ = 0.96

6.3. Numerical Assessment

In Figure 3a–j, simulations of the various classes comparing the two-step Lagrange polynomial method and the fractional Euler method are presented. We take the fractional order equal to

ζ = 0.96

and the step size

h = 1.0

. It is observed that the behaviours of both methods are close, but not the same. However, by decreasing the step size h to

0.1

, as shown in Figure 4a–j, both methods behave well and almost the same.

In Figure 5a–f, simulations of the various classes of System (1) assessing the impact of vaccination are highlighted using the two-step Lagrange polynomial method, where the fractional order is

ζ = 0.90

and the step size

h = 0.1

. In particular, it is observed that increasing the vaccination rate

γ

from

0.01

0.09

decreases the number of people infected with the Omicron SARS-CoV-2 variant, as can be seen in Figure 5a–c. Additionally, a noticeable decrease is observed in the total number of infected individuals with the Delta SARS-CoV-2 variant, as can be seen in Figure 5d–f. This is due to the high efficacy of the vaccine against both variants. Thus, a safe and effective vaccination program is crucial for the elimination of the different variants of concern.

Numerical experiments on System (1) for various fractional orders are presented in Figure 6a–j. The impact of the fractional orders are clearly observed.

Numerical assessments of the impact of various parameters on the Delta variant’s reproduction number are presented in Figure 7a–e. Specifically, Figure 7a depicts the surface plot of the Delta variant’s reproduction number as a function of the transmission rate (

β_{1}

) and vaccine efficacy (

ψ

). It is observed in this figure that, as the transmission rate is increased, the reproduction number and disease burden also increase. However, it is noted that increased vaccine efficacy decreases both the reproduction number and disease burden.

In Figure 7b, the surface plot of the Delta variant’s reproduction number as a function of transmission rate

β_{1}

and progression rate

α_{1}

is presented. It can be seen from this figure that as the transmission and progression rates increase, the reproduction number also increases.

In Figure 7c, a surface plot of the reproduction number for the Delta variant as a function of the transmission rates

β_{1}

and

θ_{1}

is presented. It can be observed that as the transmission rate and modification parameters accounting for the infectivity of symptomatic individuals increase, the reproduction number and disease burden also increases. A similar conclusion can be drawn for the surface plots presented in Figure 7d,e.

Numerical experiments to assess the impact of various parameters on the Omicron variant’s reproduction number are presented in Figure 8a–e. For instance, in Figure 8a, the surface plot for the Omicron variant’s reproduction number as a function of the transmission rate (

β_{2}

) and vaccine efficacy (

ω

) is presented. It is observed that as the transmission rate increases, the reproduction number and disease burden also increase. However, it is noted that as the vaccine efficacy increases, the reproduction number decreases. This shows that to reduce disease burden, the vaccine efficacy must be as high as possible.

In Figure 8b, a surface plot of the Omicron variant’s reproduction number as a function of the transmission rate

β_{2}

and progression rate

α_{3}

is presented. From the plot, it can be noted that, as the transmission rate and progression rates increase, both the reproduction number and the disease burden increase.

In Figure 8c, a surface plot of the Omicron variant’s reproduction number as a function of the transmission rate

β_{2}

and modification parameter

θ_{2}

is presented. From this figure, it can be observed that, as the transmission rate and the modification parameter accounting for the infectivity of individuals in the symptomatic class are increased, the reproduction number also increases. Similar conclusions can be made for the surface plots presented in Figure 8d,e.

7. Conclusions

In this paper, a new vaccination model for SARS-CoV-2 variants was developed and analyzed. The model assumes variability in transmission rates due to vaccination, where the asymptomatic infected classes were categorized based on their vaccination history. Using the Banach contraction principle and the Arzela–Ascoli theorem, the existence and uniqueness results for the fractional system were obtained. Two numerical schemes, the two-step Lagrange polynomial and the fractional Euler methods, were employed to discuss the approximate solution of the system. The model was fitted to data associated with COVID-19 deaths in Pakistan between 1 January 2022 and 10 April 2022. Assessing the fractional system numerically, the impact of vaccination and fractional order were highlighted, where it was observed that vaccination could remarkably curtail the spikes of emerging variants of SARS-CoV-2 within the population.

Important highlights from the qualitative analysis are:

(i).: The proposed fractional-order system was shown to be positively invariant using Theorem 1.
(ii).: Conditions for the existence of a solution to the fractional-order system were obtained using Theorems 2 and 3, while the uniqueness result was presented using Theorem 4.
(iii).: The designed model was also shown to be generalized Ulam–Hyers stable with the help of Theorem 6.

Highlights of the simulations include:

(iv).: The proposed model fits well to data when the order of the fractional derivative is taken as $ζ = 0.96$ , as can be observed in Figure 2.
(v).: The two numerical approaches were thus compared with different step sizes, where it was concluded that they behave closely for very small step sizes and differently for larger step sizes, as can be observed in Figure 3 and Figure 4.
(vi).: It was observed that increasing the vaccination rate $γ$ from $0.01$ to $0.09$ caused a significant decrease in the number of people infected with the Omicron SARS-CoV-2 variant, as can be seen in Figure 5a–c.
(vii).: Surface plots showing the important parameters driving the dynamics of the disease are presented in Figure 7 and Figure 8.

The model discussed in this study has some limitations. The model did not capture any interaction between the variants. Thus, future work in this direction may consider infection-acquired cross-immunity, cross-enhancement, super-infection, and dual co-infection with both variants. One could also study the in-host dynamics between both variants of SARS-CoV-2. In addition, one could also consider stochastic equivalence of the current model for a possible research problem. Moreover, one could also establish existence, uniqueness, and stability results using some novel fixed-point theorems.

Author Contributions

Conceptualization, M.A. and A.O.; methodology, M.U. and M.A.; software, M.U. and A.O.; formal analysis, M.U., M.A. and A.O.; investigation, A.O.; supervision, M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are openly available via the link: https://ourworldindata.org/coronavirus/country/pakistan, accessed on 1 May 2023.

Conflicts of Interest

The authors declare no conflict of interest.

References

Hu, B.; Guo, H.; Zhou, P.; Shi, Z.L. Characteristics of SARS-CoV-2 and COVID-19. Nat. Rev. Microbiol. 2021, 19, 141–154. [Google Scholar] [CrossRef]
United States Food and Drug Administration. FDA Takes Key Action in Fight Against COVID-19 By Issuing Emergency Use Authorization for First COVID-19 Vaccine. 2020. Available online: https://www.fda.gov/news-events/press-announcements/fda-takes-key-action-fight-against-covid-19-issuing-emergency-use-authorization-first-covid-19 (accessed on 17 June 2021).
Interim Clinical Considerations for Use of COVID-19 Vaccines Currently Authorized in the United States. Available online: https://www.cdc.gov/vaccines/covid-19/clinical-considerations/covid-19-vaccines-us.html (accessed on 14 July 2021).
Tang, P.; Hasan, M.R.; Chemaitelly, H.; Yassine, H.M.; Benslimane, F.M.; Al Khatib, H.A.; AlMukdad, S.; Coyle, P.; Ayoub, H.H.; Al Kanaani, Z.; et al. BNT162b2 and mRNA-1273 COVID-19 vaccine effectiveness against the SARS-CoV-2 Delta variant in Qatar. Nat. Med. 2021; Epub ahead of print. [Google Scholar] [CrossRef]
Andrews, N.; Tessier, E.; Stowe, J.; Gower, C.; Kirsebom, F.; Simmons, R.; Gallagher, E.; Chand, M.; Brown, K.; Ladhani, S.N.; et al. Vaccine effectiveness and duration of protection of Comirnaty, Vaxzevria and Spikevax against mild and severe COVID-19 in the UK. medRxiv 2021. [Google Scholar] [CrossRef]
Nasreen, S.; Chung, H.; He, S.; Brown, K.A.; Gubbay, J.B.; Buchan, S.A.; Fell, D.B.; Austin, P.C.; Schwartz, K.L.; Sundaram, M.E.; et al. Effectiveness of COVID-19 vaccines against symptomatic SARS-CoV-2 infection and severe outcomes with variants of concern in Ontario. Nat. Microbiol. 2022, 7, 379–385. [Google Scholar] [CrossRef] [PubMed]
Puranik, A.; Lenehan, P.J.; Silvert, E.; Niesen, M.J.M.; Corchado-Garcia, J.; O’Horo, J.C.; Virk, A.; Swift, M.D.; Halamka, J.; Badley, A.D.; et al. Comparison of two highly-effective mRNA vaccines for COVID-19 during periods of Alpha and Delta variant prevalence. medRxiv 2021. [Google Scholar] [CrossRef]
Nanduri, S.; Pilishvili, T.; Derado, G.; Soe, M.M.; Dollard, P.; Wu, H.; Li, Q.; Bagchi, S.; Dubendris, H.; Link-Gelles, R.; et al. Effectiveness of Pfizer-BioNTech and Moderna Vaccines in Preventing SARS-CoV-2 Infection Among Nursing Home Residents Before and During Widespread Circulation of the SARS-CoV-2 B.1.617.2 (Delta) Variant— National Healthcare Safety Network, March 1–August 1, 2021. Morb. Mortal. Wkly Rep. 2021, 70, 1163–1166. [Google Scholar]
Fowlkes, A.; Gaglani, M.; Groover, K.; Thiese, M.S.; Tyner, H.; Ellingson, K. Effectiveness of COVID-19 Vaccines in Preventing SARS-CoV-2 Infection Among Frontline Workers Before and During B.1.617.2 (Delta) Variant Predominance—Eight U.S. Locations, December 2020–August 2021. Morb. Mortal. Wkly Rep. 2021, 70, 1167–1169. [Google Scholar] [CrossRef]
Israel Ministry of Health. COVID-19 Vaccine Effectiveness against the Delta Variant. Israel’s Ministry of Health Report. 2021. Available online: Https://www.gov.il/BlobFolder/reports/vaccine-efficacy-safety-follow-up-committee/he/files_publications_corona_two-dose-vaccination-data.pdf (accessed on 11 November 2022).
Ali, A.; Khan, S.U.; Ali, I.; Khan, F.U. On dynamics of stochastic avian influenza model with asymptomatic carrier using spectral method. Math. Methods Appl. Sci. 2022, 45, 8230–8246. [Google Scholar] [CrossRef]
Bonyah, E.; Khan, M.A.; Okosun, K.O.; Gómez-Aguilar, J.F. On the co-infection of dengue fever and Zika virus. Optim. Contr. Appl. Meth. 2019, 40, 394–421. [Google Scholar] [CrossRef]
Rihan, F.A.; Alsakaji, H.J. Dynamics of a stochastic delay differential model for COVID-19 infection with asymptomatic infected and interacting people: Case study in the UAE. Results Phys. 2021, 28, 104658. [Google Scholar] [CrossRef]
Rihan, F.A.; Alsakaji, H.J.; Rajivganthi, C. Stochastic SIRC epidemic model with time-delay for COVID-19. Adv. Differ. Equ. 2020, 2020, 502. [Google Scholar] [CrossRef]
Omame, A.; Abbas, M. Modeling SARS-CoV-2 and HBV co-dynamics with optimal control. Physica A 2023, 615, 128607. [Google Scholar] [CrossRef] [PubMed]
Omame, A.; Abbas, M. The stability analysis of a co-circulation model for COVID-19, dengue, and zika with nonlinear incidence rates and vaccination strategies. Healthc. Anal. 2023, 3, 100151. [Google Scholar] [CrossRef] [PubMed]
Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 2000, 339, 1–77. [Google Scholar] [CrossRef]
Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; North-Holland Mathematics Studies; Elsevier Science: Amsterdam, The Netherlands, 2006; p. 204. [Google Scholar]
Lin, W. Global existence theory and chaos control of fractional differential equations. J. Math. Anal. Appl. 2007, 332, 709–726. [Google Scholar] [CrossRef]
Caputo, M.; Fabrizio, M. A new definition of fractional derivative without singular kernel. Prog. Fract. Differ. Appl. 2015, 1, 73–85. [Google Scholar]
Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and applications to heat transfer model. Therm Sci. 2016, 20, 763–769. [Google Scholar] [CrossRef]
Omame, A.; Abbas, M.; Onyenegecha, C.P. A fractional-order model for COVID-19 and tuberculosis co-infection using Atangana-Baleanu derivative. Chaos Solitons Fractals 2021, 153, 111486. [Google Scholar] [CrossRef]
Omame, A.; Abbas, M.; Nwajeri, U.K. A fractional-order control model for Diabetes COVID-19 co-dynamics with Mittag-Leffler function. Alex. Eng. J. 2022, 61, 7619–7635. [Google Scholar] [CrossRef]
Ogunrinde, R.B.; Nwajeri, U.K.; Fadugba, S.E.; Ogunrinde, R.R.; Oshinubi, K.I. Dynamic model of COVID-19 and citizens reaction using fractional derivative. Alex. Eng. J. 2021, 60, 2001–2012. [Google Scholar] [CrossRef]
Asamoah, J.K.K. Fractal-fractional model and numerical scheme based on Newton polynomial for Q fever disease under Atangana-Baleanu derivative. Results Phys. 2022, 34, 105189. [Google Scholar] [CrossRef]
Asamoah, J.K.K.; Okyere, E.; Yankson, E.; Opoku, A.A.; Adom-Konadu, A.; Acheampong, E.; Arthur, Y.D. Non-fractional and fractional mathematical analysis and simulations for Q fever. Chaos Solitons Fractals 2022, 156, 111821. [Google Scholar] [CrossRef]
Safari, F.; Chen, W. Numerical approximations for space-time fractional Burgers equations via a new semi-analytical method. Comput. Math. Appl. 2021, 96, 55–66. [Google Scholar] [CrossRef]
Jafari, H.; Khalique, C.M.; Nazari, M. Application of the Laplace decomposition method for solving linear and nonlinear fractional diffusion-wave equations. Appl. Math. Lett. 2011, 24, 1799–1805. [Google Scholar] [CrossRef]
Omame, A.; Zaman, F.D. Solution of the modified time fractional coupled Burgers equation using Laplace Adomian decomposition method. Acta Mech. Autom. 2023, 17, 124–132. [Google Scholar] [CrossRef]
Podlubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
Rothana, H.A.; Byrareddy, S.N. The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak. J. Autoimmun. 2020, 109, 102433. [Google Scholar] [CrossRef] [PubMed]
Chen, T.-M.; Rui, J.; Wang, Q.-P.; Zhao, Z.-Y.; Cui, J.-A.; Yin, L. A mathematical model for simulating the phase-based transmissibility of a novel coronavirus. Infect. Dis. Poverty 2020, 9, 24. [Google Scholar] [CrossRef]
Pakistan: Coronavirus Pandemic Country Profile. Available online: https://ourworldindata.org/coronavirus/country/pakistan (accessed on 11 November 2022).
Yong, Z.; Jinrong, W.; Lu, Z. Basic Theory of Fractional Differential Equations; World Scientific: Singapore, 2016. [Google Scholar]
van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef] [PubMed]
Ulam, S.M. A Collection of Mathematical Problems; Interscience Publishers: New York, NY, USA, 1960. [Google Scholar]
Ulam, S.M. Problem in Modern Mathematics; Dover Publications: Mineola, NY, USA, 2004. [Google Scholar]
Toufik, M.; Atangana, A. New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models. Eur. Phys. J. Plus 2017, 132, 444. [Google Scholar] [CrossRef]
Li, C.; Zeng, F. The finite difference methods for fractional ordinary differential equations. Numer. Funct. Anal. Optim. 2013, 34, 149–179. [Google Scholar] [CrossRef]
Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52. [Google Scholar] [CrossRef]
Baleanu, D.; Jajarmi, A.; Hajipour, M. On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag-Leffler kernel. Nonlinear Dyn. 2018, 94, 397–414. [Google Scholar] [CrossRef]

Figure 1. Schematic diagram of Model (1), with

λ_{1} = (\frac{β_{1} (I_{1} + ϕ_{1 v} I_{1}^{v} + θ_{1} H_{1})}{N})

λ_{2} = (\frac{β_{2} (I_{2} + ϕ_{2 v} I_{2}^{v} + θ_{2} H_{2})}{N})

Figure 1. Schematic diagram of Model (1), with

λ_{1} = (\frac{β_{1} (I_{1} + ϕ_{1 v} I_{1}^{v} + θ_{1} H_{1})}{N})

λ_{2} = (\frac{β_{2} (I_{2} + ϕ_{2 v} I_{2}^{v} + θ_{2} H_{2})}{N})

Figure 2. Fitting the model to data.

Figure 3. Simulations of the various classes comparing the two-step Lagrange polynomial method and the fractional Euler method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 1.0

Figure 3. Simulations of the various classes comparing the two-step Lagrange polynomial method and the fractional Euler method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 1.0

Figure 4. Simulations of the various classes comparing the two-step Lagrange polynomial method and the fractional Euler method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 0.1

Figure 4. Simulations of the various classes comparing the two-step Lagrange polynomial method and the fractional Euler method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 0.1

Figure 5. Simulations of the various classes of System (1) assessing the impact of vaccination using the two-step Lagrange polynomial method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 0.1

Figure 5. Simulations of the various classes of System (1) assessing the impact of vaccination using the two-step Lagrange polynomial method. Here, the fractional order is

ζ = 0.96

, and the step size

h = 0.1

Figure 6. Simulations of the various classes of System (1) assessing the impact of fractional order using the two-step Lagrange polynomial method. Here, the step size

h = 0.1

Figure 6. Simulations of the various classes of System (1) assessing the impact of fractional order using the two-step Lagrange polynomial method. Here, the step size

h = 0.1

Figure 7. Simulations of the reproduction number for Delta variant for System (1) assessing the impact of different parameters.

Figure 8. Simulations of the reproduction number for Omicron variant for System (1) assessing the impact of different parameters.

Table 1. Parameters describing the flows in Model (1).

Parameter	Description	Value	Source
$ψ$	Vaccine efficacy against the Delta SARS-CoV-2 variant	0.80	[4]
$ω$	Vaccine efficacy against the Omicron SARS-CoV-2 variant	0.57	Assumed
$η_{1}, η_{2}$	Loss of infection-acquired immunity for Delta
	and Omicron variants, respectively	0.01	Assumed
$ϕ_{1 v}, ϕ_{2 v}$	Modification parameters accounting for the reduced
	transmission rate of individuals in $I_{1}^{v}$ and $I_{2}^{v}$ classes	0.6–1.0	Assumed
$θ_{1}, θ_{2}$	Modification parameters accounting for increased
	transmission rate of symptomatic individuals
	in $H_{1}$ and $H_{1}$ , respectively	1.2	Inferred from [3]
$α_{1}, α_{2}, α_{3}, α_{4}$	Progression rates	$\frac{1}{14}$	[31]
$τ_{1}, τ_{2}, q_{1}, q_{2}, q_{3}, q_{4}$	Recovery rates	$\frac{1}{15}$	[32]
$β_{1}$	Contact rate for Delta variant transmission	0.1940	Fitted
$β_{2}$	Contact rate for Omicron variant transmission	0.1043	Fitted
$δ_{1}$	Delta-variant-induced death rate	$6.0461 \times 10^{- 5}$	Fitted
$δ_{2}$	Omicron-variant-induced death rate	$1.2295 \times 10^{- 7}$	Fitted
$Ψ$	Recruitment rate for individuals	$\frac{231, 402, 116}{67.27 \times 365}$	[33]
$ξ$	Fraction of vaccinated individuals	0.5	Assumed
$γ$	Vaccination rate	0.05	Assumed
$ϑ$	Natural death rate	$\frac{1}{67.27 \times 365}$	[33]
$R_{0 D}$	Delta-variant-associated reproduction number	1.9795	Fitted
$R_{0 M}$	Omicron-variant-associated reproduction number	1.0846	Fitted

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 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/).

Share and Cite

MDPI and ACS Style

Usman, M.; Abbas, M.; Omame, A. Analysis of the Solution of a Model of SARS-CoV-2 Variants and Its Approximation Using Two-Step Lagrange Polynomial and Euler Techniques. Axioms 2023, 12, 480. https://doi.org/10.3390/axioms12050480

AMA Style

Usman M, Abbas M, Omame A. Analysis of the Solution of a Model of SARS-CoV-2 Variants and Its Approximation Using Two-Step Lagrange Polynomial and Euler Techniques. Axioms. 2023; 12(5):480. https://doi.org/10.3390/axioms12050480

Chicago/Turabian Style

Usman, Muhammad, Mujahid Abbas, and Andrew Omame. 2023. "Analysis of the Solution of a Model of SARS-CoV-2 Variants and Its Approximation Using Two-Step Lagrange Polynomial and Euler Techniques" Axioms 12, no. 5: 480. https://doi.org/10.3390/axioms12050480

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Analysis of the Solution of a Model of SARS-CoV-2 Variants and Its Approximation Using Two-Step Lagrange Polynomial and Euler Techniques

Abstract

1. Introduction

2. Model Formulation

Model’s Basic Properties

3. Existence and Uniqueness of The Solution

3.1. Existence

3.2. Uniqueness

3.3. The Basic Reproduction Number of the Model

3.4. Local Asymptotic Stability of the Disease-Free Equilibrium (DFE) of the Model

4. Ulam–Hyers Stability

5. Numerical Scheme

5.1. Two-Step Lagrange Polynomial Method

5.2. Fractional Euler Method

6. Simulations of the SARS-CoV-2 Model (1)

6.1. Baseline Values and Initial Conditions

6.2. Model Fitting

6.3. Numerical Assessment

7. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI