HIV/AIDS-Pneumonia Coinfection Model with Treatment at Each Infection Stage: Mathematical Analysis and Numerical Simulation

In the paper, we have considered a nonlinear compartmental mathematical model that assesses the effect of treatment on the dynamics of HIV/AIDS and pneumonia coinfection in a human population at different infection stages. Our model revealed that the disease-free equilibrium points of the HIV/AIDS and pneumonia submodels are both locally and globally asymptotically stable whenever the associated basic reproduction numbers (RH and RP) are less than unity. Both the submodel endemic equilibrium points are locally and globally asymptotically stable whenever the associated basic reproduction numbers (RP and RH) are greater than unity. The full HIV/AIDS-pneumonia coinfection model has both locally and globally asymptotically stable disease-free equilibrium points whenever the basic reproduction number of the coinfection model ðRHPÞ is less than unity. Using standard values of parameters collected from different kinds of literature, we found that the numerical values of the basic reproduction numbers of the HIV/AIDS-only submodel and pneumonia-only submodel are 17 and 7, respectively, and the basic reproduction number of the HIV/AIDS-pneumonia coinfection model is max f7, 17g = 17. Applying sensitive analysis, we identified the most influential parameters to change the behavior of the solution of the considered coinfection dynamical system are the HIV/AIDS and pneumonia transmission rates β1 and β2, respectively. The coinfection model was numerically simulated to investigate the stability of the coinfection endemic equilibrium point, the impacts of transmission rates, and treatment strategies for HIV/AIDS-only, pneumonia-only, and HIV/AIDS-pneumonia coinfected individuals. Finally, we observed that numerical simulations indicate that treatment against infection at every stage lowers the rate of infection or disease prevalence.


