Fractional Order Model for Co-Infection of Streptococcus pneumoniae and COVID-19 with Saturated Incidence Force of Infection

Chiekezi Zephaniah Opara, Ephraim Ngozi Erumaka, Nneka Onyinyechi Iheonu, Nnamdi Nonso Araka

Federal University of Technology, Owerri, Nigeria.

DOI: 10.4236/oalib.1112569                

Abstract

This study presented and analyzed the fractional order model for co-infection of Streptococcus pneumoniae and COVID-19 by adopting the Atangana-Baleanu derivative in Caputo sense. The solution’s boundedness and non-negativity are derived by applying the Laplace transform of the Atangana-Baleanu derivative in Caputo sense. We established the existence and uniqueness of the solutions of the proposed model using Atangana-Baleanu Caputo Integral, and Banach fixed point theorem. The disease-free equilibrium and basic reproduction number of the model were obtained. The stability of the model was investigated by applying Ulam-Hyers-Rassias stability. The deduced and presented fractional model is fitted to Nigeria’s confirmed daily COVID-19 cases as at February 2024. The numerical solution of the system is derived through the Adams-Bashforth linear multi-step method. The simulation of the overall patients co-infected with COVID-19 and Streptococcus pneumoniae, at different fractional orders was carried out.

Keywords

COVID-19 Pneumonia Co-Infection, Banach Fixed-Point Theorem, Atangana-Baleunu in Caputo Sense, Fractional Derivative

1. Background

COVID-19 pandemic has spread across over 150 countries, and led to headlines in the world in the year 2020. The COVID-19 is a contagious disease caused by the virus SARS-CoV-2. Many people thought the virus was first identified in December 2019 in the city of Wuhan, Hubei province, China. It was initially referred to as the “2019-n-CoV”, the virus was later officially named SARS-CoV-2. The COVID-19 disease is an infectious disease caused by the SARS-CoV-2 virus. It is primarily transmitted from an infected to a healthy individual through droplets of saliva or discharge from the nose, eyes, and mouth when an infected person coughs or sneezes [1].

Available data recorded over 775.6 million infection cases, with over 7 million deaths, and 675.6 million recoveries across the globe [2] [3]. Also, over 13.59 billion doses of the COVID-19 vaccine have been administered, with above 69.4% (5.5 billion individuals) of the world population vaccinated (received at least one dose of a COVID-19 vaccine) [1]. Similarly, the record also shows that Nigeria has confirmed 267,188 COVID-19 cases, discharged 259,953, and has recorded 3155 deaths, with 4080 active cases, with 17,371 total cases per one million populations, and 374 total deaths per one million populations across the 36 states, and has administered over 133.05 million doses of the several COVID-19 vaccine with 39% of her citizens fully vaccinated, and 46% received at least one dose of the vaccine as at February 2024 [1]-[3].

On the other hand, Streptococcus pneumoniae (S. pneumoniae), also known as pneumococcus, is a Gram-positive bacterium that commonly colonizes the human upper respiratory tract. While it is typically asymptomatic, it is also a leading cause of various respiratory and invasive infections, particularly in young children, the elderly, and immunocompromised individuals. Characterizing the mechanisms of pathogenesis, S. pneumoniae utilizes a variety of virulence factors to establish and maintain infections, including adhesins, pneumolysin, and various enzymes that contribute to host tissue damage and immune evasion. Pneumococcal infections are typically initiated by the colonization of the upper respiratory tract, followed by invasion and dissemination to other sites of the body, such as the lungs, bloodstream, and central nervous system [4]. Over the past thirty years, pneumonia has consistently been the primary factor contributing to the death of children worldwide, resulting in the loss of 671,928 lives among those under the age of five in 2019. Nigeria has the highest number of deaths, with 129,444 children under five dying from pneumonia in 2019 at a death rate of 386.15 per 100,000 [5]. Streptococcus pneumoniae is primarily transmitted through respiratory droplets produced when an infected person coughs or sneezes. The transmission can occur through direct contact with respiratory secretions from an infected individual or indirectly by touching surfaces contaminated with these secretions and then touching the face, particularly the mouth or nose [6]. Streptococcus pneumoniae has been detected in at least 63% of severe COVID-19 cases, and has been confirmed to be the common cause of complications, severe illness, and death of COVID-19 patients across the globe. Many studies have shown that co-infected patients have a higher risk of severe illness and death compared to those who only had COVID-19. The co-infection dynamics between COVID-19 and Streptococcus pneumoniae pose significant challenges to public health systems worldwide. Therefore, there is need to develop a mathematical model that will proffer an effective solution which will drastically reduce or eliminate the burden on he global population by the COVID-19 and its co-infections with other diseases at various stages. Hence the challenge of developing a fractional order model for co-infection of streptococcus pneumoniae and COVID-19 with saturated incidence force of infection.

Several mathematical model research have been conducted for COVID-19, and the co-infection of COVID-19 with other diseases. For instance, Omame in his fractional order model for the co-interaction of COVID-19 and Hepatitis B virus, fitted and estimated real data from Wuhan, China, and important parameters relating to each disease and their co-infection, and concluded that HBV and COVID-19 transmission rates can greatly impact the dynamics of the co-infection of both diseases [7]. His research submitted that to control the co-circulation of both diseases in a population, efforts must be geared toward preventing incident infection with either or both diseases. Baleanu analyzed and presented a fractional-order model for COVID-19 transmission with the Caputo-Fabrizio derivative [8]. They used the homotopy analysis transform method (HATM) to solve the model and also provided a solution in convergent series. Also, Omame in his fractional order model for dual variants of COVID-19 and HIV co-infection applied the Atangana-Baleanu derivative and concluded that the fractional derivatives had a great influence on the dynamics of the diseases across the epidemiological states over time which is as a result of the memory effect from the application of fractional derivative [9].

Aba in their study developed a mathematical model to depict the dynamics of COVID-19 with the consideration of quarantine, isolation, and environmental viral load [10]. The research showed that reducing the number of effective contacts with infected individuals leads to a substantial decline in the total number of infected cases. The model also highlighted that reducing the environmental viral load significantly reduces the pandemic’s peak. Additionally, the study examined the effectiveness of quarantine or contact-tracing policies and hospitalization or self-isolation measures. The results demonstrated that contact tracing of exposed individuals is more effective in eliminating COVID-19 infections.

Olumuyiwa in his research analyzed the fractional order of the pneumococcal pneumonia infection model using the Caputo Fabrizio operator [11]. He used fixed-point theory to prove the existence of a solution and computed the iterative solution by applying the fractional Adams-Bashforth method. In conclusion, he suggested that the dynamics of the disease were more apparent with the fractional order model than with the integer order model. Furthermore, Omame, Abbas, & Onyenegecha (2021) in their fractional-order model for COVID-19 and tuberculosis co-infection, applied Atangana-Baleauna derivative [12]. The study revealed that both the burden of COVID-19 and the co-infection of both diseases in the population will be reduced if the risk of COVID-19 infection by latently-infected TB individuals is reduced. Aga, Keno, Terfasa, & Berhe (2024) in their study, presented a nonlinear deterministic mathematical model for co-infection of pneumonia and COVID-19 transmission dynamics and applied the optimal control theory to describe the optimal control model that incorporates four controls, namely, prevention of pneumonia, prevention of COVID-19, treatment of infected pneumonia and treatment of infected COVID-19 and further concluded that the combination of several control variables is most effective in the reduction or elimination of the co-infection [13].

Research Objectives

The objective of this study is to analyze the fractional order model for COVID-19 and Streptococcus pneumoniae co-infection, using the Atangana-Baleanu fractional order derivative in Caputo sense, establish the solution’s positivity and Boundedness of the model using Laplace transform of the Atangana-Baleanu, to establish the existence and uniqueness of the model solutions by applying the Atangana-Baleanu Caputo Integral, and the fixed-point theorem, to obtain the disease free equilibrium points and the basic reproductive number, and prove the Ulam-Hyers stability of the solutions. We also analyzed the numerical solution of the model.

2. Method

We focus particularly on the Atangana-Baleanu fractional derivatives and integrals in Caputo Sense (ABC), which are characterized by their non-local and non-singular properties. These derivatives provide a more accurate representation of our model compared to other common fractional derivatives. The ABC fractional derivatives not only facilitate easier computation but also allow for more comprehensive analysis within our framework.

Definition 2.1 (Atangana and, Baleanu, 2016). Consider a function f belonging to the Sobolev space $H^1\left(a_1\mathrm{,}{a}_2\right)$ with $a_2>a_1$ and a fractional order $\vartheta \in \left[\mathrm{0,1}\right]$ . The Atangana-Baleanu (AB) fractional derivative of $f\left(t\right)$ of order $\vartheta \in {ℝ}^{+}$ is defined as:

${}_{a_1}^{abc}{D}_t^{\vartheta }f\left(t\right)=\frac{\beta \left(\vartheta \right)}{1-\vartheta }{\displaystyle {\int }_{a_1}^t}f\left(t\right)E_{\vartheta }\left[-\vartheta \frac{{\left(t-ϵ\right)}^{\vartheta }}{1-\vartheta }\right]\text{d}ϵ$ (1)

where $\beta \left(\vartheta \right)$ is a normalization function that satisfies $\beta \left(0\right)=\beta \left(1\right)=1$ , and $E_{\vartheta }(\cdot )$ denotes the Mittag-Leffler function, which generalizes the exponential function in the context of fractional calculus.

Definition 2.2 (Atangana and, Baleanu, 2016) For the same function $f\in H^1\left(a_1\mathrm{,}{a}_2\right)$ , the Atangana-Baleanu (AB) fractional derivative of order $\vartheta \in {ℝ}^{+}$ can also be represented as:

${}_{a_1}^{abc}{D}_t^{\vartheta }f\left(t\right)=\frac{1-\vartheta }{\beta \left(\vartheta \right)f\left(t\right)}+\frac{\vartheta }{\beta \left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_{a_1}^t}f\left(t\right){\left(t-ϵ\right)}^{\vartheta -1}\text{d}ϵ$ (2)

This alternative form highlights the dual nature of the AB fractional derivative, combining both a local term and a non-local integral term [14] [15].

Lemma 2.1 (Atangana and, Baleanu, 2016) For a function $f\in H^1\left(a_1\mathrm{,}{a}_2\right)$ , the Newton-Leibniz formula for the AB fractional integral is given by:

${}_{a_1}^{abc}{I}_1^{\vartheta }\left({}_{a_1}^{abc}{D}_t^{\vartheta }f\left(t\right)\right)=f\left(t\right)-f\left(a_1\right)$ (3)

This lemma establishes a fundamental relationship between the fractional derivative of AB and its corresponding integral.

Theorem 2.1 (Atangana and, Baleanu, 2016) Consider a continuous function $f\in C\left[a\mathrm{,}b\right]$ . The AB fractional derivatives satisfy the inequality:

$\Vert {}_{a_1}^{abc}{D}_t^{\vartheta }{f}_1\left(t\right)-{}_{a_1}^{abc}{D}_t^{\vartheta }{f}_2\left(t\right)\Vert <\delta \Vert f_1\left(t\right)-f_2\left(t\right)\Vert$ (4)

where $\delta$ is a positive constant. This theorem demonstrates the stability of the AB fractional derivatives under small perturbations of the function.

