# HIV/AIDS-Pneumonia Codynamics Model Analysis with Vaccination and Treatment

Authors: Shewafera Wondimagegnhu Teklu; Koya Purnachandra Rao
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

In this paper, we proposed and analyzed a realistic compartmental mathematical model on the spread and control of HIV/AIDS-pneumonia coepidemic incorporating pneumonia vaccination and treatment for both infections at each infection stage in a population. The model exhibits six equilibriums: HIV/AIDS only disease-free, pneumonia only disease-free, HIV/AIDS-pneumonia coepidemic disease-free, HIV/AIDS only endemic, pneumonia only endemic, and HIV/AIDS-pneumonia coepidemic endemic equilibriums. The HIV/AIDS only submodel has a globally asymptotically stable disease-free equilibrium ifR1<1. Using center manifold theory, we have verified that both the pneumonia only submodel and the HIV/AIDS-pneumonia coepidemic model undergo backward bifurcations whenever R2<1 and R3=maxR1,R2<1, respectively. Thus, for pneumonia infection and HIV/AIDS-pneumonia coinfection, the requirement of the basic reproduction numbers to be less than one, even though necessary, may not be sufficient to completely eliminate the disease. Our sensitivity analysis results demonstrate that the pneumonia disease transmission rate β2 and the HIV/AIDS transmission rate β1 play an important role to change the qualitative dynamics of HIV/AIDS and pneumonia coinfection. The pneumonia infection transmission rate β2 gives rises to the possibility of backward bifurcation for HIV/AIDS and pneumonia coinfection if R3=maxR1,R2<1, and hence, the existence of multiple endemic equilibria some of which are stable and others are unstable. Using standard data from different literatures, our results show that the complete HIV/AIDS and pneumonia coinfection model reproduction number is R3=maxR1,R2=max1.386,9.69=9.69 at β1=2 and β2=0.2 which shows that the disease spreads throughout the community. Finally, our numerical simulations show that pneumonia vaccination and treatment against disease have the effect of decreasing pneumonia and coepidemic disease expansion and reducing the progression rate of HIV infection to the AIDS stage.


---

## Body

## 1. Introduction

HIV/AIDS remains a major global health problem affecting approximately 70 million people worldwide causing significant morbidity and mortality (WHO, 2018) [1]. Over two-thirds of HIV/AIDS-infected population throughout the world is living in the sub-Saharan African Region [1–6]. AIDS is a common individual immune system disease caused by human immunodeficiency virus (HIV), i.e., RNA retrovirus which has developed into a global pandemic since the first patient was identified in 1981, making it one of the most destructive epidemics in history. HIV attacks human white blood cells and is transmitted through open sex, needle sharing, infected blood, and at childbirth [3, 6–9].Pneumonia is one of the leading airborne infectious diseases caused by microorganisms such as bacteria, viruses, or fungi. It has been the common cause of morbidity and mortality in adults, children under five years of age, and HIV-mediated immunosuppression worldwide, and it is a treatable respiratory lung infectious disease [5, 10–14]. In most prospective microbiology-based studies, bacteria especially Streptococcus bacteria are identified in 30-50% of pneumonia cases which are a leading cause of pneumonia in developing countries [13, 15–17]. However, over the past, our understanding about transmission of pneumonia is basically based on research from high-income western countries but the WHO, 2018 report assessed that from 9.5 million annual death worldwide, pneumonia and other respiratory infections cause about 2 million child deaths yearly in developing countries [14, 18].A coepidemic is the coexistence of two or more infections on a single individual at the population level [19]. HIV/AIDS-associated opportunistic infectious diseases are more common or more dangerous because of HIV immunosuppression [10].Mathematical and statistical models of infectious diseases have, historically, provided useful insight into the transmission dynamics and control of infectious diseases [14]. Mathematical models have been used to investigate the dynamics of single infections and coepidemics, and HIV/AIDS-pneumonia is among the diseases that infect a large number of individuals worldwide [10, 17, 20, 21].Babaei et al. [8] developed and analyzed a simple mathematical model for the interaction between drug addiction and the contagion of HIV/AIDS in Iranian prisons. They analyzed the stability of drug addiction and HIV/AIDS models separately with no medical treatment and investigated the impact of rehabilitating treatments on the control of HIV/AIDS spread in prisons, and finally, the reproduction numbers are compared in cases where there is no cure or some treatment methods are available. From their analysis, we have shown that their treatment methods for addiction withdrawal have a direct impact on the decrement and control of HIV/AIDS infection in prisons. Kizito et al. [13] constructed and discussed a mathematical model of treatment and vaccination impacts on pneumococcal pneumonia transmission dynamics. They found that, with treatment and vaccination combined, pneumonia can be eradicated; however, with treatment intervention alone, pneumonia remains in the population. Bakare and Nwozo [22] construct and analyzed a mathematical model for malaria–schistosomiasis coinfection. They have calculated the basic reproduction numbers and discussed the stability of equilibrium points of the model. They have shown the region where their model state variables become both mathematically and epidemiologically well-posed. They showed the model did not undergo backward bifurcation. Their mathematical modeling analysis result shows that intervention strategy suppresses the human-mosquito contact rate and human-snail contact rate to achieve malaria–schistosomiasis coepidemic free community. Shah et al. [3] formulated and analyzed a mathematical model for HIV/AIDS-TB coinfection considering HIV-infected population, and they found that medication plays a vital role in controlling the spread of the disease.Limited mathematical modeling research analysis has been conducted on HIV/AIDS-pneumonia coepidemics, for prevention and controlling of the disease transmission with controlling and prevention mechanisms; however, theoretical sources such as [10, 15, 20, 21] show the coexistence of HIV/AIDS-pneumonia. For our new research article, we reviewed only two published HIV/AIDS-pneumonia coepidemic model articles. Nthiiri et al. [5] constructed mathematical modeling on HIV/AIDS-pneumonia coinfection with maximum protection against single HIV/AIDS, and pneumonia infections were their basic concern. They did not consider maximum protection against coinfection. In their model analysis, we have found that when protection is maximum, the number of HIV/AIDS and pneumonia cases is going down. Teklu and Mekonnen [6] constructed a deterministic mathematical model and analyzed it both mathematically and numerically. Our model considered treatment at each infection stage of the coinfection model, and we found that when the treatment rate increases, the number of infectious population at each infection stage decreases. Our model did not consider pneumonia vaccination.We are motivated by the above studies especially the HIV/AIDS-pneumonia coexistence in the community; therefore, in this study, we considered the three center for disease control and prevention (CDC) stages of the HIV infection which are acute HIV infection, chronic HIV infection, and AIDS stage; we presented and analyzed a mathematical model describing the transmission dynamics of HIV/AIDS and pneumonia coinfection in a population where treatment for HIV/AIDS and both vaccination and treatment for pneumonia are available, respectively, in the community. Our model will be used to evaluate the effect of treatment at every infection stage of the HIV/AIDS only model, pneumonia only model, HIV/AIDS-pneumonia coinfection model, and effect of vaccination for pneumonia only model as control strategies for minimizing incidences of coinfections in the target population. The paper is organized as follows. The model is formulated in Section2 and is analyzed in Section 3. Sensitivity analysis and numerical simulation were carried out in Section 4. Finally, discussion, conclusion, and recommendation of the study are carried out in Sections 5 and 6, respectively.


## 2. Mathematical Model Formulation

### 2.1. Assumptions and Descriptions

According to CDC the three HIV/AIDS infection stages, we divide the human populationNt into twelve distinct classes as susceptible class to both HIV and pneumonia infections Y1t, pneumonia vaccine class Y2t, pneumonia-infected class Y3t, acute HIV-infected class Y4t, chronic HIV-infected class Y5t, AIDS patient class Y6t, acute HIV-pneumonia coepidemic class Y7t, chronic HIV-pneumonia coepidemic class Y8t, AIDS-pneumonia coepidemic class Y9t, pneumonia treatment class Y10t, HIV/AIDS treatment class Y11t entered from the three infection stages Y4t, Y5t, and Y6t, and HIV/AIDS-pneumonia coepidemic treatment class Y12t entered from Y7t, Y8t, and Y9t cases such that
(1)Nt=Y1t+Y2t+Y3t+Y4t+Y5t+Y6t+Y7t+Y8t+Y9t+Y10t+Y11t+Y12t.The susceptible class acquires HIV at the standard incidence rate given by(2)λHCt=β1NY4t+ρ1Y5t+ρ2Y7t+ρ3Y8t,where ρ3≥ρ2≥ρ1≥1 is the modification rate that increases infectivity and β1 is the HIV/AIDS contagion rate.The susceptible class acquires pneumonia at the mass action incidence rate(3)λPCt=β2Y3t+ω1Y7t+ω2Y8t+ω3Y9t,where ω3>ω2>ω1 is the modification rate that increases infectivity and β2 is the pneumonia contagion rate.To construct the complete coepidemic dynamical system, we have assumed the following:(i)
A fraction of the population has been vaccinated before the disease outbreak at the portion ofπ and 1−π fraction of population entered to the vulnerable class(ii)
The susceptible class is increased from the vaccinated class in which those individuals who are vaccinated but did not respond to vaccination with the waning rate ofτ and from pneumonia-treated class in which those individuals who lose their temporary immunity by the rateθ(iii)
Assume vaccination is not 100% effective, so vaccinated individuals also have a chance of being infected with proportionϵ of the serotype not covered by the vaccine where 0≤ϵ<1(iv)
Individuals in a given compartment are homogeneous(v)
Assume no HIV transmission fromY6t and Y9t classes due to their reduced daily activities(vi)
Individuals in each class are subject to natural mortality rated(vii)
The human population is variable(viii)
We assumed there is no dual-infection transmission simultaneously(ix)
Assume HIV has no vertical transmission and pneumonia is not naturally recovered(x)
No permanent immunity for pneumonia-infected individuals and become susceptible again after treatment


### 2.2. Schematic Diagram of the HIV/AIDS-Pneumonia Coepidemic Model

In this subsection using parameters in Table1, variable definitions in Table 2, and the model assumptions and descriptions given in (2.1), the schematic diagram for the transmission of HIV/AIDS-pneumonia coepidemic is given by the diagram.Table 1
Descriptions of model parameters.
ParameterInterpretationsdNatural mortality rateΛHuman recruitment rateδ1Development rate from acute HIV to chronic HIV infectionδ2Development rate from chronic HIV to AIDS stageϵThe proportion of the serotype not covered by the vaccineθImmunity loss rateψ1Alteration rate indicating acute HIV infection is more vulnerable to pneumoniaψ2Alteration rate indicating chronic HIV infection is more vulnerable to pneumoniaψ3Alteration rate indicating AIDS patient is more vulnerable to pneumoniaλHCHIV/AIDS standard incidence rateλPCPneumonia mass action incidence rateδ3Development rate from acute HIV-pneumonia to chronic HIV-pneumonia coepidemicsδ4Development rate from chronic HIV-pneumonia to AIDS-pneumonia coepidemicsdPPneumonia death ratedAAIDS death ratedAPAIDS-pneumonia death rateκPneumonia infection treatment rateκ1Acute HIV infection treatment rateτ1Vaccination waning rateκ2Chronic HIV infection treatment rateκ3AIDS patients treatment rateσ1Acute HIV-pneumonia coepidemic treatment rateσ2Chronic HIV-pneumonia coepidemic treatment rateσ3AIDS-pneumonia coepidemic treatment rateβ1Transmission rate of HIVβ2Transmission rate of pneumoniaTable 2
Definitions of variables.
VariablesDefinitionsY1Vulnerable to both HIV and pneumonia classY2Pneumonia-vaccinated classY3Pneumonia-infected classY4Acute HIV-infected classY5Chronic HIV-infected classY6AIDS patients classY7Acute HIV-pneumonia coepidemic classY8Chronic HIV-pneumonia coepidemic classY9AIDS-pneumonia coepidemic classY10HIV/AIDS treatment classY11Pneumonia treatment classY12Coepidemics treatment class


### 2.3. The HIV/AIDS-Pneumonia Coepidemic Dynamical System

From Figure1, the HIV/AIDS and pneumonia coinfection dynamical system is given by
(4)dY1dt=1−πΛ+τ1Y2+θY10−d+λHC+λPCY1,dY2dt=πΛ−ϵλPCY2−d+τ1+λHCY2,dY3dt=ϵλPCY2+λPCY1−νλHC+d+κ+dPY3,dY4dt=λHCY1+λHCY2−d+κ1+δ1+ψ1λPCY4,dY5dt=δ1Y4−d+κ2+δ2+ψ2λPCY5,dY6dt=δ2H2−d+κ3+dA+ψ3λPCY6,dY7dt=ψ1λPCY4+νλHCY3−d+dP+σ1+δ3Y7,dY8dt=ψ2λPCY5+δ3Y7−d+dP+σ2+δ4Y8,dY9dt=ψ3λPCY6+δ4Y8−d+dAP+σ3Y9,dY10dt=κY3−d+θY10,dY11dt=κ1Y4+κ2Y5+κ3Y6−dY11,dY12dt=σ1Y7+σ2Y8+σ3Y9−dY12.Figure 1
Flowchart of the HIV/AIDS-pneumonia coinfection model (4) where λHC and λPC are given in (2) and (3), respectively.With initial conditions,(5)Y10>0,Y20≥0,Y30≥0,Y40≥0,Y50≥0,Y60≥0,Y70≥0,Y80≥0,Y90≥0,Y10>0,Y11>0,Y120≥0.The sum of all the differential equations in (4) is
(6)dNdt=Λ−dN−dPY3+dAY6+dPY7+dPY8+dAPY9,


### 2.4. Positivity and Boundedness of the Solutions of the Model (4)

The model is mathematically analyzed by proving various theorems and algebraic computation 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 nonnegative well as the solutions of system (4) remain positive with positive initial conditions (5) in the bounded region
(7)Ω=Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,Y9,Y10,Y11,Y12∈ℝ+12,N≤Λd.Here, in order for the model (4) to be epidemiologically well-posed, it is important to show that each state variable defined in Table 2 with positive initial conditions (5) is nonnegative for all time t>0 in the bounded region given in (7).Theorem 1.
At the initial conditions (5), the solutions Y1t, Y2t,Y3t,Y4t,Y5t,Y6t,Y7t,Y8t, Y9t, Y10t,Y11t, and Y12t of system (4) are nonnegative for all time t>0.Proof.
AssumeY10>0, Y20>0, Y30>0, Y40>0, Y50>0, Y60>0, Y70>0, Y80>0,Y90>0, Y100>0, Y110, and Y120>0; then, for all t>0, we have to prove that Y1t>0, Y2t>0, Y3t>0,Y4t>0, Y5t>0, Y6t>0, Y7t>0, Y8t>0, Y9t>0, Y10t>0, Y11t>0, and Y12t>0.
Define:τ=supt>0:Y1t>0,Y2t>0,Y3t>0,Y4t>0,Y5t>0,Y6t>0,Y7t>0,Y8t>0,Y9t>0,Y10t>0,Y11t>0,andY12t>0.
From the continuity ofY1t,Y2t,Y3t,Y4t,Y5t,Y6t,Y7t,Y8t,Y9t,Y10t,Y11t, and Y12tt, we deduce that τ>0. If τ=+∞, then positivity holds. But, if 0<τ<+∞, Y1τ=0 or Y2τ=0 or Y3τ=0 or Y4τ=0 or Y5τ=0 or Y6τ=0 or Y7τ=0 or Y8τ=0 or Y9τ=0 or Y10τ=0 or Y11τ=0 or Y120=0.
Here, from the first equation of the model (4), we have dY1/dt=1−πΛ+θY10+τY2−d+λHC+λPCY1.
Using the method of integrating factor, we obtainedY1τ=M1Y10+M1∫0τexp∫d+λHct+λPctdt1−πΛ+θY10t+τY2tdt>0whereM1=exp−dτ+∫0τλHCw+λPCwdw>0,Y10>0, and from the definition ofτ, we see that Y2t>0, Y10t>0, and also the exponential function is always positive; then, the solution Y1τ>0; hence, Y1τ≠0. From the second equation of the model (4), we have dY2/dt=πΛ−d+τ1+ϵλPc+λHcY2 and we have got Y2τ=M1Y20+M1∫0τexp∫d+τ1+ϵλPct+λHctdtπΛdt>0, where M1=exp−dτ+τ1τ+∫0τλHcw+ϵλPcwdw>0,Y20>0, and also, the exponential function always is positive; then, the solutionY2τ>0; hence, Y2τ≠0. Similarly, all the remaining state variables Y3τ>0; hence,Y3τ≠0 and Y4τ>0; hence, Y4τ≠0 and Y5τ>0; hence, Y5τ≠0 and Y6τ>0; hence, Y6τ≠0 and Y7τ>0; hence, Y7τ≠0 and Y8τ>0; hence, Y8τ≠0 and Y9τ>0; hence, Y9τ≠0 and Y10τ>0; hence, Y10τ≠0 and Y11τ>0; hence, Y11τ≠0 and Y12τ>0; hence Y12τ≠0. Thus, based on the definition ofτ, it is not finite which means τ=+∞, and hence, all the solutions of system (2) are nonnegative.☐Theorem 2.
The regionΩ given by (7) is bounded in ℝ+12.Proof.
Using equation (6), since all the state variables are nonnegative by Theorem 1, in the absence of infections, we have got dN/dt≤Λ−dN. By applying standard comparison theorem, we have got ∫dN/Λ−dN≤∫dt and integrating both sides gives −1/dlnΛ−dN≤t+c where c is some constant, and after some steps of calculations, we have got 0≤Nt≤Λ/d which means all possible solutions of system (4) with positive initial conditions given in (5) enter in the bounded region (6).☐


