A dynamical analysis and numerical simulation of COVID-19 and HIV/AIDS co-infection with intervention strategies

Abstract

HIV/AIDS-COVID-19 co-infection is a major public health concern especially in developing countries of the world. This paper presents HIV/AIDS-COVID-19 co-infection to investigate the impact of interventions on its transmission using ordinary differential equation. In the analysis of the model, the solutions are shown to be non-negative and bounded, using next-generation matrix approach the basic reproduction numbers are computed, sufficient conditions for stabilities of equilibrium points are established. The sensitivity analysis showed that transmission rates are the most sensitive parameters that have direct impact on the basic reproduction numbers and protection and treatment rates are more sensitive and have indirect impact to the basic reproduction numbers. Numerical simulations shown that some parameter effects on the transmission of single infections as well as co-infection, and applying the protection rates and treatment rates have effective roles to minimize and also to eradicate the HIV/AIDS-COVID-19 co-infection spreading in the community.

KEYWORDS:

• COVID-19
• co-infection
• protection
• treatment
• stability analysis

1. Introduction

Infectious diseases are clinically verified illnesses caused by pathogenic microbial agents’ such as viruses, bacteria, fungi, and parasites [14]. COVID-19 and HIV/AIDS are most common infectious diseases caused by viruses [12,27].

The human immunodeficiency virus (HIV) discovered in 1981, has been among the leading life-affecting viruses expanding throughout nations in the world. Acquired immunodeficiency syndrome abbreviated as AIDS is one of the most community-affecting epidemics in history caused by HIV and it has been a global pandemic since 1981 [28,29,12,22,18]. HIV/AIDs remain a significant global health problem that affects almost 70 million people worldwide and has been a significant cause of morbidity and mortality [12,2,21]. It mainly attacks human white blood cells and is transmitted through sexual intercourse, needle sharing, and direct contact with blood or other body fluids containing the virus and from mother to child during childbirth [19,29,22].

COVID-19 a communicable disease caused by the novel coronavirus (SARS-CoV-2) was discovered for the first time in the middle of December 2019 in China and it has been declared an extremely contagious acute respiratory disease [31,24,26,7,20,22,27]. Since it has rapidly spread throughout nations in the world and causing major health concern the World Health Organization (WHO) declared it as a global pandemic on 11 March 2020 [32,3,20]. It has a direct transmission from human being to human being with direct touch of contaminated materials and indirect transmission with inhalation of respiratory air droplets by coughing and sneezing [23,25,9,17]. The pandemic has an extra health system burden and a great negative impact on nation’s population and economies throughout the whole world [16,10,26]. Recently, WHO approved vaccination, quarantine, isolation, using face masks, washing hands with alcohol, and social distancing are possible prevention and control strategies [32,23,33,10].

Different researchers have investigated that COVID-19 infection could be high in people living with other infections like TB, HIV and cholera who have compromised immunity [24,11,26,1,22,4]. HIV/AIDS infection is the main cause of individuals’ weak immunity and individuals living with HIV/AIDS are highly affected by opportunistic infections like TB, pneumonia and COVID-19 [2,28]. Even though there is no full evidence to justify the observation conclusively, individuals living with HIV/AIDS were more likely to have a positive medical result, and were almost 32% and 82% more likely to get hospital admitting and require ventilation, respectively, because of COVID-19 infection than those had negative HIV/AIDs result [1,22]. Different scholars have formulated and analysed mathematical models to investigate the transmission dynamics of infectious diseases using ordinary differential equations approach like [31,19,6,11,32,3,23,25,16,28,29,10,17,26,7] and using fractional order derivative approach like [1].

Mathematical modelling approaches have been crucial to provide basic frameworks in order to understand the transmission dynamics of infectious diseases [26]. In the development of this paper, we have reviewed research papers that have been done on the transmission dynamics of different infectious diseases. Tchoumi et al. [26] proposed and investigated the co-dynamics of malaria and COVID-19 co-dynamics: with optimal control strategies. The numerical simulation results verifies the theoretical optimal control analysis and illustrates that using malaria and COVID-19 protection measures concurrently can help mitigate there transmission compared with applying single infections protection measures. Teklu, Shewafera Wondimagegnhu, and Koya Purnachandra Rao, 2022 [29] constructed and examined HIV/AIDS and pneumonia co-infection model with control measures such as pneumonia vaccination and treatments of pneumonia and HIV/AIDS infections. H.F. Huo and R. Chen [12] formulated a stage structure HIV/AIDS model to study the transmission of HIV/AIDS with treatment. The model analyses show that ART at the asymptomatic initial stages of the HIV incidence or before-AIDS stage is effective to decrease its spreading rate. Omondi et al. (2018) [21], analysed a sex-structured community infection model and discuss male and female HIV infection trends with heterosexual activities. The analysis of the study deduced that ART treatment has a considerable impact on the HIV/AIDS spreading rate. Hezam et al. [11], formulated a mathematical model for cholera and COVID-19 co-infection which describes the transmission dynamics of COVID-19 and cholera in Yemen. The model analysis, examined four controlling measures such as social distancing, lockdown, the number of test kits to control the COVID-19 outbreak, and the number of susceptible individuals who can get CWTs for water purification. Anwar et al. [33], constructed a mathematical model on COVID-19 with the isolation controlling measure on the COVID-19-infected individuals throughout the community. Ahmed et al. [1] formulated and analysed HIV and COVID-19 co-infection model with ABC-fractional operator approach to investigate an epidemic prediction of a combined HIV-COVID-19 co-infection model. Numerical simulations were carried out to justify that the disease will stabilize at a later stage when enough protection strategies are taken. Ringa et al. [22] developed and examined a mathematical model on HIV and COVID-19 co-infection dynamics with optimal control strategies. The results suggested that HIV prevention strategies can significantly reduce the burden of co-infections with COVID-19, while effective treatment of COVID-19 could reduce co-infections with opportunistic infections such as HIV/AIDS.

Researches show that HIV/AIDS and COVID-19 co-infection is a major public health concern especially in developing nations of the world. The main purpose of this paper is to investigate the impacts of COVID-19 protection, COVID-19 treatment, HIV prevention and HIV treatment prevention and controlling strategies simultaneously on the transmission dynamics of HIV/AIDS and COVID-19 co-infection in the community. Some researchers invested much effort in studying HIV/AIDS and COVID-19 co-infection, none of them has considered COVID-19 protection, COVID-19 treatment, HIV protection, and HIV treatment as prevention and control strategies simultaneously in a single model formulation. Based on the findings of different literatures we have realized the gaps and highly motivated to undertake this study and fill the gap we have mentioned above. The remaining part of the paper is organized in the following sequence; the model is constructed in section 2 and is analysed in section 3. Sensitivity analysis and numerical simulation, and discussions and conclusions of the study are carried out in sections 4, and 5, respectively.

2. Model formulation

Motivated by observation of different scholars epidemiological models in a real world situations have played a fundamental role in the better understanding of the transmission dynamics and controlling and preventive strategies against infectious diseases, in this paper, we have proposed a mathematical model with system of ordinary differential equations. The total human population at a given time $t$ which is denoted by $N(t)$ is classified into eight mutually exclusive classes of susceptible individuals denoted by$S(t)$, COVID-9 protected individuals denoted by$P_1(t)$, HIV protected individuals denoted by $P_2(t)$, COVID-19 mono-infected individuals denoted by $I_C(t)$, HIV mono-infected individuals denoted by $I_H(t)$, HIV/AIDS and COVID-19 co-infected individuals denoted by $C(t)$, COVID-19 recovered individuals denoted by $R(t)$, and HIV/AIDS treated individuals denoted by $T_H(t)$ such that the total population is given by $N(t)=S(t)+P_1(t)+P_2(t)+I_H(t)+I_C(t)+R(t)+C(t)+T_H(t).$The susceptible individuals acquires HIV/AIDS at the standard incidence rate given by (1) ${\lambda }_H(t)=\frac{{\beta }_1}{N}(I_H(t)+{\rho }_{1}C(t)).$(1) where $1\le {\rho }_1<\mathrm{\infty }$ is the modification parameter that increases infectivity and ${\beta }_1$ is the HIV/AIDS transmission rate.

The susceptible individuals acquires COVID-19 at the standard incidence rate given by (2) ${\lambda }_{\text{C}}(\text{t})=\frac{{\beta }_2}{N}(I_C(t)+\omega C(t))$(2) where $1\le \omega <\mathrm{\infty }$ is the modification parameter that increases infectivity and ${\beta }_2$ is the COVID-19 transmission rate.

To formulate the HIV/AIDS and COVID-19 co-infection model we have the following assumptions: ${\pi }_1$, ${\pi }_2$ and $(1-{\pi }_1-{\pi }_2)$ are portions of the recruited individuals entered the COVID-19 protected class, the HIV protected class, and the susceptible class respectively, the susceptible class is increased by individuals entered from COVID-19 recovered class in which those individuals lost their immunity by the rate $\eta$, human populations are homogeneous and variable, no HIV transmission from HIV treated class and no vertical transmission, there is no simultaneous dual-infection transmission, HIV/AIDS and COVID-19 co-infected individuals are treated with the rate $\theta$ and entered to the HIV/AIDS treated class.

In this section using parameters given in Table 1, model variables in Table 2, and the model assumptions given, the schematic diagram for the transmission dynamics of HIV/AIDS and COVID co-infection is given in Figure 1.

Figure 1. Schematic diagram of HIV/AIDS and COVID-19 transmission dynamics with ${\lambda }_H(t)$ and ${\lambda }_C(t)$ given in (1) and (2), respectively.

Table 1. Model parameters symbols and definitions.