Theorem 2.2 (Atangana and, Baleanu, 2016) The unique solution to the fractional differential equation of order $\vartheta$ ,

${}_{a_1}^{abc}{D}_t^{\vartheta }f\left(t\right)=r\left(t\right)$

is given by:

${}_{a_1}^{abc}{D}_t^{\vartheta }r\left(t\right)=\frac{1-\vartheta }{\beta \left(\vartheta \right)r\left(t\right)}+\frac{\vartheta }{\beta \left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_{a_1}^t}f\left(t\right){\left(t-ϵ\right)}^{\vartheta -1}\text{d}ϵ$ (5)

This theorem provides the explicit form of the AB fractional derivative for the given differential equation.

Lemma 2.2. Let $\omega \in {ℝ}^m$ be a differentiable vector function such that for all $t\ge 0$ ,

${}_0^{abc}{D}_t^{\vartheta }\left({\left[\omega -{\omega }^{\mathrm{<mo>*</mo>}}\right]}^{\text{T}}Q\left[\omega -{\omega }^{\mathrm{<mo>*</mo>}}\right]\right)\le {\left[\omega -{\omega }^{\mathrm{<mo>*</mo>}}\right]}^{\text{T}}Q{}_0^{abc}{D}_t^{\vartheta }\left[\omega -{\omega }^{\mathrm{<mo>*</mo>}}\right]$ (6)

where $Q\in {ℝ}^{m\times m}$ is a constant, symmetric, positive definite matrix. This lemma provides a condition for the fractional derivative of a quadratic form involving $\omega$ .

3. Formulation of Model

3.1. Assumptions of the Model

The assumptions of the model are stated below:

1) The population is homogeneously mixed (A population that interacts with one another to the same degree and fixed).

2) Individuals infected with COVID-19 are susceptible to infection with Streptococcus pneumoniae and vice versa.

3) Co-infected individuals can transmit either COVID-19 or Streptococcus pneumoniae but not the mixed infections at the same time.

4) Co-infected infected individuals can recover either from COVID-19 or Streptococcus pneumoniae but not from the mixed infections at the same time.

5) Rate of transmissibility for singly infected and co-infected individuals are assumed same.

6) We assume that the saturated incidence force of infection rate for both disease is

${\pi }_i=\frac{{\beta }_i}{1+{\alpha }_1{I}_i+{\alpha }_2{I}_{ij}}\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}i\mathrm{,}j=\text{coinfected}\text{diseases}$

7) People in the Vaccinated compartment can be infected for both diseases.

8) That for both diseases, people can be infected only through contacts with infectious people.

9) Recovered individuals from both diseases may go back to the susceptible class.

The epidemiological states of the model at a given time is divided into ten compartments: unvaccinated susceptible humans $S_h\left(t\right)$ , vaccinated individuals against COVID-19 $V_c\left(t\right)$ , vaccinated individuals against S. pneumoniae $V_p\left(t\right)$ , COVID-19 infected only individuals $I_c\left(t\right)$ , isolated individuals with COVID-19 $Q_c\left(t\right)$ , COVID-19 recovered individuals (non-S. pneumoniae infection) $R_c\left(t\right)$ S. pneumoniae infected (non-COVID-19 infection) $I_p\left(t\right)$ , S. pneumoniae recovered individuals (non-COVID-19 infection) $R_p\left(t\right)$ , Individuals with co-infections $I_{cp}\left(t\right)$ , Individuals with co-infections (hospitalized) $H_{cp}\left(t\right)$ .

Individuals are recruited into the population at the rate, ${\Theta }_h$ . Susceptible humans $S\left(t\right)$ , acquire COVID-19 disease at the rate

$\frac{{\beta }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)}{1+{\alpha }_1{I}_c+{\alpha }_2{I}_{cp}}={\pi }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)$ where, ${\pi }_c=\frac{{\beta }_c}{1+{\alpha }_1{I}_c+{\alpha }_2{I}_{cp}}$

denotes the saturated incidence force of S. pneumoniae infection. Also is the face-mask compliance rate and denotes the efficacy of the face-mask for both diseases.

Furthermore, it is assumed that natural mortality is the same and occurs across all the model compartment, at the rate $\mu$ . COVID-19 Vaccinated susceptible, $V_c\left(t\right)$ , acquire COVID-19 at a reduced rate $\frac{{\beta }_c\left(1-{\chi }_c\right)\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)}{1+{\alpha }_1{I}_c+{\alpha }_2{I}_{cp}}={\pi }_c\left(1-{\chi }_c\right)\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)$ , where ${\chi }_c$ is the COVID-19 vaccine efficacy.

While, S. pneumoniae vaccinated susceptible $V_p\left(t\right)$ , acquires S. pneumoniae at a reduced rate of $\frac{{\beta }_p\left(1-{\chi }_p\right)\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)}{1+{\alpha }_1{I}_p+{\alpha }_2{I}_{cp}}\mathrm{<mo>=</mo>}{\pi }_p\left(1-{\chi }_p\right)\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)$ , where ${\chi }_p$ is the S. pneumoniae vaccine efficacy.

Also, we have assumed that the transmissibility rate for singly infected and co-infected persons is the same, as there is no evidence yet to justify otherwise. All other transitions of the model that are given could be found in the model Equation (1). As all the model parameters are well defined in Table 1 below.

3.2. Symbols and Parameters of the Model

Below in Table 1 are the symbols and parameters of the proposed model.

Table 1. Description of the model variables and parameters.

Variable	Description
$S_h$	Unvaccinated susceptible humans
$I_c\left(t\right)$	Individuals infected with COVID-19 only
$V_c\left(t\right)$	Vaccinated individuals against COVID-19
$Q_c$	Isolated individuals with COVID-19
$R_c$	COVID-19 recovered individuals
$I_p$	Individuals infected with S. pneumoniae only
$V_p$	Individuals vaccinated against pneumonia
$R_p$	S. pneumoniae recovered individuals
$I_{cp}$	Individuals with co-infection
$H_{cp}$	Individuals with co-infection (hospitalized)
Parameter	Description
${\beta }_c$	COVID-19 transmission rate
${\beta }_p$	pneumonia transmission rate
$\eta$	Face-mask compliance
$\vartheta$	Face-mask efficacy
$\Theta$	Recruitment rate into susceptible compartment
${\omega }_{vc}$	COVID-19 vaccination rate
${\chi }_c$	COVID-19 vaccine efficacy
${\gamma }_c$	Recovery rate of COVID-19 infected individuals (natural recovery)
${\theta }_c$	Per capita rate of COVID recovered individuals going back to the susceptible class
${\omega }_{vp}$	S. pneumoniae average vaccination rate
${\chi }_p$	S. pneumoniae vaccine efficacy
${\gamma }_p$	Recovery rate of S. pneumoniae infected individuals
${\theta }_p$	Per capita rate of pneumonia recovered individuals going back to the susceptible class
$\mu$	Natural death rate
${\nu }_1$	Rate for individuals in isolation compartment $Q_c$
${\nu }_2$	Rate for individuals in hospitalization compartment $H_{cp}$ class
${\rho }_1$ , ${\omega }_1$	COVID-19 treatment rate for individuals in $I_c$ and $Q_c$ respectively
${\rho }_2$ , ${\omega }_2$	COVID-19 treatment rate for individuals in $I_{cp}$ and $H_{cp}$ respectively
${\sigma }_1$	COVID-19 related death rate
${\sigma }_2$	S. pneumoniae-related death rate
${\kappa }_c$	Modification parameter accounting for increased susceptibility to COVID-19 infection
${\kappa }_p$	Modification parameter accounting for increased susceptibility to S. pneumoniae infection
$\pi =\frac{{\beta }_i}{1+{\alpha }_1{I}_i+{\alpha }_2{I}_{i\mathrm{,}j}}$	Saturated incidence force of infection; _1 and _2 are the saturation factor that measures the inhibitory effect of COVID-19 and S. pneumoniae respectively
${\gamma }_{cd}$	COVID-19 recovery rate for $I_{cp}$ and $H_{cp}$
${\gamma }_{pd}$	pneumonia recovery rate for $I_{cp}$ and $H_{cp}$ classes respectively

The total population at time t is given by

$N_h=S+V_c+I_c+Q_c+R_c+V_p+I_p+R_p+I_{cp}+H_{cp}+R_{cp}$

3.3. Flow Diagram of the Model

Figure 1. Model flow diagram.

Applying the symbols and parameters in Table 1, assumptions, and the model flow diagram in Figure 1, we now formulate the model equations. The model employs fractional calculus, specifically the Atangana-Baleanu (AB) fractional derivatives in Caputo sense, which offer a more comprehensive framework for capturing the memory and hereditary properties of biological processes. This approach allows us to better simulate the complex interactions and delayed effects inherent in the spread and progression of infectious diseases. The model is given thus:

(7)

with corresponding initial conditions

$\begin{array}{l}{S}_h\left(0\right)=S_h^0,V_c\left(0\right)=V_c^0,I_c\left(0\right)=I_c^0,Q_c\left(0\right)=Q_c^0,R_c\left(0\right)=R_c^0,V_p\left(0\right)=V_p^0\\ I_p\left(0\right)=I_p^0,R_p\left(0\right)=R_p^0,I_{cp}\left(0\right)=I_{cp}^0,H_{cp}\left(0\right)=H_{cp}^0,R_{cp}\left(0\right)=R_{cp}^0\end{array}$

4. Analysis of the Model

4.1. Boundedness of the Solution

The boundedness and non-negativity of the solutions which shows that the system (7) is both mathematically and biologically well posed is presented

$\mathcal{D}=\left\{\left(S_H+V_c+I_c+Q_c+R_c+V_p+I_p+R_p+I_{cp}+H_{cp}\right)\le \frac{\Theta }{{\mu }_h}\right\}$

Proof. Adding all the equations of the system 7 gives

$\frac{\text{d}{N}_h}{\text{d}t}=\Theta -{\mu }_h{N}_h\left(t\right)-\left({\sigma }_1{I}_c+{\sigma }_1{Q}_c+{\sigma }_2{I}_p+\left({\sigma }_1+{\sigma }_2\right)I_{cp}+\left({\sigma }_1+{\sigma }_2\right)H_{cp}\right)$ (8)

From the inequality given in Equation (8),

$\Theta -\left({\mu }_h+2\delta \right)N_h\le \frac{\text{d}{N}_h}{\text{d}t}\le \Theta -{\mu }_h{N}_h\mathrm{,}$ (9)

where $\delta =\min \left\{{\sigma }_1\mathrm{,}{\sigma }_2\right\}$ , we proceed by applying the Laplace transform.

Applying the Laplace transform to each term in the inequality, we have:

$ℒ\left\{\Theta -\left({\mu }_h+2\delta \right)N_h\right\}\le ℒ\left\{\frac{\text{d}{N}_h}{\text{d}t}\right\}\le ℒ\left\{\Theta -{\mu }_h{N}_h\right\}.$ (10)

Using the properties of the Laplace transform, we obtain

$\frac{\Theta }{s}-\left({\mu }_h+2\delta \right)ℒ\left\{N_h\left(t\right)\right\}\le sℒ\left\{N_h\left(t\right)\right\}-N_h\left(0\right)\le \frac{\Theta }{s}-{\mu }_hℒ\left\{N_h\left(t\right)\right\}\mathrm{.}$ (11)