## 2.1. Assumptions and Descriptions

According to CDC the three HIV/AIDS infection stages, we divide the human populationNt into twelve distinct classes as susceptible class to both HIV and pneumonia infections Y1t, pneumonia vaccine class Y2t, pneumonia-infected class Y3t, acute HIV-infected class Y4t, chronic HIV-infected class Y5t, AIDS patient class Y6t, acute HIV-pneumonia coepidemic class Y7t, chronic HIV-pneumonia coepidemic class Y8t, AIDS-pneumonia coepidemic class Y9t, pneumonia treatment class Y10t, HIV/AIDS treatment class Y11t entered from the three infection stages Y4t, Y5t, and Y6t, and HIV/AIDS-pneumonia coepidemic treatment class Y12t entered from Y7t, Y8t, and Y9t cases such that
(1)Nt=Y1t+Y2t+Y3t+Y4t+Y5t+Y6t+Y7t+Y8t+Y9t+Y10t+Y11t+Y12t.The susceptible class acquires HIV at the standard incidence rate given by(2)λHCt=β1NY4t+ρ1Y5t+ρ2Y7t+ρ3Y8t,where ρ3≥ρ2≥ρ1≥1 is the modification rate that increases infectivity and β1 is the HIV/AIDS contagion rate.The susceptible class acquires pneumonia at the mass action incidence rate(3)λPCt=β2Y3t+ω1Y7t+ω2Y8t+ω3Y9t,where ω3>ω2>ω1 is the modification rate that increases infectivity and β2 is the pneumonia contagion rate.To construct the complete coepidemic dynamical system, we have assumed the following:(i)
A fraction of the population has been vaccinated before the disease outbreak at the portion ofπ and 1−π fraction of population entered to the vulnerable class(ii)
The susceptible class is increased from the vaccinated class in which those individuals who are vaccinated but did not respond to vaccination with the waning rate ofτ and from pneumonia-treated class in which those individuals who lose their temporary immunity by the rateθ(iii)
Assume vaccination is not 100% effective, so vaccinated individuals also have a chance of being infected with proportionϵ of the serotype not covered by the vaccine where 0≤ϵ<1(iv)
Individuals in a given compartment are homogeneous(v)
Assume no HIV transmission fromY6t and Y9t classes due to their reduced daily activities(vi)
Individuals in each class are subject to natural mortality rated(vii)
The human population is variable(viii)
We assumed there is no dual-infection transmission simultaneously(ix)
Assume HIV has no vertical transmission and pneumonia is not naturally recovered(x)
No permanent immunity for pneumonia-infected individuals and become susceptible again after treatment


## 2.2. Schematic Diagram of the HIV/AIDS-Pneumonia Coepidemic Model

In this subsection using parameters in Table1, variable definitions in Table 2, and the model assumptions and descriptions given in (2.1), the schematic diagram for the transmission of HIV/AIDS-pneumonia coepidemic is given by the diagram.Table 1
Descriptions of model parameters.
ParameterInterpretationsdNatural mortality rateΛHuman recruitment rateδ1Development rate from acute HIV to chronic HIV infectionδ2Development rate from chronic HIV to AIDS stageϵThe proportion of the serotype not covered by the vaccineθImmunity loss rateψ1Alteration rate indicating acute HIV infection is more vulnerable to pneumoniaψ2Alteration rate indicating chronic HIV infection is more vulnerable to pneumoniaψ3Alteration rate indicating AIDS patient is more vulnerable to pneumoniaλHCHIV/AIDS standard incidence rateλPCPneumonia mass action incidence rateδ3Development rate from acute HIV-pneumonia to chronic HIV-pneumonia coepidemicsδ4Development rate from chronic HIV-pneumonia to AIDS-pneumonia coepidemicsdPPneumonia death ratedAAIDS death ratedAPAIDS-pneumonia death rateκPneumonia infection treatment rateκ1Acute HIV infection treatment rateτ1Vaccination waning rateκ2Chronic HIV infection treatment rateκ3AIDS patients treatment rateσ1Acute HIV-pneumonia coepidemic treatment rateσ2Chronic HIV-pneumonia coepidemic treatment rateσ3AIDS-pneumonia coepidemic treatment rateβ1Transmission rate of HIVβ2Transmission rate of pneumoniaTable 2
Definitions of variables.
VariablesDefinitionsY1Vulnerable to both HIV and pneumonia classY2Pneumonia-vaccinated classY3Pneumonia-infected classY4Acute HIV-infected classY5Chronic HIV-infected classY6AIDS patients classY7Acute HIV-pneumonia coepidemic classY8Chronic HIV-pneumonia coepidemic classY9AIDS-pneumonia coepidemic classY10HIV/AIDS treatment classY11Pneumonia treatment classY12Coepidemics treatment class


## 2.3. The HIV/AIDS-Pneumonia Coepidemic Dynamical System

From Figure1, the HIV/AIDS and pneumonia coinfection dynamical system is given by
(4)dY1dt=1−πΛ+τ1Y2+θY10−d+λHC+λPCY1,dY2dt=πΛ−ϵλPCY2−d+τ1+λHCY2,dY3dt=ϵλPCY2+λPCY1−νλHC+d+κ+dPY3,dY4dt=λHCY1+λHCY2−d+κ1+δ1+ψ1λPCY4,dY5dt=δ1Y4−d+κ2+δ2+ψ2λPCY5,dY6dt=δ2H2−d+κ3+dA+ψ3λPCY6,dY7dt=ψ1λPCY4+νλHCY3−d+dP+σ1+δ3Y7,dY8dt=ψ2λPCY5+δ3Y7−d+dP+σ2+δ4Y8,dY9dt=ψ3λPCY6+δ4Y8−d+dAP+σ3Y9,dY10dt=κY3−d+θY10,dY11dt=κ1Y4+κ2Y5+κ3Y6−dY11,dY12dt=σ1Y7+σ2Y8+σ3Y9−dY12.Figure 1
Flowchart of the HIV/AIDS-pneumonia coinfection model (4) where λHC and λPC are given in (2) and (3), respectively.With initial conditions,(5)Y10>0,Y20≥0,Y30≥0,Y40≥0,Y50≥0,Y60≥0,Y70≥0,Y80≥0,Y90≥0,Y10>0,Y11>0,Y120≥0.The sum of all the differential equations in (4) is
(6)dNdt=Λ−dN−dPY3+dAY6+dPY7+dPY8+dAPY9,


## 2.4. Positivity and Boundedness of the Solutions of the Model (4)

The model is mathematically analyzed by proving various theorems and algebraic computation 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 nonnegative well as the solutions of system (4) remain positive with positive initial conditions (5) in the bounded region
(7)Ω=Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,Y9,Y10,Y11,Y12∈ℝ+12,N≤Λd.Here, in order for the model (4) to be epidemiologically well-posed, it is important to show that each state variable defined in Table 2 with positive initial conditions (5) is nonnegative for all time t>0 in the bounded region given in (7).Theorem 1.
At the initial conditions (5), the solutions Y1t, Y2t,Y3t,Y4t,Y5t,Y6t,Y7t,Y8t, Y9t, Y10t,Y11t, and Y12t of system (4) are nonnegative for all time t>0.Proof.
AssumeY10>0, Y20>0, Y30>0, Y40>0, Y50>0, Y60>0, Y70>0, Y80>0,Y90>0, Y100>0, Y110, and Y120>0; then, for all t>0, we have to prove that Y1t>0, Y2t>0, Y3t>0,Y4t>0, Y5t>0, Y6t>0, Y7t>0, Y8t>0, Y9t>0, Y10t>0, Y11t>0, and Y12t>0.
Define:τ=supt>0:Y1t>0,Y2t>0,Y3t>0,Y4t>0,Y5t>0,Y6t>0,Y7t>0,Y8t>0,Y9t>0,Y10t>0,Y11t>0,andY12t>0.
From the continuity ofY1t,Y2t,Y3t,Y4t,Y5t,Y6t,Y7t,Y8t,Y9t,Y10t,Y11t, and Y12tt, we deduce that τ>0. If τ=+∞, then positivity holds. But, if 0<τ<+∞, Y1τ=0 or Y2τ=0 or Y3τ=0 or Y4τ=0 or Y5τ=0 or Y6τ=0 or Y7τ=0 or Y8τ=0 or Y9τ=0 or Y10τ=0 or Y11τ=0 or Y120=0.
Here, from the first equation of the model (4), we have dY1/dt=1−πΛ+θY10+τY2−d+λHC+λPCY1.
Using the method of integrating factor, we obtainedY1τ=M1Y10+M1∫0τexp∫d+λHct+λPctdt1−πΛ+θY10t+τY2tdt>0whereM1=exp−dτ+∫0τλHCw+λPCwdw>0,Y10>0, and from the definition ofτ, we see that Y2t>0, Y10t>0, and also the exponential function is always positive; then, the solution Y1τ>0; hence, Y1τ≠0. From the second equation of the model (4), we have dY2/dt=πΛ−d+τ1+ϵλPc+λHcY2 and we have got Y2τ=M1Y20+M1∫0τexp∫d+τ1+ϵλPct+λHctdtπΛdt>0, where M1=exp−dτ+τ1τ+∫0τλHcw+ϵλPcwdw>0,Y20>0, and also, the exponential function always is positive; then, the solutionY2τ>0; hence, Y2τ≠0. Similarly, all the remaining state variables Y3τ>0; hence,Y3τ≠0 and Y4τ>0; hence, Y4τ≠0 and Y5τ>0; hence, Y5τ≠0 and Y6τ>0; hence, Y6τ≠0 and Y7τ>0; hence, Y7τ≠0 and Y8τ>0; hence, Y8τ≠0 and Y9τ>0; hence, Y9τ≠0 and Y10τ>0; hence, Y10τ≠0 and Y11τ>0; hence, Y11τ≠0 and Y12τ>0; hence Y12τ≠0. Thus, based on the definition ofτ, it is not finite which means τ=+∞, and hence, all the solutions of system (2) are nonnegative.☐Theorem 2.
The regionΩ given by (7) is bounded in ℝ+12.Proof.
Using equation (6), since all the state variables are nonnegative by Theorem 1, in the absence of infections, we have got dN/dt≤Λ−dN. By applying standard comparison theorem, we have got ∫dN/Λ−dN≤∫dt and integrating both sides gives −1/dlnΛ−dN≤t+c where c is some constant, and after some steps of calculations, we have got 0≤Nt≤Λ/d which means all possible solutions of system (4) with positive initial conditions given in (5) enter in the bounded region (6).☐


## 3. The Mathematical Model Analysis

Before we analyze the HIV/AIDS-pneumonia coinfection model (4), we need to gain some background about the HIV/AIDS-only submodel and pneumonia-only submodel transmission dynamics.


### 3.1. HIV/AIDS Submodel Analysis

We have the HIV/AIDS submodel of (4) whenY2=Y3=Y7=Y8=Y9=Y10=Y12=0 which is given by
(8)dY1dt=Λ−d+λHY1,dY4dt=λHY1−d+κ1+δ1Y4,dY5dt=δ1Y4−d+κ2+δ2Y5,dY6dt=δ2Y5−d+κ3+dAY6,dY11dt=κ1Y4+κ2Y5+κ3Y6−dY11,where the total population isN1t=Y1t+Y4t+Y5t+Y6t+Y11t and the HIV/AIDS single infection force of infection is given by λH=β1/N1Y4+ρ1Y5 with initial conditions Y10>0, Y40≥0, Y50≥0, Y60≥0, andY110≥0.Here, the detailed HIV/AIDS submodel model analysis is given in [6].


### 3.2. Pneumonia Submodel Analysis

From model (4), we have got the pneumonia submodel at Y4=Y5=Y6=Y7=Y8=Y9=Y11=Y12=0, which is given by
(9)dY1dt=1−πΛ+τ1Y2+θY10−d+λPY1,dY2dt=πΛ−ϵλPY2−d+τ1Y2,dY3dt=ϵλPY2+λ4Y1−d+κ+dPY3,dY10dt=κY3−d+θY10.With initial conditions,Y10>0, Y20≥0, Y30≥0, Y100≥0, total population N2t=Y1t+Y2t+Y3t+Y10t, and pneumonia force of infection λP=β2Y3t.In the regionΩ2=Y1,Y2,Y3,Y10∈ℝ4+,N2≤Λ/d, it is easy to show that the set Ω2 is positively invariant and a global attractor of all positive solutions of submodel (9). Hence, it is sufficient to consider the dynamics of model (9) in Ω2 as epidemiologically and mathematically well-posed.


#### 3.2.1. Disease-Free Equilibrium Point of the Pneumonia Submodel

The disease-free equilibrium point of the pneumonia submodel is obtained by making the right-hand side of the system (15) as zero and setting the infectious class and treatment class to zero asY3=Y10=0 we have gotY10=Λd+τ1−Λπd/dd+τ1 and Y20=Λπ/d+τ1 such that E20=Y10,Y20,Y30,Y100=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0.


#### 3.2.2. The Effective Reproduction Number of the Pneumonia Submodel

The effective reproduction number measures the average number of new infections generated by a typically infectious individual in a community when some strategies are in place, like vaccination or treatment. We calculate the effective reproduction numberR2 using the van den Driesch and Warmouth next-generation matrix approach [23]. The Effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE20/∂xj∂νiE20/∂xj−1 where Fi is the rate of appearance of new infection in compartment i, νi is the transfer of infections from one compartment i to another, and E20 is the disease-free equilibrium point. Then, after a long calculation, we have got
(10)F=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1000,V=d+κ+dP0−κd+θ.Then, using Mathematica, we have got(11)V−1=1d+κ+dP0γθ+dd+κ+dP1θ+d,FV−1=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP000.The characteristic equation of the matrixFV−1is
(12)β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP−λ000−λ=0.Then, the spectral radius (effective reproduction numberR2) of FV−1 of the pneumonia submodel (9) is R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP. Here, R2 is the effective reproduction number for pneumonia infection.


#### 3.2.3. Local and Global Stability of the Disease-Free Equilibrium Point

Theorem 3.
The disease-free equilibrium point (DFE)E20 of the pneumonia submodel (9) is locally asymptotically stable if R2<1, otherwise unstable.Proof.
The local stability of the disease-free equilibrium of the system (9) can be studied from its Jacobian matrix at the disease-free equilibrium point E20=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0 and Routh Hurwitz stability criteria. The Jacobian matrix of a dynamical system (9) at the disease-free equilibrium point is given by
(13)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπd+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1−d+κ+dP000κ−d+θ.
Then, the characteristic equation of the above Jacobian matrix is given by(14)−d−λτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−λ−ϵβ2Λπd+τ1000M−λ000κ−d+θ−λ=0,where M=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1−d+κ+dP.
After some steps, we have gotλ1=−d<0 orλ2=−d+τ1<0 orλ3=d+κ+dPR2−1<0 ifR2<1 or λ4=−d+θ<0. Therefore, since all the eigenvalues of the characteristics polynomial of the system (9) are negative if R2<1, the disease-free equilibrium point of the pneumonia submodel is locally asymptotically stable.☐


#### 3.2.4. Existence of EEP for the Pneumonia Submodel