Introduction
Infectious diseases are a clinically evident illness, and commonly, they have a great influence on the human population. They are induced by a pathogenic microbial agent such as bacterial, viral, fungal, parasitic, or it can be toxic proteins, called prions. The most common ones are tuberculosis and pneumonia caused by bacteria, HIV, and influenza caused by the virus [1,2].
HIV/AIDS is a global health issue affecting approximately 70 million people worldwide causing significant morbidity and mortality [3]. Over two-thirds of people living with HIV live in the Sub-Saharan African region [4]. Human immunodeficiency virus (HIV) is a retrovirus virus which attacks and weakens the human body immunity and the central nervous system and if untreated it continues to multiply into the host until it reaches the peak leading into a very serious disease called AIDS, the stage where the symptoms of the disease occur frequently [5][6][7]. HIV is transmitted through sexual intercourse, needle sharing, and direct contact of blood or other body fluids containing the virus and mother to child during childbirth [1,8]. According to the center for disease control and prevention (CDC), when individuals get HIV and do not receive treatment known as antiretroviral therapy (ART), they will typically progress through three stages of disease: acute HIV infection, clinical latency (HIV inactivity or dormancy), and acquired immunodeficiency syndrome (AIDS).
Pneumonia is one of airborne infectious disease caused by caused by bacteria, viruses, fungi, or parasites which attacks the human lungs or alveoli [9][10][11]. Most of the time pneumonia affects older adults, babies, and people with other diseases or impaired immune systems worldwide. Its most common cause is the Streptococcus pneumoniae, also known as pneumococcus [10,12,13]. The basic controlling strategies of pneumonia infection are treatment and vaccination interventions [11].
A coinfection is the infection of a host with two or more different pathogens or different strains of the same pathogens, leading to coexistence of strains (pathogens) at population level [14]. Mathematical modeling of infectious diseases such as coinfection of HIV/AIDS and opportunistic infection is regarded as a fundamental tool in understanding the dynamics and helpful in the decision-making processes regarding intervention strategies and measures required for disease elimination and/or control [15,16]. Some mathematical models have been used to investigate the transmission dynamics of coinfection of two or more diseases where HIV/AIDS and pneumonia coinfection is among the diseases that infect a large number of individuals worldwide. HIVinfected persons are particularly susceptible to the development of severe pneumococcal disease, even in the setting of combination antiretroviral therapy (ART) [4,7,12].
Over the past, mathematical model shave been developed to analyze the population dynamics of HIV/AIDS and pneumonia single infections. Authors in [3] developed and analyzed only a sex-structured population model and studied the HIV infection trends in males and females. Their model assumed that the main mode of HIV transmission is heterosexual. They showed that prevention of HIV infection still remains the most important way of controlling further spread in the community. HIV/AIDS patients under ART treatment are possibly capable of helping the eradication of HIV by convincing their sexual partners of the need to adhere to protection via use of preexposure prophylaxis (PrEP) or any other protection means and ART treatment. Huo and Chen [5] developed and analyzed a mathematical analysis to study the spread of HIV/AIDS with treatment at different stages. Their results show that early treatment for individuals in asymptomatic stage of HIV infection or the pre-AIDS stage is very important. Rahman [2] have analyzed a seven dimension in the living organism HIV model with optimal control of in-host HIV dynamics using different control strategies such as three drug combinations, that is, FIs, RTIs, and PIs to determine the optimal treatment regime. From their results, they recommend that RTIs be used as initial therapy for HIV and FI should be introduced to the patient after the RTIs but should never be used alone.
Mbabazi et al. [13] formulated a mathematical model to study the global stability of pneumococcal pneumonia with awareness, and saturated treatment is presented. Their results showed that the family of decaying curves could help in providing mechanisms to design awareness strategies for containing pneumococcal pneumoniae threshold parameter could be reduced to less than unity if antibiotic resistance awareness and treatment are implemented simultaneously to ensure eradication of pneumococcus bacteria; thus, spread of pneumococcus pneumonia in the population will die out. Tilahun et al. [17] proposed and analyzed a nonlinear mathematical model for the transmission dynamics of pneumonia disease in a population of varying size with optimal control of pneumonia disease and cost-effective strategies. Their cost-effectiveness analysis of the adopted control strategies showed that the combination of prevention and treatment is the most cost-effective intervention strategies to fight the pneumonia pandemic. Ndelwa et al. [18] developed a mathematical model and analyzed treatment and screening strategies on pneumonia infection, and from their numerical results, they concluded screening and treating at the same time can eradicate the pneumonia epidemic from the community.
Nwankwo and Okuonghae [19] formulated and analyzed a mathematical model for the transmission dynamics of syphilis and HIV coinfection in a community to assess the impact of treatment of syphilis on the coendemicity of both diseases in a population where treatment for HIV is not readily available (or easily accessible) but with syphilis treatment sufficiently available. Their syphilis-only model and the full coinfection models undergo the phenomenon of backward bifurcation due to syphilis reinfection after recovery from a previous syphilis infection. They have got the treatment of primary and secondary syphilis in both singly and dually infected individuals; especially with high treatment rates for primary syphilis, this will result in a reduction in the incidence of HIV and its coinfection with syphilis in the community. Kaur et al. [6] formulated and analyzed a mathematical model for HIV/AIDS-TB coinfection with screening and treatment of both HIV and TB infective. Their numerical results suggested that the rates of transmission of both TB and HIV should be decreased, as an increase causes a rise in the number of infective at the equilibrium level.
However, most microbiology, epidemiology, and medical sources like [7,12] shows the coexistence of HIV/AIDS and pneumonia infection but coinfection mathematical models of HIV/AIDS and pneumonia are rare in literature yet the coexistence between the two infections exist. In our review of literatures, we have got one mathematical model of HIV/AIDS and pneumonia infection in a population and we used it as initial literature reviewed as follows. According to Nthiiri et al. [4] maximum protection against the HIV/AIDS-pneumonia coinfection was analyzed where the maximum protection against HIV/AIDS and the maximum protection against pneumonia was the main concern of their project. In their model, they did not considered maximum protection against the coinfection rather considered maximum protection against single infections. Also, they did not considered treatment on either the submodels or the coinfection. Their analysis found that when protection is high, the number of HIV/AIDS and pneumonia cases is low.
Our paper, therefore, presents a mathematical model describing the transmission dynamics of HIV/AIDS and pneumonia coinfection in a population where treatment 2 Journal of Applied Mathematics for both HIV/AIDS and pneumonia are available in the community. Basically, the model will be used to evaluate the effect of treatment at every infection stage of either the single infected individuals with HIV/AIDS or pneumonia or HIV/AIDS and pneumonia coinfection as a control strategy for minimizing incidences of coinfections in the target population. In this work, we applied the center for disease control and prevention (CDC) human immunodeficiency virus (HIV) infection stages and the control measure treatment at each stage of the single infections and coinfection model. We have checked this case has never been done before. We discussed the effects of treatment for single infected individuals with either HIV/AIDS or pneumonia and the coinfected patient with HIV/AIDS and pneumonia at each infection stage. The paper is organized as follows. The model is formulated in Section 2 and is analyzed in Section 3. Sensitivity analysis, numerical results, and discussion are carried out in Section 4. Finally, conclusion and limitations of the study are carried out in Sections 5 and 6, respectively.