Reorganizing to isolate $ℒ\left\{N_h\left(t\right)\right\}$ , we get the bounds:

$ℒ\left\{N_h\left(t\right)\right\}\le \frac{\frac{\Theta }{s}+N_h\left(0\right)}{s+{\mu }_h}\mathrm{,}$

$ℒ\left\{N_h\left(t\right)\right\}\ge \frac{\frac{\Theta }{s}+N_h\left(0\right)}{s+{\mu }_h+2\delta }\mathrm{.}$

Focusing on the upper bound and applying the inverse Laplace transform, we recognize that the solution involves the Mittag-Leffler function $E_{\vartheta \mathrm{,1}}\left(z\right)$ and $E_{\vartheta \mathrm{,}\vartheta +1}\left(z\right)$ . Therefore, the upper bound solution is given by:

$\begin{array}{c}{N}_h\left(t\right)\le \left(\frac{ℬ\left(\vartheta \right)}{ℬ\left(\vartheta \right)+\left(1-\vartheta \right){\mu }_h}{N}_h\left(0\right)+\frac{\left(1-\vartheta \right)\Theta }{ℬ\left(\vartheta \right)+\left(1-\vartheta \right){\mu }_h}\right)\\ +E_{\vartheta \mathrm{,1}}\left(-\frac{\vartheta {\mu }_h}{ℬ\left(\vartheta \right)+\left(1-\vartheta \right){\mu }_h}{t}^{\vartheta }\right)\\ +\frac{\vartheta \Theta }{ℬ\left(\vartheta \right){\mu }_h}{E}_{\vartheta \mathrm{,}\vartheta +1}\left(-\frac{\vartheta {\mu }_h}{ℬ\left(\vartheta \right)+\left(1-\vartheta \right){\mu }_h}{t}^{\vartheta }\right)\mathrm{.}\end{array}$ (12)

Given the asymptotic properties of the Mittag-Leffler function $E_{\alpha ,\beta }\left(z\right)$ , as $t\to \infty$ , we have:

$N_h\left(t\right)\le \frac{\Theta }{{\mu }_h}\mathrm{.}$ (13)

Thus, the system described by Equation (7) is positively invariant, meaning the solutions remain bounded within a region $\mathcal{D}$ as $t\to \infty$ .□

4.2. Non-Negativity of the Solution

Following the pattern in the work of (Ogunrinde et al. 2021), by contradiction we assume that equation two of the model is false. Then let $t_1=\min \left\{t:S\left(h\right),V_c\left(t\right),I_c\left(t\right),Q_c\left(t\right),R_c\left(t\right),V_p\left(t\right),I_p\left(t\right),R_p\left(t\right),I_c{}_p\left(t\right),H_c{}_p\left(t\right)=0\right\}$ . Suppose $V_c\left(t_1\right)=0$ , it implies that $S\left(t\right)>0$ , $I_c\left(t\right)>0$ , $Q_c\left(t\right)>0$ , $R_c\left(t\right)>0$ , $V_p\left(t\right)>0$ , $I_c\left(t\right)>0$ , $R_p\left(t\right)>0$ , $I_{cp}\left(t\right)>0$ , $H_{cp}\left(t\right)>0$ , for all $\left[\mathrm{0,}{t}_1\right]$ . We can assume that there exists the following expression,

${\theta }_1=\underset{0\le t\le t_t}{\min }\left\{{\omega }_{vc}{S}_h-\left[{\mu }_h+{\pi }_c\left(1-{\chi }_c\right)\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)\right]V_c-{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)V_c\right\}\mathrm{.}$

It follows that

$D_t^{\varrho }{V}_c-{\theta }_1{V}_c>0.$ (14)

We can also determine a continuous function ${\Phi }_1$ to ascertain the following equation

$D_t^{\varrho }{V}_c-{\theta }_1{V}_c=-{\Phi }_1\left(t\right)\mathrm{.}$

By Laplace transform, the above inequality becomes

$s^{\varrho }{\tilde{V}}_c\left(s\right)-s^{\varrho -1}{V}_c\left(0\right)-{\theta }_1{\tilde{V}}_c\left(s\right)=-{\tilde{\Phi }}_1\left(s\right)\mathrm{,}$

from which

$\begin{array}{c}{\tilde{V}}_c\left(s\right)=V_c\left(0\right)\frac{s^{\varrho -1}}{s^{\varrho }-{\theta }_1}-\frac{{\Phi }_1\left(s\right)}{s^{\varrho }-{\theta }_1}\mathrm{,}\\ =\frac{V_c\left(0\right)}{s}{\left(1-\frac{{\theta }_1}{s^{\varrho }}\right)}^{-1}-\frac{{\Phi }_1\left(s\right)}{s^{\varrho }}{\left(1-\frac{{\theta }_1}{s^{\varrho }}\right)}^{-1}\\ =V_c\left(0\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\theta }_1^k}{s^{\varrho k+1}}-{\Phi }_1\left(s\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\theta }_1^k}{s^{\varrho k+\varrho }}\mathrm{.}\end{array}$ (15)

Ignoring the non-positive term, the inverse Laplace transform gives the solution of (14) (using Mittag-Leffler function), which satisfies the following expression

$V_c>V_c\left(0\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\left({\theta }_1{t}^{\varrho }\right)}^k}{\Gamma \left(\omega k+1\right)}=V_c\left(0\right)E_{\varrho }\left({\theta }_1{t}^{\varrho }\right)\mathrm{.}$ (16)

such that the positivity of the solution $V_c$ is given by

$V_c>V_c\left(0\right)E_{\varrho }\left({\theta }_1{t}^{\varrho }\right)>0.$

which contradicts $V_c\left(t_1\right)=0$ .

Similarly, suppose $R_c\left(t_1\right)=0$ which implies that $S\left(t\right)>0$ , $V_c\left(t\right)>0$ , $I_c\left(t\right)>0$ , ${\varphi }_c\left(t\right)>0$ , $V_p\left(t\right)>0$ , $I_p\left(t\right)>0$ , $R_p\left(t\right)>0$ , $I_{cp}>0$ , $H_{cp}>0$ , for all $0\le t\le t_1$ . We assume that there exists the following expression

${\theta }_2=\underset{0\le t\le t_t}{\min }\left\{\left[\left(Y_c+p_1\right)I_c+\left(Y_c+{\omega }_1\right){\varphi }_c\right]+\left(\mu +{\theta }_c\right)R_c\right\}\mathrm{.}$

So that

$D_t^{\varrho }{R}_c\left(t\right)>{\theta }_2{R}_c\left(t\right)\mathrm{.}$ (17)

We can still determine a continuous function ${\Phi }_2\left(t\right)$ to ascertain the following equation

$D_t^{\varrho }{R}_c\left(t\right)-{\theta }_2{R}_c\left(t\right)=-{\Phi }_2\left(t\right)\mathrm{.}$ (18)

By Laplace transform the above inequality becomes

$s^{\varrho }{\tilde{R}}_c\left(s\right)-s^{\varrho -1}{R}_c\left(0\right)-{\theta }_2{\tilde{R}}_c\left(s\right)=-{\tilde{\Phi }}_2\left(s\right)\mathrm{,}$

from which

${\tilde{R}}_c\left(s\right)=R_c\left(0\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\theta }_2^k}{s^{\varrho k+1}}-{\Phi }_2\left(s\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\theta }_2^k}{s^{\varrho k+\varrho }}\mathrm{.}$ (19)

Ignoring the non-positive term, the inverse Laplace transform gives the solution of Equation (17) (using Mittag-Leffler function), satisfies the following expression

$R_c\left(t\right)>R_c\left(0\right){\displaystyle \sum _{k=0}^{\infty }}\frac{{\left({\theta }_2{t}^{\varrho }\right)}^k}{\Gamma \left(\varrho k+1\right)}=R_c\left(0\right)E_{\varrho }\left({\theta }_2{t}^{\varrho }\right)\mathrm{.}$ (20)

Hence the positivity of this other solution $R_c$ is given by

$R_c\left(t\right)>R_c\left(0\right)E_{\varrho }\left({\theta }_2{t}^{\varrho }\right)>0.$

and contradicts $R_c\left(t_1\right)=0$ . Moreso, since the above have similar results, the same pattern will show that the positivity of the solutions $S$ , $I_c$ , ${\varphi }_c$ , $V_p$ , $I_p$ , $R_p$ , $I_{cp}$ , and $H_{cp}$ respectively are given by

$S>S\left(0\right)E_{\varrho }\left({\theta }_3{t}^{\varrho }\right)>\mathrm{0,}$

$I_c>I_c\left(0\right)E_{\varrho }\left({\theta }_4{t}^{\varrho }\right)>\mathrm{0,}{\varphi }_c>{\varphi }_c\left(0\right)E_{\varrho }\left({\theta }_5{t}^{\varrho }\right)>\mathrm{0,}$

$V_p\left(t\right)>V_p\left(0\right)E_{\varrho }\left({\theta }_6{t}^{\varrho }\right)>\mathrm{0,}{I}_p\left(t\right)>I_p\left(0\right)E_{\varrho }\left({\theta }_7{t}^{\varrho }\right)>\mathrm{0,}$

$R_p\left(t\right)>R_p\left(0\right)E_{\varrho }\left({\theta }_8{t}^{\varrho }\right)>\mathrm{0,}{I}_{cp}>I_{cp}\left(0\right)E_{\varrho }\left({\theta }_9{t}^{\varrho }\right)>\mathrm{0,}$

$H_{cp}>H_{cp}\left(0\right)E_{\varrho }\left({\theta }_{10}{t}^{\varrho }\right)>\mathrm{0,}{R}_{cp}>R_{cp}\left(0\right)E_{\varrho }\left({\theta }_{11}{t}^{\varrho }\right)>0.$

4.3. Existence and Uniqueness of Solution

In this case we consider the shortened expression for our model 7 so as to establish the existence and uniqueness of solution