Let an arbitrary endemic equilibrium point of pneumonia-only dynamical system (9) be denoted by E2∗=Y1∗,Y2∗,Y3∗,Y10∗. Moreover, let λP∗=β2Y3∗ be the associated pneumonia mass action incidence rate (“force of infection”) at an equilibrium point. To find conditions for the existence of an arbitrary equilibrium point(s) for which pneumonia infection is endemic in the population, the equations of model (9) are solved in terms of the force of infection rate λP∗=β2Y3∗ at an endemic equilibrium point. Setting the right-hand sides of the equations of model (9) to zero and we have got Y2∗=πΛ/ϵλP∗+d+τ1, Y10∗=κY3∗/d+θ and substitute Y2∗ and Y10∗ in to Y1∗, we obtain Y1∗=1−πΛ+τ1Y2∗+θTP∗/d+λP∗=1−πΛ/d+λP∗d+λP∗+πΛτ1/ϵλP∗+d+τ1d+λP∗+θγY3∗/d+θd+λP∗ and substitute Y2∗ and Y1∗ inY3∗, we obtain
(15)Y3∗=πΛϵλP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+1−πΛλP∗d+κ+dPd+θd+λP∗d+κ+dPd+λP∗d+κ+dPd+θd+λP∗−θκλP∗+πΛτ1λP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+λP∗d+γ+dPd+θd+λP∗−θκλP∗.Finally, substituteY3∗ in to pneumonia submodel (9) force of infection λP∗=β2Y3∗ as
(16)λP∗=β2Y3∗=β2πΛϵλP∗d+θd+λP∗ϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+β21−πΛλP∗d+θd+κ+dPd+θd+λP∗−θκλP∗+β2πΛτ1λP∗d+θϵλP∗+κ+τ1d+κ+dPd+θd+λP∗−θκλP∗,and lettingm1=d+κ+dP,m2=d+τ1, and m3=d+θ, we have got a2λP∗2+a1λP∗+a0=0 where a2=m1m3ϵ−θκϵ>0, a1=m1m3dϵ+m1m2m3−m2θκ−β2Λm3ϵ, and a0=m1m2m3μ1−R2>0 if R2<1.Here, the nonzero equilibrium(s) of the model (9) satisfies fλP∗=a2λP∗2+a1λP∗+a0=0 so that the quadratic equation can be analyzed for the possibility of multiple equilibriums. It is worth noting that the coefficienta2is always positive and a0 is positive (negative) if RP is less than (greater than) unity, respectively. Hence, we have established the following result.Theorem 4.
The pneumonia submodel (9) has the following:
(i)
Exactly one unique endemic equilibrium ifa0<0 (i.e.,R2>1)(ii)
Exactly one unique endemic equilibrium ifa1<0, and a0=0 or a12−4a2a0=0(iii)
Exactly two endemic equilibriums ifa0>0(i.e., R2<1), a1<0, anda12−4a2a0>0(iv)
No endemic equilibrium otherwise
Here, item (iii) shows the happening of the backward bifurcation in pneumonia submodel (9), i.e., the locally asymptotically stable disease-free equilibrium point coexists with a locally asymptotically stable endemic equilibrium point ifR2<1; examples of the existence of backward bifurcation phenomenon in mathematical epidemiological models, and the causes, can be seen in [2, 9, 22, 24–26]. The epidemiological consequence is that the classical epidemiological requirement of having the reproduction number R2to be less than one, even though necessary, is not sufficient for the effective control of the disease. The existence of the backward bifurcation phenomenon in submodel (9) is now explored.


#### 3.2.5. Bifurcation Analysis

It is instructive to explore the possibility of backward bifurcation in model (9).Theorem 5.
Model (9) exhibits backward bifurcation at R2=1 whenever the inequality D1>D2holds.Here, we apply the center manifold theory in [27]; however, to apply this theory, the following simplification and change of variables are made.LetY1=x1,Y2=x2,Y3=x3, and Y10=x4 such thatN2=x1+x2+x3+x4. Furthermore, by using vector notationX=x1,x2,x3,x4T, pneumonia submodel (9) can be written in the form dX/dt=FX withF=f1,f2,f3,f4T, as follows:
(17)dx1dt=f1=1−πΛ+τ1x2+θx4−d+λPx1,dx2dt=f2=πΛ−ϵλP+d+τ1x2,dx3dt=f3=ϵλPx2+λPx1−d+κ+dPx3,dx4dt=f4=κx3−d+θx4,with λP=β2x3, then the method entails evaluating the Jacobian of system (17) at the DFE pointE20, denoted by JE20, and this gives us
(18)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπμ+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπμdd+τ1−d+κ+dP000κ−d+θ.Consider, next, the case whenRP=1. Suppose, further, that β2=β∗ is chosen as a bifurcation parameter.Solving forβ2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1 and we have gotβ2=β∗=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd and
(19)Jβ∗=−dτ1−β∗Λd+τ1+β∗Λπddd+τ1θ0−d+τ1−ϵβ∗Λπd+τ1000β∗ϵΛπd+β∗Λd+τ1−β∗Λπddd+τ1−d+κ+dP000γ−μ+θ.After some steps of the calculation, we have got the eigenvalues ofJβ∗ as λ1=−d or λ2=−d+τ1 or λ3=0 or λ4=−d+θ.It follows that the JacobianJE20 of (17) at the DFE, with β2=β∗, denoted by Jβ∗, has a simple zero eigenvalue with all the remaining eigenvalues having a negative real part. Hence, the center manifold theory [27] can be used to analyze the dynamics of model (9). In particular, Theorem 2 of Castillo-Chavez and Song [28] will be used to show that model (9) undergoes backward bifurcation atR2=1Eigenvectors ofJβ∗: for the caseR2=1, it can be shown that the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) has a right eigenvectors associated with the zero eigenvalue given by u=u1,u2,u3,u4T with values
(20)u1=−ϵβ∗Λπdτ1d+θ−β∗Λd+τ12d+θ+β∗Λπdd+τ1d+θ+θκdd+τ12d2d+τ12u3,u2=−ϵβ∗Λπd+τ12u3,u3=u3>0,u4=κd+θu3.Similarly, the left eigenvector associated with the zero eigenvalues atβ2=β∗ given by v=v1,v2,v3,v4T arev1=v2=v4=0,v3=v3>0.After long calculations, the bifurcation coefficientsa and b are obtained as a=D1−D2 where D1=β∗Λπdd+τ1d+θ+θκdd+τ12/d2d+τ12, and D2=ϵβ∗Λπdτ1d+θ+β∗Λd+τ12d+θ/d2d+τ12+ϵϵβ∗Λπ/d+τ12.Thus, the bifurcation coefficienta is positive ifD1>D2. Furthermore, b=v3u2u3Λd+τ1−Λπd/dd+τ1>0.Hence, from in Castillo-Chavez and Song [28], model (9) exhibits a backward bifurcation at R2=1 wheneverD1>D2.


### 3.3. Analysis of the Full HIV/AIDS-Pneumonia Coinfection

Having analyzed the dynamics of the two submodels, that is, HIV/AIDS submodel (8) and the pneumonia submodel (9), the complete HIV/AIDS-pneumonia coinfection model (4) is now considered (the analysis is done in the positively invariant regionΩ given in (7)).


#### 3.3.1. Disease-Free Equilibrium Point of the HIV/AIDS-Pneumonia Coinfection

The disease-free equilibrium point of model (4) is obtained by setting all the infectious classes and treatment classes to zero such that Y3=Y4=Y5=Y6=Y7=Y8=Y9=Y10=Y11=Y12=0and hence E30== Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.


#### 3.3.2. Effective Reproduction Number of the HIV/AIDS-Pneumonia Coinfection

The basic reproduction number, denoted byR0, is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [6, 23, 28]. For simple classical models if R0<1, then it means that on average, an infected individual infects less than one susceptible over the course of its infectious period and the disease cannot grow. If however, R0>1, then an infected individual infects more than one susceptible over the course of its infectious period and the disease will persist. For more complicated models with several infected compartments, this simple heuristic definition of R0 is insufficient [23]. Due to its importance, researchers have sought to find ways of determiningR0. Two important concepts in modeling outbreaks of infectious diseases are the basic reproduction number, universally denoted by R0 and the generation time (the average time from symptom onset in a primary case to symptom onset in a secondary case), which jointly determine the likelihood and speed of epidemic outbreaks [29].Here, we calculated the HIV/AIDS-pneumonia coinfection effective reproduction numberR3 of model (4) using the van den Driesch and Warmouth next-generation matrix approach [23]. The effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE30/∂xj∂νiE30/∂xj−1 whereFiis the rate of appearance of new infection in compartmenti, νiis the transfer of infections from one compartment ito another, and E30 is the disease-free equilibrium point E30=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.After long detailed calculations, the transition matrixF is given by
(21)F=A0000000000β1β1ρ1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,and the transmission matrix V is given by
(22)V=D10000000000D2000000000−δ1D3000000000−δ2D40000000000D5000000000−δ3D6000000000−δ4D7000−κ000000θ000−κ1−κ2−κ30000d00000−σ1−σ2−σ300d,whereD1=d+κ+dP,D2=d+κ1+δ1,D3=d+κ2+δ2,D4=d+κ3+dA,D5=d+dP+σ1+δ3,D6=d+dP+σ2+δ4,andD7=d+dAP+ε3.Then, by using Mathematica, we have got(23)FV−1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1D10000000000β1D2+δ1β1ρ1D2D3β1ρ1D3000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.The characteristic equation of the matrixFV−1 is given by
(24)A1−λ0000000000B−λβ1ρ1D30000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ=0,whereA1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1/D1,B=β1/D2+δ1β1ρ1/D2D3; then, the eigenvalues of FV−1 are λ1=β2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1orλ2=β1/D2+δ1β1ρ1/D2D3 or λ3=λ4=λ5=λ6=λ7=λ8=λ9=λ10=0.Thus, the effective reproduction number of the HIV/AIDS-pneumonia coinfection dynamical system (4) is the dominant eigenvalue of the matrixFV−1 which is given byR3=maxλ1,λ2=maxβ2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1,β1/D2+δ1β1ρ1/D2D3. Here, R2=β2ϵπdΛ+β2Λd+τ1−β2πΛd/d+κ+dPd+τ1 is the effective reproduction number for pneumonia-only infected individual and R1=β1/d+κ1+δ1+β1ρ1δ1/d+κ1+δ1d+κ2+κ2 is the basic reproduction for HIV/AIDS-only infected individual.Here,R1 represent the basic reproduction number for HIV/AIDS submodel, R2 and R3 are the effective reproduction numbers for the pneumonia submodel and HIV/AIDS-pneumonia coinfection model, respectively. The following three outcomes are possible: (i) for R1<1, the HIV/AIDS submodel disease-free steady state E1 is globally stable in the region Ω1, and HIV is not spreading in the community; (ii) for R2<1, then E2 is not globally stable in the region Ω2, and pneumonia may spread through the community; (iii) for R3<1, the steady state E3is not globally stable in the region Ω, and HIV/AIDS-pneumonia coinfection may spread through the community.Note that none of the parameters corresponding to coinfection treatment (i.e.,σ1orσ2orσ3) are present in the expression forR3, indicating no impact of treating coinfected population onR3.


#### 3.3.3. Locally Asymptotically Stability of the Disease-Free Equilibrium (DFE)

Theorem 6.
The disease-free equilibrium of model (4) above is locally asymptotically stable if R3<1, otherwise unstable.Proof.
The Jacobian matrixJE30 of model (4) at E30 is given by
(25)JE30=−dτ1−β2Y10−β1N0Y10−β1N0ρ1Y10000θ0000−d−τ1−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1000000000000Z2β1N0ρ1Y100000000000δ1Z300000000000δ2Z4000000000000Z500000000000δ3Z600000000000δ4Z700000κ000000Z800000κ1κ2κ30000−d0000000σ1σ2σ300−d,
whereZ1=β2ϵY20+β2Y10−d+κ+dP,Z2=β1/N0Y10−d+κ1+δ1, Z3=−d+κ2+δ2, Z4=−d+κ3+dA,Z5=−d+dP+σ1+δ3,Z6=−d+dP+σ2+δ4,Z7=−d+dAP+σ3, andZ8=−d+θ.
Then, the characteristic equation of the Jacobian matrixJE30 is given by
(26)−d−λτ1−β2Y10−β1N0Y10−β1N0ρ1Y100000θ000−d−τ1−λ−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1−λ000000000000Z2−λβ1N0ρ1Y100000000000δ1Z3−λ00000000000δ2Z4−λ000000000000Z5−λ00000000000δ3Z6−λ00000000000δ4Z7−λ00000κ000000Z8−λ00000κ1κ2κ30000−d−λ0000000σ1σ2σ300−d−λ=0.
After detailed calculations, we have got that
λ1=λ2=λ3=−d<0 or λ4=d+κ+dPR2−1<0 if R2<1 or λ5=−d+τ1<0 or λ6=−d+κ3+dA<0 or λ7=−d+dP+σ1+δ3<0 or λ8=−d+dP+σ2+δ4<0 or λ9=−d+dAP+σ3<0 or λ10=−d+θ<0or a2λ2+a1λ+a0=0 for a2=1,a1=d+κ2+δ2+d+κ1+δ11−Y10/N0RY4>0 if RY4<1 anda0=d+κ2+δ2d+κ1+δ11−Y10/N0RY5>0 if RY5<1.
Then, by applying Routh-Hurwitz stability criteria sincea2=1>0, a1>0, and a0>0, all the eigenvalues of the Jacobian matrix are negative if R1<1 andR2<1, i.e.,R3=maxR1,R2<1. Thus, the disease-free equilibrium point (DFE) of HIV/AIDS-pneumonia coinfection model (4) is locally asymptotically stable if
(27)R3=maxR1,R2<1.☐


#### 3.3.4. Existence of Endemic Equilibrium Point (EEP) for the Full Model

The endemic equilibrium point (EEP) of full model (4) is denoted by E3∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ which occurs when the disease persists in the community. From the analysis of HIV/AIDS-only submodel (8) and the pneumonia-only submodel from (9), we have shown that there is no endemic equilibrium point if R1<1and there is/are an endemic equilibrium point(s) if R2<1 implies that there is/are endemic equilibrium point(s) if R3<1 for the coinfection model and hence there is a bifurcation point for the full model. The endemic equilibrium of system (4) is obtained as
(28)Y1∗=1−πΛ+τ1Y2∗+θY10∗d+λHC∗+λPC∗,Y2∗=πΛϵλPC∗+d+τ1+λHC∗,Y3∗=ϵλPC∗Y2∗+λPC∗Y1∗νλHC∗+d+κ+dP,Y4∗=λHC∗Y1∗d+κ1+δ1+ψ1λPC∗,Y5∗=δ1H1∗d+κ2+κ2+ψ2λPC∗,Y6∗=δ2Y5∗d+κ3+dA+δ3λPC∗,Y7∗=ψ1λPC∗Y4∗+νλHC∗Y3∗d+dP+σ1+δ3,Y8∗=δ2λPC∗Y5∗+δ3Y7∗d+dP+σ2+δ4,Y9∗=ψ3λPC∗Y6∗+δ4Y8∗d+dAP+κ3,Y10∗=κY3∗d+θ,Y11∗=κ1Y4∗+κ2Y5∗+κ3Y6∗d,andY12∗=κ1Y7∗+κ2Y8∗+κ3Y9∗d.Summary of endemic equilibrium points: the explicit computation of the endemic equilibrium of coinfection model (4) given in (28) in terms of model parameters is difficult analytically; however, model (4) endemic equilibriums correspond to the following:
(1)
E4∗=Y1∗,0,Y4∗,Y5∗,Y6∗,0,0,0,0,0,Y11∗,0,ifR1>1 is the pneumonia free (HIV) endemic equilibrium point. The analysis of the equilibrium E1∗ is similar to the endemic equilibrium E1∗ in model (7)(2)
E5∗=Y1∗,Y2∗,Y3∗,0,0,0,0,0,0,Y10∗,0,0, if R2>1 is the HIV/AIDS free (pneumonia) endemic equilibrium point. The analysis of the equilibrium E5∗ is similar to the endemic equilibrium E2∗ in equation (9)(3)
E6∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ is the HIV/AIDS-pneumonia coinfection endemic equilibrium point. It exists when each component of E6∗ in equation (28) is positive and summarizes the existence of the endemic equilibrium points in the following theorem


#### 3.3.5. Bifurcation Analysis