The Mathematical Model
According to the three center for disease control and prevention (CDC) HIV infection stages, we have divide the total population NðtÞ in to eleven mutually exclusive compartments stated in Table 1, so that NðtÞ = SðtÞ + I P ðtÞ + H 1 ðtÞ + H 2 ðtÞ + H 3 ðtÞ + C 1 ðtÞ + C 2 ðtÞ + C 3 ðtÞ + T P ðtÞ + T H ðtÞ + TðtÞ. We assumed that coinfected individuals can only transmit either pneumonia or HIV but not both infections simultaneously. Individuals acquire HIV infection following effective contacts with people infected with HIV only, (H 1 and H 2 classes), at the rate given by where ρ > 1 is the modification parameter accounting for the assumed increased infectivity due to chronic HIV infected than acute HIV-infected one and β 1 is the HIV transmission rate. Also, individuals acquire pneumonia infection from those in the I P , C 1 , C 2 , and C 3 infectious at the rate where ω 3 > ω 2 > ω 1 are modification parameters accounting for the assumed increased infectivity due to coinfections and β 2 is the pneumonia transmission rate. The derivation of the model differential equations is given in "Appendix A." 2.1. Flow Chart of the Dynamical System. Here, parameter descriptions in Table 2, state variable descriptions in Table 1 above, and based on the model assumptions that led to the formulation of the model (as stated in"Appendix A"), the flow diagram for the transmission dynamics of HIV/AIDS and pneumonia coinfection is given by Figure 1.

Dynamical System of HIV/AIDS-Pneumonia
Coinfection. Based on Figure 1 above, the dynamical system of HIV/AIDS-pneumonia coinfection becomes With initial conditions,    Journal of Applied Mathematics The sum of all the differential equations in (3) is given by Since model (3) monitors human population, it is assumed that all variables and parameters are nonnegative. The dynamics of model (3) will be analyzed in the following invariant region: Then, we have proved the positivity and boundedness of solutions of (3) in Ω in "Appendix B."