$\{\begin{array}{l}{}_0^{ABC}{D}_t^{m}G\left(t\right)=\mathcal{G}\left(t\mathrm{,}G\left(t\right)\right)\mathrm{,}\\ G\left(0\right)=G_0\mathrm{,}\end{array}$

where the vector $G=\left(S_h\mathrm{,}{V}_c\mathrm{,}{I}_c\mathrm{,}{Q}_c\mathrm{,}{R}_c\mathrm{,}{V}_p\mathrm{,}{I}_p\mathrm{,}{R}_p\mathrm{,}{I}_{cp}\mathrm{,}{H}_{cp}\right)$ denotes the model compartments and $\mathcal{G}=\left({\mathcal{G}}_1\mathrm{,}{\mathcal{G}}_2\mathrm{,}{\mathcal{G}}_3\mathrm{,}{\mathcal{G}}_4\mathrm{,}{\mathcal{G}}_5\mathrm{,}{\mathcal{G}}_6\mathrm{,}{\mathcal{G}}_7\mathrm{,}{\mathcal{G}}_8\mathrm{,}{\mathcal{G}}_9\mathrm{,}{\mathcal{G}}_{10}\right)$ denotes the continuous vector function such that

(21)

Definition 4.1. A Lipschitz condition for the second argument of a function $f\left(t\right)$ is defined as $\left|f\left(t\mathrm{,}{y}_1\right)-f\left(t\mathrm{,}{y}_2\right)\right|\le L\left|y_1-y_2\right|$ , where L is the positive constant, and $\left(t\mathrm{,}{y}_1\right)\mathrm{,}\left(t\mathrm{,}{y}_2\right)$ are a set $D\in {ℝ}^2$ .

To prove our result, we consider the following assumptions $\begin{array}{l}{S}_h\left(t\right),{\tilde{S}}_h\left(t\right),V_c\left(t\right),{\tilde{V}}_c\left(t\right),I_c\left(t\right),{\tilde{I}}_c\left(t\right),Q_c\left(t\right),{\tilde{Q}}_c\left(t\right),R_c\left(t\right),{\tilde{R}}_c\left(t\right),V_p\left(t\right),{\tilde{V}}_p\left(t\right),\\ I_p\left(t\right),{\tilde{I}}_p\left(t\right),R_p\left(t\right),{\tilde{R}}_p\left(t\right),I_c{}_p\left(t\right),{\tilde{I}}_c{}_p\left(t\right),H_c{}_p\left(t\right),{\tilde{H}}_c{}_p\left(t\right)\in z\left[0,1\right]\end{array}$ be continuous function, such that

$\begin{array}{l}\Vert S_h\left(t\right)\Vert \le z_1\mathrm{,}\Vert V_c\left(t\right)\Vert \le z_2\mathrm{,}\Vert I_c\left(t\right)\Vert \le z_3\mathrm{,}\Vert Q_c\left(t\right)\Vert \le z_4\mathrm{,}\Vert R_c\left(t\right)\Vert \le z_5\mathrm{,}\\ \Vert V_p\left(t\right)\Vert \le z_6,\Vert I_p\left(t\right)\Vert \le z_7\mathrm{,}\Vert R_p\left(t\right)\Vert \le z_8\mathrm{,}\Vert I_{cp}\left(t\right)\Vert \le z_9\mathrm{,}\Vert H_{cp}\left(t\right)\Vert \le z_1{}_0\end{array}$

for non-negative constant $z_1,z_2,z_3,z_4,z_5,z_6,z_7,z_8,z_9,z_{10}>0$ .

Theorem 4.1. The Lipschitz condition satisfy the model Equation (21) for $i\in {ℝ}_{10}$ , if the above assumptions holds true and fullfils an $\psi <1$ for $i\in {ℝ}^{10}$ .

Proof. We prove that $z\left(t\mathrm{,}{S}_h\right)$ fulfills the Lipschitz condition above for $S_h\left(t\right)\mathrm{,}{\tilde{S}}_h\left(t\right)$ , we get

$\begin{array}{l}\Vert z_1\left(t\mathrm{,}{S}_h\right)-z_1\left(t\mathrm{,}{\tilde{S}}_h\right)\Vert \\ =\Vert [\Theta -\left[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)+{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]S_h\\ +{\theta }_c{R}_c+{\theta }_p{R}_p]S_h-[\Theta -[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)\\ +{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)]S_h+{\theta }_c{R}_c+{\theta }_p{R}_p]{\tilde{S}}_h\Vert \\ \le \Theta -\left[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)+{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]S_h\\ +{\theta }_c{R}_c+{\theta }_p{R}_p\Vert S_h-{\tilde{S}}_h\Vert \\ \le {\psi }_1\Vert S_h-{\tilde{S}}_h\Vert \end{array}$(22)

where $\begin{array}{l}{\psi }_1=\Theta -\left[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_c+I_{cp}\right)+{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]S_h+{\theta }_c{R}_c\\ +{\theta }_p{R}_p\end{array}$ .

Hence, $z_1$ satisfies the Lipschitz condition with Lipschitz constant ${\psi }_1$ . Similarly, the other kernels satisfy the Lipschitz condition as follows

$\begin{array}{l}\Vert z_2\left(t\mathrm{,}{V}_c\right)-z_2\left(t\mathrm{,}{\tilde{V}}_c\right)\Vert \le {\psi }_2\Vert V_c-{\tilde{V}}_c\Vert \\ \Vert z_3\left(t\mathrm{,}{I}_c\right)-z_3\left(t\mathrm{,}{\tilde{I}}_c\right)\Vert \le {\psi }_3\Vert I_c-{\tilde{I}}_c\Vert \\ \Vert z_4\left(t\mathrm{,}{Q}_c\right)-z_4\left(t\mathrm{,}{\tilde{Q}}_c\right)\Vert \le {\psi }_4\Vert Q_c-{\tilde{Q}}_c\Vert \\ \Vert z_5\left(t\mathrm{,}{R}_c\right)-z_5\left(t\mathrm{,}{\tilde{R}}_c\right)\Vert \le {\psi }_5\Vert R_c-{\tilde{R}}_c\Vert \\ \Vert z_6\left(t\mathrm{,}{V}_p\right)-z_6\left(t\mathrm{,}{\tilde{V}}_p\right)\Vert \le {\psi }_6\Vert V_p-{\tilde{V}}_p\Vert \\ \Vert z_7\left(t\mathrm{,}{I}_p\right)-z_7\left(t\mathrm{,}{\tilde{I}}_p\right)\Vert \le {\psi }_7\Vert I_p-{\tilde{I}}_p\Vert \\ \Vert z_8\left(t\mathrm{,}{R}_p\right)-z_8\left(t\mathrm{,}{\tilde{R}}_p\right)\Vert \le {\psi }_8\Vert R_p-{\tilde{R}}_p\Vert \\ \Vert z_9\left(t\mathrm{,}{I}_{cp}\right)-z_9\left(t\mathrm{,}{\tilde{I}}_{cp}\right)\Vert \le {\psi }_9\Vert I_{cp}-{\tilde{I}}_{cp}\Vert \\ \Vert z_1{}_0\left(t\mathrm{,}{H}_{cp}\right)-z_1{}_0\left(t\mathrm{,}{\tilde{H}}_{cp}\right)\Vert \le {\psi }_1{}_0\Vert H_{cp}-{\tilde{H}}_{cp}\Vert \end{array}$ (23)

Hence, all the kernels $z_i$ for $i\in {ℝ}_{10}$ satisfies the Lipschitz properties with Lipschitz constant ${\psi }_1<1$ for $i\in {ℝ}_{10}$ .

Applying the Atangana-Baleanu fractional integral operator in Caputo sense on the model 7 and utilizing the initial conditions, we obtain the following Voltera integral equations

$\begin{array}{l}{S}_h\left(t\right)=S\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}S\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}S\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ V_c\left(t\right)=V_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{V}_c\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}V\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ I_c\left(t\right)=I_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{I}_c\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{I}_c\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ Q_c\left(t\right)=Q_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{Q}_c\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{Q}_c\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ R_c\left(t\right)=S\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{R}_c\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{R}_c\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ V_p\left(t\right)=V_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{V}_p\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{V}_p\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ I_p\left(t\right)=I_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{I}_p\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{I}_p\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ R_p\left(t\right)=R_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{R}_p\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{R}_p\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ I_{cp}\left(t\right)=I_{cp}\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{I}_{cp}\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{I}_{cp}\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\\ H_{cp}\left(t\right)=H_{cp}\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{H}_{cp}\left(t\right)\right)+\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{H}_{cp}\left(ϵ\right)\right)\text{d}ϵ\mathrm{,}\end{array}$ (24)

Theorem 4.2 The continuous function ${\mathcal{G}}_k$ , for $k=\mathrm{1,2,}\cdots \mathrm{,10}$ satisfies the Lipschitz condition for the second argument if and only if

$\begin{array}{l}{\sup }_{0<t\le 1}\Vert S_h-{\tilde{S}}_h\Vert \le z_1,\mathrm{<mi>\; </mi>}{\sup }_{0<t\le 1}\Vert V_c-{\tilde{V}}_c\Vert \le z_2,\mathrm{<mi>\; </mi>}{\sup }_{0<t\le 1}\Vert I_c-{\tilde{I}}_c\Vert \le z_3,\mathrm{<mi>\; </mi>}\\ {\sup }_{0<t\le 1}\Vert Q_c-{\tilde{Q}}_c\Vert \le z_4,\mathrm{<mi>\; </mi>}{z}_j>0,\mathrm{<mi>\; </mi>}z=1,2,3,4.\end{array}$

Proof. Let $C^1\left(\left[\mathrm{0,}T\right]\right)\mathrm{<mo>:</mo>}ℝ\to ℝ$ be the Banach space of all the continuous and differentiable functions from $\left[\mathrm{0,}T\right]\to ℝ$ so that using the fixed point theorem, we show that the function ${\mathcal{G}}_k\mathrm{,}k=\mathrm{1,2,}\cdots \mathrm{,10}$ is Lipschitz continuous with Lipschitz constants. Let ${\overline{S}}_h$ be the second solution and using triangle inequality endowed with Chebychev norm we have the following relationships