Symbol definitions
$\mu$	Natural death rate
$\mathrm{\Delta }$	Recruitment rate of human being
${\alpha }_1$	COVID-19 protection lose rate
${\alpha }_2$	HIV protection lose rate
$\theta$	Co-infected treatment rate
$\mathrm{\Phi }$	Modification parameter
$\phi$	Modification parameter
$d_2$	COVID-19 death rate
$d_1$	HIV/AIDS death rate
$\kappa$	COVID-19 recovery rate
$\gamma$	HIV infection treatment rate
${\beta }_1$	HIV/AIDS transmission rate
${\beta }_2$	COVID-19 transmission rate
${\pi }_1$	Portion of recruitment entered to COVID-19 protected class
${\pi }_2$	Portion of recruitment entered to HIV protected class
$\eta$	COVID-19 immunity lose rate
$d_3$	Death rate by HIV/AIDS and COVID-19 co-infection

Table 2. Biological meaning of variables.

Variables	Biological meaning
$S$	Susceptible individuals
$P_1$	Individuals who are protected from COVID-19
$P_2$	Individuals who are protected from HIV
$I_H$	HIV/AIDS infected individuals
$I_C$	COVID-19 infected individuals
$C$	Individuals co-infected with HIV and COVID-19
$R$	Individuals recovered from COVID-19 infection
$T_H$	HIV treated individuals

Using Figure 1 we derive the following system of non-linear differential equations of the co-infection (3) $\begin{array}{rl}\frac{dS}{dt}& =(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }+{\alpha }_1{P}_1+{\alpha }_2{P}_2+\eta \mathrm{R}-({\lambda }_H+{\lambda }_C+\mu )S,\\ \frac{d{P}_1}{dt}& ={\pi }_1\mathrm{\Delta }-(\delta {\lambda }_H+{\alpha }_1+\mu )P_1,\\ \frac{d{P}_2}{dt}& ={\pi }_2\mathrm{\Delta }-({\alpha }_2+\mu +\sigma {\lambda }_C)P_2,\\ \frac{d{I}_H}{dt}& ={\lambda }_{H}S+\delta {\lambda }_H{P}_1-(\mu +d_1+\gamma +\varphi {\lambda }_C)I_H,\\ \frac{d{I}_C}{dt}& ={\lambda }_{C}S+\sigma {\lambda }_C{P}_2-(\kappa +\mu +d_2+\phi {\lambda }_H)I_C,\\ \frac{dC}{dt}& =\phi {\lambda }_H{I}_C+\varphi {\lambda }_C{I}_H+\rho {\lambda }_C{T}_H-(\mu +d_3+\theta )C\\ \frac{dR}{dt}& =\kappa I_C-(\mu +\eta )R\\ \frac{d{T}_H}{dt}& =\gamma I_H+\theta C-\rho {\lambda }_C{T}_H-\mu T_H.\end{array}$(3) With initial conditions$S(0)>0,P_1(0)\ge 0,P_2(0)\ge 0,I_H(0)\ge 0,I_C(0)\ge 0,C(0)\ge 0,R(0)\ge 0$ and (4) $T_H(0)\ge 0.$(4) The sum of all the differential equations in (3) is (5) $\frac{dN}{dt}=\mathrm{\Delta }-\mu N-(d_1{I}_H+d_2{I}_C+d_{3}C).$(5)

2.1. Basic properties of the COVID-19 and HIV/AIDS co-infection model (3)

The model can be analysed by dealing with different quantitative and qualitative attributes. Since the system deals with human populations which cannot be negative, we need to show that all the state variables are always non-negative as well as the solutions of the system (3) remain positive with positive initial conditions (4) in the bounded region (6) $\mathrm{\Omega }=\left\{\begin{array}{c}(S,P_1,P_2,I_H,I_C,C,R,T_H)\in {\mathbb{R}}_{+}^8,N\le \frac{\mathrm{\Delta }}{\mu }\end{array}\right\}.$(6) It is important to show that each state variable defined in Table 2 with positive initial conditions (4) is non-negative for all-time $t>0$ in the bounded region given in (6) and the model (3) is mathematically and epidemiologically meaning full.

Theorem 1

(Non-negativity of the model solutions): At the initial conditions (4) the solutions $S(t)$, $P_1(t)$, $P_2(t)$, $I_H(t),I_C(t)$,$C(t),R(t)$ and $T_H(t)$ of the system (3) are non-negative for all time $t>0$.

Proof: Assume $S(0)>0$, $P_1(0)>0,P_2(0)>0$, $I_H(0)>0$, $I_C(0)>0$, $C(0)>0$, $R(0)>0,$ and $T_H(0)>0$ then for all t > 0, we have to prove that $S$ (t) > 0,$P_1(t)>0,P_2(t)>0$, $I_H(t)$ > 0, $I_C(t)$ > 0, $C(t)$ > 0,$R(t)$ > 0, and $T_H(t)$ > 0.

Define:$\tau$ = sup$\{t>0:S(\text{t})>0,P_1(t)>0,P_2(t)>0$, $I_H(t)>0,I_C(t)>0,C(t)>0$,$R(t)>0$ and $T_H(t)>0\}$.

Since $S(t),P_1(t),P_2(t),I_H(t),I_C(t),C(t),R(t)\text{ and }{T}_H(t)$ are continuous we deduce that $\tau >0$. If $\tau$  = ∞, then positivity holds.

But, if 0 < $\tau$<∞, $S(\tau )=0$ or $P_1(\tau )=0$ or $P_2(\tau )=0$ or $I_H(\tau )=0$ or $I_C(\tau )=0$ or $C(\tau )=$ or $R(\tau )=0$ or $T_H(0)=0.$

