Deterministic and Fractional-Order Co-Infection Model of Omicron and Delta Variants of Asymptomatic SARS-CoV-2 Carriers

Abstract

The Delta and Omicron variants’ system was used in this research study to replicate the complex process of the SARS-CoV-2 outbreak. The generalised fractional system was designed and rigorously analysed in order to gain a comprehensive understanding of the transmission dynamics of both variants. The proposed dynamical system has heredity and memory effects, which greatly improved our ability to perceive the disease propagation dynamics. The non-singular Atangana–Baleanu fractional operator was used to forecast the current pandemic in order to meet this challenge. The Picard recursions approach can be used to ensure that the designed fractional system has at least one solution occupying the growth condition and memory function regardless of the initial conditions. The Hyers–Ulam–Rassias stability criteria were used to carry out the stability analysis of the fractional governing system of equations, and the fixed-point theory ensured the uniqueness of the solution. Additionally, the model exhibited global asymptotically stable behaviour in some conditions. The approximate behaviour of the fatal virus was investigated using an efficient and reliable fractional numerical Adams–Bashforth approach. The outcome demonstrated that there will be a significant decline in the population of those infected with the Omicron and Delta SARS-CoV-2 variants if the vaccination rate is increased (in both the symptomatic and symptomatic stages).

1. Introduction

In the early Twentieth Century, the continuous advancements in mathematical modelling were used in computational biology [1,2] to investigate the spread of diseases. Analysts and investigators acquired valuable knowledge about numerous transmissible diseases by utilising deterministic and stochastic models. In 1927, Kermack developed a useful system for evolving and executing epidemic models, which has been treated as a fundamental model in epidemiology to this day [3]. Many infectious diseases, such as Hepatitis B, HIV/AIDS, Herpes Simplex Virus, Rubella, etc., can transmit horizontally and vertically. Such types of diseases are horizontally spread in animals and humans through proximity among the hosts or through disease carriers, e.g., mosquitoes, flies, etc. In 2003, Severe Acute Respiratory Syndrome (SARS) was discovered in China [4,5], and the outbreak of MERS appeared in South Korea in 2015 [6,7].

The recent outbreak of the coronavirus pandemic (COVID-19) firstly appeared in Wuhan, China, in 2019 and has precipitously spread throughout the entire world due to its high rate of transmission. COVID-19 was announced as a global pandemic by the World Health Organization (WHO) on 8 March 2020, and the ICTV named itSevere Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [8]. There have been 0.1 billion cases reported, and deaths surpassed 2 million on 1 January 2021. Numerous studies have declared that it could have originated from pangolins or bats [9], and the origination of the virus in the human populationcould have been associated with a seafood market, but the intermediary host has not been proven [10,11]. The infected individuals initially had symptoms similar to SARS-COVID and MERS-COVID infections, such as fatigue, cough, fever, conjunctivitis, sore throat, and, in many cases, bilateral lung penetration [12]. Moreover, many patients may lose their sense of taste and/or smell, may have diarrhoea, and/or may have signs of a breathing disorder [13,14].

Fractional calculus is one of the most-promising avenues in modern science. There has been a huge amount of research performed in this area, which has contributed to remarkable applications. The reason for this is that the fractional operator is a global operator that predicts the behaviour of physical phenomena at infinity tails, whereas the local operator ignores the influence of larger neighbours. The fractional derivative can also be used to solve real-world problems where the integer-order derivative is singular. Numerous significant fractional operators have been developed such as the Caputo fractional derivative [15], the conformal fractional differential operator [16], the $\beta -$fractional operator [17], the truncated M-fractional operator [18], the Riemann–Liouville fractional derivative [19], the Caputo–Fabrizio differential operator [20], and the fractional Atangana–Baleanu derivative [21]. In the last few decades, remarkable developments have emerged within the field of fractional partial differential equations because of their applicability to various areas of technology and science [22,23,24]. Imran et al. [25,26] discussed an MHD viscous fluid and bi-convection using fractional differential operators and provided applications to fluid dynamics. Maryam et al. [27] established the application of fractional calculus to the field of nano-fluids. Naik et al. [28,29,30] studied the fractional-order model of the HIV transmission epidemic with optimal control. Mishra et al. [31] proposed a COVID-19 model with asymptotic and quarantine classes. Owolabi et al. [32,33] analysed breast cancer and brain tumour disease using a type of fractional calculus. A dynamical analysis of river blindness disease was presented by Atangana [34]. Salman examined the HBV infection disease model under the influence of the fractional-order derivative [35]. Area investigated the fractional Ebola epidemic model [36]. Faridi developed a mathematical fractional cardiovascular disease model applying the Atangana–Baleanu fractional operator with the memory function and growth condition [37]. The Rubella disease model was examined by the Caputo–Fabrizio fractional operator [38]. Addai et al. [39] proposed a new Caputo–Fabrizio fractional SARS-CoV-2 epidemiological model with Alzheimer’s disease, applied Lagrange interpolation for the approximation, as well as discussed the stability analysis. Omame et al. [40,41,42] developed many co-infectious disease models, analysed them through different techniques, and presented a numerical simulation.

There are many epidemiological models such as the SIR [43,44,45], SIS [46], SEIQR [47,48], and SIRC models [49,50]. The SIR model can be used to develop the elementary governing model for COVID-19 disease, which presents the population as three different classes, which are susceptible, infected, and recovered individuals [51]. The divisionof these categories has been examined in a type of classical-order derivative by many analysts to obtain better simulation results, such as the mathematical model of COVID-19’s spread in India, as presented by Biswas [52]. A dynamical investigation of the COVID-19 model in Hubei Province, China, was presented using an SEIQR model [53]. The COVID-19 model was studied for forecasting in Spain along with the study of the impacts of various parameters [54]. A new governing system was proposed by researchers to check the new cases and deaths from the virus’s outbreak in India [55]. The investigation of the virus’s impact in Brazil using a type of multiple delay model and control schemes was provided in [56].

Different variants of Severe Acute Respiratory Syndrome CoronaVirus 2 (SARS-CoV-2) have emerged in recent times, including the highly transmissible and deadlier SARS-CoV-2 Delta variant [57]. This variant has been declared as a Variant Of Concern (VOC) by the World Health Organization (WHO). In almost every country on the planet, this variant has been detected [57]. The Delta variant has also been linked with increased SARS-CoV-2 cases in places where it is the dominant variant in the population [57]. Very recently (precisely, 24 November 2021), South Africa reported the discovery of another SARS-CoV-2 variant, Omicron B.1.1.529, to the WHO. It has already been declared as a VOC as well, even though studies are still ongoing to determine its transmissibility rate and severity [58]. The WHO has authorised a number of vaccines for emergency use, notwithstanding the fact that there are multiple control measures against SARS-CoV-2. Some of these are: the BNT162b2 Pfizer-BioNtech vaccine, the mRNA-1273 Moderna vaccine, the Johnson and Johnson vaccine, and several others [59,60]. These vaccines are incredibly effective against several SARS-CoV-2 variants, including the Delta variant, according to studies [60]. Recently, much epidemiological research has been performed on the efficacy of current SARS-CoV-2 vaccines against Variants Of Concern (VOCs) [61,62,63,64,65]. Similar results were obtained in studies conducted in Israel [61], the United Kingdom [62], Canada [63], and the USA [64,65], where the Pfizer and Moderna vaccines showed 39–75% efficacy against the Delta variant. In a recent study conducted by Tang et al. [66], to measure the effectiveness of the Pfizer and Moderna vaccines against the SARS-CoV-2 Delta variant in Qatar, they reported that those fully vaccinated with any of these two vaccines have 51.9% protection against the Delta variant.

Some of the novelties of this work are as follows: we considered and analysed a model of SARS-CoV-2 with two variants and vaccination-induced cross-immunity, proved the stability of the equilibrium point, as well as analysed how the fractional-order model impacts the dynamics of SARS-CoV-2 variants.

1.1. Preliminaries

Definition 1

([21]). Let $\daleth \in H^1(a,b),b>a$, then the Atangana–Baleanu fractional integral with $0<\varsigma \le 1$ for the function $\daleth \left(t\right)$ is

$$\begin{array}{c}\hfill {}_b^{ABC}{I}^{\varsigma }\left({\daleth }_t\left(t\right)\right)=\frac{(1-\varsigma )\daleth \left(t\right)}{AB\left(\varsigma \right)}+\frac{\varsigma }{\mathsf{\Gamma }\left(\varsigma \right)AB\left(\varsigma \right)}{\int }_b^t{\daleth }^{\prime }\left(\zeta \right){(t-\zeta )}^{\varsigma -1}d\zeta \end{array}$$

(1)

Definition 2

([21]). Let $\daleth \in H^1(a,b),b>a$, then the Atangana–Baleanu time fractional derivative with fractional order $0<\varsigma \le 1$ for the function $\daleth \left(t\right)$ is

$$\begin{array}{c}\hfill {}_b^{ABC}{D}^{\varsigma }{\daleth }_t\left(\left(t\right)\right)=\frac{AB\left(\varsigma \right)}{1-\varsigma }{\int }_b^t{\daleth }^{\prime }\left(\zeta \right)E_{\varsigma }\left[\frac{\varsigma }{1-\varsigma }{(t-\zeta )}^{\varsigma }\right]d\zeta ,\end{array}$$

(2)

where the normalisation function $AB\left(0\right)=AB\left(1\right)=1$.

1.2. Linear Growth Condition

Suppose ∃$ϵ>o$, as in [67]:

$$\begin{array}{c}\hfill {\left|(t,y)\right|}^2\le ϵ\left(\right|y|{|}^2+1),\forall (t,y)\in \mathfrak{R}\times [t_0,\beth ].\end{array}$$

(3)

1.3. Data Operator

The integral operator conserves the initial condition, while the derivative loses the initial condition. This implies that the derivative has no memory, while the integral does, i.e., it can remember the initial state. The integral operator is well-posed, while the derivative is ill-posed, as is common knowledge. This has implications for inverse problems. Assume an unknown function g, with h representing the data of the inverse problem. Obviously, an operator equation gives the inverse problem [67]:

$$\begin{array}{c}\hfill \mathsf{\Delta }g=h,\end{array}$$

(4)

where $\mathsf{\Delta }\in \mathfrak{U}$ for the data.

1.4. Memory Function

We discuss a memory function, which is associated with the Caputo–Fabrizio (CF) operator. Angel et al. [68] observed that the (CF) operator shows a cross-over effect in terms of the probability distribution and the mean-squared displacement. Moreover, the (CF) operator performs the usual diffusion and confined diffusion before and after the time $t_1=a$, respectively. Therefore, the memory function can be defined as [67]

$${\gamma }^{\prime }\left(t\right)=\left\{\begin{array}{cc}0,\hfill & \hfill \mathrm{for}a\ge t,\\ \frac{(1-\varsigma )f\left(0\right)}{aM\left(\varsigma \right)\mathsf{\Gamma }(2-\varsigma )\mathsf{\Gamma }\left(\varsigma \right)}exp(-\frac{\varsigma }{1-\varsigma }t){(a-t)}^{\varsigma -1}{t}^{1-\varsigma },\hfill & \hfill \mathrm{for}a<t.\end{array}\right\}$$

This is a piecewise function, which is zero at the beginning of the process and zero after the cross-over since the history has been captured.