$\begin{array}{l}\Vert {\mathcal{G}}_1\left(S_h\right)-{\mathcal{G}}_1\left(\overline{S}\right)\Vert \\ \le \Vert (\Theta -\left[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)+{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]S_h\\ +{\theta }_c{R}_c+{\theta }_p{R}_p)-(\Theta -[{\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\\ +{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)]{\overline{S}}_h+{\theta }_c{R}_c+{\theta }_p{R}_p)\Vert \end{array}$

$\begin{array}{l}\le \Vert {\mu }_h+{\omega }_{vp}+{\omega }_{vc}+{\pi }_c\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)+{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\Vert \Vert S_h-{\overline{S}}_h\Vert \\ \le \Vert {\mu }_h+{\omega }_{vp}+{\omega }_{vc}+\left(1-\eta \rho \right)\left(z_1+z_2\right)+\left(1-\eta \rho \right)\left(z_3+z_4\right)\Vert \mathrm{<mi>\; </mi>}\Vert S_h-{\overline{S}}_h\Vert \\ =Z_{{\mathcal{G}}_1}\Vert S_h-{\overline{S}}_h\Vert \end{array}$

where $Z_{{\mathcal{G}}_1}={\mu }_h+{\omega }_{vp}+{\omega }_{vc}+\left(1-\eta \rho \right)\left(z_1+z_2\right)+\left(1-\eta \rho \right)\left(z_3+z_4\right)$ .

Similarly, we applied same for ${\mathcal{G}}_2\left(t\mathrm{,}{V}_c\left(t\right)\right)$ up to ${\mathcal{G}}_1{}_0\left(t\mathrm{,}{V}_c\left(t\right)\right)$ , having $Z_{{\mathcal{G}}_2}={\mu }_h+\left(1-{\chi }_c\right)\left(1-\eta \rho \right)\left(z_1+z_2\right)-\left(1-\eta \rho \right)\left(z_3+z_4\right)$

$Z_{{\mathcal{G}}_3}={\kappa }_p{\pi }_p\left(1-\eta \rho \right)\left(z_3+z_4\right)+{\mu }_h+{\sigma }_1+{\gamma }_c+{\rho }_1+v_1$

$Z_{{\mathcal{G}}_4}={\gamma }_c+{\omega }_1+{\mu }_h+{\sigma }_1$

$Z_{{\mathcal{G}}_5}={\mu }_h+{\theta }_c$

$Z_{{\mathcal{G}}_6}=\left[{\mu }_h+\left(1-{\chi }_p\right)\left(1-\eta \rho \right)\left(z_3+z_4\right)\right]+\left(1-\eta \rho \right)\left(z_1+z_2\right)$

$Z_{{\mathcal{G}}_7}={\kappa }_c\left(1-\eta \rho \right)\left(z_1+z_2\right)+{\mu }_h+{\sigma }_2+{\gamma }_p$

$Z_{{\mathcal{G}}_8}={\mu }_h+{\theta }_p$

$Z_{{\mathcal{G}}_9}={\mu }_h+{\sigma }_1+{\sigma }_2+{\gamma }_{cd}+{\gamma }_{pd}+v_2+{\rho }_2$

$Z_{{\mathcal{G}}_{10}}={\mu }_h+{\sigma }_1+{\sigma }_2+{\omega }_2+{\gamma }_{pd}$

Thus the argument in the above equations ${\mathcal{G}}_z\mathrm{,}z=\mathrm{1,2,3},4$ satisfies the Lipschitz condition in its second argument with the Lipschitz constant

$L_{{\mathcal{G}}_k}<{\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^{-1}$ □

We apply the Neumann series on the model 21 to test its convergence. Consider the folowing Neumann series,

$S_{hn}\left(t\right)=S_h\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_1\left(t\mathrm{,}{S}_{hn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_1\left(ϵ\mathrm{,}{S}_{hn-1}\left(ϵ\right)\right)\text{d}ϵ$ (25)

$V_{cn}\left(t\right)=V_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_2\left(t\mathrm{,}{V}_{cn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_2\left(ϵ\mathrm{,}{V}_{cn-1}\left(ϵ\right)\right)\text{d}ϵ$ (26)

$I_{cn}\left(t\right)=I_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_3\left(t\mathrm{,}{I}_{cn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_3\left(ϵ\mathrm{,}{I}_{cn-1}\left(ϵ\right)\right)\text{d}ϵ$ (27)

$Q_{cn}\left(t\right)=Q_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_4\left(t\mathrm{,}{Q}_{cn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_4\left(ϵ\mathrm{,}{Q}_{cn-1}\left(ϵ\right)\right)\text{d}ϵ$ (28)

$R_{cn}\left(t\right)=R_c\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_5\left(t\mathrm{,}{R}_{cn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_5\left(ϵ\mathrm{,}{S}_{hn-1}\left(ϵ\right)\right)\text{d}ϵ$ (29)

$V_{pn}\left(t\right)=V_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_6\left(t\mathrm{,}{V}_{pn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_6\left(ϵ\mathrm{,}{V}_{pn-1}\left(ϵ\right)\right)\text{d}ϵ$ (30)

$I_{pn}\left(t\right)=I_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_7\left(t\mathrm{,}{I}_{pn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_7\left(ϵ\mathrm{,}{I}_{pn-1}\left(ϵ\right)\right)\text{d}ϵ$ (31)

$R_{pn}\left(t\right)=R_p\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_8\left(t\mathrm{,}{R}_{pn-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_7\left(ϵ\mathrm{,}{R}_{pn-1}\left(ϵ\right)\right)\text{d}ϵ$ (32)

$I_{cpn}\left(t\right)=I_{cp}\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_9\left(t\mathrm{,}{S}_{cp-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_9\left(ϵ\mathrm{,}{S}_{cpn-1}\left(ϵ\right)\right)\text{d}ϵ$ (33)

$H_{cpn}\left(t\right)=H_{cp}\left(0\right)+\frac{1-\vartheta }{ℬ\left(\vartheta \right)}{\mathcal{G}}_{10}\left(t\mathrm{,}{H}_{cp-1}\left(t\right)\right)+\frac{\vartheta }{ℬ\vartheta \Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}{\mathcal{G}}_{10}\left(ϵ\mathrm{,}{H}_{cpn-1}\left(ϵ\right)\right)\text{d}ϵ$ (34)

It follows that the Neumann series is convergent by considering

$\begin{array}{l}\Vert S_{hn+1}\left(t\right)-S_{hn}\left(t\right)\Vert \\ \le \frac{1-\vartheta }{ℬ}\Vert {\mathcal{G}}_1\left(t\mathrm{,}{S}_{hn}\left(t\right)\right)-{\mathcal{G}}_1\left(t\mathrm{,}{S}_{hn-1}\left(t\right)\right)\Vert \\ -\frac{\vartheta }{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}\Vert {\mathcal{G}}_1\left(ϵ\mathrm{,}{S}_{hn}\right)-{\mathcal{G}}_1\left(ϵ\mathrm{,}{S}_{hn-1}\left(ϵ\right)\right)\Vert \text{d}ϵ\end{array}$ (35)

$\le \frac{\left(1-\vartheta \right)L_{{\mathcal{G}}_1}}{ℬ\left(\vartheta \right)}\Vert S_{hn}\left(t\right)-S_{hn-1}\left(t\right)\Vert +\frac{\vartheta L_{{\mathcal{G}}_1}}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{\vartheta -1}\text{d}ϵ\Vert S_{hn}\left(t\right)-S_{hn-1}\left(t\right)\Vert$ (36)

$\le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_1}{\Vert S_{hn}\left(t\right)-S_{hn-1}\left(t\right)\Vert }_{\infty }$ (37)

Applying similar process yields

$\Vert V_{cn+1}\left(t\right)-V_{cn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_2}{\Vert V_{cn}\left(t\right)-V_{cn-1}\left(t\right)\Vert }_{\infty }$ (38)

$\Vert I_{cn+1}\left(t\right)-I_{cn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_3}{\Vert I_{cn}\left(t\right)-I_{cn-1}\left(t\right)\Vert }_{\infty }$ (39)

$\Vert Q_{cn+1}\left(t\right)-Q_{cn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_4}{\Vert Q_{cn}\left(t\right)-Q_{cn-1}\left(t\right)\Vert }_{\infty }$ (40)

$\Vert R_{cn+1}\left(t\right)-R_{cn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_5}{\Vert R_{cn}\left(t\right)-R_{cn-1}\left(t\right)\Vert }_{\infty }$ (41)

$\Vert V_{pn+1}\left(t\right)-V_{pn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_6}{\Vert V_{pn}\left(t\right)-V_{pn-1}\left(t\right)\Vert }_{\infty }$ (42)

$\Vert I_{pn+1}\left(t\right)-I_{pn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_7}{\Vert I_{pn}\left(t\right)-I_{pn-1}\left(t\right)\Vert }_{\infty }$ (43)

$\Vert R_{pn+1}\left(t\right)-R_{pn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_8}{\Vert R_{pn}\left(t\right)-R_{pn-1}\left(t\right)\Vert }_{\infty }$ (44)

$\Vert I_{cpn+1}\left(t\right)-I_{cpn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_9}|{\Vert I_{cpn}\left(t\right)-I_{cpn-1}\left(t\right)\Vert }_{\infty }$ (45)

$\Vert H_{cpn+1}\left(t\right)-H_{cpn}\left(t\right)\Vert \le \left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)L_{{\mathcal{G}}_{10}}{\Vert H_{cpn}\left(t\right)-H_{cpn-1}\left(t\right)\Vert }_{\infty }$ (46)

It follows that

$\Vert S_{hn+1}\left(t\right)-S_{hn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_1^n}{\Vert S_{h1}\left(t\right)-S_{h0}\left(t\right)\Vert }_{\infty }$ (47)

$\Vert V_{cn+1}\left(t\right)-V_{cn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_2^n}{\Vert V_{c1}\left(t\right)-V_{c0}\left(t\right)\Vert }_{\infty }$

$\Vert I_{cn+1}\left(t\right)-I_{cn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_3^n}{\Vert I_{c1}\left(t\right)-I_{c0}\left(t\right)\Vert }_{\infty }$

$\Vert Q_{cn+1}\left(t\right)-Q_{cn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_4^n}{\Vert Q_{c1}\left(t\right)-Q_{c0}\left(t\right)\Vert }_{\infty }$

$\Vert R_{cn+1}\left(t\right)-R_{cn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_5^n}{\Vert R_{c1}\left(t\right)-R_{c0}\left(t\right)\Vert }_{\infty }$

$\Vert V_{pn+1}\left(t\right)-V_{pn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_6^n}{\Vert V_{p1}\left(t\right)-V_{p0}\left(t\right)\Vert }_{\infty }$

$\Vert I_{pn+1}\left(t\right)-I_{pn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_7^n}{\Vert I_{p1}\left(t\right)-I_{p0}\left(t\right)\Vert }_{\infty }$

$\Vert R_{pn+1}\left(t\right)-R_{pn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_8^n}{\Vert R_{p1}\left(t\right)-R_{p0}\left(t\right)\Vert }_{\infty }$

$\Vert I_{cpn+1}\left(t\right)-I_{cpn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_9^n}{\Vert I_{cp1}\left(t\right)-I_{cp0}\left(t\right)\Vert }_{\infty }$

$\Vert H_{cpn+1}\left(t\right)-H_{cpn}\left(t\right)\Vert \le {\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^n{L}_{{\mathcal{G}}_{10}^n}{\Vert H_{cp1}\left(t\right)-H_{cp0}\left(t\right)\Vert }_{\infty }$

Since

$L_{{\mathcal{G}}_k}<{\left(\frac{1-\vartheta }{ℬ\left(\vartheta \right)}+\frac{T_{\max }^{\vartheta }}{ℬ\left(\vartheta \right)\Gamma \left(\vartheta \right)}\right)}^{-1}$ (48)

Evidently, the right hand sides of inequalities in (48) tends to zero as $n\to \infty$ uniformly on $\left[\mathrm{0,}T\right]$ . By taking limits in (47), we see that $S_h\mathrm{,}{V}_c\mathrm{,}{I}_c\mathrm{,}{Q}_c\mathrm{,}{R}_c\mathrm{,}{V}_p\mathrm{,}{I}_p\mathrm{,}{R}_p\mathrm{,}{I}_{cp}\mathrm{,}{H}_{cp}$ satisfies the original integral, thereby proving the existence and the continuous functions.

4.4. Disease Free Equilibrium and Basic Reproduction Number

The model (7) has a DFE, obtained by setting the right hand sides of the equations in model (7) to zero, given by

${ℰ}_0=\left(S_h^{\mathrm{<mo>*</mo>}}\mathrm{,}{V}_c^{\mathrm{<mo>*</mo>}}\mathrm{,}{I}_c^{\mathrm{<mo>*</mo>}}\mathrm{,}{Q}_c^{\mathrm{<mo>*</mo>}}\mathrm{,}{R}_c^{\mathrm{<mo>*</mo>}}\mathrm{,}{V}_p\mathrm{,}{I}_p\mathrm{,}{R}_p\mathrm{,}{I}_{cp}\mathrm{,}{H}_{cp}\right)$ (49)

$=\left(\frac{\Theta }{{\mu }_h+{\omega }_{vp}+{\omega }_{vc}}\mathrm{,}\frac{{\omega }_{vc}{S}_h^{\mathrm{<mo>*</mo>}}}{{\mu }_h}\mathrm{,0,0,0,}\frac{{\omega }_{vp}{S}_h^{\mathrm{<mo>*</mo>}}}{{\mu }_h}\mathrm{,0,0,0,0}\right)$ (50)

The local stability of the disease free equilibrium, E₀ can be established using the next generation operator method [16] on system (7). Using the notation in lemma (3.2.1).

We first obtain $ℱ$ given as

$\begin{array}{l}{}_0^{abc}{D}_t^{\vartheta }{I}_c=\left[{\pi }_c\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]\left(S_h+V_p\right)+\left[{\pi }_c\left(1-{\chi }_c\right)\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]V_c\\ {}_0^{abc}{D}_t^{\vartheta }{Q}_c=0\\ {}_0^{abc}{D}_t^{\vartheta }{I}_p=\left[{\pi }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]\left(S_h+V_c\right)+\left[{\pi }_p\left(1-{\chi }_p\right)\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)\right]V_p\\ {}_0^{abc}{D}_t^{\vartheta }{I}_{cp}={\pi }_p{\kappa }_p\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)I_c+{\pi }_c{\kappa }_c\left(1-\eta \rho \right)\left(I_p+I_{cp}\right)I_p\\ {}_0^{abc}{D}_t^{\vartheta }{H}_{cp}=0\end{array}$

and $\mathcal{V}$ given as

$\begin{array}{l}{}_0^{abc}{D}_t^{\vartheta }{I}_c=\left[{\mu }_h+{\sigma }_1+{\gamma }_c+{\rho }_1+v_1\right]I_c+{\gamma }_{pd}{I}_{cp}+{\gamma }_{pd}{H}_{cp}\\ {}_0^{abc}{D}_t^{\vartheta }{Q}_c=v_1{I}_c-\left[{\gamma }_c+{\omega }_1+{\mu }_h+{\sigma }_1\right]Q_c\\ {}_0^{abc}{D}_t^{\vartheta }{I}_p=\left[{\mu }_h+{\sigma }_2+{\gamma }_p\right]I_p+\left({\gamma }_{cd}+{\rho }_2\right)I_{cp}+{\omega }_2{H}_{cp}\\ {}_0^{abc}{D}_t^{\vartheta }{I}_{cp}=\left({\mu }_h+{\sigma }_1+{\sigma }_2+{\gamma }_{cd}+{\gamma }_{pd}+v_2+{\rho }_2\right)I_{cp}\\ {}_0^{abc}{D}_t^{\vartheta }{H}_{cp}=v_2{I}_{cp}-\left({\mu }_h+{\sigma }_1+{\sigma }_2+{\omega }_2+{\gamma }_{pd}\right)H_{cp}\end{array}$

the matrix F and V are respectively given as

$F=\left[\begin{array}{ccccc}\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right){\pi }_c+{\pi }_c\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_c^{\mathrm{<mo>*</mo>}}& 0& 0& {\pi }_c\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right)+{\pi }_c\left(1-{\chi }_c\right)\left(1-\eta \rho \right)V_c^{\mathrm{<mo>*</mo>}}& 0\\ 0& 0& 0& 0& 0\\ 0& 0& {\pi }_p\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right)+{\pi }_p\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_p^{\mathrm{<mo>*</mo>}}& {\pi }_p\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_c^{\mathrm{<mo>*</mo>}}\right)+{\pi }_p\left(1-{\chi }_p\right)\left(1-\eta \rho \right)V_p^{\mathrm{<mo>*</mo>}}& 0\\ {\kappa }_p{\pi }_p\left(1-\eta \rho \right)\left(I_p^{\mathrm{<mo>*</mo>}}+I_{cp}^{\mathrm{<mo>*</mo>}}\right)+{\kappa }_c{\pi }_c\left(1-\eta \rho \right)I_p^{\mathrm{<mo>*</mo>}}& 0& {\kappa }_p{\pi }_p\left(1-\eta \rho \right)I_c^{\mathrm{<mo>*</mo>}}+{\kappa }_c{\pi }_c\left(1-\eta \rho \right)\left(I_p^{\mathrm{<mo>*</mo>}}+I_{cp}^{\mathrm{<mo>*</mo>}}\right)& {\kappa }_p{\pi }_p\left(1-\eta \rho \right)I_c^{\mathrm{<mo>*</mo>}}+{\kappa }_c{\pi }_c\left(1-\eta \rho \right)I_p^{\mathrm{<mo>*</mo>}}& 0\\ 0& 0& 0& 0& 0\end{array}\right]$ (51)

$V=\left[\begin{array}{ccccc}{\mu }_h+{\sigma }_1+{\gamma }_c+{\rho }_1+v_1& {\gamma }_c+{\omega }_1+{\mu }_h+{\sigma }_1& 0& {\gamma }_{pd}& {\gamma }_{pd}\\ v_1& {\mu }_h+{\sigma }_2+{\gamma }_p& 0& 0& 0\\ 0& 0& {\mu }_h+{\gamma }_p+{\sigma }_2& {\gamma }_{cd}+{\rho }_2& {\omega }_2\\ 0& 0& 0& {\mu }_h+{\sigma }_1+{\sigma }_2+{\gamma }_{cd}+{\gamma }_{pd}+v_2+{\rho }_2& 0\\ 0& 0& 0& v_2& -\left({\mu }_h+{\sigma }_1+{\sigma }_2+{\omega }_2+{\gamma }_{pd}\right)\end{array}\right]$ (52)

$F{V}^{-1}=\left[\begin{array}{c}0\\ 0\\ \frac{\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right){\pi }_c+{\pi }_c\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_c^{\mathrm{<mo>*</mo>}}}{{\mu }_h+{\sigma }_1+{\gamma }_c+{\rho }_1+v_1}\\ \frac{\left(1-\eta \rho \right){\pi }_p\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right)+{\pi }_p\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_p^{\mathrm{<mo>*</mo>}}}{{\mu }_h+{\sigma }_1+{\sigma }_2+{\gamma }_{cd}+{\gamma }_{pd}+v_2+{\rho }_2}\end{array}\right]$ (53)

the basic reproduction number of the model (7), is given by ${ℛ}_0=\max \left\{{ℛ}_{01},{ℛ}_{02}\right\}$ , where ${ℛ}_{01}$ and ${ℛ}_{02}$ are respective reproduction numbers given by

$\begin{array}{l}{ℛ}_{01}=\frac{{\pi }_c\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right)+{\pi }_c\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_c^{\mathrm{<mo>*</mo>}}}{{\mu }_h+{\sigma }_1+{\gamma }_c+{\rho }_1+v_1}\mathrm{,}\\ {ℛ}_{02}=\frac{{\pi }_c\left(1-\eta \rho \right)\left(S_h^{\mathrm{<mo>*</mo>}}+V_p^{\mathrm{<mo>*</mo>}}\right)+{\pi }_c\left(1-\eta \rho \right)\left(1-{\chi }_c\right)V_c^{\mathrm{<mo>*</mo>}}}{{\mu }_h+{\sigma }_1+{\sigma }_2+{\gamma }_{cd}+{\gamma }_{pd}+v_2+{\rho }_2}\mathrm{.}\end{array}$

4.5. Generalized Ulam-Hyers-Rassias Stability

In this section, we establish the generalized Ulam-Hyers-Rassias (UHR) stability for the fractional model 7. This type of stability assesses whether approximate solutions stay close to true solutions under certain conditions. We state the required definition.

Definition 4.2. The model 7 has Ulam-Hyers stability if there exist constant ${\vartheta }_i>0$ , $i\in {ℝ}^{10}$ , satisfying for every ${ϵ}_i>0$ , $i\in {ℝ}^{10}$ . if

$\left|{}_0^{abc}{D}_t^{\vartheta }{S}_h\left(t\right)-Z_1\left(t\mathrm{,}s\right)\right|\le {ϵ}_1$ (54)

$\left|{}_0^{abc}{D}_t^{\vartheta }{V}_c\left(t\right)-Z_2\left(t\mathrm{,}{V}_c\right)\right|\le {ϵ}_2$ (55)

$\left|{}_0^{abc}{D}_t^{\vartheta }{I}_c\left(t\right)-Z_3\left(t\mathrm{,}{I}_c\right)\right|\le {ϵ}_3$ (56)

$\left|{}_0^{abc}{D}_t^{\vartheta }{Q}_c\left(t\right)-Z_4\left(t\mathrm{,}{Q}_c\right)\right|\le {ϵ}_4$ (57)

$\left|{}_0^{abc}{D}_t^{\vartheta }{R}_c\left(t\right)-Z_5\left(t\mathrm{,}{R}_c\right)\right|\le {ϵ}_5$ (58)

$\left|{}_0^{abc}{D}_t^{\vartheta }{V}_p\left(t\right)-Z_6\left(t\mathrm{,}{V}_p\right)\right|\le {ϵ}_6$ (59)

$\left|{}_0^{abc}{D}_t^{\vartheta }{I}_p\left(t\right)-Z_7\left(t\mathrm{,}{I}_p\right)\right|\le {ϵ}_7$ (60)

$\left|{}_0^{abc}{D}_t^{\vartheta }{R}_p\left(t\right)-Z_8\left(t\mathrm{,}{R}_p\right)\right|\le {ϵ}_8$ (61)

$\left|{}_0^{abc}{D}_t^{\vartheta }{I}_{cp}\left(t\right)-Z_9\left(t\mathrm{,}{I}_{cp}\right)\right|\le {ϵ}_9$ (62)

$\left|{}_0^{abc}{D}_t^{\vartheta }{H}_{cp}\left(t\right)-Z_{10}\left(t\mathrm{,}{H}_{cp}\right)\right|\le {ϵ}_{10}$ (63)

and there exist a solution of the COVID-19 and Streptococcus pneumoniae Model 7 given by

$\tilde{S}\left(t\right)\mathrm{,}{\tilde{V}}_c\left(t\right)\mathrm{,}{\tilde{I}}_c\left(t\right)\mathrm{,}{\tilde{Q}}_c\left(t\right)\mathrm{,}{\tilde{R}}_c\left(t\right)\mathrm{,}{\tilde{V}}_p\left(t\right)\mathrm{,}{\tilde{I}}_p\left(t\right)\mathrm{,}{\tilde{R}}_p\left(t\right)\mathrm{,}{\tilde{I}}_{cp}\left(t\right)\mathrm{,<mi>\; </mi>}{H}_{cp}\left(t\right)$ (64)

that satisfying the given model such that

$\begin{array}{l}\Vert S_h-{\tilde{S}}_h\Vert \le {\vartheta }_1{ϵ}_1\mathrm{,}\Vert V_c-{\tilde{V}}_c\Vert \le {\vartheta }_2{ϵ}_2\mathrm{,}\Vert I_c-{\tilde{I}}_c\Vert \le {\vartheta }_3{ϵ}_3\mathrm{,}\Vert Q_c-{\tilde{Q}}_c\Vert \le {\vartheta }_4{ϵ}_4\mathrm{,}\\ \Vert R_c-{\tilde{R}}_c\Vert \le {\vartheta }_5{ϵ}_5,\Vert V_p-{\tilde{V}}_p\Vert \le {\vartheta }_6{ϵ}_6\mathrm{,}\Vert I_p-{\tilde{I}}_p\Vert \le {\vartheta }_7{ϵ}_7\mathrm{,}\Vert R_p-{\tilde{R}}_p\Vert \le {\vartheta }_8{ϵ}_8\mathrm{,}\\ \Vert I_{cp}-{\tilde{I}}_{cp}\Vert \le {\vartheta }_9{ϵ}_9\mathrm{,}\Vert H_{cp}-{\tilde{H}}_{cp}\Vert \le {\vartheta }_1{}_0{ϵ}_1{}_0\end{array}$

Remark 1. Consider that the function $S_h$ is a solution of the first inequality in 4.2 above, if a continuous function $h_1$ exist so that

$\left|h_1\left(t\right)\right|\le {ϵ}_1$

${}_0^{abc}{D}_t^{\vartheta }{S}_h\left(t\right)=Z_1\left(t\mathrm{,}{S}_h\right)+h_1\left(t\right)$ (65)

Similarly, we apply this to the other inequalities by finding $h_i$ for $i\in {ℝ}^1{}^0$

$\left|h_2\left(t\right)\right|\le {ϵ}_2$ (66)

${}_0^{abc}{D}_t^{\vartheta }{V}_c\left(t\right)=Z_2\left(t\mathrm{,}{V}_c\right)+h_2\left(t\right)$ (67)

$\left|h_3\left(t\right)\right|\le {ϵ}_3$ (68)

${}_0^{abc}{D}_t^{\vartheta }{I}_c\left(t\right)=Z_3\left(t\mathrm{,}{I}_c\right)+h_3\left(t\right)$ (69)

$\left|h_4\left(t\right)\right|\le {ϵ}_4$ (70)

${}_0^{abc}{D}_t^{\vartheta }{Q}_c\left(t\right)=Z_4\left(t\mathrm{,}{Q}_c\right)+h_4\left(t\right)$ (71)

$\left|h_5\left(t\right)\right|\le {ϵ}_5$ (72)

${}_0^{abc}{D}_t^{\vartheta }{R}_c\left(t\right)=Z_5\left(t\mathrm{,}{R}_c\right)+h_5\left(t\right)$ (73)

$\left|h_6\left(t\right)\right|\le {ϵ}_6$ (74)

${}_0^{abc}{D}_t^{\vartheta }{V}_p\left(t\right)=Z_6\left(t\mathrm{,}{V}_p\right)+h_6\left(t\right)$ (75)

$\left|h_7\left(t\right)\right|\le {ϵ}_7$ (76)

${}_0^{abc}{D}_t^{\vartheta }{I}_p\left(t\right)=Z_7\left(t\mathrm{,}{I}_p\right)+h_7\left(t\right)$ (77)

$\left|h_8\left(t\right)\right|\le {ϵ}_8$ (78)

${}_0^{abc}{D}_t^{\vartheta }{R}_p\left(t\right)=Z_8\left(t\mathrm{,}{R}_p\right)+h_8\left(t\right)$ (79)

$\left|h_9\left(t\right)\right|\le {ϵ}_9$ (80)

${}_0^{abc}{D}_t^{\vartheta }{I}_{cp}\left(t\right)=Z_9\left(t\mathrm{,}{I}_{cp}\right)+h_9\left(t\right)$ (81)

$\left|h_{10}\left(t\right)\right|\le {ϵ}_{10}$ (82)

${}_0^{abc}{D}_t^{\vartheta }{H}_{cp}\left(t\right)=Z_1{}_0\left(t\mathrm{,}{H}_{cp}\right)+h_1{}_0\left(t\right)$ (83)

Let us assume that the model 7 is Ulam-Hyers Stable.

Proof. To prove the above assumption, let ${ϵ}_1>0$ and the function $S_h$ be arbitrary so that

$\left|{}_0^{abc}{D}_t^{{\vartheta }_1,{\vartheta }_2}{S}_h\left(t\right)-Z_1\left(t\mathrm{,}{S}_h\right)\right|\le {ϵ}_1$ (84)

Referencing remark 1, we have a function $h_1$ with $\left|h_1\left(t\right)\right|<{ϵ}_1$ satisfies

${}_0^{abc}{D}_t^{{\vartheta }_1,{\vartheta }_2}{S}_h\left(t\right)=Z_1\left(t\mathrm{,}{S}_h\right)+h_1\left(t\right)$

Consequently,

$\begin{array}{c}{S}_h\left(t\right)=S_h\left(0\right)+\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1\right)}{\displaystyle {\int }_0^t}{\left(t-ϵ\right)}^{{\vartheta }_1-1}{Z}_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\text{d}ϵ\\ +\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)+{\vartheta }_2-1}{ℬ\left({\vartheta }_1\right)}{Z}_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\\ +\frac{{\vartheta }_1{\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1\right)}{\displaystyle {\int }_0^t}{t}^{{\vartheta }_2-1}{\left(t-ϵ\right)}^{{\vartheta }_1-1}{h}_1\left(t\right)\text{d}ϵ+\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)t^{{\vartheta }_2-1}}{ℬ\left({\vartheta }_1\right)}{h}_1\left(s\right)\end{array}$

Let ${\tilde{S}}_h$ be the solution of the given model, then

$\begin{array}{c}{\tilde{S}}_h\left(t\right)=S_h\left(0\right)+\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1\right)}{\displaystyle {\int }_0^t}{t}^{{\vartheta }_2-1}{\left(t-ϵ\right)}^{{\vartheta }_1-1}{Z}_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\text{d}ϵ\\ +\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)t^{{\vartheta }_2-1}}{ℬ\left({\vartheta }_1\right)}{Z}_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\end{array}$