The threshold quantityR3=maxR1,R2 is the effective reproduction number of the system (4) where R1 and R2 are defined as above.Theorem 7.
Model (4) exhibits the phenomenon of backward bifurcation at R3=1 whenever the inequality G1>G2 holds.The phenomenon of backward bifurcation can be proved with the concept of the center manifold theory [10, 27] on coepidemic model (4). To apply this theory, the following simplification and change of variables are made.LetY1=x1, Y2=x2, Y3=x3,Y4=x4,Y5=x5, Y6=x6,Y7=x7,Y8=x8,Y9=x9,Y10=x10,Y11=x11,and Y12=x12 so that N=x1+x2+x3+x4+x5,+x6+x7,+x8+x9+x10+x11+x12.Further, by using vector notationX=x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12T, complete model (4) can be written in the form dX/dt=FX with F=f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12T, as follows:
(29)dx1dt=f1=1−πΛ+τ1x2+θx10−d+λHC+λPCx1,dx2dt=f2=πΛ−ϵλPCx2−d+τ1+λHCx2,dx3dt=f3=ϵλPCx2+λPCx1−νλHC+d+κ+dPx3,dx4dt=f4=λHCx1−d+κ1+δ1+ψλPCx4,dx5dt=f5=δ1x4−d+κ2+κ2+ψ2λPCx5,dx6dt=f6=δ2x5−d+κ3+dA+ψ3λPCx6,dx7dt=f7=ψ1λPCx4+νλHCx3−d+dP+σ1+δ3x7,dx8dt=f8=ψ2λPCx5+δ3x7−d+dP+σ2+δ4x8,dx9dt=f9=ψ3λPCx6+δ4x8−d+dAP+κ3x9,dx10dt=f10=κx3−d+θx10,dx11dt=f11=κ1x4+κ2x5+κ3x6−dx11,dx12dt=f12=κ1x7+κ2x8+κ3x9−dx12,with λHC=β1/Nx4+ρ1x5+ρ2x7+ρ2x8whereρ3≥ρ2≥ρ1≥1 andλPC=β2x3+ω1x7+ω2x8+ω3x9 whereω3≥ω2≥ω1≥1; then, the method entails evaluating the Jacobian of system (29) at the DFEE30, denoted by JE30, and this gives us
(30)JE30=−dτ1F1F2F30F4F5F6θ000F7F8F9F100F11F12F1300000F14000F15F16F17000000F18F190F20F210000000δ1F2200000000000δ2F23000000000000F2400000000000δ3F2500000000000d4F2600000κ000000F2700000κ1κ2κ30000−d0000000κ1κ2κ300−d,whereF1=−β2Y10,F2=−β1Y10/Y10+Y20,F3=−β1ρ1Y10/Y10+Y20, F4=−β1ρ2Y10/Y10+Y20−β2ω1Y10,F5=−β1ρ3Y10/Y10+Y20−β2ω2Y10,F6=−β2ω3Y10,F7=−d+τ1,F8=−ϵβ2Y20,F9=−β1Y20/Y10+Y20,F10=−β1ρ1Y20/Y10+Y20,F11=−ϵβ2ω1VP0−β1ρ2YP0/Y10+YP0,F12=−ϵβ2ω2VP0−β1ρ3Y20/Y10+Y20,F13=−ϵβ2ω3Y20,F14=ϵβ2Y20+β2Y10−d+κ+dP,F15=ϵβ2ω1Y20+β2ω1Y10,F16=ϵβ2ω2Y20+β2ω2Y10,F17=ϵβ2ω3Y20+β2ω3Y10, F18=β1Y10/Y10+Y20−d+κ1+δ1,F19=β1ρ1Y10/Y10+Y20,F20=β1ρ2Y10/Y10+Y20,F21=β1ρ3Y10/Y10+Y20,F22=−d+κ2+δ2,F23=−d+κ3+dA,F24=−d+dP+σ1+δ3,F25=−d+dP+σ2+δ4,F26=−d+dAP+σ3, and F27=−d+θ.Without loss of generality, consider the case whenR2>R1, and R3=1, so that R2=1. Furthermore, let β2=β∗ is chosen as a bifurcation parameter. Solving for β2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1, we have got the value β∗=β2=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd.After solving the JacobianJE30 of the system (29) at the DFE, with β2=β∗, we obtained the eigenvalues as λ1=−d<0 or λ2=−d+τ1<0 or λ3=0 or λ4=−d<0 orλ5=−d<0 orλ6=−d+θ<0 or λ7=−d+dAP+κ3<0 orλ8=−d+κ3+dA<0 or λ9=−d+dP+σ1+δ3<0 or λ10=−d+dP+σ1+δ4<0or
(31)a2λ2+a1λ+a0=0,wherea2=1>0, a1=d+κ1+δ11−Y10/Y10+Y40RY4+d+κ2+δ2>0 if RY4<1, and a0=a0=d+κ1+δ1d+κ2+δ21−Y10/Y10+Y20R1>0 if R1<1.Equation (31) has/have no positive root/s whenever R1<1, and hence, both eigenvalues are negative. It follows that the Jacobian JE30 of (29) at the DFE, withβ2=β∗, denoted by Jβ∗, has a simple zero eigenvalue (with all other eigenvalues having negative real part). Hence, the center manifold theory [27] can be used to analyze the dynamics of model (4). In particular, the Castillo-Chavez and Song theorem [28] will be used to show that model (4) undergoes backward bifurcation at RP=1.Eigenvectors ofJβ∗: for the case when RP=1, the right eigenvectors of the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) associated with the zero eigenvalue given by u=u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12T areu1=−τ1F8F27+F1F7F27−θγF7/dF7F27u3,u2=−F8/F7u3,u3=u3>0,u10=−κ/F27u3, and u4=u5=u6=u7=u8=u9=u11=u12=0.The left eigenvectors associated with the zero eigenvalue atβ2=β2∗ satisfyingv.w=1 given by v=v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12are v1=v2=v4=v5=v6=v10=v11=v12=0, v3=v3>0, v7=−δ3δ4F17+δ3F16F26−F15F25F26/F24F25F26v3,v8=δ4F17−F16F26/F25F26v3, and v9=−F17/F26v3.After going through detailed computations and simplification, we have the following bifurcation coefficientsa and b as
(32)a=2v3u1u3∂2f30,0∂x1∂x3+2v3u2u3∂2f30,0∂x2∂x3=2β2∗v3u3u1+ϵu2.⟹a=2β2∗v3u32G1−G2 where G1=θκ/dd+θ and G2=ϵβ∗τ1VP0+β∗d+τ1Y10+ϵβ∗dY20/dd+τ1. Thus, the bifurcation coefficient a is positive wheneverG1>G2. Furthermore, b=v3u3∂2f20,0/∂x3∂β2=v3u3ϵY20+Y10>0.Hence, it follows from in Castillo-Chavez and Song [28] that model (4) exhibits a backward bifurcation at R3=R2=1 wheneverG1>G2.Theorem 8.
(i)
Model (4) will undergo backward bifurcation ifa=G1>G2>0(ii)
Model (4) will undergo forward bifurcation if a=G1>G2<0


## 3.1. HIV/AIDS Submodel Analysis

We have the HIV/AIDS submodel of (4) whenY2=Y3=Y7=Y8=Y9=Y10=Y12=0 which is given by
(8)dY1dt=Λ−d+λHY1,dY4dt=λHY1−d+κ1+δ1Y4,dY5dt=δ1Y4−d+κ2+δ2Y5,dY6dt=δ2Y5−d+κ3+dAY6,dY11dt=κ1Y4+κ2Y5+κ3Y6−dY11,where the total population isN1t=Y1t+Y4t+Y5t+Y6t+Y11t and the HIV/AIDS single infection force of infection is given by λH=β1/N1Y4+ρ1Y5 with initial conditions Y10>0, Y40≥0, Y50≥0, Y60≥0, andY110≥0.Here, the detailed HIV/AIDS submodel model analysis is given in [6].


## 3.2. Pneumonia Submodel Analysis

From model (4), we have got the pneumonia submodel at Y4=Y5=Y6=Y7=Y8=Y9=Y11=Y12=0, which is given by
(9)dY1dt=1−πΛ+τ1Y2+θY10−d+λPY1,dY2dt=πΛ−ϵλPY2−d+τ1Y2,dY3dt=ϵλPY2+λ4Y1−d+κ+dPY3,dY10dt=κY3−d+θY10.With initial conditions,Y10>0, Y20≥0, Y30≥0, Y100≥0, total population N2t=Y1t+Y2t+Y3t+Y10t, and pneumonia force of infection λP=β2Y3t.In the regionΩ2=Y1,Y2,Y3,Y10∈ℝ4+,N2≤Λ/d, it is easy to show that the set Ω2 is positively invariant and a global attractor of all positive solutions of submodel (9). Hence, it is sufficient to consider the dynamics of model (9) in Ω2 as epidemiologically and mathematically well-posed.


### 3.2.1. Disease-Free Equilibrium Point of the Pneumonia Submodel

The disease-free equilibrium point of the pneumonia submodel is obtained by making the right-hand side of the system (15) as zero and setting the infectious class and treatment class to zero asY3=Y10=0 we have gotY10=Λd+τ1−Λπd/dd+τ1 and Y20=Λπ/d+τ1 such that E20=Y10,Y20,Y30,Y100=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0.


### 3.2.2. The Effective Reproduction Number of the Pneumonia Submodel

The effective reproduction number measures the average number of new infections generated by a typically infectious individual in a community when some strategies are in place, like vaccination or treatment. We calculate the effective reproduction numberR2 using the van den Driesch and Warmouth next-generation matrix approach [23]. The Effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE20/∂xj∂νiE20/∂xj−1 where Fi is the rate of appearance of new infection in compartment i, νi is the transfer of infections from one compartment i to another, and E20 is the disease-free equilibrium point. Then, after a long calculation, we have got
(10)F=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1000,V=d+κ+dP0−κd+θ.Then, using Mathematica, we have got(11)V−1=1d+κ+dP0γθ+dd+κ+dP1θ+d,FV−1=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP000.The characteristic equation of the matrixFV−1is
(12)β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP−λ000−λ=0.Then, the spectral radius (effective reproduction numberR2) of FV−1 of the pneumonia submodel (9) is R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP. Here, R2 is the effective reproduction number for pneumonia infection.


### 3.2.3. Local and Global Stability of the Disease-Free Equilibrium Point

Theorem 3.
The disease-free equilibrium point (DFE)E20 of the pneumonia submodel (9) is locally asymptotically stable if R2<1, otherwise unstable.Proof.
The local stability of the disease-free equilibrium of the system (9) can be studied from its Jacobian matrix at the disease-free equilibrium point E20=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0 and Routh Hurwitz stability criteria. The Jacobian matrix of a dynamical system (9) at the disease-free equilibrium point is given by
(13)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπd+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1−d+κ+dP000κ−d+θ.
Then, the characteristic equation of the above Jacobian matrix is given by(14)−d−λτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−λ−ϵβ2Λπd+τ1000M−λ000κ−d+θ−λ=0,where M=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1−d+κ+dP.
After some steps, we have gotλ1=−d<0 orλ2=−d+τ1<0 orλ3=d+κ+dPR2−1<0 ifR2<1 or λ4=−d+θ<0. Therefore, since all the eigenvalues of the characteristics polynomial of the system (9) are negative if R2<1, the disease-free equilibrium point of the pneumonia submodel is locally asymptotically stable.☐


### 3.2.4. Existence of EEP for the Pneumonia Submodel

Let an arbitrary endemic equilibrium point of pneumonia-only dynamical system (9) be denoted by E2∗=Y1∗,Y2∗,Y3∗,Y10∗. Moreover, let λP∗=β2Y3∗ be the associated pneumonia mass action incidence rate (“force of infection”) at an equilibrium point. To find conditions for the existence of an arbitrary equilibrium point(s) for which pneumonia infection is endemic in the population, the equations of model (9) are solved in terms of the force of infection rate λP∗=β2Y3∗ at an endemic equilibrium point. Setting the right-hand sides of the equations of model (9) to zero and we have got Y2∗=πΛ/ϵλP∗+d+τ1, Y10∗=κY3∗/d+θ and substitute Y2∗ and Y10∗ in to Y1∗, we obtain Y1∗=1−πΛ+τ1Y2∗+θTP∗/d+λP∗=1−πΛ/d+λP∗d+λP∗+πΛτ1/ϵλP∗+d+τ1d+λP∗+θγY3∗/d+θd+λP∗ and substitute Y2∗ and Y1∗ inY3∗, we obtain
(15)Y3∗=πΛϵλP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+1−πΛλP∗d+κ+dPd+θd+λP∗d+κ+dPd+λP∗d+κ+dPd+θd+λP∗−θκλP∗+πΛτ1λP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+λP∗d+γ+dPd+θd+λP∗−θκλP∗.Finally, substituteY3∗ in to pneumonia submodel (9) force of infection λP∗=β2Y3∗ as
(16)λP∗=β2Y3∗=β2πΛϵλP∗d+θd+λP∗ϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+β21−πΛλP∗d+θd+κ+dPd+θd+λP∗−θκλP∗+β2πΛτ1λP∗d+θϵλP∗+κ+τ1d+κ+dPd+θd+λP∗−θκλP∗,and lettingm1=d+κ+dP,m2=d+τ1, and m3=d+θ, we have got a2λP∗2+a1λP∗+a0=0 where a2=m1m3ϵ−θκϵ>0, a1=m1m3dϵ+m1m2m3−m2θκ−β2Λm3ϵ, and a0=m1m2m3μ1−R2>0 if R2<1.Here, the nonzero equilibrium(s) of the model (9) satisfies fλP∗=a2λP∗2+a1λP∗+a0=0 so that the quadratic equation can be analyzed for the possibility of multiple equilibriums. It is worth noting that the coefficienta2is always positive and a0 is positive (negative) if RP is less than (greater than) unity, respectively. Hence, we have established the following result.Theorem 4.
The pneumonia submodel (9) has the following:
(i)
Exactly one unique endemic equilibrium ifa0<0 (i.e.,R2>1)(ii)
Exactly one unique endemic equilibrium ifa1<0, and a0=0 or a12−4a2a0=0(iii)
Exactly two endemic equilibriums ifa0>0(i.e., R2<1), a1<0, anda12−4a2a0>0(iv)
No endemic equilibrium otherwise
Here, item (iii) shows the happening of the backward bifurcation in pneumonia submodel (9), i.e., the locally asymptotically stable disease-free equilibrium point coexists with a locally asymptotically stable endemic equilibrium point ifR2<1; examples of the existence of backward bifurcation phenomenon in mathematical epidemiological models, and the causes, can be seen in [2, 9, 22, 24–26]. The epidemiological consequence is that the classical epidemiological requirement of having the reproduction number R2to be less than one, even though necessary, is not sufficient for the effective control of the disease. The existence of the backward bifurcation phenomenon in submodel (9) is now explored.


### 3.2.5. Bifurcation Analysis

It is instructive to explore the possibility of backward bifurcation in model (9).Theorem 5.
Model (9) exhibits backward bifurcation at R2=1 whenever the inequality D1>D2holds.Here, we apply the center manifold theory in [27]; however, to apply this theory, the following simplification and change of variables are made.LetY1=x1,Y2=x2,Y3=x3, and Y10=x4 such thatN2=x1+x2+x3+x4. Furthermore, by using vector notationX=x1,x2,x3,x4T, pneumonia submodel (9) can be written in the form dX/dt=FX withF=f1,f2,f3,f4T, as follows:
(17)dx1dt=f1=1−πΛ+τ1x2+θx4−d+λPx1,dx2dt=f2=πΛ−ϵλP+d+τ1x2,dx3dt=f3=ϵλPx2+λPx1−d+κ+dPx3,dx4dt=f4=κx3−d+θx4,with λP=β2x3, then the method entails evaluating the Jacobian of system (17) at the DFE pointE20, denoted by JE20, and this gives us
(18)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπμ+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπμdd+τ1−d+κ+dP000κ−d+θ.Consider, next, the case whenRP=1. Suppose, further, that β2=β∗ is chosen as a bifurcation parameter.Solving forβ2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1 and we have gotβ2=β∗=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd and
(19)Jβ∗=−dτ1−β∗Λd+τ1+β∗Λπddd+τ1θ0−d+τ1−ϵβ∗Λπd+τ1000β∗ϵΛπd+β∗Λd+τ1−β∗Λπddd+τ1−d+κ+dP000γ−μ+θ.After some steps of the calculation, we have got the eigenvalues ofJβ∗ as λ1=−d or λ2=−d+τ1 or λ3=0 or λ4=−d+θ.It follows that the JacobianJE20 of (17) at the DFE, with β2=β∗, denoted by Jβ∗, has a simple zero eigenvalue with all the remaining eigenvalues having a negative real part. Hence, the center manifold theory [27] can be used to analyze the dynamics of model (9). In particular, Theorem 2 of Castillo-Chavez and Song [28] will be used to show that model (9) undergoes backward bifurcation atR2=1Eigenvectors ofJβ∗: for the caseR2=1, it can be shown that the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) has a right eigenvectors associated with the zero eigenvalue given by u=u1,u2,u3,u4T with values
(20)u1=−ϵβ∗Λπdτ1d+θ−β∗Λd+τ12d+θ+β∗Λπdd+τ1d+θ+θκdd+τ12d2d+τ12u3,u2=−ϵβ∗Λπd+τ12u3,u3=u3>0,u4=κd+θu3.Similarly, the left eigenvector associated with the zero eigenvalues atβ2=β∗ given by v=v1,v2,v3,v4T arev1=v2=v4=0,v3=v3>0.After long calculations, the bifurcation coefficientsa and b are obtained as a=D1−D2 where D1=β∗Λπdd+τ1d+θ+θκdd+τ12/d2d+τ12, and D2=ϵβ∗Λπdτ1d+θ+β∗Λd+τ12d+θ/d2d+τ12+ϵϵβ∗Λπ/d+τ12.Thus, the bifurcation coefficienta is positive ifD1>D2. Furthermore, b=v3u2u3Λd+τ1−Λπd/dd+τ1>0.Hence, from in Castillo-Chavez and Song [28], model (9) exhibits a backward bifurcation at R2=1 wheneverD1>D2.