Mathematical Model Analysis
Before we analyzed the HIV/AIDS-pneumonia coinfection model (3), it is useful to gain some background about the HIV-only submodel and pneumonia-only submodel transmission dynamics.
3.1. HIV/AIDS-Only Submodel Analysis. The HIV submodel of (3) (obtained by setting I P = C 1 = C 2 = C 3 = T P = T = 0Þ is given by where the total population is N 1 ðtÞ = SðtÞ +H 1 ðtÞ + H 2 ðtÞ + H 3 ðtÞ + T H ðtÞ and the HIV/AIDS force of infection is given by λ H = β 1 ðH 1 + ρH 2 Þ with initial conditions: The sum of all the differential equations in (7) is obtained as Consider the region It is easy to show that the set Ω 1 is positively invariant and a global attractor of all positive solution of submodel. Hence, it is sufficient to consider the dynamics of model (7) in Ω 1 . In this region, the model is considered epidemiologically and mathematically well posed.

Disease-Free Equilibrium Point of the HIV/AIDS-Only
Submodel. The disease-free equilibrium point (DFE) of the HIV/AIDS-only submodel (7) is denoted by E 0 H = ðS 0 , H 0 1 , H 0 2 , H 0 3 , T 0 H Þ and obtained by making the right hand side of the system as zero and setting all the infectious classes and treatment class to zero as H 1 3.1.2. The Basic Reproduction Number of the Submodel. The basic reproduction number of HIV/AIDS-infected individuals denoted by R H is defined as the average number of secondary infections produced by a single HIV/AIDS infectious individual introduced in a wholly susceptible population during his or her entire infectious period [14,20]. This definition is given for the dynamical system that represents the spread of infection in a population. We calculate the basic reproduction number by using the nextgeneration operator method on the dynamical system (7). The basic reproduction number is obtained by taking the largest (dominant) eigenvalue (spectral radius) of the matrix: [21,22] where F i is the rate of appearance of new infection in compartment i, ν i is the transfer of infections from one compartment i to another, and E H 0 is the disease-free equilibrium point.
After some calculations we have got F = Then, the spectral radius (reproduction number R H ) of FV −1 of the HIV/AIDS subdynamical system (7) (7). Let an arbitrary equilibrium point of a HIV/AIDS-only dynamical system (7) is denoted by Þ be the associated infection rate ("force of infection") at endemic equilibrium point. After some calculations, we have got λ * Theorem 2. The submodel (7) has a unique endemic equilibrium point if and only if R H > 1.

Global Asymptotic Stability (GAS) of the Disease-Free Equilibrium Point
Theorem 3. The disease-free equilibrium point E 0 H = ðΛ/μ, 0 , 0, 0, 0Þ of the dynamical system (7) is globally asymptotically Proof. To show global stability of the DFE applied Lyapunov function method as [13,23].
Thus, by LaSalle's invariance principle [24], it implies that the disease-free equilibrium point Proof. Let the Lyapunov function V : Here, we have dL/dt = a 1 ððdS/dtÞ − ððS * /SÞðdS/dtÞÞÞ + a 2 ððd At the endemic equilibrium point E * H = ðS * , H * 1 , H * 2 , H * 3 , T * H Þ, we obtain from the system (7): Journal of Applied Mathematics Then, we obtained Choose a 1 , a 2 , and a 3 such that the expressions in the brackets vanish and take a 1 = a 2 and solve for a 3 , by making expressions in the closed bracket zero, we obtain as Grouping some terms in the expression above yields Using the arithmetic-geometric mean inequality property, we have 2 − ðS/S * Þ − ðS * /SÞ ≤ 0 Therefore, we conclude by LaSalle's invariance principle [24]

Pneumonia Submodel Analysis.
We have the pneumonia submodel of (3) when where the total population is N 2 ðtÞ = SðtÞ + I p ðtÞ + T p ðtÞ and the pneumonia force of infection is given by λ P = β 2 I P with initial conditions The sum of all the differential equations in (18) above is obtained as Consider the region Ω 2 = ðS, I P , T P Þ ∈ ℝ 3 + , N 2 ≤ Λ/μ È É . It is easy to show that the set Ω 2 is positively invariant and a global attractor of all positive solution of submodel (18). Hence, it is sufficient to consider the dynamics of model (18) in Ω 1 . In this region, the model is epidemiologically and mathematically well posed.  (18) is obtained by making the right hand side of the system as zero and setting all the infectious classes and treatment class to zero as I P = T P = 0 we have got S 0 = Λ/μ such that E 0 P = ðΛ/μ, 0, 0Þ.

The Reproduction Number of the Pneumonia-Only
Submodel. We calculate the basic reproduction number denoted by R P using the van den Driesch and Warmouth next-generation matrix approach from [22]. The basic reproduction number is obtained by taking the largest (dominant) eigenvalue (spectral radius) of the matrix: where F i is the rate of appearance of new infection in compartment i , ν i is the transfer of infections from one compartment i to another, and E 0 P is 7 Journal of Applied Mathematics the disease-free equilibrium point. The reproduction number R P of the pneumonia-only dynamical system (18) is obtained by rearrange the differential equation of the dynamical system (18)  Proof. The local stability of the disease-free equilibrium of the system (18) can be studied from its Jacobian matrix at the disease-free equilibrium point E 0 P = ðS 0 , I P 0 , T P 0 Þ = ðΛ/ μ, 0, 0Þ and Routh-Hurwitz stability criteria. Then, the Jacobian matrix of the dynamical system (18) at E 0 P = ðΛ/ μ, 0, 0Þ is given by Then, the characteristic equation of the above Jacobian matrix is given by Therefore, since all the eigenvalues of the characteristics polynomial of the system (18) are negative for R P < 1, the disease-free equilibrium point of the system (18) is locally asymptotically stable.

Existence of Endemic Equilibrium Point of the Pneumonia Submodel.
Before investigating the global asymptotic stability of the DFE, it is instructive to determine the number of endemic equilibrium solutions of the model (18). Let an arbitrary endemic equilibrium point of the pneumonia-only dynamical system (18) be denoted by E P * = ðS * , I P * , T P * Þ. Now, after some calculations, we have got a unique endemic equilibrium point E P Theorem 7. The model (18) has a unique endemic equilibrium point if R P > 1.

Local and Global Stabilities of Endemic Equilibrium
Point. The following theorem studies the local stability of 8 Journal of Applied Mathematics the endemic equilibrium. The result is obtained by means of the Routh-Hurwitz stability criteria.

Theorem 9.
The endemic equilibrium point of the system (18) is locally asymptotically stable if R P > 1.
Proof. To show the local stability of the endemic equilibrium point, we use the method of the Jacobian matrix and Routh Hurwitz stability criteria. Then, the Jacobian matrix of the dynamical system (18) at the endemic equilibrium point and T * P = ðD 2 ½R P − 1Þ/ðd 2 m 2 + d 2 m 3 D 1 ½R P − 1Þ, D 2 = d 2 γΛD 1 for the system (18) whenever R P > 1 is given by Then, the characteristic equation of the above Jacobian matrix is given by Here, we apply the necessary condition of Routh-Hurwitz stability criteria since a 3 = 1 is positive in sign, indicating all a 2 , a 1 , and a 0 should be positive.
Hence, a 2 = −ðA 1 Thus, all the coefficients of the characteristic's polynomial a 0 , a 1 , a 2 , and a 3 have the same sign. Then, applying the Routh-Hurwitz criteria to determine the sign of the root without calculating their values of the root of the characteristics equation a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 and we have got the Routh-Hurwitz array that has no sign change, indicating the roots of the characteristics equation of the dynamical system (18) are negatives. Hence, the endemic equilibrium point E * P = ðS * , I * P , T * P Þ of the dynamical system (18) is locally asymptotically stable.
9 Journal of Applied Mathematics

Global Stability of Endemic Equilibrium of Pneumonia Submodel
Theorem 10. If R P > 1 the unique endemic equilibrium point E * P = ðS * , I * P , T * P Þ is globally asymptotically stable in the interior of Ω 2 .
Proof. Let us take the Lyapunov function L : ℝ 3 + ⟶ ℝ by LðS, I P , T P Þ = A 1 ½ðS − S * Þ + ðI P − I * P Þ+ ðT P − T * P Þ 2 + A 2 ½I P − I * P − I * P ln ðI P /I * P Þ + A 3 ½T P − T * P 2 for positive constants A 1 = 1/2, A 2 = ðδ P + 2μÞ/β 2 and A 3 = ðδ P + 2μÞ/2γ. Then, L is C 1 on the interior of Ω 2 , E * P = ðS * , I * P , T * P Þ is the global minimum of L on Ω 2 , and LðS * , I * P , T * P Þ = 0. Then, the time derivative of L is given by Using expressions at the endemic equilibrium, we have got Hence, dL/dt is negative. Note that dL/dt = 0 if and only if S = S * , I P = I * P and T P = T * P . Therefore, the largest compact invariant set in fðS, I P , T P Þ ∈ Ω 2 : dL/dt = 0g is the singleton {E * P } set where E * P the endemic equilibrium point. Then, by LaSalle's invariant principle, [24] then implies that E * P is globally asymptotically stable in the interior of Ω 2 if R P > 1.

Analysis of the HIV-Pneumonia Coinfection Model (3).
In this section, we analyze the main dynamical system (3). Epidemiologically, the HIV/AIDS-pneumonia coinfection dynamical system (3) (3) is obtained by setting all the infectious classes and treatment classes to zero and making the right hand side of the system as zero. Then, the disease-free equilibrium of the HIV/AIDS coinfection model is E 0 HP = ð Λ/μ, 0, 0, 0, 0, 0, 0, 0, 0, 0Þ.

The Basic Reproduction Number ðR HP Þ of the HIV-Pneumonia Coinfection
Model. The basic reproduction number of the dynamical system (3) by applying the nextgeneration operator method [21,22] is the largest (dominant) eigenvalue (spectral radius) of the matrix: where F i is the rate of appearance of new infection in compartment i, ν i is the transfer of infections from one compartment i to another, and E 0 HP is the disease-free equilibrium point. Here, we obtained the following matrices: Then, using Mathematica, we have got the spectral radius (dominant eigenvalue) of the matrix FV −1 is R HP = max fðΛβ 1 /μðμ + γ 1 + α 1 ÞÞ + ðΛρα 1 β 1 /ðμðμ + γ 1 + α 1 Þðμ + γ 2 + α 2 ÞÞÞ, Λβ 2 /μðμ + γ + δ P Þg where R P = Λβ 2 /μðμ + γ + δ P Þ is the basic reproduction number for pneumonia infection and R H = ðΛβ 1 /μðμ + γ 1 + α 1 ÞÞ + ðΛρα 1 β 1 /ðμðμ + γ 1 + α 1 Þðμ + γ 2 + α 2 ÞÞÞ is the basic reproduction number for HIV/AIDS infection, and hence, R HP = max fR P , R H g is the basic reproduction number of the HIV and pneumonia co-infection.

Local Stability of the Disease-Free Equilibrium Point
Theorem 11. The disease-free equilibrium point of the model (3) above is locally asymptotically stable if R HP < 1 and unstable if R HP > 1.

Sensitivity Analysis of the Model Parameters and Numerical Simulations
Results of sensitivity analysis and the numerical simulation are given in this section, and the set of parameters used are given in Table 3 below. The full model (3) is now simulated, using the parameter estimates in Table 3 (unless otherwise stated), to assess the potential impact of treatment strategies against pneumonia and HIV/AIDS coinfection, as follows.

Sensitivity Analysis of the Model Parameters.
Sensitivity is performed to identify the most dominant parameters for the spreading out as well as control of infection in the community. To go through sensitivity analysis, we follow the technique described in [9,17]. Results of sensitivity analysis are given in this section, and the set of parameters used are given in Table 3 below where N 0 is the total number of the assumed initial population under consideration in the model numerical simulation part.
We are able to know how important each parameter is to the spread of the disease through sensitivity indices of R HP to all different parameters. We carried out the sensitivity analysis to determine the model robustness to parameter values.
The normalized forward sensitivity index of a variable R HP that depends differentiably on a parameter p is defined as SIðpÞ = ð∂R HP /∂pÞ * ðp/R HP Þ [9].
These sensitivity indices allow us to determine the relative importance of different parameters in pneumonia and HIV/AIDS transmission 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 of R H and R P since The sensitivity indices in terms of R HP = R H = ðΛβ 1 /μ ðμ + γ 1 + α 1 ÞÞ + ðΛρα 1 β 1 /ðμðμ + γ 1 + α 1 Þðμ + γ 2 + α 2 ÞÞÞand the sensitivity indices in terms of R HP = R P = Λβ 2 /μðμ + γ + δ P Þ are given in Tables 4 and 5, respectively.   Table 5: Sensitivity indices of R HP = R P .

Sensitivity index Value
Journal of Applied Mathematics Using the parameter values in Table 3, the sensitivity indexes are computed in Tables 4 and 5 as above.

Numerical Results and Discussion.
Using the data provided in Table 3 and different initial conditions, the numer-ical results are generated for the dynamics of model (3) using MATLAB numerical solver (ode45) which generate results for ode45. The ode45 was chosen because of its computational speed and increased level of accuracy for solving nonstiff ordinary differential equations and we simulated the  Table 3 above whose sources are mainly calculated from literatures in order to have realistic simulation results. In this analysis, we also discussed the effect of parameters change on the basic reproduction number graphically using MATLAB ode45 software.  (3) for different values of the associated reproduction thresholds R H and R P and we plot the graphs of the time verses infected population for different values of reproduction numbers. Figure 2 above was plotted using MATLAB ode45 program under consideration of the basic reproduction numbers being less than a unity and shows the behavior of the infectious classes of the HIV/AIDS-pneumonia coinfection model (3) at R H < 1 and R P < 1 (i.e., R HP < 1). The simulation given in Figure 2 above shows that each of infectious classes (I P , H 1 ,H 2 , H 3 , C 1 , C 2 , C 3 ) is converging to the disease-free equilibrium point of the model. This was obtained when R P =0.1445 at β 2 = 0:0000079 and R H = 0:1374 at β 1 = 0:00000289 with all other parameters are given as in Table 3. This indicates that the disease-free equilibrium point of the full HIV/AIDS-pneumonia coinfection model is globally asymptotically stable.     Table 3. The simulation in Figure 3 above shows the stability of the HIV/AIDS-pneumonia coinfection model endemic equilibrium point. From Figure 3 above, in the long run, the convergence of the solutions is observed at the values greater than 7 years. The plot shows that the HIV/AIDS-pneumonia coinfection model (3) endemic equilibrium point is locally asymptotically stable.

Reproduction Number Simulations with Variable
Parameter Values. Here, we have taken parameters from Table 3, and we have done numerical simulations of repro-duction numbers with variable parameter values using ode45 method, and we obtained figures from Figures 4-9. Thus, Figure 4 shows the pneumonia-only submodel reproduction number simulation at variable treatment rate, and from the graph, we see that pneumonia infection dies out at pneumonia treatment rate γ > 0:89. Figure 6 showed that pneumonia infection dies out whenever β 2 < 0:00003. Figure 7 shows that in the long run the HIV/AIDS    transmission decreases whenever β 1 < 0:00002. Figure 8 shows the HIV/AIDS reproduction number simulation at variable acute-infected treatment rate, and from the graph, we see that the HIV/AIDS-only submodel reproduction number is less than unity whenever γ 1 > 0:93 at β 1 = 0:00015, and Figure 9 shows the HIV/AIDS reproduction number simulation at variable chronic-infected treatment rate, and from the graph, we see that the HIV/AIDS-only submodel reproduction number is less than unity whenever the γ 2 > 0:75 at β 1 = 0:00015. Similarly, using parameter values from Table 3, the numerical simulation in Figure 5 shows that comparison of the basic reproduction numbers R H and R P at different values of β 1 and β 2 , respectively; also, from the simulations, we see that the basic reproduction number for HIV/AIDS is greater than the basic reproduction number of pneumoniaonly infection where HIV/AIDS-only submodel reproduction number is also the basic reproduction number of the coinfection model. Also, R H < 1 whenever β 1 < 0:000006 and R P < 1 if β 2 < 0:000014. Antiretroviral therapy (ART) is used to suppress the HIV virus and stop the progression to AIDS stage,     Figure 10 shows that when pneumonia treatment rate γ increases from 0.035 to 0.072 the pneumonia infection decreases. Figure 11 shows that when acute HIV infection treatment rate γ 1 increases from 0.81 to 0.91 the acute HIV infection decreases. Figure 12 shows that when chronic HIV infection treatment rate γ 2 increases from 0.59 to 0.89 the chronic HIV infection decreases. Figure 13 shows that when AIDS patients treatment rate γ 3 increases from 0.061 to 0.192 the AIDS patients decreases. Figure 14 shows that when acute HIV/AIDSpneumonia coinfection treatment rate ε 1 increases from 0.0.1 to 0.6 the acute HIV/AIDS-pneumonia coinfection decreases. Figure 15 shows that when chronic HIV/AIDSpneumonia coinfection treatment rate ε 2 increases from 0.0.34 to 0.69 chronic HIV/AIDS-pneumonia coinfection decreases. Figure 16 shows that when AIDS-pneumonia coinfection treatment rate ε 3 increases from 0.32 to 0.75 AIDS-pneumonia coinfection decreases. Thus, our model considers treatment at every infection stages of the full HIV/AIDS-pneumonia coinfection model; all numerical simulation graphs from Figures 10-16 show those effects of treatment on the infected population in the corresponding compartments, and also, all simulated curves show that the infected population in the compartment decreases whenever the corresponding treatment rate increases which is the finding of this work.

Conclusion
In this work, we formulated a mathematical model of eleven nonlinear differential equations on HIV/AIDS and pneumonia coinfection with the assumption of mass action incidence and incorporating treatment at each stage of the infection. We have shown the positivity and boundedness of the solutions of the model. The threshold parameter R HP was calculated and used to determine the conditions under which the HIV/AIDS and pneumonia could be transmitted and remained endemic in the population. We thus showed that three disease-free equilibrium points E 0 H , E 0 P , and E 0 HP , respectively, for the HIV submodel, pneumonia submodel, and full model are locally asymptotically stable when R HP < 1, i.e., R H < 1 and R P < 1. We also showed that the population with both HIV/AIDS and pneumonia infection have three endemic equilibrium points E * H , E * P , and E * HP respectively, for the HIVsub model, pneumonia submodel, and full model which were locally asymptotically stable when R HP > 1, i.e., R H > 1 or R P > 1. Global stability analysis of the submodels disease-free equilibrium points was established whenever R HP < 1, i.e., R H < 1 and R P < 1 . To investigate, the effect of treatment at each infected compartment was considered for both the submodels and the full model, namely, treatment of pneumonia infection, treatment of acute HIV infection, treatment of chronic HIV infection, treatment of AIDS patients, treatment of acute HIV/AIDSpneumonia coinfection, treatment of chronic HIV/AIDSpneumonia coinfection, and treatment of AIDS-pneumonia coinfection. We have showed the most sensitive parameters of our model which can be epidemiologically controlled are the HIV/AIDS transmission rate β 1 and pneumonia transmission rate β 2 so it is reasonable to recommend the use of intervention strategy for HIV/AIDS transmission in making β 1 less than 0:00014 and treatment for pneumonia transmission in making β 2 less than 0:00006:  where A 1 = −μR H , A 2 = β 1 S * − ðμ + γ 1 + α 1 Þ = ðβ 1 Λ/μR H Þ − ðμ + γ 1 + α 1 Þ = ðβ 1 Λ/μR H Þ − d 1 with d 1 = μ + γ 1 + α 1 , A 3 = −ðμ + γ 2 + α 2 Þ, A 4 = −ðμ + γ 3 + δ 1 Þ, A 5 = −β 1 Λ/μR H , and A 6 = −β 1 ρΛ/μR H .
Then, the characteristic equation of the above Jacobian matrix is given by