Hence,

$\begin{array}{c}\left|S_h\left(t\right)-{\tilde{S}}_h\left(t\right)\right|\le \frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1\right)}{\displaystyle {\int }_0^t}{t}^{{\vartheta }_2-1}{\left(t-ϵ\right)}^{{\vartheta }_1-1}\left|Z_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\right|\text{d}ϵ\\ +\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)t^{{\vartheta }_2-1}}{ℬ\left({\vartheta }_1\right)}\left|Z_1\left(t\mathrm{,}{S}_h\left(t\right)\right)-Z_1\left(t\mathrm{,}{\tilde{S}}_h\left(t\right)\right)\right|+\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}\\ +\frac{{\vartheta }_1{\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1\right)}{\displaystyle {\int }_0^t}{t}^{{\vartheta }_2-1}{\left(t-ϵ\right)}^{{\vartheta }_1-1}\left|h_1\left(t\mathrm{,}{S}_h\left(t\right)\right)\right|\text{d}ϵ\\ +\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)t^{{\vartheta }_2-1}}{ℬ\left({\vartheta }_1\right)}\left|h_1\left(t\right)\right|\\ \le \left[\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}+\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)}{ℬ\left({\vartheta }_1\right)}\right]{\psi }_1\left|S_h-{\tilde{S}}_h\right|\\ +\left[\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}+\frac{{\vartheta }_2\left(1-{\vartheta }_1\right)}{ℬ\left({\vartheta }_1\right)}\right]{ϵ}_1\end{array}$