## 3.2.1. Disease-Free Equilibrium Point of the Pneumonia Submodel

The disease-free equilibrium point of the pneumonia submodel is obtained by making the right-hand side of the system (15) as zero and setting the infectious class and treatment class to zero asY3=Y10=0 we have gotY10=Λd+τ1−Λπd/dd+τ1 and Y20=Λπ/d+τ1 such that E20=Y10,Y20,Y30,Y100=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0.


## 3.2.2. The Effective Reproduction Number of the Pneumonia Submodel

The effective reproduction number measures the average number of new infections generated by a typically infectious individual in a community when some strategies are in place, like vaccination or treatment. We calculate the effective reproduction numberR2 using the van den Driesch and Warmouth next-generation matrix approach [23]. The Effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE20/∂xj∂νiE20/∂xj−1 where Fi is the rate of appearance of new infection in compartment i, νi is the transfer of infections from one compartment i to another, and E20 is the disease-free equilibrium point. Then, after a long calculation, we have got
(10)F=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1000,V=d+κ+dP0−κd+θ.Then, using Mathematica, we have got(11)V−1=1d+κ+dP0γθ+dd+κ+dP1θ+d,FV−1=β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP000.The characteristic equation of the matrixFV−1is
(12)β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1d+κ+dP−λ000−λ=0.Then, the spectral radius (effective reproduction numberR2) of FV−1 of the pneumonia submodel (9) is R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP. Here, R2 is the effective reproduction number for pneumonia infection.


## 3.2.3. Local and Global Stability of the Disease-Free Equilibrium Point

Theorem 3.
The disease-free equilibrium point (DFE)E20 of the pneumonia submodel (9) is locally asymptotically stable if R2<1, otherwise unstable.Proof.
The local stability of the disease-free equilibrium of the system (9) can be studied from its Jacobian matrix at the disease-free equilibrium point E20=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0 and Routh Hurwitz stability criteria. The Jacobian matrix of a dynamical system (9) at the disease-free equilibrium point is given by
(13)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπd+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπddd+τ1−d+κ+dP000κ−d+θ.
Then, the characteristic equation of the above Jacobian matrix is given by(14)−d−λτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−λ−ϵβ2Λπd+τ1000M−λ000κ−d+θ−λ=0,where M=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1−d+κ+dP.
After some steps, we have gotλ1=−d<0 orλ2=−d+τ1<0 orλ3=d+κ+dPR2−1<0 ifR2<1 or λ4=−d+θ<0. Therefore, since all the eigenvalues of the characteristics polynomial of the system (9) are negative if R2<1, the disease-free equilibrium point of the pneumonia submodel is locally asymptotically stable.☐


## 3.2.4. Existence of EEP for the Pneumonia Submodel

Let an arbitrary endemic equilibrium point of pneumonia-only dynamical system (9) be denoted by E2∗=Y1∗,Y2∗,Y3∗,Y10∗. Moreover, let λP∗=β2Y3∗ be the associated pneumonia mass action incidence rate (“force of infection”) at an equilibrium point. To find conditions for the existence of an arbitrary equilibrium point(s) for which pneumonia infection is endemic in the population, the equations of model (9) are solved in terms of the force of infection rate λP∗=β2Y3∗ at an endemic equilibrium point. Setting the right-hand sides of the equations of model (9) to zero and we have got Y2∗=πΛ/ϵλP∗+d+τ1, Y10∗=κY3∗/d+θ and substitute Y2∗ and Y10∗ in to Y1∗, we obtain Y1∗=1−πΛ+τ1Y2∗+θTP∗/d+λP∗=1−πΛ/d+λP∗d+λP∗+πΛτ1/ϵλP∗+d+τ1d+λP∗+θγY3∗/d+θd+λP∗ and substitute Y2∗ and Y1∗ inY3∗, we obtain
(15)Y3∗=πΛϵλP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+1−πΛλP∗d+κ+dPd+θd+λP∗d+κ+dPd+λP∗d+κ+dPd+θd+λP∗−θκλP∗+πΛτ1λP∗d+κ+dPd+θd+λP∗d+κ+dPϵλP∗+d+τ1d+λP∗d+γ+dPd+θd+λP∗−θκλP∗.Finally, substituteY3∗ in to pneumonia submodel (9) force of infection λP∗=β2Y3∗ as
(16)λP∗=β2Y3∗=β2πΛϵλP∗d+θd+λP∗ϵλP∗+d+τ1d+κ+dPd+θd+λP∗−θκλP∗+β21−πΛλP∗d+θd+κ+dPd+θd+λP∗−θκλP∗+β2πΛτ1λP∗d+θϵλP∗+κ+τ1d+κ+dPd+θd+λP∗−θκλP∗,and lettingm1=d+κ+dP,m2=d+τ1, and m3=d+θ, we have got a2λP∗2+a1λP∗+a0=0 where a2=m1m3ϵ−θκϵ>0, a1=m1m3dϵ+m1m2m3−m2θκ−β2Λm3ϵ, and a0=m1m2m3μ1−R2>0 if R2<1.Here, the nonzero equilibrium(s) of the model (9) satisfies fλP∗=a2λP∗2+a1λP∗+a0=0 so that the quadratic equation can be analyzed for the possibility of multiple equilibriums. It is worth noting that the coefficienta2is always positive and a0 is positive (negative) if RP is less than (greater than) unity, respectively. Hence, we have established the following result.Theorem 4.
The pneumonia submodel (9) has the following:
(i)
Exactly one unique endemic equilibrium ifa0<0 (i.e.,R2>1)(ii)
Exactly one unique endemic equilibrium ifa1<0, and a0=0 or a12−4a2a0=0(iii)
Exactly two endemic equilibriums ifa0>0(i.e., R2<1), a1<0, anda12−4a2a0>0(iv)
No endemic equilibrium otherwise
Here, item (iii) shows the happening of the backward bifurcation in pneumonia submodel (9), i.e., the locally asymptotically stable disease-free equilibrium point coexists with a locally asymptotically stable endemic equilibrium point ifR2<1; examples of the existence of backward bifurcation phenomenon in mathematical epidemiological models, and the causes, can be seen in [2, 9, 22, 24–26]. The epidemiological consequence is that the classical epidemiological requirement of having the reproduction number R2to be less than one, even though necessary, is not sufficient for the effective control of the disease. The existence of the backward bifurcation phenomenon in submodel (9) is now explored.


## 3.2.5. Bifurcation Analysis

It is instructive to explore the possibility of backward bifurcation in model (9).Theorem 5.
Model (9) exhibits backward bifurcation at R2=1 whenever the inequality D1>D2holds.Here, we apply the center manifold theory in [27]; however, to apply this theory, the following simplification and change of variables are made.LetY1=x1,Y2=x2,Y3=x3, and Y10=x4 such thatN2=x1+x2+x3+x4. Furthermore, by using vector notationX=x1,x2,x3,x4T, pneumonia submodel (9) can be written in the form dX/dt=FX withF=f1,f2,f3,f4T, as follows:
(17)dx1dt=f1=1−πΛ+τ1x2+θx4−d+λPx1,dx2dt=f2=πΛ−ϵλP+d+τ1x2,dx3dt=f3=ϵλPx2+λPx1−d+κ+dPx3,dx4dt=f4=κx3−d+θx4,with λP=β2x3, then the method entails evaluating the Jacobian of system (17) at the DFE pointE20, denoted by JE20, and this gives us
(18)JE20=−dτ1−β2Λd+τ1+β2Λπddd+τ1θ0−d+τ1−ϵβ2Λπμ+τ1000β2ϵΛπd+β2Λd+τ1−β2Λπμdd+τ1−d+κ+dP000κ−d+θ.Consider, next, the case whenRP=1. Suppose, further, that β2=β∗ is chosen as a bifurcation parameter.Solving forβ2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1 and we have gotβ2=β∗=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd and
(19)Jβ∗=−dτ1−β∗Λd+τ1+β∗Λπddd+τ1θ0−d+τ1−ϵβ∗Λπd+τ1000β∗ϵΛπd+β∗Λd+τ1−β∗Λπddd+τ1−d+κ+dP000γ−μ+θ.After some steps of the calculation, we have got the eigenvalues ofJβ∗ as λ1=−d or λ2=−d+τ1 or λ3=0 or λ4=−d+θ.It follows that the JacobianJE20 of (17) at the DFE, with β2=β∗, denoted by Jβ∗, has a simple zero eigenvalue with all the remaining eigenvalues having a negative real part. Hence, the center manifold theory [27] can be used to analyze the dynamics of model (9). In particular, Theorem 2 of Castillo-Chavez and Song [28] will be used to show that model (9) undergoes backward bifurcation atR2=1Eigenvectors ofJβ∗: for the caseR2=1, it can be shown that the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) has a right eigenvectors associated with the zero eigenvalue given by u=u1,u2,u3,u4T with values
(20)u1=−ϵβ∗Λπdτ1d+θ−β∗Λd+τ12d+θ+β∗Λπdd+τ1d+θ+θκdd+τ12d2d+τ12u3,u2=−ϵβ∗Λπd+τ12u3,u3=u3>0,u4=κd+θu3.Similarly, the left eigenvector associated with the zero eigenvalues atβ2=β∗ given by v=v1,v2,v3,v4T arev1=v2=v4=0,v3=v3>0.After long calculations, the bifurcation coefficientsa and b are obtained as a=D1−D2 where D1=β∗Λπdd+τ1d+θ+θκdd+τ12/d2d+τ12, and D2=ϵβ∗Λπdτ1d+θ+β∗Λd+τ12d+θ/d2d+τ12+ϵϵβ∗Λπ/d+τ12.Thus, the bifurcation coefficienta is positive ifD1>D2. Furthermore, b=v3u2u3Λd+τ1−Λπd/dd+τ1>0.Hence, from in Castillo-Chavez and Song [28], model (9) exhibits a backward bifurcation at R2=1 wheneverD1>D2.


## 3.3. Analysis of the Full HIV/AIDS-Pneumonia Coinfection

Having analyzed the dynamics of the two submodels, that is, HIV/AIDS submodel (8) and the pneumonia submodel (9), the complete HIV/AIDS-pneumonia coinfection model (4) is now considered (the analysis is done in the positively invariant regionΩ given in (7)).


### 3.3.1. Disease-Free Equilibrium Point of the HIV/AIDS-Pneumonia Coinfection

The disease-free equilibrium point of model (4) is obtained by setting all the infectious classes and treatment classes to zero such that Y3=Y4=Y5=Y6=Y7=Y8=Y9=Y10=Y11=Y12=0and hence E30== Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.


### 3.3.2. Effective Reproduction Number of the HIV/AIDS-Pneumonia Coinfection

The basic reproduction number, denoted byR0, is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [6, 23, 28]. For simple classical models if R0<1, then it means that on average, an infected individual infects less than one susceptible over the course of its infectious period and the disease cannot grow. If however, R0>1, then an infected individual infects more than one susceptible over the course of its infectious period and the disease will persist. For more complicated models with several infected compartments, this simple heuristic definition of R0 is insufficient [23]. Due to its importance, researchers have sought to find ways of determiningR0. Two important concepts in modeling outbreaks of infectious diseases are the basic reproduction number, universally denoted by R0 and the generation time (the average time from symptom onset in a primary case to symptom onset in a secondary case), which jointly determine the likelihood and speed of epidemic outbreaks [29].Here, we calculated the HIV/AIDS-pneumonia coinfection effective reproduction numberR3 of model (4) using the van den Driesch and Warmouth next-generation matrix approach [23]. The effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE30/∂xj∂νiE30/∂xj−1 whereFiis the rate of appearance of new infection in compartmenti, νiis the transfer of infections from one compartment ito another, and E30 is the disease-free equilibrium point E30=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.After long detailed calculations, the transition matrixF is given by
(21)F=A0000000000β1β1ρ1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,and the transmission matrix V is given by
(22)V=D10000000000D2000000000−δ1D3000000000−δ2D40000000000D5000000000−δ3D6000000000−δ4D7000−κ000000θ000−κ1−κ2−κ30000d00000−σ1−σ2−σ300d,whereD1=d+κ+dP,D2=d+κ1+δ1,D3=d+κ2+δ2,D4=d+κ3+dA,D5=d+dP+σ1+δ3,D6=d+dP+σ2+δ4,andD7=d+dAP+ε3.Then, by using Mathematica, we have got(23)FV−1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1D10000000000β1D2+δ1β1ρ1D2D3β1ρ1D3000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.The characteristic equation of the matrixFV−1 is given by
(24)A1−λ0000000000B−λβ1ρ1D30000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ=0,whereA1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1/D1,B=β1/D2+δ1β1ρ1/D2D3; then, the eigenvalues of FV−1 are λ1=β2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1orλ2=β1/D2+δ1β1ρ1/D2D3 or λ3=λ4=λ5=λ6=λ7=λ8=λ9=λ10=0.Thus, the effective reproduction number of the HIV/AIDS-pneumonia coinfection dynamical system (4) is the dominant eigenvalue of the matrixFV−1 which is given byR3=maxλ1,λ2=maxβ2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1,β1/D2+δ1β1ρ1/D2D3. Here, R2=β2ϵπdΛ+β2Λd+τ1−β2πΛd/d+κ+dPd+τ1 is the effective reproduction number for pneumonia-only infected individual and R1=β1/d+κ1+δ1+β1ρ1δ1/d+κ1+δ1d+κ2+κ2 is the basic reproduction for HIV/AIDS-only infected individual.Here,R1 represent the basic reproduction number for HIV/AIDS submodel, R2 and R3 are the effective reproduction numbers for the pneumonia submodel and HIV/AIDS-pneumonia coinfection model, respectively. The following three outcomes are possible: (i) for R1<1, the HIV/AIDS submodel disease-free steady state E1 is globally stable in the region Ω1, and HIV is not spreading in the community; (ii) for R2<1, then E2 is not globally stable in the region Ω2, and pneumonia may spread through the community; (iii) for R3<1, the steady state E3is not globally stable in the region Ω, and HIV/AIDS-pneumonia coinfection may spread through the community.Note that none of the parameters corresponding to coinfection treatment (i.e.,σ1orσ2orσ3) are present in the expression forR3, indicating no impact of treating coinfected population onR3.


### 3.3.3. Locally Asymptotically Stability of the Disease-Free Equilibrium (DFE)

Theorem 6.
The disease-free equilibrium of model (4) above is locally asymptotically stable if R3<1, otherwise unstable.Proof.
The Jacobian matrixJE30 of model (4) at E30 is given by
(25)JE30=−dτ1−β2Y10−β1N0Y10−β1N0ρ1Y10000θ0000−d−τ1−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1000000000000Z2β1N0ρ1Y100000000000δ1Z300000000000δ2Z4000000000000Z500000000000δ3Z600000000000δ4Z700000κ000000Z800000κ1κ2κ30000−d0000000σ1σ2σ300−d,
whereZ1=β2ϵY20+β2Y10−d+κ+dP,Z2=β1/N0Y10−d+κ1+δ1, Z3=−d+κ2+δ2, Z4=−d+κ3+dA,Z5=−d+dP+σ1+δ3,Z6=−d+dP+σ2+δ4,Z7=−d+dAP+σ3, andZ8=−d+θ.
Then, the characteristic equation of the Jacobian matrixJE30 is given by
(26)−d−λτ1−β2Y10−β1N0Y10−β1N0ρ1Y100000θ000−d−τ1−λ−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1−λ000000000000Z2−λβ1N0ρ1Y100000000000δ1Z3−λ00000000000δ2Z4−λ000000000000Z5−λ00000000000δ3Z6−λ00000000000δ4Z7−λ00000κ000000Z8−λ00000κ1κ2κ30000−d−λ0000000σ1σ2σ300−d−λ=0.
After detailed calculations, we have got that
λ1=λ2=λ3=−d<0 or λ4=d+κ+dPR2−1<0 if R2<1 or λ5=−d+τ1<0 or λ6=−d+κ3+dA<0 or λ7=−d+dP+σ1+δ3<0 or λ8=−d+dP+σ2+δ4<0 or λ9=−d+dAP+σ3<0 or λ10=−d+θ<0or a2λ2+a1λ+a0=0 for a2=1,a1=d+κ2+δ2+d+κ1+δ11−Y10/N0RY4>0 if RY4<1 anda0=d+κ2+δ2d+κ1+δ11−Y10/N0RY5>0 if RY5<1.
Then, by applying Routh-Hurwitz stability criteria sincea2=1>0, a1>0, and a0>0, all the eigenvalues of the Jacobian matrix are negative if R1<1 andR2<1, i.e.,R3=maxR1,R2<1. Thus, the disease-free equilibrium point (DFE) of HIV/AIDS-pneumonia coinfection model (4) is locally asymptotically stable if
(27)R3=maxR1,R2<1.☐