2. Model Formulation

At any given time, t is the total population, and $\mathcal{N}$ consists of nine mutually exclusive compartments: $\mathcal{S}\left(t\right)$: susceptible individuals; $\mathcal{V}\left(t\right)$: vaccinated individuals; ${\mathcal{A}}_1\left(t\right),{\mathcal{A}}_2\left(t\right),$${\mathcal{A}}_{12}\left(t\right)$: infectious individuals with the Delta or Omicron variant or co-infection with both variants (asymptomatic stage), respectively; ${\mathcal{I}}_1\left(t\right),{\mathcal{I}}_2\left(t\right),{\mathcal{I}}_{12}\left(t\right)$: infectious individuals with the Delta or Omicron variant or co-infection with both variants (in the symptomatic stage), respectively; $\mathcal{R}\left(t\right)$: individuals who have recovered from either the Delta or Omicron variant or from co-infection. Unvaccinated susceptible individuals are recruited into the population at the rate $(1-\chi )\mathsf{\Pi }$, with $\chi$ as the vaccination rate. Individuals in the unvaccinated susceptible class acquire the Omicron SARS-CoV-2 variant at the rate ${\beta }_1({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)$. They also acquire the Delta SARS-CoV-2 variant at the rate ${\beta }_2({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)$. Infection with incident co-infection is at the rate ${\beta }_{12}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})$. Individuals in this state are vaccinated at the rate $\varpi$. Natural death occurs in this state, as well as in all other epidemiological states at the rate $\mu$. Vaccinated susceptible individuals are recruited into the population at the rate $\chi \mathsf{\Pi }$. Individuals in the vaccinated class are infected by the Delta SARS-CoV-2 variant at the rate $(1-{\psi }_1){\beta }_1({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)$, where ${\psi }_1$ is the efficacy of the vaccine against the Delta SARS-CoV-2 variant. This epidemiological state also is infected with the Omicron SARS-CoV-2 variant at the rate $(1-{\psi }_2){\beta }_2({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)$, where ${\psi }_2$ is the efficacy of the vaccine against the Omicron SARS-CoV-2 variant. Vaccinated individuals can also be co-infected at the reduced rate $(1-{\psi }_{12}){\beta }_{22}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})$, where, again, ${\psi }_{12}$ is the vaccine’s efficacy against co-infection with both strains. Individuals infected with either strain can be additionally infected with the other strain at the rates ${\xi }_1{\beta }_2({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)$ and ${\xi }_2{\beta }_1({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)$, respectively. The parameters that describe the flows from one epidemiological state to the other are defined in

Table 1. Parameters describing the flows in the model (7) (note that the values for some parameters fit by the model using the data from India and Pakistan are given in Table 2).

Parameter	Description	Value	Source
$\mathsf{\Pi }$	Recruitment rate for individuals	See Table 2
${\beta }_1,{\beta }_2$	Contact rates for Omicron and Delta variants’ transmission	See Table 2	Fit
${\beta }_{12}$	Contact rate for co-infection transmission	See Table 2	Fit
$\chi$	Fraction of vaccinated individuals	See Table 2
$\varpi$	Vaccination rate	See Table 2
${\psi }_1$	Vaccine efficacy against the Omicron SARS-CoV-2 variant	0.80	[66]
${\psi }_2$	Vaccine efficacy against the Delta SARS-CoV-2 variant	0.57	[66]
$\varrho$	Cross-immunity parameter	1.0	Assumed
$\mu$	Natural death rate	See Table 2
${\mathsf{\Delta }}_1,{\mathsf{\Delta }}_2$	Omicron and Delta variants induced death rates	See Table 2	Fit
${\mathsf{\Delta }}_{12}$	Co-infection induced death rates	See Table 2	Fit
${\xi }_{12}$	Modification parameter accounting for the infectivity
	of individuals in the $I_{12}$ class	1.0	Assumed
${\xi }_1,{\xi }_2$	Modification parameters for the infectiousness of
	symptomatic individuals in ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$, respectively	1.5	[60]
${\eta }_1,{\eta }_2$	Omicron, Delta SARS-CoV-2 variant’s progression rates	$\frac{1}{14}$	[69]
${\eta }_{12}$	Progression rate for co-infection of SARS-CoV-2 variants	$\frac{1}{14}$	[70]
${\zeta }_1,{\zeta }_2,{\zeta }_{12},q_1,q_2,q_{12}$	Recovery rates	$\frac{1}{15}$	[70]

. Other assumptions in the model are as follows:

• Incident co-infection with both strains is assumed, and also, that the vaccine has some efficacy against incident co-infection.
• The model considers the SARS-CoV-2 Omicron variant (denoted by “1”) and the SARS-CoV-2 Delta variant (denoted by “2”).
• The transmissibility of the Omicron variant is assumed higher than that of the Delta variant [57].
• Vaccinated susceptible individuals (assumed also to have completed two doses of any of the available vaccines) have a reduced rate of infection by both variants.
• It is further assumed that immigrants in the population have completed their vaccination dosage.

The classical governing system of differential equations to simulate the infectious Delta and Omicron variants is

$$\begin{array}{cc}\hfill \frac{d\mathcal{S}}{dt}& =(1-\chi )\mathsf{\Pi }-(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu )\mathcal{S},\hfill \\ \hfill \frac{d\mathcal{V}}{dt}& =\chi \mathsf{\Pi }+\varpi \mathcal{S}-[(1-{\psi }_1)\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)+(1-{\psi }_2)\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\hfill \\ \hfill & +(1-{\psi }_{12})\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\mu ]\mathcal{V},\hfill \\ \hfill \frac{d{\mathcal{A}}_1}{dt}& =\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)\mathcal{S}+(1-{\psi }_1)\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)\mathcal{V}-\left({\eta }_1+{\zeta }_1+\mu \right){\mathcal{A}}_1-{\xi }_1\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2),\hfill \\ \hfill \frac{d{\mathcal{A}}_2}{dt}& =\left(\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\right)\mathcal{S}+(1-{\psi }_2)\left(\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\right)\mathcal{V}-\left({\eta }_2+{\zeta }_2+\mu \right){\mathcal{A}}_2-{\xi }_2\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1),\hfill \\ \hfill \frac{d{\mathcal{A}}_{12}}{dt}& =\left(\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})\right)\mathcal{S}+(1-{\psi }_{12})\left(\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})\right)\mathcal{V}-\left({\eta }_{12}+{\zeta }_{12}+\mu \right){\mathcal{A}}_{12}\hfill \\ \hfill & +{\xi }_1\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+{\xi }_2\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1),\hfill \\ \hfill \frac{d{\mathcal{I}}_1}{dt}& ={\eta }_1{\mathcal{A}}_1-\left(q_1+{\mathsf{\Delta }}_1+\mu \right){\mathcal{I}}_1,\hfill \\ \hfill \frac{d{\mathcal{I}}_2}{dt}& ={\eta }_2{\mathcal{A}}_2-\left(q_2+{\mathsf{\Delta }}_2+\mu \right){\mathcal{I}}_2,\hfill \\ \hfill \frac{d{\mathcal{I}}_{12}}{dt}& ={\eta }_{12}{\mathcal{A}}_{12}-\left(q_{12}+{\mathsf{\Delta }}_{12}+\mu \right){\mathcal{I}}_{12},\hfill \\ \hfill \frac{d\mathcal{R}}{dt}& ={\zeta }_1{\mathcal{A}}_1+{\zeta }_2{\mathcal{A}}_2+{\zeta }_{12}{\mathcal{A}}_{12}+q_1{\mathcal{I}}_1+q_2{\mathcal{I}}_2+q_{12}{\mathcal{I}}_{12}-\mu \mathcal{R},\hfill \end{array}$$

(5)

with the initial conditions:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathcal{S}\left(0\right)={\mathcal{S}}_0\ge 0,\hfill & \hfill \\ \mathcal{V}\left(0\right)={\mathcal{V}}_0\ge 0,\hfill & \hfill \\ {\mathcal{A}}_1\left(0\right)={\mathcal{A}}_{1,0}\ge 0\hfill \\ {\mathcal{A}}_2\left(0\right)={\mathcal{A}}_{2,0}\ge 0,\hfill & \hfill \\ {\mathcal{A}}_{12}\left(0\right)={\mathcal{A}}_{12,0}\ge 0,\hfill & \hfill \\ {\mathcal{I}}_1\left(0\right)={\mathcal{I}}_{1,0}\ge 0,\hfill & \hfill \\ {\mathcal{I}}_2\left(0\right)={\mathcal{I}}_{2,0}\ge 0,\hfill & \hfill \\ {\mathcal{I}}_{12}\left(0\right)={\mathcal{I}}_{12,0}\ge 0,\hfill & \hfill \\ \mathcal{R}\left(0\right)={\mathcal{R}}_0\ge 0.\hfill & \hfill \end{array}\right.\end{array}$$

(6)

The generalised time fractional governing model with the non-integer order by applying the Atangana–Baleanu operator (2) to the model (5) is

$$\begin{array}{cc}\hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{S}=(1-\chi )\mathsf{\Pi }-\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu \right)\mathcal{S}\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{V}=\chi \mathsf{\Pi }+\varpi \mathcal{S}-[(1-{\psi }_1)\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)+(1-{\psi }_2)\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+\hfill \\ \hfill & (1-{\psi }_{12})\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\mu ]\mathcal{V},\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_1=\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)\mathcal{S}+(1-{\psi }_1)\left(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)\right)\mathcal{V}-\left({\varsigma }_1+{\zeta }_1+\mu \right){\mathcal{A}}_1-{\xi }_1{\beta }_2({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2),\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_2=\left(\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\right)\mathcal{S}+(1-{\psi }_2)\left(\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\right)\mathcal{V}-\left({\varsigma }_2+{\zeta }_2+\mu \right){\mathcal{A}}_2-{\xi }_2\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1),\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_{12}=\left(\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})\right)\mathcal{S}+(1-{\psi }_{12})\left(\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})\right)\mathcal{V}\hfill \\ \hfill & -\left({\eta }_{12}+{\zeta }_{12}+\mu \right){\mathcal{A}}_{12}+{\xi }_1\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+{\xi }_2\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1),\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_1={\eta }_1{\mathcal{A}}_1-\left(q_1+{\mathsf{\Delta }}_1+\mu \right){\mathcal{I}}_1,\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_2={\eta }_2{\mathcal{A}}_2-\left(q_2+{\mathsf{\Delta }}_2+\mu \right){\mathcal{I}}_2,\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_{12}={\eta }_{12}{\mathcal{A}}_{12}-\left(q_{12}+{\mathsf{\Delta }}_{12}+\mu \right){\mathcal{I}}_{12},\hfill \\ \hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{R}={\zeta }_1{\mathcal{A}}_1+{\zeta }_2{\mathcal{A}}_2+{\zeta }_{12}{\mathcal{A}}_{12}+q_1{\mathcal{I}}_1+q_2{\mathcal{I}}_2+q_{12}{\mathcal{I}}_{12}-\mu \mathcal{R}.\hfill \end{array}$$

(7)