$\Vert S_h-{\tilde{S}}_h\Vert \le \frac{\left[\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}\right]{ϵ}_1}{1-\left[\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}\right]{\psi }_1}$

Then

$\Vert S-h-{\tilde{S}}_h\Vert \le {\vartheta }_1{ϵ}_1$ (85)

Therefore,

${\vartheta }_1=\frac{\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}}{1-\left[\frac{1-{\vartheta }_1{\vartheta }_2}{ℬ\left({\vartheta }_1\right)}+\frac{{\vartheta }_1{\vartheta }_2\Gamma {\vartheta }_2}{ℬ{\vartheta }_1\Gamma \left({\vartheta }_1+{\vartheta }_2\right)}\right]{\psi }_1}$ (86)

Now applying similar approach to the other functions, we obtain

$\Vert V_c-{\tilde{V}}_c\Vert \le {\vartheta }_2{ϵ}_2$ (87)

$\Vert I_c-{\tilde{I}}_c\Vert \le {\vartheta }_3{ϵ}_3$ (88)

$\Vert Q_c-{\tilde{Q}}_c\Vert \le {\vartheta }_4{ϵ}_4$ (89)

$\Vert R_c-{\tilde{R}}_c\Vert \le {\vartheta }_5{ϵ}_5$ (90)

$\Vert V_p-{\tilde{V}}_p\Vert \le {\vartheta }_6{ϵ}_6$ (91)

$\Vert I_p-{\tilde{I}}_p\Vert \le {\vartheta }_7{ϵ}_7$ (92)

$\Vert R_p-{\tilde{R}}_p\Vert \le {\vartheta }_8{ϵ}_8$ (93)

$\Vert I_{cp}-{\tilde{I}}_{cp}\Vert \le {\vartheta }_9{ϵ}_9$ (94)

$\Vert H_{cp}-{\tilde{H}}_{cp}\Vert \le {\vartheta }_1{}_0{ϵ}_1{}_0$ (95)

Hence, we conclude that the function model 7 is Ulam-Hyers stable

4.6. Numerical Scheme: Fractional Adams-Bashforth-Moulton Method

We adopted the well-known fractional Predictor-Corrector method proposed by Diethelm et al in solving the numerical scheme of the fractional COVID and pneumonia model. The method combines the Predict-Evaluate-Correctional Euler or one-step Adams-Bashforth method for prediction with the fractional one-step Adams-Moulton method for correction, collectively known as the fractional Adams-Bashforth-Moulton method (ABMM). We consider the initial value problem endowed with the Caputo fractional derivative:

$\{\begin{array}{l}{}^C{D}_t^{\phi }\Psi =\Phi \left(t,\Psi \left(t\right)\right),0<\phi \le 1,\mathrm{<mi>\; </mi>}t\in \mathcal{J}\in \left[0,T\right],\\ \Psi \left(0\right)={\Psi }_0\ge 0.\end{array}$

This problem is equivalent to the following Volterra integral equation:

$\Psi \left(t\right)={\Psi }_0+\frac{1}{\Gamma \left(\phi \right)}{\displaystyle {\int }_0^t}{\left(t-\xi \right)}^{\phi -1}\Phi \left(\xi ,\Psi \left(\xi \right)\right)\text{d}\xi ,0<\phi \le 1.$

To approximate the integral, we divide the time T into N equal parts with a step-size of $h=T/N$ . The points in the grid are defined by the nodes $t_j=jh$ for $j=\mathrm{0,1,2,}\cdots \mathrm{,}k+1$ . We then replace the integral part of the equation with:

$\Psi \left(t\right)={\Psi }_0+\frac{h^{\phi }}{\Gamma \left(\phi +1\right)}{\displaystyle \sum _{j=0}^{k+1}}{\vartheta }_{j,k+1}\Phi \left(t_j\right),$

where ${\vartheta }_{j\mathrm{,}k+1}$ are weights defined as follows:

${\vartheta }_{j,k+1}=\{\begin{array}{ll}{k}^{\phi +1}-\left(k-\phi \right){\left(k+1\right)}^{\phi },\hfill & j=0,\hfill \\ {\left(k-j+2\right)}^{\phi +1}+{\left(k-j\right)}^{\phi +1}-2\left(k-j+1\right),\hfill & 1\le j\le k,\hfill \\ 1,\hfill & j=k+1.\hfill \end{array}$

The fractional Adams-Bashforth predictor method is obtained by substituting the above approximation into the equation:

$S\left(t_{t+1}\right)=S\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,S\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,S^p\left(t_{k+1}\right)\right)\right)$