### 3.3.4. Existence of Endemic Equilibrium Point (EEP) for the Full Model

The endemic equilibrium point (EEP) of full model (4) is denoted by E3∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ which occurs when the disease persists in the community. From the analysis of HIV/AIDS-only submodel (8) and the pneumonia-only submodel from (9), we have shown that there is no endemic equilibrium point if R1<1and there is/are an endemic equilibrium point(s) if R2<1 implies that there is/are endemic equilibrium point(s) if R3<1 for the coinfection model and hence there is a bifurcation point for the full model. The endemic equilibrium of system (4) is obtained as
(28)Y1∗=1−πΛ+τ1Y2∗+θY10∗d+λHC∗+λPC∗,Y2∗=πΛϵλPC∗+d+τ1+λHC∗,Y3∗=ϵλPC∗Y2∗+λPC∗Y1∗νλHC∗+d+κ+dP,Y4∗=λHC∗Y1∗d+κ1+δ1+ψ1λPC∗,Y5∗=δ1H1∗d+κ2+κ2+ψ2λPC∗,Y6∗=δ2Y5∗d+κ3+dA+δ3λPC∗,Y7∗=ψ1λPC∗Y4∗+νλHC∗Y3∗d+dP+σ1+δ3,Y8∗=δ2λPC∗Y5∗+δ3Y7∗d+dP+σ2+δ4,Y9∗=ψ3λPC∗Y6∗+δ4Y8∗d+dAP+κ3,Y10∗=κY3∗d+θ,Y11∗=κ1Y4∗+κ2Y5∗+κ3Y6∗d,andY12∗=κ1Y7∗+κ2Y8∗+κ3Y9∗d.Summary of endemic equilibrium points: the explicit computation of the endemic equilibrium of coinfection model (4) given in (28) in terms of model parameters is difficult analytically; however, model (4) endemic equilibriums correspond to the following:
(1)
E4∗=Y1∗,0,Y4∗,Y5∗,Y6∗,0,0,0,0,0,Y11∗,0,ifR1>1 is the pneumonia free (HIV) endemic equilibrium point. The analysis of the equilibrium E1∗ is similar to the endemic equilibrium E1∗ in model (7)(2)
E5∗=Y1∗,Y2∗,Y3∗,0,0,0,0,0,0,Y10∗,0,0, if R2>1 is the HIV/AIDS free (pneumonia) endemic equilibrium point. The analysis of the equilibrium E5∗ is similar to the endemic equilibrium E2∗ in equation (9)(3)
E6∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ is the HIV/AIDS-pneumonia coinfection endemic equilibrium point. It exists when each component of E6∗ in equation (28) is positive and summarizes the existence of the endemic equilibrium points in the following theorem


### 3.3.5. Bifurcation Analysis

The threshold quantityR3=maxR1,R2 is the effective reproduction number of the system (4) where R1 and R2 are defined as above.Theorem 7.
Model (4) exhibits the phenomenon of backward bifurcation at R3=1 whenever the inequality G1>G2 holds.The phenomenon of backward bifurcation can be proved with the concept of the center manifold theory [10, 27] on coepidemic model (4). To apply this theory, the following simplification and change of variables are made.LetY1=x1, Y2=x2, Y3=x3,Y4=x4,Y5=x5, Y6=x6,Y7=x7,Y8=x8,Y9=x9,Y10=x10,Y11=x11,and Y12=x12 so that N=x1+x2+x3+x4+x5,+x6+x7,+x8+x9+x10+x11+x12.Further, by using vector notationX=x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12T, complete model (4) can be written in the form dX/dt=FX with F=f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12T, as follows:
(29)dx1dt=f1=1−πΛ+τ1x2+θx10−d+λHC+λPCx1,dx2dt=f2=πΛ−ϵλPCx2−d+τ1+λHCx2,dx3dt=f3=ϵλPCx2+λPCx1−νλHC+d+κ+dPx3,dx4dt=f4=λHCx1−d+κ1+δ1+ψλPCx4,dx5dt=f5=δ1x4−d+κ2+κ2+ψ2λPCx5,dx6dt=f6=δ2x5−d+κ3+dA+ψ3λPCx6,dx7dt=f7=ψ1λPCx4+νλHCx3−d+dP+σ1+δ3x7,dx8dt=f8=ψ2λPCx5+δ3x7−d+dP+σ2+δ4x8,dx9dt=f9=ψ3λPCx6+δ4x8−d+dAP+κ3x9,dx10dt=f10=κx3−d+θx10,dx11dt=f11=κ1x4+κ2x5+κ3x6−dx11,dx12dt=f12=κ1x7+κ2x8+κ3x9−dx12,with λHC=β1/Nx4+ρ1x5+ρ2x7+ρ2x8whereρ3≥ρ2≥ρ1≥1 andλPC=β2x3+ω1x7+ω2x8+ω3x9 whereω3≥ω2≥ω1≥1; then, the method entails evaluating the Jacobian of system (29) at the DFEE30, denoted by JE30, and this gives us
(30)JE30=−dτ1F1F2F30F4F5F6θ000F7F8F9F100F11F12F1300000F14000F15F16F17000000F18F190F20F210000000δ1F2200000000000δ2F23000000000000F2400000000000δ3F2500000000000d4F2600000κ000000F2700000κ1κ2κ30000−d0000000κ1κ2κ300−d,whereF1=−β2Y10,F2=−β1Y10/Y10+Y20,F3=−β1ρ1Y10/Y10+Y20, F4=−β1ρ2Y10/Y10+Y20−β2ω1Y10,F5=−β1ρ3Y10/Y10+Y20−β2ω2Y10,F6=−β2ω3Y10,F7=−d+τ1,F8=−ϵβ2Y20,F9=−β1Y20/Y10+Y20,F10=−β1ρ1Y20/Y10+Y20,F11=−ϵβ2ω1VP0−β1ρ2YP0/Y10+YP0,F12=−ϵβ2ω2VP0−β1ρ3Y20/Y10+Y20,F13=−ϵβ2ω3Y20,F14=ϵβ2Y20+β2Y10−d+κ+dP,F15=ϵβ2ω1Y20+β2ω1Y10,F16=ϵβ2ω2Y20+β2ω2Y10,F17=ϵβ2ω3Y20+β2ω3Y10, F18=β1Y10/Y10+Y20−d+κ1+δ1,F19=β1ρ1Y10/Y10+Y20,F20=β1ρ2Y10/Y10+Y20,F21=β1ρ3Y10/Y10+Y20,F22=−d+κ2+δ2,F23=−d+κ3+dA,F24=−d+dP+σ1+δ3,F25=−d+dP+σ2+δ4,F26=−d+dAP+σ3, and F27=−d+θ.Without loss of generality, consider the case whenR2>R1, and R3=1, so that R2=1. Furthermore, let β2=β∗ is chosen as a bifurcation parameter. Solving for β2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1, we have got the value β∗=β2=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd.After solving the JacobianJE30 of the system (29) at the DFE, with β2=β∗, we obtained the eigenvalues as λ1=−d<0 or λ2=−d+τ1<0 or λ3=0 or λ4=−d<0 orλ5=−d<0 orλ6=−d+θ<0 or λ7=−d+dAP+κ3<0 orλ8=−d+κ3+dA<0 or λ9=−d+dP+σ1+δ3<0 or λ10=−d+dP+σ1+δ4<0or
(31)a2λ2+a1λ+a0=0,wherea2=1>0, a1=d+κ1+δ11−Y10/Y10+Y40RY4+d+κ2+δ2>0 if RY4<1, and a0=a0=d+κ1+δ1d+κ2+δ21−Y10/Y10+Y20R1>0 if R1<1.Equation (31) has/have no positive root/s whenever R1<1, and hence, both eigenvalues are negative. It follows that the Jacobian JE30 of (29) at the DFE, withβ2=β∗, denoted by Jβ∗, has a simple zero eigenvalue (with all other eigenvalues having negative real part). Hence, the center manifold theory [27] can be used to analyze the dynamics of model (4). In particular, the Castillo-Chavez and Song theorem [28] will be used to show that model (4) undergoes backward bifurcation at RP=1.Eigenvectors ofJβ∗: for the case when RP=1, the right eigenvectors of the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) associated with the zero eigenvalue given by u=u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12T areu1=−τ1F8F27+F1F7F27−θγF7/dF7F27u3,u2=−F8/F7u3,u3=u3>0,u10=−κ/F27u3, and u4=u5=u6=u7=u8=u9=u11=u12=0.The left eigenvectors associated with the zero eigenvalue atβ2=β2∗ satisfyingv.w=1 given by v=v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12are v1=v2=v4=v5=v6=v10=v11=v12=0, v3=v3>0, v7=−δ3δ4F17+δ3F16F26−F15F25F26/F24F25F26v3,v8=δ4F17−F16F26/F25F26v3, and v9=−F17/F26v3.After going through detailed computations and simplification, we have the following bifurcation coefficientsa and b as
(32)a=2v3u1u3∂2f30,0∂x1∂x3+2v3u2u3∂2f30,0∂x2∂x3=2β2∗v3u3u1+ϵu2.⟹a=2β2∗v3u32G1−G2 where G1=θκ/dd+θ and G2=ϵβ∗τ1VP0+β∗d+τ1Y10+ϵβ∗dY20/dd+τ1. Thus, the bifurcation coefficient a is positive wheneverG1>G2. Furthermore, b=v3u3∂2f20,0/∂x3∂β2=v3u3ϵY20+Y10>0.Hence, it follows from in Castillo-Chavez and Song [28] that model (4) exhibits a backward bifurcation at R3=R2=1 wheneverG1>G2.Theorem 8.
(i)
Model (4) will undergo backward bifurcation ifa=G1>G2>0(ii)
Model (4) will undergo forward bifurcation if a=G1>G2<0


## 3.3.1. Disease-Free Equilibrium Point of the HIV/AIDS-Pneumonia Coinfection

The disease-free equilibrium point of model (4) is obtained by setting all the infectious classes and treatment classes to zero such that Y3=Y4=Y5=Y6=Y7=Y8=Y9=Y10=Y11=Y12=0and hence E30== Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.


## 3.3.2. Effective Reproduction Number of the HIV/AIDS-Pneumonia Coinfection

The basic reproduction number, denoted byR0, is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [6, 23, 28]. For simple classical models if R0<1, then it means that on average, an infected individual infects less than one susceptible over the course of its infectious period and the disease cannot grow. If however, R0>1, then an infected individual infects more than one susceptible over the course of its infectious period and the disease will persist. For more complicated models with several infected compartments, this simple heuristic definition of R0 is insufficient [23]. Due to its importance, researchers have sought to find ways of determiningR0. Two important concepts in modeling outbreaks of infectious diseases are the basic reproduction number, universally denoted by R0 and the generation time (the average time from symptom onset in a primary case to symptom onset in a secondary case), which jointly determine the likelihood and speed of epidemic outbreaks [29].Here, we calculated the HIV/AIDS-pneumonia coinfection effective reproduction numberR3 of model (4) using the van den Driesch and Warmouth next-generation matrix approach [23]. The effective reproduction number is the largest (dominant) eigenvalue (spectral radius) of the matrix FV−1=∂FiE30/∂xj∂νiE30/∂xj−1 whereFiis the rate of appearance of new infection in compartmenti, νiis the transfer of infections from one compartment ito another, and E30 is the disease-free equilibrium point E30=Λd+τ1−Λπd/dd+τ1,Λπ/d+τ1,0,0,0,0,0,0,0,0,0,0.After long detailed calculations, the transition matrixF is given by
(21)F=A0000000000β1β1ρ1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,and the transmission matrix V is given by
(22)V=D10000000000D2000000000−δ1D3000000000−δ2D40000000000D5000000000−δ3D6000000000−δ4D7000−κ000000θ000−κ1−κ2−κ30000d00000−σ1−σ2−σ300d,whereD1=d+κ+dP,D2=d+κ1+δ1,D3=d+κ2+δ2,D4=d+κ3+dA,D5=d+dP+σ1+δ3,D6=d+dP+σ2+δ4,andD7=d+dAP+ε3.Then, by using Mathematica, we have got(23)FV−1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1D10000000000β1D2+δ1β1ρ1D2D3β1ρ1D3000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.The characteristic equation of the matrixFV−1 is given by
(24)A1−λ0000000000B−λβ1ρ1D30000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ00000000000−λ=0,whereA1=πϵΛβ2/d+τ1+β2−πΛd+Λd+τ1/dd+τ1/D1,B=β1/D2+δ1β1ρ1/D2D3; then, the eigenvalues of FV−1 are λ1=β2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1orλ2=β1/D2+δ1β1ρ1/D2D3 or λ3=λ4=λ5=λ6=λ7=λ8=λ9=λ10=0.Thus, the effective reproduction number of the HIV/AIDS-pneumonia coinfection dynamical system (4) is the dominant eigenvalue of the matrixFV−1 which is given byR3=maxλ1,λ2=maxβ2ϵπdΛ+β2Λd+τ1−β2πΛd/D1d+τ1,β1/D2+δ1β1ρ1/D2D3. Here, R2=β2ϵπdΛ+β2Λd+τ1−β2πΛd/d+κ+dPd+τ1 is the effective reproduction number for pneumonia-only infected individual and R1=β1/d+κ1+δ1+β1ρ1δ1/d+κ1+δ1d+κ2+κ2 is the basic reproduction for HIV/AIDS-only infected individual.Here,R1 represent the basic reproduction number for HIV/AIDS submodel, R2 and R3 are the effective reproduction numbers for the pneumonia submodel and HIV/AIDS-pneumonia coinfection model, respectively. The following three outcomes are possible: (i) for R1<1, the HIV/AIDS submodel disease-free steady state E1 is globally stable in the region Ω1, and HIV is not spreading in the community; (ii) for R2<1, then E2 is not globally stable in the region Ω2, and pneumonia may spread through the community; (iii) for R3<1, the steady state E3is not globally stable in the region Ω, and HIV/AIDS-pneumonia coinfection may spread through the community.Note that none of the parameters corresponding to coinfection treatment (i.e.,σ1orσ2orσ3) are present in the expression forR3, indicating no impact of treating coinfected population onR3.


## 3.3.3. Locally Asymptotically Stability of the Disease-Free Equilibrium (DFE)

Theorem 6.
The disease-free equilibrium of model (4) above is locally asymptotically stable if R3<1, otherwise unstable.Proof.
The Jacobian matrixJE30 of model (4) at E30 is given by
(25)JE30=−dτ1−β2Y10−β1N0Y10−β1N0ρ1Y10000θ0000−d−τ1−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1000000000000Z2β1N0ρ1Y100000000000δ1Z300000000000δ2Z4000000000000Z500000000000δ3Z600000000000δ4Z700000κ000000Z800000κ1κ2κ30000−d0000000σ1σ2σ300−d,
whereZ1=β2ϵY20+β2Y10−d+κ+dP,Z2=β1/N0Y10−d+κ1+δ1, Z3=−d+κ2+δ2, Z4=−d+κ3+dA,Z5=−d+dP+σ1+δ3,Z6=−d+dP+σ2+δ4,Z7=−d+dAP+σ3, andZ8=−d+θ.
Then, the characteristic equation of the Jacobian matrixJE30 is given by
(26)−d−λτ1−β2Y10−β1N0Y10−β1N0ρ1Y100000θ000−d−τ1−λ−β2ϵY20−β1N0Y20−β1N0ρ1Y20000000000Z1−λ000000000000Z2−λβ1N0ρ1Y100000000000δ1Z3−λ00000000000δ2Z4−λ000000000000Z5−λ00000000000δ3Z6−λ00000000000δ4Z7−λ00000κ000000Z8−λ00000κ1κ2κ30000−d−λ0000000σ1σ2σ300−d−λ=0.
After detailed calculations, we have got that
λ1=λ2=λ3=−d<0 or λ4=d+κ+dPR2−1<0 if R2<1 or λ5=−d+τ1<0 or λ6=−d+κ3+dA<0 or λ7=−d+dP+σ1+δ3<0 or λ8=−d+dP+σ2+δ4<0 or λ9=−d+dAP+σ3<0 or λ10=−d+θ<0or a2λ2+a1λ+a0=0 for a2=1,a1=d+κ2+δ2+d+κ1+δ11−Y10/N0RY4>0 if RY4<1 anda0=d+κ2+δ2d+κ1+δ11−Y10/N0RY5>0 if RY5<1.
Then, by applying Routh-Hurwitz stability criteria sincea2=1>0, a1>0, and a0>0, all the eigenvalues of the Jacobian matrix are negative if R1<1 andR2<1, i.e.,R3=maxR1,R2<1. Thus, the disease-free equilibrium point (DFE) of HIV/AIDS-pneumonia coinfection model (4) is locally asymptotically stable if
(27)R3=maxR1,R2<1.☐