Let us consider nine kernel values for the sake of simplicity and clarity, then Model (7) can be written as

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[\mathcal{S}\left(t\right)\right]={\mathfrak{H}}_1(t,\mathcal{S}),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[\mathcal{V}\left(t\right)\right]={\mathfrak{H}}_2(t,\mathcal{V}),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{A}}_1\left(t\right)\right]={\mathfrak{H}}_3(t,{\mathcal{A}}_1),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{A}}_2\left(t\right)\right]={\mathfrak{H}}_4(t,{\mathcal{A}}_2),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{A}}_{12}\left(t\right)\right]={\mathfrak{H}}_5(t,{\mathcal{A}}_{12}),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{I}}_1\left(t\right)\right]={\mathfrak{H}}_6(t,{\mathcal{I}}_1),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{I}}_2\left(t\right)\right]={\mathfrak{H}}_7(t,{\mathcal{I}}_2),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[{\mathcal{I}}_{12}\left(t\right)\right]={\mathfrak{H}}_8(t,{\mathcal{I}}_{12}),\hfill & \hfill \\ {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\left[\mathcal{R}\left(t\right)\right]={\mathfrak{H}}_9(t,\mathcal{R}).\hfill & \hfill \end{array}\right.\end{array}$$

(8)

By applying the fractional integral to Model (8),

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathcal{S}\left(t\right)-\mathcal{S}\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_1(t,\mathcal{S})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_1(\top ,\mathcal{S}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ \mathcal{V}\left(t\right)-\mathcal{V}\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_2(t,\mathcal{V})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_2(\top ,\mathcal{V}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_1\left(t\right)-{\mathcal{A}}_1\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_3(t,{\mathcal{A}}_1)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_3(\top ,{\mathcal{A}}_1){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_2\left(t\right)-{\mathcal{A}}_2\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_4(t,{\mathcal{A}}_2)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_4(\top ,{\mathcal{A}}_2){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_{12}\left(t\right)-{\mathcal{A}}_{12}\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_5(t,{\mathcal{A}}_{12})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_5(\top ,{\mathcal{A}}_{12}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_1\left(t\right)-{\mathcal{I}}_1\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_6(t,{\mathcal{I}}_1)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_6(\top ,{\mathcal{I}}_1){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_2\left(t\right)-{\mathcal{I}}_2\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_7(t,{\mathcal{I}}_2)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_7(\top ,{\mathcal{I}}_2){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_{12}\left(t\right)-{\mathcal{I}}_{12}\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_8(t,{\mathcal{I}}_{12})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_8(\top ,{\mathcal{I}}_{12}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ \mathcal{R}\left(t\right)-\mathcal{R}\left(0\right)=\aleph \left(\varsigma \right){\mathfrak{H}}_9(t,\mathcal{R})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_9(\top ,\mathcal{R}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \end{array}\right.\end{array}$$

(9)

where $\aleph \left(\varsigma \right)=\frac{1-\varsigma }{AB\left(\varsigma \right)}$ and $\wp \left(\varsigma \right)=\frac{\varsigma }{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}$.

3. Mathematical Analysis

The boundedness, existence, uniqueness, and stability analyses are discussed in this section.

Assume a set U of continuous real-valued variables specified on a finite interval ℶ and U fulfilling

$$\begin{array}{cc}\hfill & \parallel (\mathcal{S},\mathcal{V},{\mathcal{A}}_1,{\mathcal{A}}_2,{\mathcal{A}}_{12},{\mathcal{I}}_1,{\mathcal{I}}_2,{\mathcal{I}}_{12},\mathcal{R})\parallel = \parallel \mathcal{S}\parallel + \parallel \mathcal{V}\parallel + \parallel {\mathcal{A}}_1\parallel + \parallel {\mathcal{A}}_2\parallel + \parallel {\mathcal{A}}_{12}\parallel + \hfill \\ & \parallel {\mathcal{I}}_1\parallel + \parallel {\mathcal{I}}_2\parallel + \parallel {\mathcal{I}}_{12}\parallel + \parallel \mathcal{R}\parallel ,\hfill \end{array}$$

(10)

in the Banach space.

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\parallel \mathcal{S}\parallel =sup\left\{\right|{\mathcal{M}}_1\left(t\right)|:t\in \beth \}\le {\omega }_1,\hfill & \hfill \\ \parallel \mathcal{V}\parallel =sup\left\{\right|{\mathcal{M}}_2\left(t\right)|:t\in \beth \}\le {\omega }_2,\hfill & \hfill \\ \parallel {\mathcal{A}}_1\parallel =sup\left\{\right|{\mathcal{M}}_3\left(t\right)|:t\in \beth \}\le {\omega }_3,\hfill & \hfill \\ \parallel {\mathcal{A}}_2\parallel =sup\left\{\right|{\mathcal{M}}_4\left(t\right)|:t\in \beth \}\le {\omega }_4,\hfill & \hfill \\ \parallel {\mathcal{A}}_{12}\parallel =sup\left\{\right|{\mathcal{M}}_5\left(t\right)|:t\in \beth \}\le {\omega }_5,\hfill & \hfill \\ \parallel {\mathcal{I}}_1\parallel =sup\left\{\right|{\mathcal{M}}_6\left(t\right)|:t\in \beth \}\le {\omega }_6,\hfill & \hfill \\ \parallel {\mathcal{I}}_2\parallel =sup\left\{\right|{\mathcal{M}}_7\left(t\right)|:t\in \beth \}\le {\omega }_7,\hfill & \hfill \\ \parallel {\mathcal{I}}_{12}\parallel =sup\left\{\right|{\mathcal{M}}_8\left(t\right)|:t\in \beth \}\le {\omega }_8,\hfill & \hfill \\ \parallel \mathcal{R}\parallel =sup\left\{\right|{\mathcal{M}}_9\left(t\right)|:t\in \beth \}\le {\omega }_9.\hfill & \hfill \end{array}\right.\end{array}$$

Moreover, in the finite domain, a set $U\left(t\right)=\left\{\upsilon \right(\beth )\times \upsilon (\beth \left)\right\}$, $\upsilon (\beth )$ of continuous real-valued functions has the supremum norm property in the Banach space.

3.1. Boundedness of the Governing Model (5)

Theorem 1.

All the solutions are bounded with positive initial conditions and $\mathcal{N}\left(t\right)\le \frac{\mathsf{\Pi }}{\mu },\forall t>0$.

Proof. 

The growth rate for the population can be expressed by (5) as

$$\begin{array}{cc}\hfill & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{N}={}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{S}+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{V}+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_1+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_2+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{A}}_{12}+\hfill \\ & {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_1+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_2+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{I}}_{12}+{}_0^{\mathcal{ABC}}{D}_t^{\varsigma }{\mathcal{R}}_2.\hfill \end{array}$$

(11)

Equation (11) yields

$$\begin{array}{c}\hfill {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{N}\left(t\right)=\mathsf{\Pi }-\mu \mathcal{N}\left(t\right).\end{array}$$

(12)

Applying the Laplace transformation to Equation (12), we obtain

$$\begin{array}{c}\hfill \frac{s^{\varsigma }\overline{\mathcal{N}}\left(s\right)-s^{\varsigma -1}\mathcal{N}\left(0\right)}{(1-\varsigma )s^{\varsigma }+\varsigma }=\frac{\mathsf{\Pi }}{s}-\mu \overline{\mathcal{N}}\left(s\right).\end{array}$$

(13)

$$\begin{array}{c}\hfill \overline{\mathcal{N}}\left(s\right)=\frac{\mathsf{\Pi }\left((1-\varsigma )s^{\varsigma }+\varsigma \right)}{s\left((1+\mu (1-\varsigma ))s^{\varsigma }+\mu \varsigma \right)}+\frac{s^{\varsigma -1}\mathcal{N}\left(0\right)}{s\left((1+\mu (1-\varsigma ))s^{\varsigma }+\mu \varsigma \right)}.\end{array}$$

(14)

After performing a simple calculation and the inverse Laplace transform to Equation (21), one can see

$$\begin{array}{c}\hfill \mathcal{N}\left(t\right)=\frac{\mathsf{\Pi }\varsigma }{\mu }+\left[\frac{\mathsf{\Pi }(1-\varsigma )-\mathsf{\Pi }\varsigma (1+\mu (1-\varsigma \left)\right)+1}{\mu (1+\mu (1-\varsigma \left)\right)}\right]E_{\varsigma }\left(-\frac{\mu \varsigma }{(1+\mu (1-\varsigma \left)\right)}{t}^{\varsigma }\right).\end{array}$$

(15)

Thus, for the limiting case $t\to 0$,

$$\begin{array}{c}\hfill \mathcal{N}\left(t\right)\le \frac{\mathsf{\Pi }}{\mu }.\end{array}$$

(16)

One can say that

$$\begin{array}{c}\hfill 0\le \mathcal{N}\left(t\right)\le \frac{\mathsf{\Pi }}{\mu }.\end{array}$$

(17)

Finally, with the evidence from (17), it is proven that each populations’ individualsare bounded. □

3.2. Positivity of Solution

Theorem 2.

The closed set:

$$\begin{array}{cc}\hfill \mathcal{D}& =\left\{(\mathcal{S},\mathcal{V},{\mathcal{A}}_1,{\mathcal{A}}_2,{\mathcal{A}}_{12},{\mathcal{I}}_1,{\mathcal{I}}_2,{\mathcal{I}}_{12},\mathcal{R})\in {\mathfrak{R}}_{+}^9:\mathcal{S}+\mathcal{V}+{\mathcal{A}}_1+{\mathcal{A}}_2+{\mathcal{A}}_{12}+{\mathcal{I}}_1+{\mathcal{I}}_2+{\mathcal{I}}_{12}+\mathcal{R}\le \frac{\mathsf{\Pi }}{\mu }\right\}.\hfill \end{array}$$

is positively invariant in the context of Model (7).

Proof. 

Adding all the equations of System (7) gives

$${}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{N}=\mathsf{\Pi }-\mu \mathcal{N}\left(t\right)-[{\mathsf{\Delta }}_1{\mathcal{I}}_1+{\mathsf{\Delta }}_2{\mathcal{I}}_2+{\mathsf{\Delta }}_{12}{\mathcal{I}}_{12}],$$

(18)

and from (18), we have that

$$\mathsf{\Pi }-(\mu +3\mathsf{\Delta })\mathcal{N}\le {}_0^{\mathcal{ABC}}{D}_t^{\varsigma }\mathcal{N}<\mathsf{\Pi }-\mu \mathcal{N},$$

where $\mathsf{\Delta }=min\{{\mathsf{\Delta }}_1,{\mathsf{\Delta }}_2,{\mathsf{\Delta }}_{12}\}$.

We obtain $\mathcal{N}\left(t\right)\le \frac{\mathsf{\Pi }}{\mu }$ as $t\to \infty$ by applying the Laplace transform of the AB derivative to the above inequality. Hence, System (7) has a solution in ℶ. Thus, the given system is positively invariant. □

3.3. Lipschitz Condition

Now, we prove the Lipschitz condition for the model.

Theorem 3.

The seven above-mentioned kernels (8):

$$\{{\mathfrak{H}}_1(t,\mathcal{S}),{\mathfrak{H}}_2(t,\mathcal{V}),{\mathfrak{H}}_3(t,{\mathcal{A}}_1),{\mathfrak{H}}_4(t,{\mathcal{A}}_2),{\mathfrak{H}}_5(t,{\mathcal{A}}_{12}),{\mathfrak{H}}_6(t,{\mathcal{I}}_1),{\mathfrak{H}}_7(t,{\mathcal{I}}_2),{\mathfrak{H}}_8(t,{\mathcal{I}}_{12}),{\mathfrak{H}}_9(t,\mathcal{R})\},$$

satisfy the Lipschitz condition.

Proof. 

Firstly, we verify Lipschitz condition for kernel ${\mathfrak{H}}_1$. Assume $\mathcal{S}$ and ${\mathcal{S}}_1$ are two functions, and their corresponding norm is

$$\begin{array}{c}\hfill \begin{array}{cc}& \parallel {\mathfrak{H}}_1(t,\mathcal{S})-{\mathfrak{H}}_1(t,{\mathcal{S}}_1)\parallel = \parallel ((1-\chi )\mathsf{\Pi }-(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\hfill \\ & +\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu \left)\mathcal{S}\right)-\hfill \\ & \left((1-\chi )\mathsf{\Pi }-(\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu ){\mathcal{S}}_1\right)\parallel .\hfill \end{array}\end{array}$$

(19)

Upon simplification and applying the norm property, we obtain

$$\begin{array}{c}\hfill \parallel {\mathfrak{H}}_1(t,S)-{\mathfrak{H}}_1(t,{\mathcal{S}}_1)\parallel \le (\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)\\ \hfill +\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu )\parallel \mathcal{S}\left(t\right)-\mathcal{S}\left(t_1\right)\parallel ,\end{array}$$

$$\begin{array}{c}\hfill \parallel {\mathfrak{H}}_1(t,S)-{\mathfrak{H}}_1(t,{\mathcal{S}}_1)\parallel \le {\phi }_1\parallel S\left(t\right)-S\left(t_1\right)\parallel ,\end{array}$$

(20)

taking ${\phi }_1=\frac{{\beta }_1}{\mathcal{N}}({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1)+\frac{{\beta }_2}{\mathcal{N}}({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2)+\frac{{\beta }_{12}}{\mathcal{N}}({\mathcal{A}}_{12}+{\theta }_{12}{\mathcal{I}}_{12})+\varpi +\mu$, where $({\mathcal{A}}_1+{\theta }_1{\mathcal{I}}_1),({\mathcal{A}}_2+{\theta }_2{\mathcal{I}}_2),$ and $({\mathcal{A}}_{12}+{\mathcal{I}}_{12})$ are bounded.

In a similar way, we can verify the Lipschitz condition for the kernel:

$${\mathfrak{H}}_2\left(\mathcal{V}\right),{\mathfrak{H}}_3\left({\mathcal{A}}_1\right),{\mathfrak{H}}_4\left({\mathcal{A}}_2\right),{\mathfrak{H}}_5\left({\mathcal{A}}_{12}\right),{\mathfrak{H}}_6\left({\mathcal{I}}_1\right),{\mathfrak{H}}_7\left({\mathcal{I}}_2\right),{\mathfrak{H}}_8\left({\mathcal{I}}_{12}\right),\mathrm{and}{\mathfrak{H}}_9\left(R\right),$$

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\parallel {\mathfrak{H}}_2(t,\mathcal{V})-{\mathfrak{H}}_2(t,{\mathcal{V}}_1)\parallel \le {\phi }_2\parallel \mathcal{V}\left(t\right)-\mathcal{V}\left(t_1\right)\parallel ,\hfill & \hfill \\ \parallel {\mathfrak{H}}_3(t,{\mathcal{A}}_1)-{\mathfrak{H}}_3(t,{\mathcal{A}}_{1,1})\parallel \le {\phi }_3\parallel {\mathcal{A}}_1\left(t\right)-{\mathcal{A}}_1\left(t_1\right)\parallel ,\hfill & \hfill \\ \parallel {\mathfrak{H}}_4(t,{\mathcal{A}}_2)-{\mathfrak{H}}_4(t,{\mathcal{A}}_{2,1})\parallel \le {\phi }_4\parallel {\mathcal{A}}_2\left(t\right)-{\mathcal{A}}_2\left(t_1\right)\parallel ,\hfill & \hfill \\ \parallel {\mathfrak{H}}_5(t,{\mathcal{A}}_{12})-{\mathfrak{H}}_5(t,{\mathcal{A}}_{12,1})\parallel \le {\phi }_5\parallel {\mathcal{A}}_{12}\left(t\right)-{\mathcal{A}}_{12}\left(t_1\right)\parallel ,\hfill & \hfill \\ \parallel {\mathfrak{H}}_6(t,{\mathcal{I}}_1)-{\mathfrak{H}}_6(t,{\mathcal{I}}_{1,1})\parallel \le {\phi }_6\parallel {\mathcal{I}}_1\left(t\right)-{\mathcal{I}}_1\left(t_1\right)\parallel ,\hfill & \hfill \\ \parallel {\mathfrak{H}}_7(t,{\mathcal{I}}_2)-{\mathfrak{H}}_7(t,{\mathcal{I}}_{2,1})\parallel \le {\phi }_7\parallel {\mathcal{I}}_2\left(t\right)-{\mathcal{I}}_2\left(t_1\right)\parallel .\hfill & \hfill \\ \parallel {\mathfrak{H}}_8(t,{\mathcal{I}}_{12})-{\mathfrak{H}}_8(t,{\mathcal{I}}_{12,1})\parallel \le {\phi }_8\parallel {\mathcal{I}}_{12}\left(t\right)-{\mathcal{I}}_{12}\left(t_1\right)\parallel .\hfill & \hfill \\ \parallel {\mathfrak{H}}_9(t,\mathcal{R})-{\mathfrak{H}}_9(t,{\mathcal{R}}_1)\parallel \le {\phi }_9\parallel \mathcal{R}\left(t\right)-\mathcal{R}\left(t_1\right)\parallel .\hfill & \hfill \end{array}\right.\end{array}$$

(21)

□

3.4. Existence of the Solution

In this part, we prove the existence of the solution of the model. Model (9) can be expressed in the presence of the memory function:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathcal{S}\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_1(t,\mathcal{S})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_1(\top ,\mathcal{S}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ \mathcal{V}\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_2(t,\mathcal{V})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_2(\top ,\mathcal{V}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_1\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_3(t,{\mathcal{A}}_1)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_3(\top ,{\mathcal{A}}_1){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_2\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_4(t,{\mathcal{A}}_2)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_4(\top ,{\mathcal{A}}_2){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{A}}_{12}\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_5(t,{\mathcal{A}}_{12})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_5(\top ,{\mathcal{A}}_{12}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_1\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_6(t,{\mathcal{I}}_1)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_6(\top ,{\mathcal{I}}_1){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_2\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_7(t,{\mathcal{I}}_2)+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_7(\top ,{\mathcal{I}}_2){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ {\mathcal{I}}_{12}\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_8(t,{\mathcal{I}}_{12})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_8(\top ,{\mathcal{I}}_{12}){(t-\top )}^{\varsigma -1}d\top ,\hfill & \hfill \\ \mathcal{R}\left(t\right)=\mathsf{{\rm Y}}\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_9(t,\mathcal{R})+\wp \left(\varsigma \right){\int }_0^t{\mathfrak{H}}_9(\top ,\mathcal{R}){(t-\top )}^{\varsigma -1}d\top .\hfill & \hfill \end{array}\right.\end{array}$$

(22)

Theorem 4.

Suppose the linear growth condition holds and $\mathcal{S}\left(t\right)$ is the solution of the governing model. Then, we have

$$\begin{array}{c}\hfill \underset{t\in [0,\beth ]}{sup}{\parallel \mathcal{S}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2{+1+3(1+\parallel \mathcal{S}\left(t\right)\parallel }^2)ϵ{\aleph }^2\left(\varsigma \right)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\end{array}$$

(23)

Proof. 

Let $\forall n\ge 1$:

$$\begin{array}{c}\hfill {\gamma }_n\left(t\right)=inf\{\beth ,inf\{t\in [0,\beth ]:\parallel \mathcal{S}\left(t\right)\parallel \ge n\left\}\right\};\end{array}$$

(24)

one can observe that ${lim}_{n\to \infty }{\gamma }_n\to A$ and also set

$$\begin{array}{c}\hfill {\mathcal{S}}_n\left(t\right)=inf\{t,{\gamma }_n\},\forall t\in [0,\beth ].\end{array}$$

(25)

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2-{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2=3{\aleph }^2\left(\varsigma \right){\parallel {\mathfrak{H}}_1(t,{\mathcal{S}}_n)\parallel }^2+3\parallel \left(\wp \left(\varsigma \right)\right){\int }_0^t\left[{\mathfrak{H}}_1(\top ,{\mathcal{S}}_n){(t-\top )}^{\varsigma -1}\right]d\top {\parallel }^2.\end{array}$$

(26)

Applying the Hölder inequality and the linear growth condition on (26),

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2-{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2=3{\aleph }^2\left(\varsigma \right)ϵ(1+\parallel {\mathcal{S}}_n\left(t\right){\parallel }^2)+3\left(\frac{t^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\left(1+\parallel {\mathcal{S}}_n{(\top )\parallel }^2\right)d\top .\end{array}$$

(27)

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2 < \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2+1\le {3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_n\left(t\right){\parallel }^2)+\\ \hfill {\int }_0^{t}3\left(\frac{t^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right)\left(1 + \parallel {\mathcal{S}}_n{(\top )\parallel }^2\right)d\top .\end{array}$$

(28)

By taking the supremum on both sides of (28),

$$\begin{array}{c}\hfill \begin{array}{cc}& \underset{\top \in [0,\beth ]}{sup}\parallel {\mathcal{S}}_n{(\top )\parallel }^2\le 1+\underset{\top \in [0,\beth ]}{sup}{\parallel {\mathcal{S}}_n(\top )\parallel }^2\le \hfill \\ & \underset{\top \in [0,\beth ]}{sup}\left\{1+{3\parallel \gamma (\top )\parallel }^2+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_n(\top ){\parallel }^2)\right\}\hfill \\ & +3\left(\frac{{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\left(1 + \underset{{\sigma }^{\varsigma }\in [0,\top ]}{sup}{\parallel {\mathcal{S}}_n\left({\sigma }^{\varsigma }\right)\parallel }^2\right)d\top .\hfill \end{array}\end{array}$$

(29)

Applying the Gronwall inequality on (29),

$$\begin{array}{c}\hfill \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{S}}_n\left(t\right)\parallel }^2+1\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_n\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\end{array}$$

(30)

$$\begin{array}{c}\hfill \underset{t\in [0,{\gamma }_n]}{sup}{\parallel {\mathcal{S}}_n\left(t\right)\parallel }^2\le \underset{t\in [0,{\gamma }_n]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_n\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\end{array}$$

(31)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{n\to \infty }{lim}\underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{S}}_n\left(t\right)\parallel }^2\le & \underset{n\to \to \infty }{lim}\underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_n\left(t\right){\parallel }^2)\right\}\times \hfill \\ & exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\hfill \end{array}\end{array}$$

(32)

In a similar way, we can verify

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{V}}_n\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{V}}_n\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{A}}_{1,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{A}}_{1,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{A}}_{2,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{A}}_{2,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{A}}_{12,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{A}}_{12,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{I}}_{1,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{I}}_{1,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{I}}_{2,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1+\parallel {\mathcal{I}}_{2,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{I}}_{12,n}\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{I}}_{12,n}\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\hfill & \hfill \\ \underset{t\in [0,\beth ]}{sup}{\parallel {\mathcal{R}}_n\left(t\right)\parallel }^2\le \underset{t\in [0,\beth ]}{sup}\left\{{3\parallel \mathsf{{\rm Y}}\left(t\right)\parallel }^2+1+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{R}}_n\left(t\right){\parallel }^2)\right\}\times exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\hfill & \hfill \end{array}\right.\end{array}$$

(33)

Now, the Picard recursions approach assists in proving the existence of the solution; thus, the recursive formula for System (9) upon setting

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathsf{{\rm Y}}\left(t\right)={\mathcal{S}}_0\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{V}}_0\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{A}}_{1,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{A}}_{2,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{A}}_{12,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{I}}_{1,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{I}}_{2,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{I}}_{12,0}\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}),\hfill & \hfill \\ \mathsf{{\rm Y}}\left(t\right)={\mathcal{R}}_0\left(t\right)\in {\beth }^2([0,\beth ],\mathfrak{R}).\hfill & \hfill \end{array}\right.\end{array}$$

Taking the first equation from System (9),

$$\begin{array}{c}\hfill {\mathcal{S}}_n\left(t\right)={\mathcal{S}}_0\left(t\right)+\aleph \left(\varsigma \right){\mathfrak{H}}_1(t,{\mathcal{S}}_{n-1}\left(t\right))+\wp \left(\varsigma \right){\int }_0^t\left[{\mathfrak{H}}_1(\top ,{\mathcal{S}}_{n-1}(\top )){(t-\top )}^{\varsigma -1}\right]d\top .\end{array}$$

(34)

The induction technique is applied, such that ${\mathcal{S}}_n\left(t\right)\in {\beth }^2[0,\beth ],\mathfrak{R},\forall n\ge 1$:

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2=3\parallel {\mathcal{S}}_0{\left(t\right)\parallel }^2+3{\aleph }^2\left(\varsigma \right)\parallel {\mathfrak{H}}_1(t,{\mathcal{S}}_{n-1}\left(t\right)){\parallel }^2+\left(\frac{3{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t{\parallel {\mathfrak{H}}_1(\top ,{\mathcal{S}}_{n-1})(\top )\parallel }^{2}d\top .\end{array}$$

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2=3\parallel {\mathcal{S}}_0{\left(t\right)\parallel }^2+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_{n-1}{\left(t\right)\parallel }^2)+\left(\frac{3{\beth }^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t(1 + \parallel {\mathcal{S}}_{n-1}(\top ){\parallel }^2)d\top .\end{array}$$

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2\le {\mu }_1+\left(\frac{3{\beth }^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t{\parallel {\mathcal{S}}_{n-1}(\top )\parallel }^{2}d\top ,\end{array}$$

(35)

where

$$\begin{array}{c}\hfill {\mu }_1=3\parallel {\mathcal{S}}_0{\left(t\right)\parallel }^2+3{\aleph }^2\left(\varsigma \right)ϵ(1 + \parallel {\mathcal{S}}_{n-1}\left(t\right){\parallel }^2)+\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\end{array}$$

and taking the max on both sides of (35),

$$\begin{array}{c}\hfill \underset{1\le n\le i}{max}\parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2\le \underset{1\le n\le i}{max}{\mu }_1+\left(\frac{3{\beth }^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\parallel {\mathcal{S}}_0{(\top )\parallel }^2+\underset{1\le n\le i}{max}{\parallel {\mathcal{S}}_n(\top )\parallel }^{2}d\top .\end{array}$$

(36)

$$\begin{array}{c}\hfill \underset{1\le n\le i}{max}\parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2\le \underset{1\le n\le i}{max}{\mu }_1+\left(\frac{3A^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\parallel {\mathcal{S}}_0{(\top )\parallel }^2+\underset{1\le n\le i}{max}{\parallel {\mathcal{S}}_n(\top )\parallel }^{2}d\top ,\end{array}$$

(37)

where

$$\begin{array}{c}\hfill {\mu }_2=\underset{1\le n\le i}{max}{\mu }_1+\left(\frac{3{\beth }^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\underset{1\le j\le \top }{sup}{\parallel {\mathcal{S}}_0\left(j\right)\parallel }^{2}d\top .\end{array}$$

Thus, the Gronwall inequality yields

$$\begin{array}{c}\hfill \underset{1\le n\le i}{max}{\parallel {\mathcal{S}}_n\left(t\right)\parallel }^2\le {\mu }_{2}exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right).\end{array}$$

(38)

Taking an arbitrary i,

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_n{\left(t\right)\parallel }^2\le {\mu }_{2}exp\left(\frac{3{\beth }^{\varsigma +1}ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right),\forall t\in [0,\beth ].\end{array}$$

(39)

Now,

$$\begin{array}{c}\hfill \begin{array}{cc}& \parallel {\mathcal{S}}_1\left(t\right)-{\mathcal{S}}_0{\left(t\right)\parallel }^2={\parallel {\mathcal{S}}_1\left(t\right)-\mathsf{{\rm Y}}\left(t\right)\parallel }^2\hfill \\ & \le 2\parallel {\aleph }^2\left(\varsigma \right)\left({\mathfrak{H}}_1(t,{\mathcal{S}}_{n-1}\left(t\right))\right){\parallel }^2+2\parallel \left(\frac{3{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right){\int }_0^t\left({\mathfrak{H}}_1(\top ,{\mathcal{S}}_{n-1})(\top )\right)d\top {\parallel }^2\hfill \\ & \le 2\left(\aleph \left(\varsigma \right)\right)ϵ(1+\underset{t\in [0,\beth ]}{sup}\parallel {\mathcal{S}}_0\left(t\right)\parallel )+\frac{2{\beth }^{\varsigma }ϵ}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\int }_0^t(1+\underset{j\in [0,\beth ]}{sup}\parallel {\mathcal{S}}_0\left(j\right)\parallel )d\top \hfill \end{array}\end{array}$$

(40)

Similarly,

$$\begin{array}{c}\hfill \begin{array}{cc}& \parallel {\mathcal{S}}_2\left(t\right)-{\mathcal{S}}_1{\left(t\right)\parallel }^2\le 2{\aleph }^2\left(\varsigma \right){\phi }_1\parallel {\mathcal{S}}_1\left(t\right)-{\mathcal{S}}_0{\left(t\right)\parallel }^2+\frac{2{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\phi }_1{\int }_0^t{\parallel {\mathcal{S}}_1(\top )-{\mathcal{S}}_0(\top )\parallel }^{2}d\top \hfill \\ & \le \mathsf{\Delta }\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\phi }_1\right].\hfill \end{array}\end{array}$$

(41)

In the same way,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel {\mathcal{S}}_3\left(t\right)-{\mathcal{S}}_2{\left(t\right)\parallel }^2\le & 2{\aleph }^2\left(\varsigma \right){\phi }_1\underset{1\le t\le \beth }{sup}{\parallel {\mathcal{S}}_2\left(t\right)-{\mathcal{S}}_1\left(t\right)\parallel }^2+\hfill \\ & \frac{2{\beth }^{\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\int }_0^t\underset{1\le j\le \top }{sup}{\parallel {\mathcal{S}}_2\left(j\right)-{\mathcal{S}}_1\left(j\right)\parallel }^{2}d\top ,\hfill \end{array}\\ \hfill \le \mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^2.\end{array}$$

(42)

We claim that

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_{n+1}\left(t\right)-{\mathcal{S}}_n{\left(t\right)\parallel }^2\le \mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n,\end{array}$$

(43)

is valid $\forall n\ge 1$. Clearly, one can obtain the requested result upon plugging in $n=0$. We need to show its validity for $n+1$.

$$\begin{array}{c}\hfill \begin{array}{cc}& \parallel {\mathcal{S}}_{n+2}\left(t\right)-{\mathcal{S}}_{n+1}{\left(t\right)\parallel }^2\le 2{\aleph }^2\left(\varsigma \right){\phi }_1\parallel {\mathcal{S}}_{n+1}\left(t\right)-{\mathcal{S}}_n{\left(t\right)\parallel }^2+\frac{2{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\phi }_1{\int }_0^t{\parallel {\mathcal{S}}_{n+1}\left(j\right)-{\mathcal{S}}_n\left(j\right)\parallel }^{2}d\top \hfill \\ & \le 2{\aleph }^2\left(\varsigma \right){\phi }_1\underset{1\le t\le A}{sup}\parallel {\mathcal{S}}_{n+1}\left(t\right)-{\mathcal{S}}_n{\left(t\right)\parallel }^2+\frac{2{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\phi }_1{\int }_0^t\underset{1\le j\le \top }{sup}{\parallel {\mathcal{S}}_{n+1}\left(j\right)-{\mathcal{S}}_n\left(j\right)\parallel }^{2}d\top \hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \le \mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^{n+1}.\end{array}$$

(44)

Finally, we acquire upon the induction hypothesis

$$\begin{array}{c}\hfill \parallel {\mathcal{S}}_{n+2}\left(t\right)-{\mathcal{S}}_{n+1}{\left(t\right)\parallel }^2\le \mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\phi }_1\right]}^{n+1}.\end{array}$$

(45)

A feasible choice of ${\phi }_1$ yields

$$\begin{array}{c}\hfill \sum _{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_1+\frac{2{\beth }^{2\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty .\end{array}$$

(46)

Similarly,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_2+\frac{2{\beth }^{2\varsigma }{\phi }_2}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_3+\frac{2{\beth }^{2\varsigma }{\phi }_3}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_4+\frac{2{\beth }^{2\varsigma }{\phi }_4}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_5+\frac{2{\beth }^{2\varsigma }{\phi }_5}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_6+\frac{2{\beth }^{2\varsigma }{\phi }_6}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_7+\frac{2{\beth }^{2\varsigma }{\phi }_7}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_8+\frac{2{\beth }^{2\varsigma }{\phi }_8}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty ,\hfill & \hfill \\ {\sum }_{n=0}^{\infty }\mathsf{\Delta }{\left[2{\aleph }^2\left(\varsigma \right){\phi }_9+\frac{2{\beth }^{2\varsigma }{\phi }_9}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\right]}^n<\infty .\hfill & \hfill \end{array}\right.\end{array}$$

(47)

□

3.5. Uniqueness of the Solution

Let $\mathcal{S}\left(t\right)$ and ${\mathcal{S}}_1$ be two solutions of the fractional governing model, then

$$\begin{array}{c}\hfill \mathcal{S}\left(t\right)-{\mathcal{S}}_1\left(t\right)=\aleph \left(\varsigma \right)({\mathfrak{H}}_1(t,\mathcal{S})-{\mathfrak{H}}_1(t,{\mathcal{S}}_1))+\wp \left(\varsigma \right){\int }_0^t\left[({\mathfrak{H}}_1(\top ,\mathcal{S})-{\mathfrak{H}}_1(\top ,{\mathcal{S}}_1)){(t-\top )}^{\varsigma -1}\right]d\top .\end{array}$$

Evaluating the norm, such as

$$\begin{array}{c}\hfill \parallel \mathcal{S}\left(t\right)-{\mathcal{S}}_1{\left(t\right)\parallel }^2=2\parallel \aleph \left(\varsigma \right)({\mathfrak{H}}_1(t,\mathcal{S})-{\mathfrak{H}}_1(t,{\mathcal{S}}_1)){\parallel }^2+2\\ \hfill \times \parallel \wp \left(\varsigma \right){\int }_0^t\left[({\mathfrak{H}}_1(\top ,\mathcal{S})-{\mathfrak{H}}_1(\top ,{\mathcal{S}}_1)){(t-\top )}^{\varsigma -1}\right]d\top {\parallel }^2,\end{array}$$

$$\begin{array}{c}\hfill \parallel \mathcal{S}\left(t\right)-{\mathcal{S}}_1{\left(t\right)\parallel }^2<2{\aleph }^2\left(\varsigma \right){\parallel ({\mathfrak{H}}_1(t,\mathcal{S})-{\mathfrak{H}}_1(t,{\mathcal{S}}_1))\parallel }^2+\\ \hfill \frac{2{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\int }_0^t{\parallel ({\mathfrak{H}}_1(\top ,\mathcal{S})-{\mathfrak{H}}_1(\top ,{\mathcal{S}}_1))\parallel }^{2}d\top .\end{array}$$

Applying the Lipschitz condition,

$$\begin{array}{c}\hfill \parallel \mathcal{S}\left(t\right)-{\mathcal{S}}_1{\left(t\right)\parallel }^2\le 2{\aleph }^2\left(\varsigma \right){\phi }_1\parallel \mathcal{S}-{\mathcal{S}}_1{\parallel }^2+\frac{2{\beth }^{\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\int }_0^t{\parallel \mathcal{S}-{\mathcal{S}}_1\parallel }^{2}d\top ,\end{array}$$

$$\begin{array}{c}\hfill \parallel \mathcal{S}\left(t\right)-{\mathcal{S}}_1{\left(t\right)\parallel }^2\le \frac{2{\beth }^{\varsigma }{\phi }_1}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}\parallel \frac{1}{1-2{\left(\aleph \left(\varsigma \right)\right)}^2{\phi }_1}\parallel \times {\int }_0^t{\parallel \mathcal{S}-{\mathcal{S}}_1\parallel }^{2}d\top ,\end{array}$$

with

$$\begin{array}{c}\hfill \left(\frac{1}{1-2{\left(\aleph \left(\varsigma \right)\right)}^2{\phi }_1}\right)\ne 0,\end{array}$$

the suggested Gronwall inequality is

$$\begin{array}{c}\hfill \parallel \mathcal{S}-{\mathcal{S}}_1{\parallel }^2=0.\end{array}$$

Therefore, we see that

$$\begin{array}{c}\hfill \mathcal{S}={\mathcal{S}}_1.\end{array}$$

(48)

Similarly,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathcal{V}={\mathcal{V}}_1,\hfill & \hfill \\ {\mathcal{A}}_1={\mathcal{A}}_{1,1},\hfill & \hfill \\ {\mathcal{A}}_2={\mathcal{A}}_{2,1},\hfill & \hfill \\ {\mathcal{A}}_{12}={\mathcal{A}}_{12,1},\hfill & \hfill \\ {\mathcal{I}}_1={\mathcal{I}}_{1,1},\hfill & \hfill \\ {\mathcal{I}}_2={\mathcal{I}}_{2,1},\hfill & \hfill \\ {\mathcal{I}}_{12}={\mathcal{I}}_{12,1},\hfill & \hfill \\ \mathcal{R}={\mathcal{R}}_1.\hfill & \hfill \end{array}\right.\end{array}$$

(49)

With the evidence of (49), we can say that the generalised fractional governing model (8) has a unique solution.

3.6. The Basic Reproduction Number of the Model

The DFE of the model (7) is

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_0& =\left({\mathcal{S}}^{*},{\mathcal{V}}^{*},{\mathcal{A}}_1^{*},{\mathcal{A}}_2^{*},{\mathcal{A}}_{12}^{*},{\mathcal{I}}_1^{*},{\mathcal{I}}_2^{*},{\mathcal{I}}_{12}^{*},{\mathcal{R}}_2^{*}\right),\hfill \\ & =\left({\mathcal{S}}^{*},{\mathcal{V}}^{*},0,0,0,0,0,0,0\right),\hfill \end{array}$$

with

$${\mathcal{S}}^{*}=\frac{(1-\chi )\mathsf{\Pi }}{\varpi +\mu },{\mathcal{V}}^{*}=\frac{(\mu \chi +\varpi )\mathsf{\Pi }}{\mu (\varpi +\mu )},$$

following the approach of [71]. The transfer matrices for the model are given, respectively,

$$F=\left(\begin{array}{cccccc}\frac{{\beta }_1[{\mathcal{S}}^{*}+(1-{\gimel }_1){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}& 0& 0& \frac{{\beta }_1{\theta }_1[{\mathcal{S}}^{*}+(1-{\gimel }_1){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}& 0& 0\\ 0& \frac{{\beta }_2[{\mathcal{S}}^{*}+(1-{\gimel }_2){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}& 0& 0& \frac{{\beta }_2{\theta }_2[{\mathcal{S}}^{*}+(1-{\gimel }_2){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}& 0\\ 0& 0& \frac{{\beta }_{12}[{\mathcal{S}}^{*}+(1-{\gimel }_{12}){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}& 0& 0& \frac{{\beta }_{12}{\theta }_{12}[{\mathcal{S}}^{*}+(1-{\gimel }_{12}){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$$

$$V=\left(\begin{array}{cccccc}{\mathcal{L}}_1& 0& 0& 0& 0& 0\\ 0& {\mathcal{L}}_2& 0& 0& 0& 0\\ 0& 0& {\mathcal{L}}_3& 0& 0& 0\\ -{\varsigma }_1& 0& 0& {\mathcal{L}}_4& 0& 0\\ 0& -{\varsigma }_2& 0& 0& {\mathcal{L}}_5& 0\\ 0& 0& -{\varsigma }_{12}& 0& 0& {\mathcal{L}}_6\end{array}\right)$$

(50)

where

$${\mathcal{L}}_1={\varsigma }_1+{\zeta }_1+\mu ,{\mathcal{L}}_2={\varsigma }_2+{\zeta }_2+\mu ,{\mathcal{L}}_3={\varsigma }_{12}+{\zeta }_{12}+\mu ,{\mathcal{L}}_4=q_1+{\mathsf{\Delta }}_1+\mu ,$$

$${\mathcal{L}}_5=q_2+{\mathsf{\Delta }}_2+\mu ,{\mathcal{L}}_6=q_{12}+{\mathsf{\Delta }}_{12}+\mu .$$

The basic reproduction number of the model (7) is given by ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\max \{{\mathcal{R}}_{0D},{\mathcal{R}}_{0M},{\mathcal{R}}_{0DM}\}$, where ${\mathcal{R}}_{0D}$, ${\mathcal{R}}_{0M}$, and ${\mathcal{R}}_{0DM}$ are the associated reproduction numbers for the Delta and Omicron SARS-CoV-2 variants and their co-infection, respectively, given by

$${\mathcal{R}}_{0D}=\frac{{\beta }_1\left({\mathcal{L}}_4+-{\varsigma }_1{\theta }_1\right)[(1-{\gimel }_1)(\varpi +\mu )+{\gimel }_1\varpi \mu ]}{\mu (\varpi +\mu ){\mathcal{L}}_1{\mathcal{L}}_4{\mathcal{N}}^{*}},$$

$${\mathcal{R}}_{0M}=\frac{{\beta }_2\left({\mathcal{L}}_5+-{\varsigma }_2{\theta }_2\right)[(1-{\gimel }_2)(\varpi +\mu )+{\gimel }_2\varpi \mu ]}{\mu (\varpi +\mu ){\mathcal{L}}_2{\mathcal{L}}_5{\mathcal{N}}^{*}},$$

$${\mathcal{R}}_{0DM}=\frac{{\beta }_{12}\left({\mathcal{L}}_6+-{\varsigma }_{12}{\theta }_{12}\right)[(1-{\gimel }_{12})(\varpi +\mu )+{\gimel }_{12}\varpi \mu ]}{\mu (\varpi +\mu ){\mathcal{L}}_3{\mathcal{L}}_6{\mathcal{N}}^{*}}.$$

3.7. Local Asymptotic Stability of the Disease-Free Equilibrium of the Model

Theorem 5.

The DFE ${\mathcal{Z}}_0$ of the model (7) is locally asymptotically stable (LAS) if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. 

The local stability of the model (7) was analysed by the Jacobian matrix of System (7) evaluated at the disease-free equilibrium ${\mathcal{Z}}_0$, given by

$$\left(\begin{array}{ccccccccc}-(\varpi +\mu )& 0& -\frac{{\beta }_1}{{\mathcal{N}}^{*}}{\mathcal{S}}^{*}& -\frac{{\beta }_2}{{\mathcal{N}}^{*}}{\mathcal{S}}^{*}& -\frac{{\beta }_{12}}{{\mathcal{N}}^{*}}{\mathcal{S}}^{*}& -\frac{{\beta }_1}{{\mathcal{N}}^{*}}{\theta }_1{\mathcal{S}}^{*}& -\frac{{\beta }_2}{{\mathcal{N}}^{*}}{\theta }_2{\mathcal{S}}^{*}& -\frac{{\beta }_{12}}{{\mathcal{N}}^{*}}{\theta }_{12}{\mathcal{S}}^{*}& 0\\ \varpi & -\mu & -\frac{{\beta }_1^{\prime }}{{\mathcal{N}}^{*}}{\mathcal{V}}^{*}& -\frac{{\beta }_2^{\prime }}{{\mathcal{N}}^{*}}{\mathcal{V}}^{*}& -\frac{{\beta }_{12}^{\prime }}{{\mathcal{N}}^{*}}{\mathcal{V}}^{*}& -\frac{{\beta }_1^{\prime }}{{\mathcal{N}}^{*}}{\theta }_1{\mathcal{V}}^{*}& -\frac{{\beta }_2^{\prime }}{{\mathcal{N}}^{*}}{\theta }_2{\mathcal{V}}^{*}& -\frac{{\beta }_{12}^{\prime }}{{\mathcal{N}}^{*}}{\theta }_{12}{\mathcal{V}}^{*}& 0\\ 0& 0& \frac{{\beta }_1{\mathcal{K}}_1^{*}}{{\mathcal{N}}^{*}}-{\mathcal{L}}_1& 0& 0& \frac{{\beta }_1{\theta }_1{\mathcal{K}}_1^{*}}{{\mathcal{N}}^{*}}& 0& 0& 0\\ 0& 0& 0& \frac{{\beta }_2{\mathcal{K}}_2^{*}}{{\mathcal{N}}^{*}}-{\mathcal{L}}_2& 0& 0& \frac{{\beta }_2{\theta }_2{\mathcal{K}}_2^{*}}{{\mathcal{N}}^{*}}& 0& 0\\ 0& 0& 0& 0& \frac{{\beta }_{12}{\mathcal{K}}_{12}^{*}}{{\mathcal{N}}^{*}}-{\mathcal{L}}_3& 0& 0& \frac{{\beta }_{12}{\theta }_{12}{\mathcal{K}}_{12}^{*}}{{\mathcal{N}}^{*}}& 0\\ 0& 0& {\varsigma }_1& 0& 0& 0& -{\mathcal{L}}_4& 0& 0\\ 0& 0& 0& {\varsigma }_2& 0& 0& 0& -{\mathcal{L}}_5& 0\\ 0& 0& 0& 0& {\varsigma }_{12}& 0& 0& 0& -{\mathcal{L}}_6\\ 0& 0& {\zeta }_1& {\zeta }_2& {\zeta }_{12}& q_1& q_2& q_{12}& -\mu \end{array}\right)$$

where ${\mathcal{K}}_i^{*}={\mathcal{S}}^{*}+(1-{\gimel }_i){\mathcal{V}}^{*}$ and ${\beta }_i^{{}^{\prime }}=(1-{\gimel }_i){\beta }_i,(i=1,2,12)$. The eigenvalues are given by

$$\begin{array}{c}\hfill {\phi }_1=-\mu ,{\phi }_2=-\mu ,{\phi }_3=-\mu ,{\phi }_4=-(\mu +\varpi ),\end{array}$$

and the solutions of the characteristic polynomial equations:

$${\phi }^2+\left({\mathcal{L}}_1+{\mathcal{L}}_2-\frac{{\beta }_1[{\mathcal{S}}^{*}+(1-\gimel ){\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}\right)\phi +{\mathcal{L}}_1{\mathcal{L}}_2(1-{\mathcal{R}}_{0W})=0,$$

(51)

and

$${\phi }^3+{\chi }_{11}{\phi }^2+{\chi }_{22}\phi +{\chi }_{33}=0,$$

(52)

where

$$\begin{array}{cc}\hfill {\chi }_{11}& =\left({\mathcal{L}}_3+{\mathcal{L}}_4+{\mathcal{L}}_5-\frac{{\beta }_2[{\mathcal{S}}^{*}+\varrho (1-\varsigma ){\varphi }_v{\mathcal{V}}^{*}]}{{\mathcal{N}}^{*}}\right)\hfill \\ \hfill {\chi }_{22}& =\left({\mathcal{L}}_3{\mathcal{L}}_4+{\mathcal{L}}_3{\mathcal{L}}_5+{\mathcal{L}}_4{\mathcal{L}}_5-\frac{{\beta }_2({\mathcal{L}}_4+{\mathcal{L}}_5+{\varsigma }_2{\theta }_2){\mathcal{S}}^{*}}{{\mathcal{N}}^{*}}-\frac{{\beta }_2\varrho (1-\varsigma )({\varsigma }_3{\theta }_2+{\mathcal{L}}_3{\varphi }_v+{\mathcal{L}}_5{\varphi }_v){\mathcal{V}}^{*}}{{\mathcal{N}}^{*}}\right)\hfill \\ \hfill {\chi }_{33}& ={\mathcal{L}}_3{\mathcal{L}}_4{\mathcal{L}}_5(1-{\mathcal{R}}_{0D})=0.\hfill \end{array}$$

From the Routh–Hurwitz criterion, Equation (51) has roots with negative real parts if and only if ${\mathcal{R}}_{0W}<1$. Furthermore, Equation (52) has roots with negative real parts if and only if ${\chi }_{11}>0$, ${\chi }_{33}>0$, and ${\chi }_{11}{\chi }_{22}>{\chi }_{33}$. If ${\mathcal{R}}_{0D}<1$, ${\chi }_{11}>0$, and ${\chi }_{33}>0$, hence, the DFE, ${\mathcal{Z}}_0$, is locally asymptotically stable if ${\mathcal{R}}_0=\max \{{\mathcal{R}}_{0W},{\mathcal{R}}_{0D}\}<1$.

□

Since Im $\left({\phi }_k\right)=0$, for $k=0,1,2,3,...,9$,

$$|\arg \left({\phi }_k\right)|=\pi >\frac{\varsigma \pi }{2},\mathrm{for}0<\varsigma <1.$$

4. Stability Analysis

Here, we discuss the stability analysis of the designed fractional governing model.

Hyers–Ulam–Rassias Stability

The stability analysis of the Atangana–Baleanu fractional COVID-19 model is presented here.

Let us re-write Model (9):

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{}^{ABC}{D}_t^{\varsigma }[\gimel \left(t\right)]=\mathfrak{H}(t,\gimel \left(t\right)),\hfill & \hfill \\ \gimel \left(0\right)={\gimel }_0,0<t<\beth <\infty .\hfill \end{array}\right.\end{array}$$

(53)

where the vectors $\gimel \left(t\right)=(\mathcal{S}\left(t\right),\mathcal{V}\left(t\right),{\mathcal{A}}_1\left(t\right),{\mathcal{A}}_2\left(t\right),{\mathcal{A}}_{12}\left(t\right),{\mathcal{I}}_1\left(t\right),{\mathcal{I}}_2\left(t\right),{\mathcal{I}}_{12}\left(t\right),\mathcal{R}\left(t\right))$ and $\mathcal{A}=({\mathfrak{H}}_1,{\mathfrak{H}}_2,{\mathfrak{H}}_3,{\mathfrak{H}}_4,{\mathfrak{H}}_5,{\mathfrak{H}}_6,{\mathfrak{H}}_7,{\mathfrak{H}}_8,{\mathfrak{H}}_9)$ are continuous vector functions.

Definition 3

([72]). Assume the fractional order is $0<\varsigma <1$ and $\mathfrak{H}:[0,\beth ]\times {\mathfrak{R}}_{+}^9\to {\mathfrak{R}}_{+}^9$ is a continuous mapping. Then, Model (53) is Hyers–Ulam-stable if $\exists \mathcal{K}>0$ and $\mathcal{N}>0$ such that, for each solution $\gimel \in V([0,\beth ],{\mathfrak{R}}_{+}^9)$, the following inequality exists:

$$\begin{array}{c}\hfill \left|\right|{}^{\mathcal{ABC}}{D}_t^{\varsigma }[\gimel \left(t\right)]-\mathfrak{H}(t,\gimel \left(t\right))\left|\right|\le \mathcal{N},\forall t\in [0,\beth ],\end{array}$$

(54)

∃ a solution ${\gimel }^{\prime }\in V([0,\beth ],{\mathfrak{R}}_{+}^9)$ of the model (53) such as

$$\begin{array}{c}\hfill \left|\right|\gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\left|\right|\le \mathcal{K}\mathcal{N},\forall t\in [0,\beth ].\end{array}$$

(55)

Definition 4

([72]). Assume the fractional order is $0<\varsigma <1$. The function $\mathfrak{H}:[0,\beth ]\times {\mathfrak{R}}_{+}^9\to {\mathfrak{R}}_{+}^9$ and $\mathsf{\Omega }:[0,\beth ]\to {\mathfrak{R}}_{+}^9$ are continuous mappings. Then, Model (53) is generalised Hyers–Ulam–Rassias-stable with respect to Ω if $\exists V_{\mathfrak{H},\omega }>0$ such that, for each solution $\gimel \in V([0,\beth ],{\mathfrak{R}}_{+}^9)$, the following inequality exists:

$$\begin{array}{c}\hfill \left|\right|{}^{ABC}{D}_t^{\varsigma }[\gimel \left(t\right)]-\mathfrak{H}(t,\gimel \left(t\right))\left|\right|\le \mathsf{\Omega }\left(t\right),\forall t\in [0,\beth ],\end{array}$$

(56)

∃ a solution ${\gimel }^{\prime }\in V([0,\beth ],{\mathfrak{R}}_{+}^9)$ of the model (53) such as

$$\begin{array}{c}\hfill \left|\right|\gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\left|\right|\le V_{\mathfrak{H},\mathsf{\Omega }}\mathsf{\Omega }\left(t\right),\forall t\in [0,\beth ].\end{array}$$

(57)

Now, to prove that Model (53) is generalised Hyers–Ulam–Rassias-stable, we assume that:

$\left[C_1\right]\mathfrak{H}:[0,\beth ]\times {\mathfrak{R}}_{+}^9\to {\mathfrak{R}}_{+}^9$ is a continuous mapping.

$\left[C_2\right]\exists V_{\mathfrak{H}}>0$ such that, for each solution $\gimel ,{\gimel }^{\prime }\in V([0,\beth ],{\mathfrak{R}}_{+}^9),$

$$\begin{array}{c}\hfill \left|\right|\gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\left|\right|\le V_{\mathfrak{H}}\left|\right|\gimel -{\gimel }^{\prime }\left|\right|,\forall t\in [0,\beth ].\end{array}$$

$\left[C_3\right]$ Let $\mathsf{\Omega }\in ([0,\beth ],{\mathfrak{R}}^{+})$ be an increasing mapping, and let $W_{\mathsf{\Omega }}>0$ such that

$${\int }_0^t\mathsf{\Omega }(\top )d\top \le W_{\mathsf{\Omega }}\mathsf{\Omega }\left(t\right),\forall \top \in [0,\beth ].$$

Theorem 6.

Assuming that $\left[C_1\right]-\left[C_3\right]$ exist, Model (53) is generalised Hyers–Ulam–Rassias-stable with respect to Ω on the interval provided that $\aleph \left(\varsigma \right)V_{\mathfrak{H}}<1$.

Proof. 

Let ${\gimel }^{\prime }\in V([0,\beth ],{\mathfrak{R}}_{+}^9)$ be a solution of Model (53). Then, the unique solution of Model (53) is as proven above.

$$\begin{array}{c}\hfill \gimel \left(t\right)=\gimel \left(0\right)+\aleph \left(\varsigma \right)\mathfrak{H}(t,\gimel \left(t\right))+\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,\gimel (\top )){(t-\top )}^{\varsigma -1}d\top .\end{array}$$

(58)

With the evidence of (56), we can say

$$\begin{array}{c}\hfill \parallel {\gimel }^{\prime }\left(t\right)-\gimel \left(0\right)+\aleph \left(\varsigma \right)\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))+\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,{\gimel }^{\prime }(\top )){(t-\top )}^{\varsigma -1}d\top \parallel \\ \hfill \le \aleph \left(\varsigma \right)\mathsf{\Omega }\left(t\right)+\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{\int }_0^t\mathsf{\Omega }(\top )d\top ,\\ \hfill \le \left(\aleph \left(\varsigma \right)+\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{W}_{\mathsf{\Omega }}\right)\mathsf{\Omega }\left(t\right).\end{array}$$

Therefore,

$$\begin{array}{c}\hfill \parallel \gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\parallel \le \parallel {\gimel }^{\prime }\left(t\right)-\gimel \left(0\right)-\aleph \left(\varsigma \right)\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))-\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,{\gimel }^{\prime }(\top )){(t-\top )}^{\varsigma -1}d\top \parallel ,\end{array}$$

$$\begin{array}{c}\hfill \le \parallel {\gimel }^{\prime }\left(t\right)-\gimel \left(0\right)-\aleph \left(\varsigma \right)\mathfrak{H}(t,\gimel \left(t\right))-\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,\gimel (\top )){(t-\top )}^{\varsigma -1}d\top \\ \hfill -\aleph \left(\varsigma \right)\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))-\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,{\gimel }^{\prime }(\top )){(t-\top )}^{\varsigma -1}d\top \\ \hfill +\aleph \left(\varsigma \right)\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))+\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,{\gimel }^{\prime }(\top )){(t-\top )}^{\varsigma -1}d\top \parallel ,\end{array}$$

$$\begin{array}{c}\hfill \le \parallel {\gimel }^{\prime }\left(t\right)-\gimel \left(0\right)-\aleph \left(\varsigma \right)\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))-\wp \left(\varsigma \right){\int }_0^t\mathfrak{H}(\top ,{\gimel }^{\prime }(\top )){(t-\top )}^{\varsigma -1}d\top \parallel +\\ \hfill \aleph \left(\varsigma \right)\parallel \mathfrak{H}(t,\gimel \left(t\right))-\mathfrak{H}(t,{\gimel }^{\prime }\left(t\right))\parallel +\wp \left(\varsigma \right){\int }_0^t\parallel (\mathfrak{H}(\top ,\gimel (\top ))-\mathfrak{H}(\top ,{\gimel }^{\prime }(\top ))){(t-\top )}^{\varsigma -1}\parallel d\top ,\end{array}$$

$$\begin{array}{c}\hfill \le \left(\aleph \left(\varsigma \right)+\frac{A{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{W}_{\mathsf{\Omega }}\right)\mathsf{\Omega }\left(t\right)+\aleph \left(\varsigma \right)V_{\mathfrak{H}}\parallel \gimel \left(t\right)-{\gimel }^{\prime }\left(t\right))\parallel \\ \hfill +\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{V}_{\mathfrak{H}}{\int }_0^t\parallel \gimel (\top )-{\gimel }^{\prime }(\top )\parallel d\top .\end{array}$$

(59)

Now, $\aleph \left(\varsigma \right)V_H<1$, so

$$\begin{array}{c}\hfill \parallel \gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\parallel \le \frac{(\aleph \left(\varsigma \right)+\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{W}_{\mathsf{\Omega }})\mathsf{\Omega }\left(t\right)}{1-\aleph \left(\varsigma \right)V_{\mathfrak{H}}}+\frac{\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{V}_{\mathfrak{H}}}{1-\aleph \left(\varsigma \right)V_{\mathfrak{H}}}{\int }_0^t\parallel \gimel (\top )-{\gimel }^{\prime }(\top )\parallel d\top .\end{array}$$

(60)

Gronwall’s inequality yields

$$\begin{array}{c}\hfill \parallel \gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\parallel \le \left[\frac{\aleph \left(\varsigma \right)+\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{W}_{\mathsf{\Omega }}}{1-\aleph \left(\varsigma \right)V_{\mathfrak{H}}}exp\left(t\right)\right]\mathsf{\Omega }\left(t\right).\end{array}$$

(61)

With the setting $V_{\mathfrak{H},\mathsf{\Omega }}=\left[\frac{\aleph \left(\varsigma \right)+\frac{{\beth }^{\varsigma }}{AB\left(\varsigma \right)\mathsf{\Gamma }\left(\varsigma \right)}{W}_{\mathsf{\Omega }}}{1-\aleph \left(\varsigma \right)V_H}exp\left(t\right)\right]$, we have

$$\begin{array}{c}\hfill \parallel \gimel \left(t\right)-{\gimel }^{\prime }\left(t\right)\parallel \le V_{\mathfrak{H},\mathsf{\Omega }}\mathsf{\Omega }\left(t\right).\end{array}$$

(62)

Inequality (62) has authenticated that Model (53) is Hyers–Ulam–Rassias-stable with respect to $\mathsf{\Omega }$. □

5. Simulations of the SARS-CoV-2 Model (7)

In this part, we give a graphical analysis of the numerical results for the cases of India and Pakistan. For this, we estimated the real values of the natural death rates and transmission rates for both countries. In India, the life expectancy is 69.66 years, and the population is 211,120,000 [73]; therefore, we set the natural death rate to be $\frac{1,393,409,033}{69.66}$ per day, while the natural death rate was obtained as $\frac{1}{69.66}$ per day. The life expectancy for Pakistan is 67.27 years, and the population is 225,199,929 [74]; therefore, we set the natural death rate to be $\frac{225,199,929}{67.27}$ per day, while the natural death rate was obtained as $\frac{1}{67.27}$ per day. The initial conditions were obtained thusly for each of the countries. For India: as at 10 May 2021, cumulative cases recorded: 22,992,517; cumulative deaths recorded: 249,992. The initial conditions were set thusly: $\mathcal{S}\left(0\right)$ = 1,300,000,000, $\mathcal{V}\left(0\right)$ = 35,906,905, ${\mathcal{A}}_1\left(0\right)={\mathcal{I}}_1\left(0\right)$ = 7,641,000, ${\mathcal{R}}_1\left(0\right)=0$, ${\mathcal{A}}_2\left(0\right)={\mathcal{A}}_2^v\left(0\right)={\mathcal{I}}_2\left(0\right)$ = 7,641,000, ${\mathcal{R}}_2\left(0\right)=0$. For Pakistan: as at 1 July 2021, cumulative cases recorded: 959,685; cumulative deaths recorded: 22,345. The initial conditions were set thusly: $\mathcal{S}\left(0\right)$ = 180,000,000,$\mathcal{V}\left(0\right)$ = 20,000,000, ${\mathcal{A}}_1\left(0\right)={\mathcal{I}}_1\left(0\right)$ = 319,895, ${\mathcal{R}}_1\left(0\right)=0,{\mathcal{A}}_2\left(0\right)={\mathcal{A}}_2^v\left(0\right)={\mathcal{I}}_2\left(0\right)$ = 319,895, ${\mathcal{R}}_2\left(0\right)=0$. Using the data available for the daily cumulative number of confirmed SARS-CoV-2 deaths in India [73], the estimated contact rates were ${\beta }_1=1.9219\times {10}^{-4}$ and ${\beta }_2=1.9219\times {10}^{-4}$, while the estimated disease-induced death rates were ${\mathsf{\Delta }}_1=0.0171$ and ${\mathsf{\Delta }}_2=4.4185\times {10}^{-5}$, respectively. Likewise, from the data available to us for daily cumulative confirmed SARS-CoV-2 deaths in Pakistan [74], the estimated contact rates were ${\beta }_1=0.0099$ and ${\beta }_2=0.3954$, while the estimated disease-induced death rates were ${\mathsf{\Delta }}_1=0.0027$ and ${\mathsf{\Delta }}_2=1.1312\times {10}^{-5}$. The model fitting was performed based on the period when the Delta variant was dominant or first reported in the selected countries. For India, the simulation period for the fitting was from 10 May 2021 to 8 November 2021. This was almost around the time when the World Health Organization (WHO) declared the Delta variant a Variant Of Concern (VOC) [57]. For Pakistan, the simulation period was from 1 July 2021 (around the time the first case of the SARS-CoV-2 Delta variant was first reported in the country) to 8 November 2021. The estimated parameter values are shown in

Table 2. Estimated parameters for both countries using the available datasets.

	Parameter	Source
${\beta }_1$	$1.9219\times {10}^{-4}$	Fit
${\beta }_2$	$0.2982$	Fit
${\beta }_{12}$	$0.2982$	Fit
${\mathsf{\Delta }}_1$	0.0171	Fit
${\mathsf{\Delta }}_2$	$4.4185\times {10}^{-5}$	Fit
$\varpi$	0.08	Estimated
$\chi$	0.26	Estimated
${\mathcal{R}}_0$	2.0904	Fit

. For both countries, we show the dynamics of the proposed co-infection model in Figure 1, Figure 2, Figure 3 and Figure 4.

The simulation for India is depicted in Figure 1 and Figure 2. The number of susceptible people decreased over time since people were exposed to the coronavirus, which spread in a particular region in Figure 1a. The vaccination effect was also studied in the proposed model. When the vaccination ratio increased in India for the Delta variant, as observed in Figure 1b, as a result, the number of asymptomatic infected and symptomatic infected people decreased, which is displayed in Figure 1c and Figure 2b. This was due to vaccination leading to a stronger immune response against the Delta variant. There was no vaccination program against the Omicron variant. Therefore, the asymptomatic infected and symptomatic people infected with the Omicron variant increased, as shown in Figure 1d and Figure 2c. Likewise, as shown in Figure 2a,d, the number of coinfected asymptomatic and symptomatic people increased since there was no strategy to develop an immune response against the Omicron and Delta variants. Figure 2e denotes the recovered people, which increased with the passage of time. This happened because some people have naturally strong immunity and can fight off some diseases without taking medicine or vaccines. Similarly, the behaviour of each class for the Pakistan data was the same as for India.

The simulation for Pakistan is depicted in Figure 3 and Figure 4. The number of susceptible people decreased over time since people were exposed to the coronavirus, which spread in a particular region in Figure 3a. The vaccination effect was also studied in the proposed model. When the vaccination ratio increased in India for the Delta variant, as observed in Figure 3b, as a result, the number of asymptomatic infected and symptomatic infected people decreased, which is displayed in Figure 3c and Figure 4b. This was due to vaccination leading to a stronger immune response against the Delta variant. There was no vaccination program against the Omicron variant. Therefore, the asymptomatic infected and symptomatic people infected with the Omicron variant increased, as shown in Figure 3d and Figure 4c. Likewise, as shown in Figure 4a,d, the number of coinfected asymptomatic and symptomatic people increased since there was no strategy to develop an immune response against the Omicron and Delta variants. Figure 4e denotes the recovered people, which increased with the passage of time. This happened because some people have naturally strong immunity and can fight off some diseases without taking medicine or vaccines.

There were some differences in the behaviours of each class regarding the maximum value of each class. This was due to the the large population of India and different natural and transmission rates for both countries.

6. Conclusions

Nature can be understood in a better way and can even be predicted more accurately if the constraints of developed approaches, theories, and techniques are updated, questioned, and revised according to modern scientific discoveries and the occurrence of unexpected physical phenomena. The global non-singular Atangana–Baleanu operator with the presence of a memory function and the growth condition was used to build a generalised fractional model for the two variants of the SARS-COVID-19 outbreak in order to capture the complex dynamics. The Picard recursions approach validated the existence of the solution of the proposed fractional model, and the Banach fixed-point theory also verified the solution’s uniqueness. The model was shown to be globally asymptotically stable in the stability analysis of the Atangana–Baleanu fractional-order dual versions utilising the Hyers–Ulam–Rassias approach. Finally, the fractional Adams–Bashforth numerical approach was applied for the numerical simulation and graphical behaviour presented. A graphic analysis was performed to show the proposed model’s asymptotically global stability. We are confident about the established results of our model, which has the capability to provide more efficient and better predictions.

Our model was put forth based on the focal points of two distinct COVID-19 variants. Those infected with either variant or both variants were called co-infection classes in the model. Further research on COVID-19 co-infections with other diseases, such as TB, influenza, malaria, and other diseases, has been prompted by the appearance of many COVID-19 variants. While this is the first study on the co-infection of the Delta and Omicron variants that we are aware of, more research on the mathematics (stochastic, agent-based modelling, within-/intra-host) and epidemiological dynamics of this co-infection is needed. We had much difficulty estimating various parameters in our study due to the inadequacies and inconsistencies of a number of components and the characteristics of the new COVID-19 variants. Therefore, future research on the interactions between various COVID-19 variants using more accurate data and in-depth information is worthwhile.