$V_c\left(t_{t+1}\right)=V_c\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,V_c\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,V_c^p\left(t_{k+1}\right)\right)\right)$

$I_c\left(t_{t+1}\right)=I_c\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,I_c\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,I_c^p\left(t_{k+1}\right)\right)\right)$

$Q_c\left(t_{t+1}\right)=Q_c\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathbb{G}}_1\left(t_j,Q_c\left(t_j\right)\right)+{\mathbb{G}}_1\left(t_k+1,Q_c^p\left(t_{k+1}\right)\right)\right)$

$R_c\left(t_{t+1}\right)=R_c\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,R_c\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,R_c^p\left(t_{k+1}\right)\right)\right)$

$V_p\left(t_{t+1}\right)=V_p\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,V_p\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,V_p^p\left(t_{k+1}\right)\right)\right)$

$I_p\left(t_{t+1}\right)=I_p\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,I_p\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,I_p^p\left(t_{k+1}\right)\right)\right)$

$R_p\left(t_{t+1}\right)=R_p\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,R_p\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,R_p^p\left(t_{k+1}\right)\right)\right)$

$I_{cp}\left(t_{t+1}\right)=I_{cp}\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,I_{cp}\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,I_{cp}^p\left(t_{k+1}\right)\right)\right)$

$H_{cp}\left(t_{t+1}\right)=H_{cp}\left(0\right)+\frac{h^{\phi }}{\Gamma \left(\phi +2\right)}\left({\displaystyle \sum }_{j=0}^{k+1}{\vartheta }_{j,k+1}{\mathcal{G}}_1\left(t_j,H_{cp}\left(t_j\right)\right)+{\mathcal{G}}_1\left(t_k+1,H_{cp}^p\left(t_{k+1}\right)\right)\right)$

where ${\iota }_{j\mathrm{,}k+1}=\frac{h^{\phi }}{\phi }\left({\left(k-j+1\right)}^{\phi }-{\left(k-j\right)}^{\phi }\right)$ are the weights for the predictor values.

The predictor-corrector method’s computational complexity, is $\mathcal{O}\left(h^{-2}\right)$ , compared to the $\mathcal{O}\left(h^{-1}\right)$ complexity of integer-order initial value problems. This increased complexity arises from the non-local nature of fractional-order derivatives. However, advanced techniques, such as the nested memory concept, offer improved computational efficiency with a complexity of $\mathcal{O}\left(h^{-1}\log \left(h^{-1}\right)\right)$ while maintaining accuracy.

4.7. Sensitivity Analysis

Sensitivity analysis is carried out in this section to analyse the influence of the different parameters involved in the reproduction numbers of model 7. We used PRCC techniques separately for the reproduction number to show the role of the parameters in the reproduction number. It can be observed from Figure 2, that the parameters that will have a high impact on $R_{01}$ are ${\beta }_c\mathrm{,<mi>\; </mi>}\Theta$ , and ${\chi }_p$ . These parameters are in direct proportion with $R_{01}$ , as any factor that increases them also raises the value of $R_{01}$ . Hence efforts target on these parameters should be ones that reduce their values. On the other hand ${\omega }_{vc}\mathrm{,<mi>\; </mi>}{\omega }_{vp}$ and ${\nu }_1$ reduces

Figure 2. Sensitivity with respect to $R_{01}$ .

the value of $R_{01}$ when we increase them. The other parameters associated with this reproduction number have minimal effect on it. As shown in Figure 3, using $R_{02}$ as the response function, $\Theta \mathrm{,}$ and ${\omega }_{vc}$ are in direct variation with $R_{02}$ , while ${\omega }_{vc}\mathrm{,<mi>\; </mi>}{\chi }_p$ and ${\sigma }_2$ play a reverse role.

Figure 3. Sensitivity with respect to $R_{02}$ .

Simulations

Applying the numerical scheme given above, we now present the plots. Figures 4(a)-(d), shows the plots of the infectious classes over time at different fractional

Figure 4. Comparison of the infectious classes at different orders.

orders of 0.95, 0.90, 0.85, 0.80. It shows that as we increase the order at DFE, we have the dynamics of the infectious classes of the model approaching zero faster, suggesting that at if we want to achieve a desirable result of reducing the disease burden, we should analyze these models at the best possible high order.

The plots in Figure 5(a) & Figure 5(b) and Figure 6(a)-6(b) below confirm the convergence of the disease classes to zero. This result confirms with what is expected at DFE, when the reproduction number is less than unity and from the definition of the reproduction number, the disease will eventually die out. The figures also confirmed the fact that the differences observed when the simulations are executed through a fractional derivative are a result of the memory effect which is absent in the classical integer order operator. The memory effect which considers past records and events of the diseases will play a great influence on the dynamics and accurate future prediction of the disease for effective control. In other words, this establishes that the result obtained through the AB derivative is more accurate than those obtained from integer order operator.

Figure 5. 3-dimensional plot of the infectious classes under different initial conditions.

Figure 6. Convergence of the infectious classes.

5. Conclusion

In this study, we presented and investigated the transmission dynamics of COVID-19 and Streptococcus pneumoniae co-infection using Atangana-Baleanu fractional derivative operator to epidemiologically understand the two diseases while reducing their spread. The model was established to be mathematically and biologically well-posed by showing that all solutions to the model are positive and remain bounded within a region. The existence and uniqueness analysis were established using Banach fixed point theorem to ensure validity of the model. Furthermore, we examine the stability of the model by applying the criteria proposed by Ulam-Hyers-Rassias Stability. The numerical scheme of the fractional model was solved using the fractional Predictor-Corrector method which combines the Predict-Evaluate-Correctional-Euler or one-step Adams-Bashforth method for prediction with the fractional one-step Adams-Moulton method.

The impact of the fractional derivative on the proposed model was presented using simulations. The simulation results as presented in the figures clearly reveal that the dynamics of both diseases in each epidemiological category over time were significantly impacted by the use of fractional derivatives. The variations observed in the simulations with fractional derivatives are attributed to the memory effect, which is absent when using the traditional integer-order operator. The existence of this memory effect, which accounts for past events, has a profound influence on the future dynamics of diseases, making them more amenable to control. Consequently, this highlights the superiority of the results obtained through the use of the AB derivative compared to those achieved with the integer-order operator.

6. Recommendations

From these work findings, we recommend the following:

1) Effective vaccination against COVID-19, and Streptococcus pneumoniae, and administering of anti-bacteria treatment will decrease the basic reproductive number of the co-infection to barest minimum.

2) Streptococcus pneumoniae clinical diagnosis should be made compulsory for every COVID-19 patient for early detection and treatment for co-infected individual to minimize the severity and mortality rate of the co-infection

Contribution to Knowledge

The study formulated an in-dept fractional mathematical model for the co-infection of COVID-19 and Streptococcus pneumoniae with saturated incidence force of infection, and integrated the isolation, and hospitalized compartments. The fractional order model integrated memory effect, and hence provided accurate predictions for the control of the diseases.

Data Availability

All data regarding the research work is clearly mentioned in the research work.

Conflicts of Interest

The authors declare no conflicts of interest.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	World Health Organization (WHO) (2023) Website Publication. https://covid19.who.int/region/afro/country/ng://covid19.ncdc.gov.ng/
[2]	Our World in Data 2023. https://ourworldindata.org/covid-vaccinations?country=OWID_WRL
[3]	WorldoMeter Data 2023. https://www.worldometers.info/coronavirus/
[4]	Mitchell, T.J. and Dalziel, C.E. (2014) The Biology of Pneumolysin. In: Anderluh, G. and Gilbert, R., Eds., MACPF/CDC Proteins-Agents of Defence, Attack and Invasion, Springer, 145-160. https://doi.org/10.1007/978-94-017-8881-6_8
[5]	Nigeria Health Watch (2023) Going Beyond Global Commitments to Local Action for Zero Childhood Pneumonia Deaths in Nigeria? https://nigeriahealthwatch.medium.com/going-beyond-global-commitments-to-local-action-for-zero-childhood-pneumonia-deaths-in-nigeria-cc32e8e59280#:~:text=For%20three%20decades%2C%20pneumonia%20has,rate%20of%20386.15%20per%20100%2C000
[6]	Centers for Disease Control and Prevention (2017) “Streptococcus pneumoniae” Published on Official Website. https://www.cdc.gov/pneumococcal/clinicians/streptococcus-pneumoniae.html
[7]	Omame, A., Abbas, M. and Onyenegecha, C.P. (2022) A Fractional Order Model for the Co-Interaction of COVID-19 and Hepatitis B Virus. Results in Physics, 37, Article 105498. https://doi.org/10.1016/j.rinp.2022.105498
[8]	Baleanu, D., Mohammadi, H. and Rezapour, S. (2020) A Fractional Differential Equation Model for the COVID-19 Transmission by Using the Caputo-Fabrizio Derivative. Advances in Difference Equations, 2020, Article No. 299. https://doi.org/10.1186/s13662-020-02762-2
[9]	Omame, A., Isah, M.E., Abbas, M., Abdel-Aty, A. and Onyenegecha, C.P. (2022) A Fractional Order Model for Dual Variants of COVID-19 and HIV Co-Infection via Atangana-Baleanu Derivative. Alexandria Engineering Journal, 61, 9715-9731. https://doi.org/10.1016/j.aej.2022.03.013
[10]	Aba Oud, M.A., Ali, A., Alrabaiah, H., Ullah, S., Khan, M.A. and Islam, S. (2021) A Fractional Order Mathematical Model for COVID-19 Dynamics with Quarantine, Isolation, and Environmental Viral Load. Advances in Difference Equations, 2021, Article No. 106. https://doi.org/10.1186/s13662-021-03265-4
[11]	Peter, O.J., Yusuf, A., Oshinubi, K., Oguntolu, F.A., Lawal, J.O., Abioye, A.I., et al. (2021) Fractional Order of Pneumococcal Pneumonia Infection Model with Caputo Fabrizio Operator. Results in Physics, 29, Article 104581. https://doi.org/10.1016/j.rinp.2021.104581
[12]	Omame, A., Abbas, M. and Onyenegecha, C.P. (2021) A Fractional-Order Model for COVID-19 and Tuberculosis Co-Infection Using Atangana-Baleanu Derivative. Chaos, Solitons & Fractals, 153, Article 111486. https://doi.org/10.1016/j.chaos.2021.111486
[13]	Aga, B.Z., Keno, T.D., Terfasa, D.E. and Berhe, H.W. (2024) Pneumonia and COVID-19 Co-Infection Modeling with Optimal Control Analysis. Frontiers in Applied Mathematics and Statistics, 9, Article 1286914. https://doi.org/10.3389/fams.2023.1286914
[14]	Atangana, A. and Baleanu, D. (2016) New Fractional Derivatives with Nonlocal and Non-Singular Kernel: Theory and Application to Heat Transfer Model. Thermal Science, 20, 763-769. https://doi.org/10.2298/tsci160111018a
[15]	Atangana, A. and Secer, A. (2013) A Note on Fractional Order Derivatives and Table of Fractional Derivatives of Some Special Functions. Abstract and Applied Analysis, 2013, Article 279681. https://doi.org/10.1155/2013/279681
[16]	van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/s0025-5564(02)00108-6