## 3.3.4. Existence of Endemic Equilibrium Point (EEP) for the Full Model

The endemic equilibrium point (EEP) of full model (4) is denoted by E3∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ which occurs when the disease persists in the community. From the analysis of HIV/AIDS-only submodel (8) and the pneumonia-only submodel from (9), we have shown that there is no endemic equilibrium point if R1<1and there is/are an endemic equilibrium point(s) if R2<1 implies that there is/are endemic equilibrium point(s) if R3<1 for the coinfection model and hence there is a bifurcation point for the full model. The endemic equilibrium of system (4) is obtained as
(28)Y1∗=1−πΛ+τ1Y2∗+θY10∗d+λHC∗+λPC∗,Y2∗=πΛϵλPC∗+d+τ1+λHC∗,Y3∗=ϵλPC∗Y2∗+λPC∗Y1∗νλHC∗+d+κ+dP,Y4∗=λHC∗Y1∗d+κ1+δ1+ψ1λPC∗,Y5∗=δ1H1∗d+κ2+κ2+ψ2λPC∗,Y6∗=δ2Y5∗d+κ3+dA+δ3λPC∗,Y7∗=ψ1λPC∗Y4∗+νλHC∗Y3∗d+dP+σ1+δ3,Y8∗=δ2λPC∗Y5∗+δ3Y7∗d+dP+σ2+δ4,Y9∗=ψ3λPC∗Y6∗+δ4Y8∗d+dAP+κ3,Y10∗=κY3∗d+θ,Y11∗=κ1Y4∗+κ2Y5∗+κ3Y6∗d,andY12∗=κ1Y7∗+κ2Y8∗+κ3Y9∗d.Summary of endemic equilibrium points: the explicit computation of the endemic equilibrium of coinfection model (4) given in (28) in terms of model parameters is difficult analytically; however, model (4) endemic equilibriums correspond to the following:
(1)
E4∗=Y1∗,0,Y4∗,Y5∗,Y6∗,0,0,0,0,0,Y11∗,0,ifR1>1 is the pneumonia free (HIV) endemic equilibrium point. The analysis of the equilibrium E1∗ is similar to the endemic equilibrium E1∗ in model (7)(2)
E5∗=Y1∗,Y2∗,Y3∗,0,0,0,0,0,0,Y10∗,0,0, if R2>1 is the HIV/AIDS free (pneumonia) endemic equilibrium point. The analysis of the equilibrium E5∗ is similar to the endemic equilibrium E2∗ in equation (9)(3)
E6∗=Y1∗,Y2∗,Y3∗,Y4∗,Y5∗,Y6∗,Y7∗,Y8∗,Y9∗,Y10∗,Y11∗,Y12∗ is the HIV/AIDS-pneumonia coinfection endemic equilibrium point. It exists when each component of E6∗ in equation (28) is positive and summarizes the existence of the endemic equilibrium points in the following theorem


## 3.3.5. Bifurcation Analysis

The threshold quantityR3=maxR1,R2 is the effective reproduction number of the system (4) where R1 and R2 are defined as above.Theorem 7.
Model (4) exhibits the phenomenon of backward bifurcation at R3=1 whenever the inequality G1>G2 holds.The phenomenon of backward bifurcation can be proved with the concept of the center manifold theory [10, 27] on coepidemic model (4). To apply this theory, the following simplification and change of variables are made.LetY1=x1, Y2=x2, Y3=x3,Y4=x4,Y5=x5, Y6=x6,Y7=x7,Y8=x8,Y9=x9,Y10=x10,Y11=x11,and Y12=x12 so that N=x1+x2+x3+x4+x5,+x6+x7,+x8+x9+x10+x11+x12.Further, by using vector notationX=x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12T, complete model (4) can be written in the form dX/dt=FX with F=f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12T, as follows:
(29)dx1dt=f1=1−πΛ+τ1x2+θx10−d+λHC+λPCx1,dx2dt=f2=πΛ−ϵλPCx2−d+τ1+λHCx2,dx3dt=f3=ϵλPCx2+λPCx1−νλHC+d+κ+dPx3,dx4dt=f4=λHCx1−d+κ1+δ1+ψλPCx4,dx5dt=f5=δ1x4−d+κ2+κ2+ψ2λPCx5,dx6dt=f6=δ2x5−d+κ3+dA+ψ3λPCx6,dx7dt=f7=ψ1λPCx4+νλHCx3−d+dP+σ1+δ3x7,dx8dt=f8=ψ2λPCx5+δ3x7−d+dP+σ2+δ4x8,dx9dt=f9=ψ3λPCx6+δ4x8−d+dAP+κ3x9,dx10dt=f10=κx3−d+θx10,dx11dt=f11=κ1x4+κ2x5+κ3x6−dx11,dx12dt=f12=κ1x7+κ2x8+κ3x9−dx12,with λHC=β1/Nx4+ρ1x5+ρ2x7+ρ2x8whereρ3≥ρ2≥ρ1≥1 andλPC=β2x3+ω1x7+ω2x8+ω3x9 whereω3≥ω2≥ω1≥1; then, the method entails evaluating the Jacobian of system (29) at the DFEE30, denoted by JE30, and this gives us
(30)JE30=−dτ1F1F2F30F4F5F6θ000F7F8F9F100F11F12F1300000F14000F15F16F17000000F18F190F20F210000000δ1F2200000000000δ2F23000000000000F2400000000000δ3F2500000000000d4F2600000κ000000F2700000κ1κ2κ30000−d0000000κ1κ2κ300−d,whereF1=−β2Y10,F2=−β1Y10/Y10+Y20,F3=−β1ρ1Y10/Y10+Y20, F4=−β1ρ2Y10/Y10+Y20−β2ω1Y10,F5=−β1ρ3Y10/Y10+Y20−β2ω2Y10,F6=−β2ω3Y10,F7=−d+τ1,F8=−ϵβ2Y20,F9=−β1Y20/Y10+Y20,F10=−β1ρ1Y20/Y10+Y20,F11=−ϵβ2ω1VP0−β1ρ2YP0/Y10+YP0,F12=−ϵβ2ω2VP0−β1ρ3Y20/Y10+Y20,F13=−ϵβ2ω3Y20,F14=ϵβ2Y20+β2Y10−d+κ+dP,F15=ϵβ2ω1Y20+β2ω1Y10,F16=ϵβ2ω2Y20+β2ω2Y10,F17=ϵβ2ω3Y20+β2ω3Y10, F18=β1Y10/Y10+Y20−d+κ1+δ1,F19=β1ρ1Y10/Y10+Y20,F20=β1ρ2Y10/Y10+Y20,F21=β1ρ3Y10/Y10+Y20,F22=−d+κ2+δ2,F23=−d+κ3+dA,F24=−d+dP+σ1+δ3,F25=−d+dP+σ2+δ4,F26=−d+dAP+σ3, and F27=−d+θ.Without loss of generality, consider the case whenR2>R1, and R3=1, so that R2=1. Furthermore, let β2=β∗ is chosen as a bifurcation parameter. Solving for β2 from R2=1 as R2=β2ϵΛπd+β2Λd+τ1−β2Λπd/dd+τ1d+κ+dP=1, we have got the value β∗=β2=dd+τ1d+κ+dP/ϵΛπd+Λd+τ1−Λπd.After solving the JacobianJE30 of the system (29) at the DFE, with β2=β∗, we obtained the eigenvalues as λ1=−d<0 or λ2=−d+τ1<0 or λ3=0 or λ4=−d<0 orλ5=−d<0 orλ6=−d+θ<0 or λ7=−d+dAP+κ3<0 orλ8=−d+κ3+dA<0 or λ9=−d+dP+σ1+δ3<0 or λ10=−d+dP+σ1+δ4<0or
(31)a2λ2+a1λ+a0=0,wherea2=1>0, a1=d+κ1+δ11−Y10/Y10+Y40RY4+d+κ2+δ2>0 if RY4<1, and a0=a0=d+κ1+δ1d+κ2+δ21−Y10/Y10+Y20R1>0 if R1<1.Equation (31) has/have no positive root/s whenever R1<1, and hence, both eigenvalues are negative. It follows that the Jacobian JE30 of (29) at the DFE, withβ2=β∗, denoted by Jβ∗, has a simple zero eigenvalue (with all other eigenvalues having negative real part). Hence, the center manifold theory [27] can be used to analyze the dynamics of model (4). In particular, the Castillo-Chavez and Song theorem [28] will be used to show that model (4) undergoes backward bifurcation at RP=1.Eigenvectors ofJβ∗: for the case when RP=1, the right eigenvectors of the Jacobian of (29) at β2=β∗ (denoted by Jβ∗) associated with the zero eigenvalue given by u=u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12T areu1=−τ1F8F27+F1F7F27−θγF7/dF7F27u3,u2=−F8/F7u3,u3=u3>0,u10=−κ/F27u3, and u4=u5=u6=u7=u8=u9=u11=u12=0.The left eigenvectors associated with the zero eigenvalue atβ2=β2∗ satisfyingv.w=1 given by v=v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12are v1=v2=v4=v5=v6=v10=v11=v12=0, v3=v3>0, v7=−δ3δ4F17+δ3F16F26−F15F25F26/F24F25F26v3,v8=δ4F17−F16F26/F25F26v3, and v9=−F17/F26v3.After going through detailed computations and simplification, we have the following bifurcation coefficientsa and b as
(32)a=2v3u1u3∂2f30,0∂x1∂x3+2v3u2u3∂2f30,0∂x2∂x3=2β2∗v3u3u1+ϵu2.⟹a=2β2∗v3u32G1−G2 where G1=θκ/dd+θ and G2=ϵβ∗τ1VP0+β∗d+τ1Y10+ϵβ∗dY20/dd+τ1. Thus, the bifurcation coefficient a is positive wheneverG1>G2. Furthermore, b=v3u3∂2f20,0/∂x3∂β2=v3u3ϵY20+Y10>0.Hence, it follows from in Castillo-Chavez and Song [28] that model (4) exhibits a backward bifurcation at R3=R2=1 wheneverG1>G2.Theorem 8.
(i)
Model (4) will undergo backward bifurcation ifa=G1>G2>0(ii)
Model (4) will undergo forward bifurcation if a=G1>G2<0


## 4. Sensitivity and Numerical Analysis

### 4.1. Sensitivity Analysis

Definition. The normalized forward sensitivity index of a variableR3that depends differentiably on a parameter p is defined as SIp=∂R3/∂p∗p/R3 [18].Sensitivity indices allow us to examine the relative importance of different parameters in pneumonia and HIV/AIDS spread and prevalence. The most sensitive parameter has the magnitude of the sensitivity index larger than that of all other parameters. We can calculate the sensitivity index in terms ofR1and R2sinceR3=maxR1,R2. Sensitivity analysis results and the numerical simulation are given in this section with parameters values given in Table 3 where N0 is the total number of the initial population of complete model (4).Table 3
Parameter values used for the full HIV/AIDS-pneumonia coepidemic model simulation.
ParameterValueSourceΛ0.0413∗N0Estimatedd0.02Estimatedδ10.498[7]δ20.08[7]δ30.2885[6]δ40.3105[6]ψ11.1Assumedψ21.2Assumedψ31.4Assumedν1AssumeddP0.1[16]θ0.1[18]dA0.333[6]π0.2[18]τ10.0025[18]ϵ0.002[18]dAP0.42Assumedκ0.2[18]κ10.2[7]κ20.15[7]κ30.13Assumedσ10.498[7]σ20.08[7]σ30.230Assumedβ1Variable[6]β2Variable[6]ρ1,ρ2,ρ3,ω1,ω2,ω31.2,1,1,1,1,1AssumedUsing the values of parameters in Table3, the sensitivity indices are calculated in Tables 4 and 5.Table 4
Sensitivity indices ofR3=R1.
Sensitivity indexValuesSIβ1+1SIρ1+0.6134SIδ1- 0.0639SId-0.3150SIκ1-0.1371SIκ2-0.0264SIδ2-0.0141Table 5
Sensitivity indices ofR3=R2.
Sensitivity indexValuesSIΛ+1SIβ2+1SId-0.4421SIκ-0.6559SIdP-0.3852SIπ-0.3852SIϵ-0.3852SIτ1-0.3852In this paper, with parameter values in Table3, we have got R1=1.386 at β1=2 implies HIV/AIDS spreads in the community and also we have got the indices as shown in Table 4. Here, sensitivity analysis shows that the human recruitment rate Λ and HIV/AIDS spreading rate β1 have the highest impact on the basic reproduction number of HIV/AIDS (R1).Similarly, with parameter values in Table3, we have got R2=9.69 at β2=0.2 imply that pneumonia spreads throughout the community and also we have got the indices as shown in Table 4. Here, sensitivity analysis shows that the foremost sensitive positive parameters are the human recruitment rate Λ and the pneumonia spreading rate β2. The foremost sensitive negative parameter is treatment rate of pneumonia κ which is inversely related to the effective reproduction number R2, i.e., a smaller amount of increase in this parameter value will lead to a greater amount of reduction in the effective reproduction number while a smaller amount of decrement will cause a big increment in the basic reproduction number. Epidemiologically, the most sensitive parameters to R1 and R2 which can be controlled through interventions and preventions are found to beβ1 andβ2, respectively.


### 4.2. Numerical Analysis

In this section, numerical simulation is performed for complete HIV/AIDS-pneumonia coepidemic model (4). With ode45, we have checked the effect of some parameters in the spreading as well as for the control of pneumonia only, HIV/AIDS only, and coepidemic of HIV/AIDS and pneumonia. The parameter values put forward in Table 3 are used for numerical simulation. In the numerical simulation part, we investigated the stability of the endemic equilibrium point of complete model (4), parameter effects on the reproduction numbers, and the impact of treatment mainly on dually infected individuals in the community.


#### 4.2.1. Local Stability of Endemic Equilibrium Point of Complete Model (4)

Figure2 shows that in the long run (after 50 years), the solutions of dynamical system (4) will be converging to its endemic equilibrium point, i.e., the endemic equilibrium point is locally asymptotically stable whenever
(33)R3=maxR1,R2=max1.386,9.69=9.69>1.Figure 2
Local stability of endemic equilibrium point of the coepidemic model (4) whenever R1=1.386 at β1=2 and R2=9.69 at β2=0.2.


#### 4.2.2. Effect of Parameters on the Threshold ParameterR2

In this subsection, as we see in Figure3, we have investigated the effect of pneumonia vaccination portion π on the pneumonia effective reproduction number R2. The figure reflects that when the value ofπ increases, the pneumonia effective reproduction number is going down, and whenever the value of π>0.64 imply R2<1. Therefore, public policymakers must concentrate on maximizing the values of pneumonia vaccination portion π to prevent and control pneumonia spreading.Figure 3
Effect of pneumonia vaccination onR2.In this subsection, as we see in Figure4, we have investigated the effect of pneumonia spreading rate β2 on the pneumonia effective reproduction number R2 by keeping the other rates as in Table 3. Figure 4 reflects that when the value ofβ2 increases, the pneumonia effective reproduction number R2 increases, and whenever the value of β2<0.022 implies R2<1. Therefore, public policymakers must concentrate on minimizing the values of pneumonia spreading rate β2 to minimize pneumonia effective reproduction numberR2.Figure 4
Effect of pneumonia transmission onR2.


#### 4.2.3. Effect of Pneumonia Treatment Rate on Infectious Population

In this subsection, as we see in Figure5, we have investigated the effect of κ in decreasing the number of pneumonia-only infectious populations. The figure reflects that when the values of κ increase, the number of pneumonia-only infectious population is going down. Therefore, public policymakers must concentrate on maximizing the values of treatment rate κ to pneumonia disease.Figure 5
Effect of treatment on pneumonia-infected population.


#### 4.2.4. Effect of Treatment Rates on HIV/AIDS Infectious Population

In this subsection, as we see in Figures6–8, respectively, we have investigated the effects of κ1,κ2, and κ3 in decreasing the number of acute HIV only, chronic HIV only, and AIDS-infected population, respectively. The figures reflect that when the values of κ1,κ2, and κ3 increase, the number of acute HIV only, chronic HIV only, and AIDS-infected population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rate of individuals to HIV/AIDS infection.Figure 6
Effect of treatment on acute HIV-infected population atβ1=0.5.Figure 7
Effect of treatment on chronic HIV-infected population atβ1=0.5.Figure 8
Effect of treatment on AIDS patients atβ1=0.5.


#### 4.2.5. Effect of HIV/AIDS Transmission Rate on Coinfectious Population

In this section, we see in Figure9 the effect of the spreading rate of HIV/AIDS β1 on the acute HIV-pneumonia coepidemic populationY7. The figure reflects that as the value of the transmission rate (β1) of HIV/AIDS is increased, the coepidemic population increases, which means the expansion of coepidemic of HIV/AIDS-pneumonia will increase. To control coepidemic of HIV/AIDS-pneumonia, decreasing the spreading rate of HIV/AIDS is important. Therefore, stakeholders must concentrate on decreasing the spreading rate of HIV/AIDS by using the treatment and appropriate method of prevention mechanism to bring down the expansion of coepidemic in the community.Figure 9
Effect ofβ1 on acute HIV-pneumonia coepidemic population.