The first equation of the COVID-19 and HIV/AIDS co-infection model (3) gives $\frac{dS}{dt}+({\lambda }_H+{\lambda }_C+\mu )S=(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }+{\alpha }_1{P}_1+{\alpha }_2{P}_2+\eta \mathrm{R}$and applying the integrating factor method and after some calculations we have got $S(\tau )=M_{1}S(0)+M_1{\int }_0^{\tau }ex{p}^{\int ({\lambda }_H+{\lambda }_C+\mu ))dt}((1-{\pi }_1-{\pi }_2)\mathrm{\Delta }+{\alpha }_1{P}_1(\text{t})+{\alpha }_2{P}_2(\text{t})+\eta \mathrm{R}(\text{t}))dt>0$where $M_1=ex{p}^{-\left(\mu \tau +\underset{0}{\overset{\tau }{\int }}({\lambda }_H(w)+{\lambda }_C(w)\right)}>0,S(0)>0,P_1(\text{t})>0,P_2(\text{t})>0,\text{R}(\text{t})>0,$ and by the definition of $\tau ,$ the solution $S(\tau )>0$ hence $S(\tau )\ne 0$.

Again the second equation of the model (3) gives us $\frac{d{P}_1}{dt}+(\delta {\lambda }_H+{\alpha }_1+\mu )P_1={\pi }_1\mathrm{\Delta }$ and we have got $P_1(\tau )=M_1{P}_1(0)+M_1{\int }_0^{\tau }ex{p}^{\int ({\alpha }_1+\mu +\delta {\lambda }_H(t))dt}{\pi }_1\mathrm{\Delta }dt>0$ where $M_1=ex{p}^{-\left({\alpha }_1\tau +\mu \tau +\underset{0}{\overset{\tau }{\int }}(\delta {\lambda }_H(w)\right)}>0,P_1(0)>0,$ and from the definition of $\tau ,$ then the solution$P_1(\tau )>0$ hence $P_1(\tau )\ne 0$.

Similarly we have the results $P_2(\tau )>0$ hence $P_2(\tau )\ne 0$, $I_H(\tau )>0$ hence $I_H(\tau )\ne 0$, $I_C(\tau )>0$ hence $I_C(\tau )\ne 0$, $C(\tau )>0$ hence $C(\tau )\ne 0$, $R(\tau )>0$ hence $Q(\tau )\ne 0$, and $T_H(\tau )>0$ hence $T_H(\tau )\ne 0$.

Thus, $\tau =\mathrm{\infty }$ by the definition of $\tau$ and hence all the solutions of the system (3) are non-negative which is the required prove of the theorem.

Theorem 2

(Boundedness of the model solutions): The region $\mathrm{\Omega }$ given by (6) is boundedin ${\mathbb{R}}_{+}^8$.

Proof:

Since all the state variables are non-negative by Theorem 1, from (5) in the absence of infections we do have $\frac{dN}{dt}\le \mathrm{\Delta }-\mu N$. By incorporating standard comparison theorem we have got $\int \frac{dN}{\mathrm{\Delta }-\mu N}\le \int dt$ and integrating both sides gives $-\frac{1}{\mu }\ln (\mathrm{\Delta }-\mu N)\le t+c$ where $c$ is some constant and after some steps of calculations we have got $0\le N(t)\le \frac{\mathrm{\Delta }}{\mu }$ which means all possible solutions of the system (3) with positive initial conditions given in (4) enter into the bounded region $(6)$.

3. Qualitative analysis of the models

Before we analysed the HIV/AIDS and COVID-19 co-infection model (3), we need to gain some background about the HIV/AIDS mono-infection and COVID-19 mono-infection transmission dynamics.

3.1. HIV/AIDS mono-infection model

The HIV/AIDS mono-infection model of the co-infection model (3) at $P_1=I_C=C=Q=0$ is given by (7) $\begin{array}{rl}\frac{dS}{dt}& =(1-{\pi }_2)\mathrm{\Delta }+{\alpha }_2{P}_2-({\lambda }_H+\mu )S,\\ \frac{d{P}_2}{dt}& ={\pi }_2\mathrm{\Delta }-({\alpha }_2+\mu )P_2,\\ \frac{d{I}_H}{dt}& ={\lambda }_{H}S-(\mu +d_1+\gamma )I_H,\\ \frac{d{T}_H}{dt}& =\gamma I_H-\mu T_H.\end{array}$(7) where the total HIV/AIDS mono-infection population is given by $N_1(t)=S(t)+P_2(t)+I_H(t)+T_H(t),$ with force of infection given by ${\lambda }_H=\frac{{\beta }_1}{N_1}{I}_H$, and initial conditions $S(0)>0$, $P_2(0)\ge 0$, $I_H(0)\ge 0$, and $I_H(0)\ge 0$. In the region ${\mathrm{\Omega }}_1=\left\{\begin{array}{c}(S,P_2,I_H,T_H)\in {{\mathbb{R}}^4}_{+},\\ N_1\le \frac{\mathrm{\Delta }}{\mu }\end{array}\right\}$, it is easy to show that the set ${\mathrm{\Omega }}_1$ is positively invariant and a global attractor of all positive solutions of the HIV/AIDS mono-infected model (7). Hence, it is sufficient to consider the dynamics of model (7) in ${\mathrm{\Omega }}_1$ as epidemiologically and mathematically well-posed.

3.1.1. Stability of the disease-free equilibrium point

The disease-free equilibrium (DFE) of the HIV/AIDS mono-infection model (7) is obtained by setting each equation of the system (7) to zero. Also, at the DFE, there are no infections, and treated. Thus, the disease free equilibrium point is $E_H^0=({\text{S}}^0,P_2^0,I_H^0{T}_H^0)=\left(\frac{\mathrm{\Delta }}{\mu }\left(\frac{{\alpha }_2+\mu (1-{\pi }_2)}{{\alpha }_2+\mu }\right),\frac{\mathrm{\Delta }{\pi }_2}{{\alpha }_2+\mu },0,0\right)$.

The linear stability of the DFE ${\mathit{E}}_{\mathit{H}}^0$ is established using the van den Driesch and Warmouth next-generation matrix approach stated in [30] on the model (7). The matrices $\mathit{F}$ and $\mathit{V}$ for the new infection terms and the remaining transfer terms are given by $F=\left(\begin{array}{cc}\frac{{\beta }_1}{N_1^0}{\text{S}}^0& 0\\ 0& 0\end{array}\right)=\left(\begin{array}{cc}\frac{{\beta }_1{\alpha }_2+{\beta }_1\mu (1-{\pi }_2)}{{\alpha }_2+\mu }& 0\\ 0& 0\end{array}\right),$and $V=\left(\begin{array}{cc}\mu +d_1+\gamma & 0\\ -\gamma & \mu \end{array}\right).$The basic reproduction number of the HIV/AIDS mono-infection model (7) is the dominant eigenvalue of the next-generation matrix $F{V}^{-1}$ given by ${\mathcal{R}}_H=\frac{{\beta }_1(\mu (1-{\pi }_2)+{\alpha }_2)}{(\gamma +\mu +d_1)(\mu +{\alpha }_2)}.$Since the basic reproduction number of HIV/AIDs denoted by ${\mathcal{R}}_H$ is defined as the expected number of secondary cases generated by one infected individual during its entire period of infectiousness in a fully susceptible population, the DFE$E_H^0=({\text{S}}^0,P_2^0,I_H^0{T}_H^0)=\left(\frac{\mathrm{\Delta }}{\mu }\left(\frac{{\alpha }_2+\mu (1-{\pi }_2)}{{\alpha }_2+\mu }\right),\frac{\mathrm{\Delta }{\pi }_2}{{\alpha }_2+\mu },0,0\right)$ is locally asymptotically stable if ${\mathcal{R}}_H$ $<1$ and unstable otherwise.

3.1.2. Existence and uniqueness of endemic equilibrium point

Now making the right-hand side of the dynamical system given in Equation (7) as zero we have determined that $\begin{array}{rl}{S}^{\ast }& =\frac{(1-{\pi }_2)\mathrm{\Delta }({\alpha }_2+\mu )+{\alpha }_2{\pi }_2\mathrm{\Delta }}{({\alpha }_2+\mu )(\mu +{\lambda }_H^{\ast })},I_H^{\ast }=\frac{{\pi }_2\mathrm{\Delta }({\alpha }_2+\mu ){\lambda }_H^{\ast }+{\alpha }_2{\pi }_2\mathrm{\Delta }{\lambda }_H^{\ast }}{(\mu +d_1+\gamma )({\alpha }_2+\mu )(\mu +{\lambda }_H^{\ast })},\\ T_H^{\ast }& =\frac{(1-{\pi }_2)\mu \mathrm{\Delta }\gamma ({\alpha }_2+\mu ){\lambda }_H^{\ast }+{\alpha }_2{\pi }_2\mu \mathrm{\Delta }\gamma {\lambda }_H^{\ast }}{(\mu +d_1+\gamma )({\alpha }_2+\mu )(\mu +{\lambda }_H^{\ast })}.\end{array}$Let us put ${\text{m}}_1={\alpha }_2+\mu$, ${\text{m}}_2=\mu +d_1+\gamma$ and substitute $I_H^{\ast }$ in the HIV/AIDS force on infection we have determined that ${\lambda }_H^{\ast }=\frac{{\beta }_1(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1{\lambda }_H^{\ast }+{\beta }_1{\alpha }_2{\pi }_2\mathrm{\Delta }{\lambda }_H^{\ast }}{\begin{array}{c}(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1{\text{m}}_2+{\alpha }_2{\pi }_2\mathrm{\Delta }{\text{m}}_2+{\pi }_2\mathrm{\Delta }{\text{m}}_2(\mu +{\lambda }_H^{\ast })+(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1{\lambda }_H^{\ast }\\ +{\alpha }_2{\pi }_2\mathrm{\Delta }{\lambda }_H^{\ast }(1+\mu \gamma )+(1-{\pi }_2)\mu \mathrm{\Delta }\gamma {\text{m}}_1{\lambda }_H^{\ast }\end{array}}$and after some simplification we have derived the non-zero linear equation as (8) ${\text{a}}_1{\lambda }_H^{\ast }+{\text{a}}_0=0.$(8) where ${\text{a}}_1=-{\pi }_2\mathrm{\Delta }{\text{m}}_2-(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1-{\alpha }_2{\pi }_2\mathrm{\Delta }-(1-{\pi }_2)\mu \mathrm{\Delta }\gamma {\text{m}}_1-{\alpha }_2{\pi }_2\mu \mathrm{\Delta }\gamma <0$and

${\text{a}}_0={\text{m}}_1{\text{m}}_2[{\mathcal{R}}_H-1]>0$ if ${\mathcal{R}}_H>1$ since all parameters are positive.

Then we do have ${\lambda }_H^{\ast }=\frac{-{\text{a}}_0}{{\text{a}}_1}=\frac{-{\text{m}}_1{\text{m}}_2[{\mathcal{R}}_H-1]}{-[{\pi }_2\mathrm{\Delta }{\text{m}}_2+(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1+{\alpha }_2{\pi }_2\mathrm{\Delta }+(1-{\pi }_2)\mu \mathrm{\Delta }\gamma {\text{m}}_1+{\alpha }_2{\pi }_2\mu \mathrm{\Delta }\gamma ]}$ and hence

${\lambda }_H^{\ast }=\frac{-{\text{a}}_0}{{\text{a}}_1}=\frac{{\text{m}}_1{\text{m}}_2[{\mathcal{R}}_H-1]}{{\pi }_2\mathrm{\Delta }{\text{m}}_2+(1-{\pi }_2)\mathrm{\Delta }{\text{m}}_1+{\alpha }_2{\pi }_2\mathrm{\Delta }+(1-{\pi }_2)\mu \mathrm{\Delta }\gamma {\text{m}}_1+{\alpha }_2{\pi }_2\mu \mathrm{\Delta }\gamma \gamma }>0$ if and only if${\mathcal{R}}_H>1$.

Thus, based on the final result ${\lambda }_H^{\ast }>0$ there is a unique positive endemic equilibrium for the HIV/AIDS sub-model given in Equation (7) if and only if ${\mathcal{R}}_H>1$.

Lemma 1:

The HIV/AIDS mono-infection model given in Equation (7) has a unique endemic equilibrium solution if and only if ${\mathcal{R}}_H>1$.

3.1.3. DFE global asymptotic stability

Lemma 2:

(The Castillo-Chavez et al. criteria stated in [5]), if the HIV/AIDS mono-infection model can be written as $\begin{array}{rl}\frac{dY}{dt}& =G(Y,W),\\ \frac{dZ}{dt}& =H(Y,W),H(Y^0,0)=0.\end{array}$where $Y\in {\mathbb{R}}^m$ be the components of non-infected individuals and $W\in {\mathbb{R}}^n$ be the components of infected individuals including treated class and $E_H^0=(Y^0,0)$ denotes the DFE point of the dynamical system (7).

Assume (i) For $\left(\frac{dY}{dt}\right)=G(Y^0,0)$, $Y^0$ is globally asymptotically stable (GAS)

(ii) $H(Y,W)=BW-\check{H}(Y,W)$, $\check{H}(Y,W)\ge 0$ for $(Y,W)\in {\mathrm{\Omega }}_1$ where $B=D_{W}H(Y^0,0)$ is an M-matrix i.e. the off diagonal elements of $B$ are non-negative and ${\mathrm{\Omega }}_1$ is the region in which the system makes biological sense. Then the fixed point $E_H^0=(Y^0,0)$ is GAS equilibrium point of the system (7) whenever ${\mathcal{R}}_H<1$.

Lemma 3:

The DFE point $E_H^0=\left(\frac{(1-{\pi }_2)\mathrm{\Delta }({\alpha }_2+\mu )+{\alpha }_2{\pi }_2\mathrm{\Delta }}{\mu ({\alpha }_2+\mu )},\frac{{\pi }_3\mathrm{\Delta }}{{\alpha }_2+\mu },0,0\right)$ of the HIV/AIDS mono-infection model (7) is GAS if ${\mathcal{R}}_H<1$ and the two sufficient conditions given in Lemma 2 are satisfied.

Proof:

Here we are applying Lemma 2 on the HIV/AIDS mono-infection model (7) and we have got the following matrices $\begin{array}{rl}\frac{dY}{dt}& =G(Y,W)=\left[\begin{array}{c}(1-{\pi }_2)\mathrm{\Delta }+{\alpha }_2{P}_2-({\lambda }_H+\mu )S\\ {\pi }_2\mathrm{\Delta }-({\alpha }_2+\mu )P_2\end{array}\right],\end{array}$ $\begin{array}{rl}\frac{dW}{dt}& =H(Y,W)=\left[\begin{array}{c}{\lambda }_{H}S-(\mu +d_1+\gamma )I_H\\ \gamma I_H-\mu T_H\end{array}\right]\end{array}$ $\begin{array}{rl}G(Y^0,0)& =\left[\begin{array}{c}(1-{\pi }_2)\mathrm{\Delta }+{\alpha }_2{P}_2^0-\mu S^0\\ {\pi }_2\mathrm{\Delta }-({\alpha }_2+\mu )P_2^0\end{array}\right]\text{■}\end{array}$

where $Y^0=(S^0,P_2^0)=\left(\frac{(1-{\pi }_2)\mathrm{\Delta }({\alpha }_2+\mu )+{\alpha }_2{\pi }_2\mathrm{\Delta }}{\mu ({\alpha }_2+\mu )},\frac{{\pi }_2\mathrm{\Delta }}{{\alpha }_2+\mu }\right)$ is globally stable which satisfies condition (i) of Lemma 2

and $B=D_{W}H(Y^{\ast },0)=\left[\begin{array}{cc}{\beta }_1-(\mu +d_1+\gamma )& 0\\ \gamma & -\mu \end{array}\right].$After some steps of calculations, we have determined that $\check{H}(Y,W)=\left[\begin{array}{c}{\check{H}}_1(Y,W)\\ {\check{H}}_2(Y,W)\end{array}\right]=\left[\begin{array}{c}{\beta }_1{I}_H-\frac{{\beta }_1{I}_H}{N_1}S\\ 0\end{array}\right]=\left[\begin{array}{c}{\beta }_1{I}_H\left(1-\frac{S}{N_1}\right)\\ 0\end{array}\right].$Now since $S\le N_1$ we have got $\frac{S}{N_1}\le 1$ and ${\check{H}}_1(Y,W)\ge 0$, which satisfies condition (ii) of Lemma 2, thus, the DFE point $E_H^0=\left(\frac{(1-{\pi }_2)\mathrm{\Delta }({\alpha }_2+\mu )+{\alpha }_2{\pi }_2\mathrm{\Delta }}{\mu ({\alpha }_2+\mu )},\frac{{\pi }_2\mathrm{\Delta }}{{\alpha }_2+\mu },0,0\right)$ of the HIV/AIDS mono-infection model (7) is GAS if ${\mathcal{R}}_H<1$.

Biologically, whenever ${\mathcal{R}}_H<1$ the HIV/AIDS mono-infection disease dies out while the total population increases [6].

3.2. COVID-19 mono-infection model

The COVID-19 mono-infection model of the system (3) at $P_2=I_H=C=T_H=0$, is given by (10) $\begin{array}{rl}\frac{dS}{dt}& =(1-{\pi }_1)\mathrm{\Delta }+{\alpha }_1{P}_1+\eta \mathrm{Q}-({\lambda }_C+\mu )S,\\ \frac{d{P}_1}{dt}& ={\pi }_1\mathrm{\Delta }-({\alpha }_1+\mu )P_1,\\ \frac{d{I}_C}{dt}& ={\lambda }_{C}S-(\kappa +\mu +d_2)I_C,\\ \frac{dR}{dt}& =\kappa I_C-(\mu +\eta )R.\end{array}$(10) with the COVID-19 mono-infection infection rate (‘force of infection’) given by (11) ${\lambda }_C=\frac{{\beta }_2}{N_2}{I}_C(t).$(11) and with COVID-19 mono-infection initial conditions $S(0)>0$,$P_1(0)\ge 0,I_C(0)\ge 0,R(0)\ge 0,$ total population $N_2(t)=S(t)+P_1(t)+I_C(t)$ + $Q(t)$.

In the region ${\mathrm{\Omega }}_2=\left\{\begin{array}{c}(S,P_1,I_C,R)\in {{\mathbb{R}}^4}_{+},N_2\le \frac{\mathbf{\Delta }}{\mathit{\mu }}\end{array}\right\}$, it is easy to show that the set ${\mathrm{\Omega }}_2$ is positively invariant and a global attractor of all positive solutions of the COVID-19 mono-infection model (10). Hence, it is sufficient to consider the dynamics of model (10) in ${\mathrm{\Omega }}_2$ as epidemiologically and mathematically well-posed.

3.2.1. Stability of DFE point

The disease-free equilibrium point of COVID-19 mono-infection model (10) denoted by $E_C^0$ is obtained by making the entire model equations equal to zero at $I_C=0$. Thus, the COVID-19 mono-infection model disease-free equilibrium point is $E_C^0=({\text{S}}^0,P_1^0,I_C^0{\text{R}}^0)=\left(\frac{\mathrm{\Delta }}{\mu }\left(\frac{{\alpha }_1+\mu (1-{\pi }_1)}{{\alpha }_1+\mu }\right),\frac{\mathrm{\Delta }{\pi }_1}{{\alpha }_1+\mu },0,0\right)$.

The basic reproduction number of COVID-19 mono-infection is the average number of new COVID-19 infection generated by a typically COVID-19 infectious individual in a community. Similarly, using the approach in [30], we find COVID-19 reproduction number ${\mathcal{R}}_C$ to be ${\mathcal{R}}_C=\frac{{\beta }_2(\mu (1-{\pi }_1)+{\alpha }_1)}{(\mu +d_2+\kappa )(\mu +{\alpha }_1)}$.

Since the basic reproduction number of COVID-19 denoted by ${\mathcal{R}}_C$ is defined as the expected number of secondary cases generated by one infected individual during its entire period of infectiousness in a fully susceptible population, the DFE$E_C^0=({\text{S}}^0,P_1^0,I_C^0{\text{R}}^0)=\left(\frac{\mathrm{\Delta }}{\mu }\left(\frac{{\alpha }_1+\mu (1-{\pi }_1)}{{\alpha }_1+\mu }\right),\frac{\mathrm{\Delta }{\pi }_1}{{\alpha }_1+\mu },0,0\right)$ is locally asymptotically stable if ${\mathcal{R}}_C<1$ and unstable otherwise.

3.2.2. Existence and uniqueness of endemic equilibrium point

The endemic equilibrium point(s) of the COVID-19 mono-infection model given in Equation (10) is/are determined by making the right-hand side of the system zero and after some steps of computations we have got that (12) $\begin{array}{rl}{S}^{\ast }& =\frac{(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m}_3+{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3}{m_1{m}_2{m}_3({\lambda }_C^{\ast }+\mu )-m_1\eta \kappa {\lambda }_C^{\ast }},\\ P_1^{\ast }& =\frac{{\pi }_1\mathrm{\Delta }}{m_1},I_C^{\ast }=\frac{(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m}_3{\lambda }_C^{\ast }+{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3{\lambda }_C^{\ast }}{m_1{m_2}^2{m}_3({\lambda }_C^{\ast }+\mu )-m_1{m}_2\eta \kappa {\lambda }_C^{\ast }},\\ R^{\ast }& =\frac{(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m}_3\kappa {\lambda }_C^{\ast }+{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3\kappa {\lambda }_C^{\ast }}{m_1{m_2}^2{m_3}^2({\lambda }_C^{\ast }+\mu )-m_1{m}_2{m}_3\eta \kappa {\lambda }_C^{\ast }}.\end{array}$(12) where $m_1={\alpha }_1+\mu$, $m_2=\kappa +\mu +d_2$ and $m_3=\mu +\eta$.

Now substitute $I_C^{\ast }$ given in Equation (12) in Equation (11) we computed as

$N_2^{\ast }{\lambda }_C^{\ast }={\beta }_2{I}_C^{\ast }$ and gives as (13) $\begin{array}{rl}& (1-{\pi }_1)\mathrm{\Delta }{m}_1{m_2}^2{m_3}^2+{\alpha }_1{\pi }_1\mathrm{\Delta }{m_2}^2{m_3}^2+{\pi }_1\mathrm{\Delta }{m_2}^2{m_3}^2{\lambda }_C^{\ast }+{\pi }_1\mathrm{\Delta }{m_2}^2{m_3}^2\mu \\ & +(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m_3}^2{\lambda }_C^{\ast }-{\pi }_1\mathrm{\Delta }{m}_2{m}_3\eta \kappa {\lambda }_C^{\ast }\\ & +{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m_3}^2{\lambda }_C^{\ast }+(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m}_3\kappa {\lambda }_C^{\ast }+{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3\kappa {\lambda }_C^{\ast }\\ & -{\beta }_2(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m_3}^2-{\beta }_2{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m_3}^2=0.\end{array}$(13) Rearranging (13) we have derived the non-zero linear equation (14) $b_1{\lambda }_C^{\ast }+b_0=0.$(14) where $\begin{array}{rl}{b}_1& ={\pi }_1\mathrm{\Delta }{m}_2{m}_3(m_2{m}_3-\eta \kappa )+(1-{\pi }_1)\mathrm{\Delta }{m}_1{m}_2{m}_3(m_3+\kappa )\\ & +{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3(m_3+\kappa )>0,\end{array}$and

$b_0=\mathrm{\Delta }{m}_2{m_3}^2((1-{\pi }_1)m_1{m}_2+{\alpha }_1{\pi }_1{m}_2)[1-R_C]<0$ if $R_C>0$ since parameters have positive values.

From Equation (14) we do have

${\lambda }_C^{\ast }=-\frac{b_0}{b_1}=\frac{\mathrm{\Delta }{m}_2{m_3}^2({\pi }_2{m}_1{m}_2+{\alpha }_1{\pi }_1{m}_2)[R_C-1]}{{\pi }_1\mathrm{\Delta }{m}_2{m}_3(m_2{m}_3-\eta \kappa )+{\pi }_2\mathrm{\Delta }{m}_1{m}_2{m}_3(m_3+\kappa )+{\alpha }_1{\pi }_1\mathrm{\Delta }{m}_2{m}_3(m_3+\kappa )}$ > 0 if and only if ${\mathcal{R}}_C>1$ since all parameters are positive.

Thus, the COVID-19 mono-infection model given in Equation (10) has a unique positive endemic equilibrium point only whenever ${\mathcal{R}}_C>1$.

Lemma 4:

The COVID-19 mono-infection model given in Equation (10) has a unique positive endemic equilibrium if and only if ${\mathcal{R}}_C>1$.

3.2.3. DFE point global asymptotic stability

Lemma 5:

The DFE point ${\mathit{E}}_{\mathit{C}}^0=\left(\frac{(1-{\mathit{\pi }}_1)\mathbf{\Delta }({\mathit{\alpha }}_1+\mathit{\mu })+{\mathit{\alpha }}_1{\mathit{\pi }}_1\mathbf{\Delta }}{\mathit{\mu }({\mathit{\alpha }}_1+\mathit{\mu })},\frac{{\mathit{\pi }}_1\mathbf{\Delta }}{{\mathit{\alpha }}_1+\mathit{\mu }},0,0\right)$ of the COVID-19 mono-infection model (10) is GAS if ${\mathcal{R}}_{\mathit{C}}<1$ whenever the two sufficient conditions given in Lemma 2 are satisfied.

Proof: Here we are applying Lemma 2 on the COVID-19 mono-infection model (10) and setting $Y\in {\mathbb{R}}^2$ be the components of non-infected individuals and $W\in {\mathbb{R}}^2$ be the components of infected individuals including recovery class. Then we have determined the following matrices $\begin{array}{rl}\frac{dY}{dt}& =G(Y,W)=\left[\begin{array}{c}(1-{\pi }_1)\mathrm{\Delta }+{\alpha }_1{P}_1+\eta \mathrm{R}-({\lambda }_C+\mu )S\\ {\pi }_1\mathrm{\Delta }-({\alpha }_1+\mu )P_1\end{array}\right],\end{array}$ $\begin{array}{rl}\frac{dW}{dt}& =H(Y,W)=\left[\begin{array}{c}{\lambda }_{C}S-(\kappa +\mu +d_2)I_C\\ \kappa I_C-(\mu +\eta )R\end{array}\right],\end{array}$ $\begin{array}{rl}G(Y,0)& =\left[\begin{array}{c}(1-{\pi }_1)\mathrm{\Delta }+{\alpha }_1{P}_1-\mu S\\ {\pi }_1\mathrm{\Delta }-({\alpha }_1+\mu )P_1\end{array}\right],\end{array}$

and $B=D_{W}H(Y^{\ast },0)=\left[\begin{array}{cc}\frac{{\beta }_2{S}^0}{S^0+P_1^0}-(\kappa +\mu +d_2)& 0\\ \kappa & -(\mu +\eta )\end{array}\right].$After some steps of calculations we have determined that $\check{H}(Y,W)=\left[\begin{array}{c}{\check{H}}_1(Y,W)\\ {\check{H}}_2(Y,W)\end{array}\right]=\left[\begin{array}{c}\frac{{\beta }_2{S}^0{I}_C}{S^0+P_1^0}-\frac{{\beta }_2{I}_{C}S}{N_2}\\ 0\end{array}\right]=\left[\begin{array}{c}{\beta }_2{I}_C\left(\frac{S^0}{S^0+P_1^0}-\frac{S}{N_2}\right)\\ 0\end{array}\right].$since $\le S^0$, $P_1<P_1^0$ we have got $S-S^0\le 1$, $P_1-P_1^0\le 1$ and ${\check{H}}_1(Y,W)\ge 0$, thus, the DFE point $E_C^0=\left(\frac{(1-{\pi }_1)\mathrm{\Delta }({\alpha }_1+\mu )+{\alpha }_1{\pi }_1\mathrm{\Delta }}{\mu ({\alpha }_1+\mu )},\frac{{\pi }_1\mathrm{\Delta }}{{\alpha }_1+\mu },0,0\right)$ of the COVID-19 mono-infection model (10) is GAS if ${\mathcal{R}}_C<1$.

Biologically, whenever ${\mathcal{R}}_C<1$ the COVID-19 mono-infection disease dies out while the total population increases [6].

3.3. COVID-19 and HIV/AIDS co-infection model

Having analysed the dynamics of the two mono-infection models, that is HIV/AIDS mono-infection model (7) and the COVID-19 mono-infection model (10) the main HIV/AIDS and COVID-19 co-infection model (3) is now considered and analysed in the bounded region $\mathrm{\Omega }$ given in Equation (6).

3.3.1. Stability of DFE point

The disease-free equilibrium point of the HIV/AIDS and COVID-19 Co-infection model (3) denoted by $E_{HC}^0$ is obtained by making all the equations equal to zero providing that $I_H=I_C=\text{C}=0$. Thus, the HIV/AIDS and COVID-19 Co-infection model (3) DFE point is $\begin{array}{rl}{E}_{HC}^0& =({\text{S}}^0,P_1^0,P_2^0,I_H^0,I_C^0,,{\text{C}}^0,{\text{R}}^0,T_H^0)\\ & \frac{(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }({\alpha }_1+\mu )({\alpha }_2+\mu )+{\alpha }_1{\pi }_1\mathrm{\Delta }+{\alpha }_2{\pi }_2\mathrm{\Delta }({\alpha }_1+\mu )}{\mu ({\alpha }_1+\mu )({\alpha }_2+\mu )},\frac{{\pi }_1\mathrm{\Delta }}{{\alpha }_1+\mu },\\ & \frac{{\pi }_2\mathrm{\Delta }}{{\alpha }_2+\mu },0,0,0,0,0).\end{array}$Similarly, using the approach in [30], the model (3) basic reproduction number denoted by ${\mathcal{R}}_{HC}$ is to be determined as $F{V}^{-1}=\left[\begin{array}{ccccc}\frac{{\beta }_1(\mu (1-{\pi }_2)+{\alpha }_2)}{(\gamma +\mu +d_1)(\mu +{\alpha }_2)}& 0& 0& 0& 0\\ 0& \frac{{\beta }_2(\mu (1-{\pi }_1)+{\alpha }_1)}{(\mu +d_2+\kappa )(\mu +{\alpha }_1)}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right].$The basic reproduction number of COVID-19 and HIV AIDS co-infected model is the dominant eigenvalue of the next-generation matrix $F.V^{-1}$ given by ${\mathcal{R}}_{HC}=max\left\{\frac{{\beta }_1(\mu (1-{\pi }_2)+{\alpha }_2)}{(\gamma +\mu +d_1)(\mu +{\alpha }_2)},\frac{{\beta }_2(\mu (1-{\pi }_1)+{\alpha }_1)}{(\mu +d_2+\kappa )(\mu +{\alpha }_1)}\right\},$ where ${\mathcal{R}}_{HC}=max\{{\mathcal{R}}_H,{\mathcal{R}}_C\}$, ${\mathcal{R}}_H$ represent the basic reproduction number for HIV/AIDS mono-infection model (7), ${\mathcal{R}}_C$ and ${\mathcal{R}}_{HC}$ are the basic reproduction numbers for COVID-19 mono-infection model (10) and COVID-19 and HIV/AIDS co-infection model (3), respectively.

Since the basic reproduction number of HIV/AIDS and COVID-19 co-infection denoted by ${\mathcal{R}}_{HC}$ is defined as the expected number of secondary cases generated by one co-infected individual during its entire period of infectiousness in a fully susceptible population, the DFE

$E_{HC}^0=\left(\frac{(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }({\alpha }_1+\mu )({\alpha }_2+\mu )+{\alpha }_1{\pi }_1\mathrm{\Delta }+{\alpha }_2{\pi }_2\mathrm{\Delta }({\alpha }_1+\mu )}{\mu ({\alpha }_1+\mu )({\alpha }_2+\mu )},\frac{{\pi }_1\mathrm{\Delta }}{{\alpha }_1+\mu },\frac{{\pi }_2\mathrm{\Delta }}{{\alpha }_2+\mu },0,0,0,0,0\right)$ islocally asymptotically stable if ${\mathcal{R}}_{HC}$ $<1$ and unstable otherwise.

3.3.2. Endemic equilibrium of the model (3)

The endemic equilibrium point of the complete model (3) obtained by making each differential equation as zero is given by $\begin{array}{rl}{S}^{\ast }& =\frac{(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }+{\alpha }_1{P}_1^{\ast }+{\alpha }_2{P}_2^{\ast }+\eta {\text{R}}^{\ast }}{({\lambda }_H^{\ast }+{\lambda }_C^{\ast }+\mu )},P_1^{\ast }=\frac{{\pi }_1\mathrm{\Delta }}{(\delta {\lambda }_H^{\ast }+{\alpha }_1+\mu )},\\ P_2^{\ast }& =\frac{{\pi }_2\mathrm{\Delta }}{(\delta {\alpha }_2+\mu +\sigma {\lambda }_C^{\ast })},I_H^{\ast }=\frac{{\lambda }_H^{\ast }{S}^{\ast }+\delta {\lambda }_H^{\ast }{P}_1^{\ast }}{(\mu +d_1+\gamma +\varphi {\lambda }_C^{\ast })},I_C^{\ast }=\frac{{\lambda }_C^{\ast }{S}^{\ast }+\delta {\lambda }_C^{\ast }{P}_2^{\ast }}{(\kappa +\mu +d_2+\phi {\lambda }_H^{\ast })},\\ C^{\ast }& =\frac{\phi {\lambda }_H^{\ast }{I}_C^{\ast }+\varphi {\lambda }_C^{\ast }{I}_H^{\ast }+\rho {\lambda }_C^{\ast }{T}_H^{\ast }}{(\mu +d_3+\theta )},R^{\ast }=\frac{\kappa I_C^{\ast }}{(\mu +\eta )},T_H^{\ast }=\frac{\gamma I_H^{\ast }+\theta C^{\ast }}{(\rho {\lambda }_C^{\ast }+\mu )}.\end{array}$Since the model is highly non-linear the explicit computation of the endemic equilibrium of the full model (3) in terms of model parameters is difficult analytically; however, based on the previous analyses of the HIV/AIDS and COVID-19 mono-infection models the model (3) endemic equilibrium point $E_{HC}^{\ast }=(S^{\ast },P_1^{\ast },P_1^{\ast },I_H^{\ast },I_C^{\ast },C^{\ast },R^{\ast },T_H^{\ast })$ exists whenever ${\mathcal{R}}_H>1$ and ${\mathcal{R}}_C>1$ i.e. ${\mathcal{R}}_{HC}>1$. We will give the stability explanation of $E_{HC}^{\ast }$ in our numerical simulation part.

3.3.3. Bifurcation analysis of the co-infection model

Let $S=y_1$, $P_1=y_2$, $P_2=y_3$, $I_H=y_4$, $I_C=y_5$, $C=y_6$, $R=y_7$, and $T_H=y_8$ so that $N=y_1+y_2+y_3+y_4$ + $y_5$, $+y_6$ + $y_7$, $+y_8$.

Further, by using vector notation $Y=(y_1,y_2,y_3,y_4,y_5,y_6,y_7,y_8{)}^T,$ the full model (3) can be written in the form $\frac{dY}{dt}=F(Y)$ with $G=(g_1,g_2,g_3,g_4,g_5,g_6,g_7,g_8{)}^T$, as follows (15) $\begin{array}{rl}\frac{d{y}_1}{dt}& =g_1=(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }+{\alpha }_1{y}_2+{\alpha }_2{y}_3+\eta y_7-({\lambda }_H+{\lambda }_C+\mu )y_1\\ \frac{d{y}_2}{dt}& =g_2={\pi }_1\mathrm{\Delta }-(\delta {\lambda }_H+{\alpha }_1+\mu )y_2,\\ \frac{d{y}_3}{dt}& =g_3={\pi }_2\mathrm{\Delta }-({\alpha }_2+\mu +\sigma {\lambda }_C)y_3,\\ \frac{d{y}_4}{dt}& =g_4={\lambda }_H{y}_1+\delta {\lambda }_H{y}_2-(\mu +d_1+\gamma +\varphi {\lambda }_C)y_4,\\ \frac{d{y}_5}{dt}& =g_5={\lambda }_C{y}_1+\sigma {\lambda }_C{y}_3-(\kappa +\mu +d_2+\phi {\lambda }_H)y_5,\\ \frac{d{y}_6}{dt}& =g_6=\phi {\lambda }_H{y}_5+\varphi {\lambda }_C{y}_4+\rho {\lambda }_C-(\mu +d_3+\theta )y_6,\\ \frac{d{y}_7}{dt}& =g_7=\kappa y_5-(\mu +\eta )y_7,\\ \frac{d{y}_8}{dt}& =g_8=\gamma y_4+\theta y_6-\rho {\lambda }_C-\mu y_8,\end{array}$(15) with ${\lambda }_H=\frac{{\beta }_1}{N}[y_4+{\rho }_1{y}_6]$ where $1\le {\rho }_1<\mathrm{\infty }$ and ${\lambda }_C={\beta }_2[y_5+\omega y_6]$ where $1\le \omega <\mathrm{\infty }$ then the method entails evaluating the Jacobian of the system (13) at the DFE $E_{HC}^0$, denoted by $J(E_{HC}^0)$ and this gives us $J(E_{HC}^0)=\left(\begin{array}{cccccccc}-\mu & {\alpha }_1& {\alpha }_2& G_1& G_2& G_3& \eta & 0\\ 0& -({\alpha }_1+\mu )& 0& G_4& 0& G_5& 0& 0\\ 0& 0& -({\alpha }_2+\mu )& 0& G_6& G_7& 0& 0\\ 0& 0& 0& G_8& 0& G_9& 0& 0\\ 0& 0& 0& 0& G_{10}& G_{11}& 0& 0\\ 0& 0& 0& 0& 0& -(\mu +d_3+\theta )& 0& 0\\ 0& 0& 0& 0& \kappa & 0& -(\mu +\eta )& 0\\ 0& 0& 0& \gamma & 0& \theta & 0& -\mu \end{array}\right).$where $G_1=-\frac{{\beta }_1}{N^0}{y}_1^0$, $G_2=-{\beta }_2{y}_1^0,G_3=-\frac{{\beta }_1}{N^0}{\rho }_1{y}_1^0-{\beta }_2\omega y_1^0$,$G_4=-\frac{{\beta }_1}{N^0}{y}_2^0,G_5=-\frac{{\beta }_1}{N^0}{\rho }_1{y}_2^0$, $G_6=-{\beta }_2{y}_3^0$, $G_7=-{\beta }_2\omega y_3^0$, $G_8=\frac{{\beta }_1}{N^0}{y}_1^0$ + $\frac{{\beta }_1}{N^0}{y}_2^0\text{ - }(\mu +d_1+\gamma )$, $G_9=\frac{{\beta }_1}{N^0}{\rho }_1{y}_1^0+\frac{{\beta }_1}{N^0}{\rho }_1{y}_2^0$, $G_{10}={\beta }_2{y}_1^0+{\beta }_2{y}_3^0-(\kappa +\mu +d_2)$, $G_{11}={\beta }_2\omega y_1^0+{\beta }_2\omega y_3^0$,

Without loss of generality, consider the case when ${\mathcal{R}}_C>{\mathcal{R}}_H,$ and ${\mathcal{R}}_{HC}=1$, so that ${\mathcal{R}}_C=1$. Furthermore, let ${\beta }_2$ $={\beta }^{\ast }$ is chosen as a bifurcation parameter. Solving for ${\beta }_2$ from ${\mathcal{R}}_2=1$ as ${\mathcal{R}}_C=\frac{{\beta }_2(\mu (1-{\pi }_1)+{\alpha }_1)}{(\mu +d_2+\kappa )(\mu +{\alpha }_1)}=1$ we have got the value ${\beta }^{\ast }={\beta }_2=\frac{(\mu +d_2+\kappa )(\mu +{\alpha }_1)}{(\mu (1-{\pi }_1)+{\alpha }_1)}$.

After solving the Jacobian $J(E_{HC}^0)$ of the system (13) at the DFE, with ${\beta }_2={\beta }^{\ast }$, we obtained the eigenvalues as

${\lambda }_1=-\mu$ < 0 or ${\lambda }_2=-({\alpha }_1+\mu )<0$ or ${\lambda }_3=-({\alpha }_2+\mu )<0$ or ${\lambda }_4=G_8=\frac{{\beta }_1}{N^0}{y}_1^0$+$\frac{{\beta }_1}{N^0}{y}_2^0\text{ - }(\mu +d_1+\gamma )=(\mu +d_1+\gamma )[{\mathcal{R}}_H-1]<0$ if ${\mathcal{R}}_H<1$ or ${\lambda }_5=0$ or ${\lambda }_6=-(\mu +d_3+\theta )<0$ or ${\lambda }_7=-(\mu +\eta )<0$ or ${\lambda }_8=-\mu <0.$

Hence all eigenvalues are negative whenever ${\mathcal{R}}_{HC}<1$. It follows that the Jacobian $J(E_{HC}^0)$ of (13) at the DFE, with ${\beta }_2={\beta }^{\ast }$, denoted by $J_{{\beta }^{\ast }}$, has a simple zero eigenvalue (with all other eigenvalues having negative real part). Hence, the Centre Manifold theory can be used to analyse the dynamics of the model (3). In particular, the Castillo-Chavez and Song theorem [6], will be used to show that the model (3) undergoes forward bifurcation at ${\mathcal{R}}_C=1$. Eigenvectors of $J_{{\beta }^{\ast }}$: For the case when ${\mathcal{R}}_C=1$, the right eigenvectors of the Jacobian of (15) at ${\beta }_2={\beta }^{\ast }$ (denoted by $J_{{\beta }^{\ast }})$ associated with the zero eigenvalue given by $u=(u_1,u_2,u_3,u_4,u_5,u_6,u_7,u_8{)}^T$. $\begin{array}{rl}{u}_1& =\frac{{\alpha }_2{G}_6(\mu +\eta )u_5+({\alpha }_2+\mu )(\mu +\eta )G_2{u}_5+({\alpha }_2+\mu )\eta \kappa u_5}{\mu ({\alpha }_2+\mu )(\mu +\eta )},u_2=0,\\ u_3& =\frac{G_6}{{\alpha }_2+\mu }{u}_5,u_4=0,u_5=u_5>0,u_6=0,u_7=\frac{\kappa }{\mu +\eta }{u}_5\text{and}{u}_8=0.\end{array}$The left eigenvectors associated with the zero eigenvalue at ${\beta }_2={\beta }_2^{\ast }$ satisfying $v.w=1,$ given by $\begin{array}{rl}v& =(v_1,v_2,v_3,v_4,v_5,v_6,v_7,v_8)\text{are}{v}_1=v_2=v_3=v_4=v_6=v_7=v_8=0\text{and}\\ v_5& =v_5>0.\end{array}$After going through detailed computations and simplification, we have the following bifurcation coefficients $a$ and $b$ as $\begin{array}{rl}a& =2v_5{u}_1{u}_5\frac{{\mathrm{\partial }}^2{f}_5(0,0)}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_5}+2v_5{u}_3{u}_5\frac{{\mathrm{\partial }}^2{f}_5(0,0)}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_5}=2{\beta }_2^{\ast }{v}_5{u}_5[u_1+u_3].\end{array}$ $\begin{array}{rl}& =2{\beta }_2^{\ast }{v}_5{u}_5^2\left[\frac{\begin{array}{c}-{\alpha }_2{\beta }_2{y}_3^0(\mu +\eta )-({\alpha }_2+\mu )(\mu +\eta ){\beta }_2{y}_1^0\\ -({\alpha }_2+\mu )\eta \kappa -\mu (\mu +\eta ){\beta }_2{y}_3^0\end{array}}{\mu ({\alpha }_2+\mu )(\mu +\eta )}\right].\end{array}$Thus (16) $\begin{array}{rl}\text{a}& =-2{\beta }_2^{\ast }{v}_5{u}_5^2\left[\frac{\begin{array}{l}{\alpha }_2{\beta }_2{y}_3^0(\mu +\eta )+({\alpha }_2+\mu )(\mu +\eta ){\beta }_2{y}_1^0\\ +({\alpha }_2+\mu )\eta \kappa +\mu (\mu +\eta ){\beta }_2{y}_3^0\end{array}}{\mu ({\alpha }_2+\mu )(\mu +\eta )}\right]<0\text{and}\\ \text{b}& =v_5{u}_5\frac{{\mathrm{\partial }}^2{\text{f}}_5(0,0)}{\mathrm{\partial }{\text{y}}_5\mathrm{\partial }{\beta }_2}=v_5{u}_5({\text{y}}_3^0+{\text{y}}_1^0)>0.\end{array}$(16) Hence, it follows from in Castillo-Chavez and Song [6] that the model (3) do not exhibits a backward bifurcation at ${\mathcal{R}}_{\text{HC}}={\mathcal{R}}_{\text{C}}=1.$ Thus, there is only a DFE point

$E_{HC}^0=\left(\frac{(1-{\pi }_1-{\pi }_2)\mathrm{\Delta }({\alpha }_1+\mu )({\alpha }_2+\mu )+{\alpha }_1{\pi }_1\mathrm{\Delta }+{\alpha }_2{\pi }_2\mathrm{\Delta }({\alpha }_1+\mu )}{\mu ({\alpha }_1+\mu )({\alpha }_2+\mu )},\frac{{\pi }_1\mathrm{\Delta }}{{\alpha }_1+\mu },\frac{{\pi }_2\mathrm{\Delta }}{{\alpha }_2+\mu },0,0,0,0,0\right)$ of theCOVID-19 and HIV/AIDS co-infection model (3) exists which means there is no positive endemic equilibrium point whenever ${\mathcal{R}}_{HC}<1$.

Note: In sections 3.1.5 and 3.2.5 we have shown both the disease-free equilibriums of the HIV/AIDS mono-infection and COVID-19 mono-infection models are GAS whenever the corresponding basic reproduction numbers values are less than unity. Thus based on this result and Equation (16) the DFE point of the COVID-19 and HIV/AIDS co-infection model (3) is GAS if ${\mathcal{R}}_{HC}=\text{max}\{{\mathcal{R}}_H,{\mathcal{R}}_C\}<1.$

4. Sensitivity and numerical analysis

To support some qualitative results in the previous sections, we have performed several sensitivity analyses and numerical simulations. To find relevant model parameters illustrated in Table 3, we have observed and reviewed different literatures that have been done by different scholars and we assume realistic values for the purpose of numerical illustration.

Table 3. Parameter values used for simulations.

Symbol	Value	Source
$\mathrm{\Delta }$	1000 humans/day	[1]
$\mu$	(1/64.5*365)/day	[22]
${\alpha }_1$	0.0015/day	Assume
${\alpha }_2$	0.0004/day	[18]
$d_1$	0.00023/day	[18]
$d_2$	0.15/day	[22]
$\theta$	0.0021/day	Assume
$\mathrm{\Phi }$	1.2 no unit	[22]
$\delta ,\sigma$	0.81 no unit	Assume
$\gamma$	0.0023/day	[13]
$d_3$	0.15/day	[22]
$\phi$	1 no unit	Assume
$\eta$	0.002/day	Assume
${\beta }_1$	0.3425/day	[22]
${\beta }_2$	0.1175/day	[22]
${\pi }_1$	0.012/day	Assume
${\pi }_2$	0.9813/day	Assume
${\pi }_3$	0.0067/day	[18]
$\kappa$	0.005/day	[1]
$\phi$	1 no unit	[22]

4.1. Investigation of sensitivity analysis

Definition: The normalized forward sensitivity index of a variable, $x$, that depends differentially on a parameter $\xi$ is defined as SI $(\xi )=\frac{\mathrm{\partial }x}{\mathrm{\partial }\xi }\ast \frac{\xi }{x}$ [8,15].

Sensitivity indices allow us to investigate the relative importance of different parameters in HIV/AIDS and COVID-19 co-infection transmission. The most sensitive parameter has a magnitude of the sensitivity index larger than that of all other parameters. We can calculate the sensitivity index in terms of basic reproduction numbers ${\mathcal{R}}_H$ and ${\mathcal{R}}_C$ since ${\mathcal{R}}_{HC}=\text{max}\{{\mathcal{R}}_H,{\mathcal{R}}_C\}$.

Using the baseline parameter values given in Table 3, we have derived Table 4 and Table 5 to show the sensitivity indices of the model parameters.

Table 4. Sensitivity indices of ${\mathcal{R}}_{HC}={\mathcal{R}}_H.$

Sensitivity Index	Values
$SI({\beta }_1)$	+1
$SI({\alpha }_2)$	0.03
$SI({\pi }_2)$	−0.71
$SI(d_1)$	0.40
$SI(\gamma )$	−0.61

Table 5. Sensitivity indices of ${\mathcal{R}}_{HC}={\mathcal{R}}_C.$

Sensitivity Index	Amanfu
$SI({\beta }_2)$	+1
$SI(d_2)$	0.01
$SI(\kappa )$	−0.6
$SI({\alpha }_1)$	0.21
$SI({\pi }_1)$	−0.72

In this paper, with baseline parameter values given in Table 3 we have computed ${\mathcal{R}}_H=1.34$ implied that HIV/AIDS spreads in the community and also we have determined the indices in Table 4. And sensitivity analysis shows that the HIV/AIDS transmission rate ${\beta }_1$ has the highest impact on the basic reproduction number of HIV/AIDS (${\mathcal{R}}_H$).

Similarly, using parameter values given in Table 3 we have computed ${\mathcal{R}}_C=3.21$ implied that COVID-19 is spreading throughout the community and also we have computed the sensitivity indices as shown in Table 5. And sensitivity analysis shows that the foremost sensitive positive parameter is the COVID-19 spreading rate ${\beta }_2$. Using Table 4 and Table 5, biologically we can conclude that the most sensitive parameters to ${\mathcal{R}}_H$ and ${\mathcal{R}}_C$ which can be controlled through prevention and control measures are found to be the HIV/AIDS transmission rate ${\beta }_1$ and the COVID-19 transmission rate ${\beta }_2,$ respectively.

Figure 2 illustrates the sensitivity indices of the model parameters where the transmission rates ${\beta }_1$ and ${\beta }_2$ are biologically the most sensitive parameters directly proportional to the basic reproduction numbers of HIV/AIDS and COVID-19 respectively. Similarly, the portion of recruitment rate entered to HIV protected group $({\pi }_2)$, the portion of recruitment rate entered to COVID-19 protected group $({\pi }_1)$, and the COVID-19 recovery rate $(\kappa )$ are the most sensitive parameters indirectly proportional to the basic reproduction numbers.

Figure 2. Sensitivity indices of model parameters on ${\mathcal{R}}_{HC}.$

4.2. Numerical simulations

In this section, we conduct numerical solutions of the system (3). For numerical solutions of the proposed model, the initial conditions have non-negative values since the number of individuals in each compartment cannot be negative. All numerical computations were performed in MATLAB using the ode45 embedded Runge–Kutta numerical method. Now, we start to investigate the stability of equilibrium points, examine the influences of various parameters, and explore the impact of intervention strategies. For simulation purposes, we specified the non-negative initial values $(\text{S}(0),P_1(0),P_2(0),I_H(0),I_C(0),\text{ C}(0),\text{R}(0),T_H(0))=(1000,300,200,100,170,80,75,75)$ and baseline parameter values in Table 3.

4.2.1. Simulation of the model whenever ${\mathcal{R}}_{\mathit{H}\mathit{C}}<1$

The scenario in Figure 3 illustrates the dynamical anticipation of the model (3) over time when the value of ${\mathcal{R}}_{HC}<1$. This numerical evaluation confirms the theoretical proof in (3.3.1). The basic reproduction number is computed as ${\mathcal{R}}_{HC}=0.38$. Figure 3 illustrated that after 20 days the solutions of the COVID-19 and HIV/AIDS co-infection dynamical system (3) are converging to the DFE point whenever ${\mathcal{R}}_{HC}=max\{{\mathcal{R}}_H,{\mathcal{R}}_C\}=max\{0.12,0.38\}=0.38<1$. Biologically it means that the outbreak of COVID-19 and HIV/AIDS co-infection throughout the community will be eradicated in the near futureto. Ultimately, the numerical simulation and the theoretical analysis in (3.3.1) are compatible.

Figure 3. Behaviour of the co-infection model (3) solutions at ${\mathcal{R}}_{HC}=0.38<1$.

4.2.2. Simulation of the model whenever ${\mathcal{R}}_{\mathit{H}\mathit{C}}>1$

Now, we turn to investigate the stability of the endemic equilibrium point of model (3) where ${\mathcal{R}}_{HC}=3.21>1.$ In Figure 4 all the numerical routines show that numerical solutions trajectories attract toward the endemic stable equilibrium of the system. Figure 4 depicts the time series evaluation of model (3) with baseline parameter values given in Table 3. Now, we computed the basic reproduction number as ${\mathcal{R}}_{HC}=max\{{\mathcal{R}}_H,{\mathcal{R}}_C\}=max\{1.34,3.21\}=3.21>1$, which is greater than one. That means that the model (3) solutions converges to the endemic equilibrium point whenever ${\mathcal{R}}_{HC}=3.21>1$ (Figure 5).

Figure 4. Behaviour of the co-infection model (3) solutions at ${\mathcal{R}}_{HC}=3.21>1$.

Figure 5. Impact of ${\beta }_1$ on $C$.

4.2.3. Effect of HIV transmission on COVID-19 infection

Figure 5 illustrates the effect of HIV transmission rate ${\beta }_1$ on the number of HIV/AIDS and COVID-19 co-infectious population $C.$ The simulation result shows that increasing the value of ${\beta }_1$ increases the number of HIV/AIDS and COVID-19 co-infectious individuals in the community. Consequently, increasing HIV transmission rate ${\beta }_1$ from 0.00001 to 0.8 leads to a highly increase of HIV/AIDS and COVID-19 co-infection $C$.

Figure 6. Impact of ${\beta }_2$ on $C$.

4.2.4. Effect of COVID-19 transmission on HIV/AIDS infection

Figure 6 illustrates the effect of COVID-19 transmission rate ${\beta }_2$ on the number of HIV/AIDS and COVID-19 co-infectious population $C.$ The simulation result shows that increasing the value of ${\beta }_2$ increases the number of HIV/AIDS and COVID-19 co-infectious individuals in the community. Consequently, increasing COVID-19 transmission rate ${\beta }_2$ from 0.00001 to 0.8 leads to a high increase of HIV/AIDS and COVID-19 co-infection $C$.

4.2.5. Effect of treatment on HIV/AIDS infection

Figure 7 illustrates the effect of treatment rate $\gamma$ on the number of HIV infectious population $I_H$. The result shows that increasing the value of $\gamma$ decreases the number of HIV/AIDS infectious individuals in the community. Consequently, increasing treatment intervention measure from 0.25 to 0.85 against HIV/AIDS highly decreases the number of HIV/AIDS infectious population $I_H$.

Figure 7. Impact of $\gamma$ on $I_H$.

4.2.6. Effect of treatment on co-infected population

Figure 8 illustrates the effect of treatment rate $\theta$ on the number of HIV/AIDS and COVID-19 co-infectious population $C$. The result shows that increasing the value of $\theta$ decreases the number of HIV/AIDS and COVID-19 co-infectious individuals in the community. Consequently, increasing treatment intervention measure from 0.2 to 0.8 against HIV/AIDS and COVID-19 co-infection highly decreases the number of co- infectious population $C$.

Figure 8. Impact of $\theta$ on $C$.

4.2.7. Effect of ${\mathit{\beta }}_1$ on ${\mathcal{R}}_{\mathit{H}}$

Figure 9 depicts the transmission rate ${\beta }_1$ have the highest direct influence on ${\mathcal{R}}_H$. The result shows that ${\beta }_1$ increase the dissemination rate of HIV/AIDS in the community. Consequently, introducing intervention measures against HIV/AIDS dissemination is essential to reduce ${\beta }_1$ value to less than 0.75.

4.2.8. Effect of ${\mathit{\pi }}_2$ on ${\mathcal{R}}_{\mathit{H}}$

Figure 10 depicts portion of recruitment rate entered to HIV protected group ${\pi }_2$ has the highest indirect influence on ${\mathcal{R}}_H$. The result shows that ${\pi }_2$ decrease the dissemination rate of HIV/AIDS in the community. Consequently, introducing the portion of recruitment rate ${\pi }_2$ more than the value 0.608 makes the value of ${\mathcal{R}}_H$ below one.

4.2.9. Effect of $\mathit{\gamma }$ on ${\mathcal{R}}_{\mathit{H}}$

Figure 11 depicts HIV/AIDS rate $\gamma$ has fundamental indirect influence on ${\mathcal{R}}_H$. The result shows that as $\gamma$ increases the dissemination rate of HIV/AIDS in the community decreases. Consequently, introducing the treatment rate $\gamma$ more than the value 0.898 makes the value of ${\mathcal{R}}_H$ below one.

4.2.10. Effect of ${\mathit{\beta }}_2$ on ${\mathcal{R}}_{\mathit{C}}$

Figure 12 investigated the impact of COVID-19 transmission rate ${\beta }_2$ on the COVID-19 basic reproduction number ${\mathcal{R}}_{\text{C}}$. The figure reflects that whenever the value of ${\beta }_2$ increases implies the COVID-19 basic reproduction number is going up, and the value of ${\beta }_2<0.105$ leads that ${\mathcal{R}}_C<1.$ Therefore public health policymakers shall concentrate on minimizing the values of COVID-19 transmission rate ${\beta }_2$ to prevent and control COVID-19 expansion in the community. Biologically, it means that the COVID-19 infection increases whenever its transmission rate ${\beta }_2$ increases and the infection will be eradicated from the community whenever ${\beta }_2<\mathrm{0.105.}$

Figure 9. Impact of ${\beta }_1$ on ${\mathcal{R}}_H$.

Figure 10. Impact of ${\pi }_2$ on ${\mathcal{R}}_H$.

Figure 11. Impact of $\gamma$ on ${\mathcal{R}}_H$.

Figure 12. Impact of ${\beta }_2$ on ${\mathcal{R}}_C$.

4.2.11. Effect of ${\mathit{\pi }}_1$ on ${\mathcal{R}}_{\mathit{C}}$

Figure 13 depicts portion of recruitment rate entered to COVID-19 protected group ${\pi }_1$ has the highest indirect influence on${\mathcal{R}}_C$. The result shows that ${\pi }_1$ decrease the dissemination rate of COVID-19 in the community. Consequently, introducing the portion of recruitment rate ${\pi }_1$ more than the value 0.898 makes the value of ${\mathcal{R}}_C$ below one.

Figure 13. Impact of ${\pi }_1$ on ${\mathcal{R}}_C$

4.2.12. Effect of $\mathit{\kappa }$ on ${\mathcal{R}}_{\mathit{C}}$

Figure 14 depicts HIV/AIDS rate $\kappa$ has fundamental indirect influence on${\mathcal{R}}_C$. The result shows that as the COVID-19 recovery rate $\kappa$ increases the dissemination rate of HIV/AIDS in the community decreases. Consequently, introducing the treatment rate $\gamma$ more than the value of 0.758 makes the value of ${\mathcal{R}}_C$ below one.

Figure 14. Impact of $\kappa$ on ${\mathcal{R}}_C$.

5. Discussions and conclusions

Recently, HIV/AIDS and COVID-19 co-infection has been reported in different literature as one of the deadly infectious diseases. It has been a common public health burden especially in developing nations of the world. Some possible control measures against HIV/AIDS and COVID-19 single infections which are considered in this study include protections and treatments. This paper presents the HIV/AIDS and COVID-19 co-infection to investigate the impact of protection measures and treatment on the HIV/AIDS and COVID-19 co-infection transmission dynamics using ordinary differential equation approach. Rigorous analysis of the model shows that, the model solutions are shown to be non-negative and bounded, using next-generation matrix approach the basic reproduction numbers are computed, sufficient conditions for the local and global stabilities of the models’ equilibrium points are determined, and each of the disease-free and endemic equilibrium points of the models are both locally and GAS whenever the corresponding basic reproduction number is below one.

The sensitivity analysis of the model shows that the most important parameters that impact the basic reproduction numbers of the model (hence, disease burden) are the HIV transmission rate, COVID-19 transmission rate, proportion of protections from HIV and COVID-19, respectively, COVID-19 recovery rate of individuals. This suggests that increasing protection rates, treatment rates and recovery rate will decrease HIV/AIDS and COVID-19 co-infection transmission in the community.

The theoretical results are verified by the numerical simulations and it can be summarized as: the complete co-infection model (3) solutions are converging to the DFE point whenever its basic reproduction number ${\mathcal{R}}_{HC}=max\{{\mathcal{R}}_H,{\mathcal{R}}_C\}=max\{0.12,0.38\}=0.38<1$ that means there is no co-existence of the DFE point together with positive endemic equilibrium point. The complete model (3) solutions are converging to its endemic equilibrium point whenever ${\mathcal{R}}_{HC}=max\{{\mathcal{R}}_H,{\mathcal{R}}_C\}=max\{1.34,3.21\}=3.21>1.$Furthermore, the HIV/AIDS basic reproduction number${\mathcal{R}}_H$ is decreasing whenever the HIV protection rate ${\pi }_2$ increases, the basic reproduction numbers ${\mathcal{R}}_C$ and ${\mathcal{R}}_{HC}$ of the COVID-19 and HIV/AIDS and COVID-19 co-infection models are decreased whenever portion of protection ${\pi }_1$ of COVID-19 increases, the transmission rates increase the corresponding basic reproduction numbers also increase, increasing COVID-19 protection, HIV/AIDS protection, treatment and COVID-19 recovery rate leads to a decrease of the corresponding basic reproduction numbers. Increasing the HIV/AIDS transmission rate ${\beta }_1$ has a direct impact on the HIV/AIDS basic reproduction number and also has a direct impact on the number of co-infectious population. Similarly, the COVID-19 transmission rate ${\beta }_2$ has a direct impact on the COVID-19 basic reproduction number and also has a direct impact on the number of co-infectious population. Increasing the co-infection treatment rate $\theta$ decreases the number of co-infected individuals. Since the results show that HIV/AIDS transmission has an impact on COVID-19 infected individuals, COVID-19 transmission has a great impact on HIV/AIDS infected individuals, and applying COVID-19 protection rate, HIV/AIDS protection rate, and treatment rates has an effective role to minimize and also to eradicate the HIV/AIDS and COVID-19 co-infection in the community. The stakeholders shall concentrate on maximizing these intervention mechanisms to prevent and control the co-infection in the community.

The study is preliminary and potential researcher can extend it in different ways, such as including the optimal control approach, stochastic approach, fractional order derivative approach, environmental impacts, age and spatial structure, HIV infection stages, validate the model by applying appropriate real infection data.

Disclosure statement

The authors declares no conflict of interest in the study.

Data availability

Data used to support the findings of this study are included in the article.