#### 4.2.6. Effect of Treatment Rates on HIV/AIDS-Pneumonia Coepidemic Population

In this subsection, as we see in Figures10–12, we have investigated the effects of treatment rates σ1,σ2, and σ3 in decreasing the number of acute HIV and pneumonia, chronic HIV and pneumonia, and AIDS and pneumonia coinfectious population, respectively. The figures reflect that when the values ofσ1,σ2, and σ3 increase, the number of acute HIV-pneumonia, chronic HIV-pneumonia, and AIDS-pneumonia coepidemic population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rates of HIV/AIDS-pneumonia coepidemic population.Figure 10
Effect of treatment on acute HIV and pneumonia coepidemic.Figure 11
Effect of treatment on chronic HIV and pneumonia.Figure 12
Effect of treatment on AIDS and pneumonia coepidemic.


## 4.1. Sensitivity Analysis

Definition. The normalized forward sensitivity index of a variableR3that depends differentiably on a parameter p is defined as SIp=∂R3/∂p∗p/R3 [18].Sensitivity indices allow us to examine the relative importance of different parameters in pneumonia and HIV/AIDS spread and prevalence. The most sensitive parameter has the magnitude of the sensitivity index larger than that of all other parameters. We can calculate the sensitivity index in terms ofR1and R2sinceR3=maxR1,R2. Sensitivity analysis results and the numerical simulation are given in this section with parameters values given in Table 3 where N0 is the total number of the initial population of complete model (4).Table 3
Parameter values used for the full HIV/AIDS-pneumonia coepidemic model simulation.
ParameterValueSourceΛ0.0413∗N0Estimatedd0.02Estimatedδ10.498[7]δ20.08[7]δ30.2885[6]δ40.3105[6]ψ11.1Assumedψ21.2Assumedψ31.4Assumedν1AssumeddP0.1[16]θ0.1[18]dA0.333[6]π0.2[18]τ10.0025[18]ϵ0.002[18]dAP0.42Assumedκ0.2[18]κ10.2[7]κ20.15[7]κ30.13Assumedσ10.498[7]σ20.08[7]σ30.230Assumedβ1Variable[6]β2Variable[6]ρ1,ρ2,ρ3,ω1,ω2,ω31.2,1,1,1,1,1AssumedUsing the values of parameters in Table3, the sensitivity indices are calculated in Tables 4 and 5.Table 4
Sensitivity indices ofR3=R1.
Sensitivity indexValuesSIβ1+1SIρ1+0.6134SIδ1- 0.0639SId-0.3150SIκ1-0.1371SIκ2-0.0264SIδ2-0.0141Table 5
Sensitivity indices ofR3=R2.
Sensitivity indexValuesSIΛ+1SIβ2+1SId-0.4421SIκ-0.6559SIdP-0.3852SIπ-0.3852SIϵ-0.3852SIτ1-0.3852In this paper, with parameter values in Table3, we have got R1=1.386 at β1=2 implies HIV/AIDS spreads in the community and also we have got the indices as shown in Table 4. Here, sensitivity analysis shows that the human recruitment rate Λ and HIV/AIDS spreading rate β1 have the highest impact on the basic reproduction number of HIV/AIDS (R1).Similarly, with parameter values in Table3, we have got R2=9.69 at β2=0.2 imply that pneumonia spreads throughout the community and also we have got the indices as shown in Table 4. Here, sensitivity analysis shows that the foremost sensitive positive parameters are the human recruitment rate Λ and the pneumonia spreading rate β2. The foremost sensitive negative parameter is treatment rate of pneumonia κ which is inversely related to the effective reproduction number R2, i.e., a smaller amount of increase in this parameter value will lead to a greater amount of reduction in the effective reproduction number while a smaller amount of decrement will cause a big increment in the basic reproduction number. Epidemiologically, the most sensitive parameters to R1 and R2 which can be controlled through interventions and preventions are found to beβ1 andβ2, respectively.


## 4.2. Numerical Analysis

In this section, numerical simulation is performed for complete HIV/AIDS-pneumonia coepidemic model (4). With ode45, we have checked the effect of some parameters in the spreading as well as for the control of pneumonia only, HIV/AIDS only, and coepidemic of HIV/AIDS and pneumonia. The parameter values put forward in Table 3 are used for numerical simulation. In the numerical simulation part, we investigated the stability of the endemic equilibrium point of complete model (4), parameter effects on the reproduction numbers, and the impact of treatment mainly on dually infected individuals in the community.


### 4.2.1. Local Stability of Endemic Equilibrium Point of Complete Model (4)

Figure2 shows that in the long run (after 50 years), the solutions of dynamical system (4) will be converging to its endemic equilibrium point, i.e., the endemic equilibrium point is locally asymptotically stable whenever
(33)R3=maxR1,R2=max1.386,9.69=9.69>1.Figure 2
Local stability of endemic equilibrium point of the coepidemic model (4) whenever R1=1.386 at β1=2 and R2=9.69 at β2=0.2.


### 4.2.2. Effect of Parameters on the Threshold ParameterR2

In this subsection, as we see in Figure3, we have investigated the effect of pneumonia vaccination portion π on the pneumonia effective reproduction number R2. The figure reflects that when the value ofπ increases, the pneumonia effective reproduction number is going down, and whenever the value of π>0.64 imply R2<1. Therefore, public policymakers must concentrate on maximizing the values of pneumonia vaccination portion π to prevent and control pneumonia spreading.Figure 3
Effect of pneumonia vaccination onR2.In this subsection, as we see in Figure4, we have investigated the effect of pneumonia spreading rate β2 on the pneumonia effective reproduction number R2 by keeping the other rates as in Table 3. Figure 4 reflects that when the value ofβ2 increases, the pneumonia effective reproduction number R2 increases, and whenever the value of β2<0.022 implies R2<1. Therefore, public policymakers must concentrate on minimizing the values of pneumonia spreading rate β2 to minimize pneumonia effective reproduction numberR2.Figure 4
Effect of pneumonia transmission onR2.


### 4.2.3. Effect of Pneumonia Treatment Rate on Infectious Population

In this subsection, as we see in Figure5, we have investigated the effect of κ in decreasing the number of pneumonia-only infectious populations. The figure reflects that when the values of κ increase, the number of pneumonia-only infectious population is going down. Therefore, public policymakers must concentrate on maximizing the values of treatment rate κ to pneumonia disease.Figure 5
Effect of treatment on pneumonia-infected population.


### 4.2.4. Effect of Treatment Rates on HIV/AIDS Infectious Population

In this subsection, as we see in Figures6–8, respectively, we have investigated the effects of κ1,κ2, and κ3 in decreasing the number of acute HIV only, chronic HIV only, and AIDS-infected population, respectively. The figures reflect that when the values of κ1,κ2, and κ3 increase, the number of acute HIV only, chronic HIV only, and AIDS-infected population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rate of individuals to HIV/AIDS infection.Figure 6
Effect of treatment on acute HIV-infected population atβ1=0.5.Figure 7
Effect of treatment on chronic HIV-infected population atβ1=0.5.Figure 8
Effect of treatment on AIDS patients atβ1=0.5.


### 4.2.5. Effect of HIV/AIDS Transmission Rate on Coinfectious Population

In this section, we see in Figure9 the effect of the spreading rate of HIV/AIDS β1 on the acute HIV-pneumonia coepidemic populationY7. The figure reflects that as the value of the transmission rate (β1) of HIV/AIDS is increased, the coepidemic population increases, which means the expansion of coepidemic of HIV/AIDS-pneumonia will increase. To control coepidemic of HIV/AIDS-pneumonia, decreasing the spreading rate of HIV/AIDS is important. Therefore, stakeholders must concentrate on decreasing the spreading rate of HIV/AIDS by using the treatment and appropriate method of prevention mechanism to bring down the expansion of coepidemic in the community.Figure 9
Effect ofβ1 on acute HIV-pneumonia coepidemic population.


### 4.2.6. Effect of Treatment Rates on HIV/AIDS-Pneumonia Coepidemic Population

In this subsection, as we see in Figures10–12, we have investigated the effects of treatment rates σ1,σ2, and σ3 in decreasing the number of acute HIV and pneumonia, chronic HIV and pneumonia, and AIDS and pneumonia coinfectious population, respectively. The figures reflect that when the values ofσ1,σ2, and σ3 increase, the number of acute HIV-pneumonia, chronic HIV-pneumonia, and AIDS-pneumonia coepidemic population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rates of HIV/AIDS-pneumonia coepidemic population.Figure 10
Effect of treatment on acute HIV and pneumonia coepidemic.Figure 11
Effect of treatment on chronic HIV and pneumonia.Figure 12
Effect of treatment on AIDS and pneumonia coepidemic.


## 4.2.1. Local Stability of Endemic Equilibrium Point of Complete Model (4)

Figure2 shows that in the long run (after 50 years), the solutions of dynamical system (4) will be converging to its endemic equilibrium point, i.e., the endemic equilibrium point is locally asymptotically stable whenever
(33)R3=maxR1,R2=max1.386,9.69=9.69>1.Figure 2
Local stability of endemic equilibrium point of the coepidemic model (4) whenever R1=1.386 at β1=2 and R2=9.69 at β2=0.2.


## 4.2.2. Effect of Parameters on the Threshold ParameterR2

In this subsection, as we see in Figure3, we have investigated the effect of pneumonia vaccination portion π on the pneumonia effective reproduction number R2. The figure reflects that when the value ofπ increases, the pneumonia effective reproduction number is going down, and whenever the value of π>0.64 imply R2<1. Therefore, public policymakers must concentrate on maximizing the values of pneumonia vaccination portion π to prevent and control pneumonia spreading.Figure 3
Effect of pneumonia vaccination onR2.In this subsection, as we see in Figure4, we have investigated the effect of pneumonia spreading rate β2 on the pneumonia effective reproduction number R2 by keeping the other rates as in Table 3. Figure 4 reflects that when the value ofβ2 increases, the pneumonia effective reproduction number R2 increases, and whenever the value of β2<0.022 implies R2<1. Therefore, public policymakers must concentrate on minimizing the values of pneumonia spreading rate β2 to minimize pneumonia effective reproduction numberR2.Figure 4
Effect of pneumonia transmission onR2.


## 4.2.3. Effect of Pneumonia Treatment Rate on Infectious Population

In this subsection, as we see in Figure5, we have investigated the effect of κ in decreasing the number of pneumonia-only infectious populations. The figure reflects that when the values of κ increase, the number of pneumonia-only infectious population is going down. Therefore, public policymakers must concentrate on maximizing the values of treatment rate κ to pneumonia disease.Figure 5
Effect of treatment on pneumonia-infected population.


## 4.2.4. Effect of Treatment Rates on HIV/AIDS Infectious Population

In this subsection, as we see in Figures6–8, respectively, we have investigated the effects of κ1,κ2, and κ3 in decreasing the number of acute HIV only, chronic HIV only, and AIDS-infected population, respectively. The figures reflect that when the values of κ1,κ2, and κ3 increase, the number of acute HIV only, chronic HIV only, and AIDS-infected population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rate of individuals to HIV/AIDS infection.Figure 6
Effect of treatment on acute HIV-infected population atβ1=0.5.Figure 7
Effect of treatment on chronic HIV-infected population atβ1=0.5.Figure 8
Effect of treatment on AIDS patients atβ1=0.5.


## 4.2.5. Effect of HIV/AIDS Transmission Rate on Coinfectious Population

In this section, we see in Figure9 the effect of the spreading rate of HIV/AIDS β1 on the acute HIV-pneumonia coepidemic populationY7. The figure reflects that as the value of the transmission rate (β1) of HIV/AIDS is increased, the coepidemic population increases, which means the expansion of coepidemic of HIV/AIDS-pneumonia will increase. To control coepidemic of HIV/AIDS-pneumonia, decreasing the spreading rate of HIV/AIDS is important. Therefore, stakeholders must concentrate on decreasing the spreading rate of HIV/AIDS by using the treatment and appropriate method of prevention mechanism to bring down the expansion of coepidemic in the community.Figure 9
Effect ofβ1 on acute HIV-pneumonia coepidemic population.


## 4.2.6. Effect of Treatment Rates on HIV/AIDS-Pneumonia Coepidemic Population

In this subsection, as we see in Figures10–12, we have investigated the effects of treatment rates σ1,σ2, and σ3 in decreasing the number of acute HIV and pneumonia, chronic HIV and pneumonia, and AIDS and pneumonia coinfectious population, respectively. The figures reflect that when the values ofσ1,σ2, and σ3 increase, the number of acute HIV-pneumonia, chronic HIV-pneumonia, and AIDS-pneumonia coepidemic population is going down, respectively. Therefore, public policymakers must concentrate on maximizing the values of treatment rates of HIV/AIDS-pneumonia coepidemic population.Figure 10
Effect of treatment on acute HIV and pneumonia coepidemic.Figure 11
Effect of treatment on chronic HIV and pneumonia.Figure 12
Effect of treatment on AIDS and pneumonia coepidemic.


## 5. Discussion

In Section1, we reviewed and introduced the epidemiology of HIV/AIDS, pneumonia, and HIV/AIDS-pneumonia coepidemic. In Section 2, we construct the compartmental HIV/AIDS-pneumonia coepidemic dynamical system using an ordinary differential equation and we partitioned it into twelve distinct compartments. In Section 3, we analyzed the model qualitatively. To study the qualitative behavior of complete model (4), first, we split the complete model into two, which are HIV/AIDS-only and pneumonia-only models. The qualitative behaviors, i.e., the positivity of future solutions of the models, boundedness of the dynamical system, disease-free equilibrium points, basic reproduction numbers, endemic equilibriums, stability analysis of disease-free equilibrium points, stability analysis of endemic equilibrium points, bifurcations analysis of pneumonia-only model and the complete HIV/AIDS-pneumonia coepidemic model, and sensitivity analysis of reproduction numbers of HIV/AIDS-only and pneumonia-only model, are analyzed in their respective order, and numerically, we experimented on the stability of endemic equilibrium point of the HIV/AIDS-pneumonia coepidemic model, effect of basic parameters in the expansion or control of pneumonia only, HIV/AIDS only, and HIV/AIDS-pneumonia coepidemic infections and parameter effects on the infected population. From the result, we conclude that increasing both the pneumonia treatment rate and pneumonia vaccination portion rate has a great contribution to bringing down pneumonia infection as well as the coepidemic in the community. Similarly, increasing the HIV/AIDS treatment rates also has a contribution to minimizing the expansion of HIV/AIDS infection. The coepidemic treatment rates also influence minimizing coepidemic population if its value is increased. The other result obtained in this section is that decreasing the transmission rates has a great influence of controlling coepidemic in the population.


## 6. Conclusion

A realistic compartmental mathematical model on the spread and control of HIV/AIDS-pneumonia coepidemic incorporating pneumonia vaccination and treatment for both infections are available at each stage of the infection in a population constructed and analyzed. We have shown the positivity and boundedness of the complete HIV/AIDS-pneumonia coepidemic model. Using center manifold theory, we have shown that the pneumonia-only infection and the complete HIV/AIDS-pneumonia coepidemic models undergo the phenomenon of backward bifurcation whenever their corresponding effective reproduction numbers are less than one. The complete model has a disease-free equilibrium that is locally asymptotically stable whenever the maximum of the reproduction numbers of the two submodels described above is less than one. Numerical simulation shows that the complete HIV/AIDS-pneumonia coepidemic model endemic equilibrium point is locally asymptotically stable when its effective reproduction number is greater than one. These results have important public health implications, as they govern the elimination and/or persistence of the two diseases in a community. By analyzing the various associated reproduction numbers, we have shown that the impact of some parameters changes on the associated reproduction numbers to give future recommendations for the stakeholders in the community. From the numerical result, we have got the complete model reproduction number isR3=maxR1,R2=max1.386,9.69=9.69 at β1=2 and β2=0.2. From our numerical result, we recommend that public policymakers must concentrate on maximizing the values of pneumonia vaccination portion and treatment rate of individuals to pneumonia disease. Finally, some of the main epidemiological findings of this study include pneumonia vaccination and treatment against disease has the effect of decreasing the pneumonia and coepidemic disease expansion and prevalence and reducing the progression rate of HIV infection to the AIDS stage and the HIV/AIDS prevalence.


### 6.1. Limitation of the Study

Due to conflict in our country Ethiopia, it is difficult to incorporate experimental data in the study.


## 6.1. Limitation of the Study

Due to conflict in our country Ethiopia, it is difficult to incorporate experimental data in the study.



---

*Source: 3105734-2022-01-11.xml*