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

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 if ℛ 1 < 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 ℛ 2 < 1  and ℛ 3 = max{ℛ 1, ℛ 2} < 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 ℛ 3 = max{ℛ 1, ℛ 2} < 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 ℛ 3 = max{ℛ 1, ℛ 2} = max{1.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.

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 Section 2 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 population N(t) into twelve distinct classes as susceptible class to both HIV and pneumonia infections Y1(t), pneumonia vaccine class Y2 (t) , pneumonia-infected class Y3(t), acute HIV-infected class Y4(t), chronic HIV-infected class Y5(t), AIDS patient class Y6(t), acute HIV-pneumonia coepidemic class Y7(t), chronic HIV-pneumonia coepidemic class Y8(t), AIDS-pneumonia coepidemic class Y9(t), pneumonia treatment class Y10(t), HIV/AIDS treatment class Y11(t) entered from the three infection stages Y4(t), Y5(t), and Y6(t), and HIV/AIDS-pneumonia coepidemic treatment class Y12(t) entered from Y7(t), Y8(t), and Y9(t) cases such that

The susceptible class acquires HIV at the standard incidence rate given by 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 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:

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

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 θ

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

Individuals in a given compartment are homogeneous

Assume no HIV transmission from Y6(t) and Y9(t) classes due to their reduced daily activities

Individuals in each class are subject to natural mortality rate d

The human population is variable

We assumed there is no dual-infection transmission simultaneously

Assume HIV has no vertical transmission and pneumonia is not naturally recovered

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 Table 1, 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.

Parameter	Interpretations
d	Natural mortality rate
Λ	Human recruitment rate
δ 1	Development rate from acute HIV to chronic HIV infection
δ 2	Development rate from chronic HIV to AIDS stage
ϵ	The proportion of the serotype not covered by the vaccine
θ	Immunity loss rate
ψ 1	Alteration rate indicating acute HIV infection is more vulnerable to pneumonia
ψ 2	Alteration rate indicating chronic HIV infection is more vulnerable to pneumonia
ψ 3	Alteration rate indicating AIDS patient is more vulnerable to pneumonia
λ HC	HIV/AIDS standard incidence rate
λ PC	Pneumonia mass action incidence rate
δ 3	Development rate from acute HIV-pneumonia to chronic HIV-pneumonia coepidemics
δ 4	Development rate from chronic HIV-pneumonia to AIDS-pneumonia coepidemics
d P	Pneumonia death rate
d A	AIDS death rate
d AP	AIDS-pneumonia death rate
κ	Pneumonia infection treatment rate
κ 1	Acute HIV infection treatment rate
τ 1	Vaccination waning rate
κ 2	Chronic HIV infection treatment rate
κ 3	AIDS patients treatment rate
σ 1	Acute HIV-pneumonia coepidemic treatment rate
σ 2	Chronic HIV-pneumonia coepidemic treatment rate
σ 3	AIDS-pneumonia coepidemic treatment rate
β 1	Transmission rate of HIV
β 2	Transmission rate of pneumonia

Table 2

Definitions of variables.

Variables	Definitions
Y 1	Vulnerable to both HIV and pneumonia class
Y 2	Pneumonia-vaccinated class
Y 3	Pneumonia-infected class
Y 4	Acute HIV-infected class
Y 5	Chronic HIV-infected class
Y 6	AIDS patients class
Y 7	Acute HIV-pneumonia coepidemic class
Y 8	Chronic HIV-pneumonia coepidemic class
Y 9	AIDS-pneumonia coepidemic class
Y 10	HIV/AIDS treatment class
Y 11	Pneumonia treatment class
Y 12	Coepidemics treatment class

2.3. The HIV/AIDS-Pneumonia Coepidemic Dynamical System

From Figure 1, the HIV/AIDS and pneumonia coinfection dynamical system is given by

Figure 1

Flowchart of the HIV/AIDS-pneumonia coinfection model (4) where λ and λ are given in (2) and (3), respectively.

With initial conditions,

The sum of all the differential equations in (4) is

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

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 Y1(t), Y2(t), Y3(t), Y4(t), Y5(t),  Y6(t), Y7(t), Y8(t), Y9(t), Y10(t), Y11(t), and Y12(t) of system (4) are nonnegative for all time  t > 0.

Proof

Assume Y1(0) > 0, Y2(0) > 0, Y3(0) > 0, Y4(0) > 0, Y5(0) > 0, Y6(0) > 0, Y7(0) > 0, Y8(0) > 0, Y9(0) > 0, Y10(0) > 0, Y11(0), and Y12(0) > 0; then, for all t > 0, we have to prove that Y1 (t) > 0, Y2(t) > 0, Y3(t) > 0, Y4(t) > 0, Y5(t) > 0, Y6(t) > 0, Y7(t) > 0, Y8(t) > 0, Y9(t) > 0,  Y10(t) > 0, Y11(t) > 0, and Y12(t) > 0.

Define: τ = sup{t > 0 : Y1 (t) > 0, Y2(t) > 0, Y3(t) > 0,  Y4(t) > 0, Y5(t) > 0, Y6(t) > 0, Y7(t) > 0, Y8(t) > 0, Y9(t) > 0, Y10(t) > 0, Y11(t) > 0, and Y12(t) > 0}.

From the continuity of Y1(t), Y2(t), Y3(t), Y4(t),  Y5(t), Y6(t), Y7(t), Y8(t), Y9(t), Y10(t), Y11(t), and Y12(t)(t), 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 Y12(0) = 0.

Here, from the first equation of the model (4), we have dY1/dt = (1 − π)Λ + θY10 + τY2 − (d + λ + λ)Y1.

Using the method of integrating factor, we obtained Y1(τ) = M1Y1(0) + M1∫0exp∫(((1 − π)Λ + θY10(t) + τY2(t))dt > 0 where M1 = exp−( > 0, Y1(0) > 0, and from the definition of τ, we see that Y2(t) > 0,  Y10(t) > 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 + ϵλ + λ)Y2 and we have got Y2(τ) = M1Y2(0) + M1∫0exp∫((πΛ)dt > 0, where M1 = exp−( > 0, Y2(0) > 0, and also, the exponential function always is positive; then, the solution Y2(τ) > 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.

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/d)ln(Λ − dN) ≤ t + c where c is some constant, and after some steps of calculations, we have got 0 ≤ N (t) ≤ Λ/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) when Y2 = Y3 = Y7 = Y8 = Y9 = Y10 = Y12 = 0 which is given by

where the total population is N1(t) = Y1(t) + Y4(t) + Y5(t) + Y6(t) + Y11(t) and the HIV/AIDS single infection force of infection is given by λ = (β1/N1)(Y4 + ρ1Y5) with initial conditions Y1(0) > 0, Y4(0) ≥ 0, Y5(0) ≥ 0, Y6(0) ≥ 0, and Y11(0) ≥ 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

With initial conditions, Y1(0) > 0, Y2(0) ≥ 0, Y3(0) ≥ 0, Y10(0) ≥ 0, total population  N2(t) = Y1(t) + Y2(t) + Y3(t) + Y10(t), and pneumonia force of infection  λ = β2Y3(t).

In the region, 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 as Y3 = Y10 = 0 we have got

Y 1 0 = Λ(d + τ1) − Λπd/d(d + τ1) and Y20 = Λπ/(d + τ1) such that E20 = (Y10, Y20, Y30, Y100) = ((Λ(d + τ1) − Λπd/d(d + τ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 number  ℛ2 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 = [∂ℱ(E20)/∂x][∂ν(E20)/∂x]−1 where ℱ is the rate of appearance of new infection in compartment 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

Then, using Mathematica, we have got

The characteristic equation of the matrix FV−1 is

Then, the spectral radius (effective reproduction number ℛ2) of FV−1 of the pneumonia submodel (9) is  ℛ2 = (β2ϵΛπd + β2Λ(d + τ1) − β2Λπd)/(d(d + τ1)(d + κ + d)). Here, ℛ2 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 ℛ2 < 1, otherwise unstable.

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)/(d(d + τ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

Then, the characteristic equation of the above Jacobian matrix is given by where M = ((β2ϵΛπd + β2Λ(d + τ1) − β2Λπd)/(d(d + τ1))) − (d + κ + d).

After some steps, we have got λ1 = −d < 0 or λ2 = −(d + τ1) < 0 or λ3 = (d + κ + d)[ℛ2 − 1] < 0 if ℛ2 < 1 or λ4 = −(d + θ) < 0. Therefore, since all the eigenvalues of the characteristics polynomial of the system (9) are negative if ℛ2 < 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 λ∗ = β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 λ∗ = β2Y3∗ at an endemic equilibrium point. Setting the right-hand sides of the equations of model (9) to zero and we have got Y2∗ = πΛ/(ϵλ∗ + d + τ1), Y10∗ = κY3∗/(d + θ) and substitute Y2∗ and Y10∗  in to Y1∗, we obtain Y1∗ = ((1 − π)Λ + τ1Y2∗ + θT∗)/(d + λ∗) = ((1 − π)Λ/(d + λ∗)d + λ∗) + (πΛτ1/(ϵλ∗ + d + τ1)(d + λ∗)) + (θγY3∗/(d + θ)(d + λ∗)) and substitute Y2∗ and Y1∗ in Y3∗, we obtain

Finally, substitute Y3∗ in to pneumonia submodel (9) force of infection λ∗ = β2Y3∗ as and letting m1 = d + κ + d, m2 = d + τ1, and m3 = d + θ, we have got a2λ∗2+a1λ∗+a0 = 0 where a2 = m1m3ϵ − θκϵ > 0,  a1 = m1m3dϵ + m1m2m3 − m2θκ − β2Λm3ϵ, and a0 = m1m2m3μ[1 − ℛ2] > 0 if ℛ2 < 1.

Here, the nonzero equilibrium(s) of the model (9) satisfies f(λ∗) = a2λ∗2 + a1λ∗ + a0 = 0 so that the quadratic equation can be analyzed for the possibility of multiple equilibriums. It is worth noting that the coefficient a2 is always positive and a0 is positive (negative) if ℛ is less than (greater than) unity, respectively. Hence, we have established the following result.

Theorem 4 .

The pneumonia submodel (9) has the following:

Exactly one unique endemic equilibrium if a0 < 0 (i.e., ℛ2>1)

Exactly one unique endemic equilibrium if a1<0, and a0 = 0 or a12 − 4a2a0 = 0

Exactly two endemic equilibriums if a0 > 0 (i.e., ℛ2 < 1), a1 < 0, and a12 − 4a2a0 > 0

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 if ℛ2 < 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 ℛ2 to 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 ℛ2 = 1 whenever the inequality D1 > D2 holds.

Here, we apply the center manifold theory in [27]; however, to apply this theory, the following simplification and change of variables are made.

Let Y1 = x1,Y2 = x2, Y3 = x3, and Y10 = x4 such that N2 = x1 + x2 + x3 + x4. Furthermore, by using vector notation X = (x1, x2, x3, x4), pneumonia submodel (9) can be written in the form dX/dt = F(X) with

F = (f1, f2, f3, f4), as follows: with λ = β2x3, then the method entails evaluating the Jacobian of system (17) at the DFE point E20, denoted by J(E20), and this gives us

Consider, next, the case when ℛ = 1. Suppose, further, that β2 = β∗  is chosen as a bifurcation parameter.

Solving for β2 from ℛ2 = 1 as ℛ2 = β2ϵΛπd + β2Λ(d + τ1) − β2Λπd/d(d + τ1)(d + κ + d) = 1 and we have got β2 = β∗ = d(d + τ1)(d + κ + d)/ϵΛπd + Λ(d + τ1) − Λπd and

After some steps of the calculation, we have got the eigenvalues of J as λ1 = −d or λ2 = −(d + τ1) or λ3 = 0 or λ4 = −(d + θ).

It follows that the Jacobian J(E20) 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 at ℛ2 = 1

Eigenvectors of J: for the case ℛ2 = 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, u4) with values

Similarly, the left eigenvector associated with the zero eigenvalues at β2 = β∗ given by v = (v1, v2, v3, v4) are v1 = v2 = v4 = 0, v3 = v3 > 0.

After long calculations, the bifurcation coefficients a and b are obtained as a = D1 − D2 where D1=β∗Λπd(d + τ1)(d + θ) + θκd(d + τ1)2/d2(d + τ1)2, and  D2 = (ϵβ∗Λπdτ1(d + θ) + β∗Λ(d + τ1)2(d + θ)/d2(d + τ1)2) + ϵ(ϵβ∗Λπ/(d + τ1)2).

Thus, the bifurcation coefficient a is positive if D1 > D2. Furthermore,  b = v3u2u3(Λ(d + τ1) − Λπd/d(d + τ1)) > 0.

Hence, from in Castillo-Chavez and Song [28], model (9) exhibits a backward bifurcation at ℛ2 = 1 whenever D1 > 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 = 0 and hence E30== (Λ(d + τ1) − Λπd/d(d + τ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 by ℛ0, 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 ℛ0 < 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, ℛ0 > 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 ℛ0 is insufficient [23]. Due to its importance, researchers have sought to find ways of determining ℛ0. Two important concepts in modeling outbreaks of infectious diseases are the basic reproduction number, universally denoted by ℛ0 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 number ℛ3 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 = [∂ℱ(E30)/∂x][∂ν(E30)/∂x]−1 where ℱ is the rate of appearance of new infection in compartment i , ν is the transfer of infections from one compartment i to another, and E30 is the disease-free equilibrium point E30 = (Λ(d + τ1) − Λπd/d(d + τ1), Λπ/d + τ1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).

After long detailed calculations, the transition matrix F is given by and the transmission matrix V is given by

where D1 = d + κ + d, D2 = d + κ1 + δ1, D3 = d + κ2 + δ2, D4 = d + κ3 + d, D5 = d + d + σ1 + δ3, D6 = d + d + σ2 + δ4, and D7 = d + d + ε3.

Then, by using Mathematica, we have got

The characteristic equation of the matrix FV−1 is given by

where A1 = (πϵΛβ2/(d + τ1)) + β2[(−πΛd + Λ[d + τ1]/d[d + τ1])]/D1, B = (β1/D2) + (δ1β1ρ1/D2D3);  then, the eigenvalues of FV−1 are λ1 = β2ϵπdΛ + β2Λ[d + τ1] − β2πΛd/D1(d + τ1) or λ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 matrix FV−1 which is given by

ℛ 3 = max{λ1, λ2} = max{β2ϵπdΛ + β2Λ[d + τ1] − β2πΛd/D1(d + τ1), (β1/D2) + (δ1β1ρ1/D2D3)}. Here,  ℛ2 = β2ϵπdΛ + β2Λ[d + τ1] − β2πΛd/(d + κ + d)(d + τ1) is the effective reproduction number for pneumonia-only infected individual and ℛ1 = (β1/(d + κ1 + δ1)) + (β1ρ1δ1/(d + κ1 + δ1)(d + κ2 + κ2)) is the basic reproduction for HIV/AIDS-only infected individual.

Here, ℛ1 represent the basic reproduction number for HIV/AIDS submodel,  ℛ2 and ℛ3 are the effective reproduction numbers for the pneumonia submodel and HIV/AIDS-pneumonia coinfection model, respectively. The following three outcomes are possible: (i) for ℛ1 < 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 ℛ2 < 1, then E2 is not globally stable in the region Ω2, and pneumonia may spread through the community; (iii) for ℛ3 < 1, the steady state E3 is 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., σ1 or σ2 or σ3) are present in the expression for ℛ3, indicating no impact of treating coinfected population on ℛ3.

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 ℛ3 < 1, otherwise unstable.

The Jacobian matrix J(E30) of model (4) at E30 is given by

where Z1 = β2ϵY20 + β2Y10 − (d + κ + d), Z2 = (β1/N0)Y10 − (d + κ1 + δ1),  Z3 = −(d + κ2 + δ2),  Z4 = −(d + κ3 + d), Z5 = −(d + d + σ1 + δ3), Z6 = −(d + d + σ2 + δ4), Z7 = −(d + d + σ3), and Z8 = −(d + θ).

Then, the characteristic equation of the Jacobian matrix J (E30) is given by

After detailed calculations, we have got that

λ 1 = λ2 = λ3 = −d < 0 or  λ4 = (d + κ + d)[ℛ2 − 1] < 0 if ℛ2 < 1 or λ5 = −(d + τ1) < 0 or λ6 = −(d + κ3 + d) < 0 or λ7 = −(d + d + σ1 + δ3) < 0 or λ8 = −(d + d + σ2 + δ4) < 0 or λ9 = −(d + d + σ3) < 0 or λ10 = −(d + θ) < 0 or a2λ2 + a1λ + a0 = 0 for a2 = 1, a1 = (d + κ2 + δ2) + (d + κ1 + δ1)[1 − (Y10/N0)ℛ] > 0 if ℛ < 1 and a0 = (d + κ2 + δ2)(d + κ1 + δ1)[1 − (Y10/N0)ℛ) > 0 if ℛ < 1.

Then, by applying Routh-Hurwitz stability criteria since a2 = 1 > 0, a1 > 0, and a0 > 0, all the eigenvalues of the Jacobian matrix are negative if ℛ1 < 1 and ℛ2 < 1, i.e., ℛ3 = max{ℛ1, ℛ2} < 1. Thus, the disease-free equilibrium point (DFE) of HIV/AIDS-pneumonia coinfection model (4) is locally asymptotically stable if ☐

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 ℛ1 < 1 and there is/are an endemic equilibrium point(s) if ℛ2 < 1 implies that there is/are endemic equilibrium point(s) if ℛ3 < 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

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:

E4∗ = (Y1∗, 0, Y4∗, Y5∗, Y6∗, 0, 0, 0, 0, 0, Y11∗, 0), if ℛ1 > 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)

E5∗ = (Y1∗, Y2∗, Y3∗, 0, 0, 0, 0, 0, 0, Y10∗, 0, 0), if ℛ2 > 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)

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 quantity  ℛ3 = max{ℛ1, ℛ2}  is the effective reproduction number of the system (4) where ℛ1 and ℛ2 are defined as above.

Theorem 7 .

Model (4) exhibits the phenomenon of backward bifurcation at ℛ3 = 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.

Let Y1 = 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 notation  X = (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12), complete model (4) can be written in the form dX/dt = F(X) with F = (f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12), as follows: with λ = β1/N[x4 + ρ1x5 + ρ2x7 + ρ2x8] where ρ3 ≥ ρ2 ≥ ρ1 ≥ 1 and λ = β2[x3 + ω1x7 + ω2x8 + ω3x9]  where ω3 ≥ ω2 ≥ ω1 ≥ 1; then, the method entails evaluating the Jacobian of system (29) at the DFE E30, denoted by J(E30), and this gives us where F1 = −β2Y10, F2 = −β1(Y10/(Y10 + Y20)), F3 = −β1ρ1(Y10/(Y10 + Y20)), F4 = −β1ρ2(Y10/(Y10 + Y20)) − β2ω1Y10, F5 = −β1ρ3(Y10/(Y10 + Y20)) − β2ω2Y10, F6 = −β2ω3Y10,  F7 = −(d + τ1), F8 = −ϵβ2Y20, F9 = −β1(Y20/(Y10 + Y20)), F10 = −β1ρ1(Y20/(Y10 + Y20)),

F 11 = −ϵβ2ω1V0 − β1ρ2(Y0/(Y10 + Y0)), F12 = −ϵβ2ω2V0 − β1ρ3(Y20/(Y10 + Y20)), F13 = −ϵβ2ω3Y20, F14 = ϵβ2Y20 + β2Y10 − (d + κ + d), F15 = ϵβ2ω1Y20 + β2ω1Y10, F16 = ϵβ2ω2Y20 + β2ω2Y10, F17 = ϵβ2ω3Y20 + β2ω3Y10, F18 = β1(Y10/(Y10 + Y20)) − (d + κ1 + δ1), F19 = β1ρ1(Y10/(Y10 + Y20)), F20 = β1ρ2(Y10/(Y10 + Y20)), F21 = β1ρ3(Y10/(Y10 + Y20)), F22 = −(d + κ2 + δ2), F23 = −(d + κ3 + d), F24 = −(d + d + σ1 + δ3), F25 = −(d + d + σ2 + δ4), F26 = −(d + d + σ3), and  F27 = −(d + θ).

Without loss of generality, consider the case when ℛ2 > ℛ1, and ℛ3 = 1, so that ℛ2 = 1. Furthermore, let β2 = β∗ is chosen as a bifurcation parameter. Solving for β2 from ℛ2 = 1 as ℛ2 = β2ϵΛπd + β2Λ(d + τ1) − β2Λπd/d(d + τ1)(d + κ + d) = 1, we have got the value β∗ = β2 = d(d + τ1)(d + κ + d)/ϵΛπd + Λ(d + τ1) − Λπd.

After solving the Jacobian J(E30) 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 + d + κ3) < 0 or λ8 = −(d + κ3 + d) < 0 or λ9 = −(d + d + σ1 + δ3) < 0 or λ10 = −(d + d + σ1 + δ4) < 0 or where a2 = 1 > 0, a1 = (d + κ1 + δ1)[1 − (Y10/(Y10 + Y40))ℛ] + (d + κ2 + δ2) > 0 if ℛ < 1, and a0 = a0 = (d + κ1 + δ1)(d + κ2 + δ2)[1 − (Y10/(Y10 + Y20))ℛ1] > 0 if ℛ1 < 1.

Equation (31) has/have no positive root/s whenever ℛ1 < 1, and hence, both eigenvalues are negative. It follows that the Jacobian J(E30) 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 ℛ = 1.

Eigenvectors of J: for the case when ℛ = 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, u12) are u1 = ((−τ1F8F27 + F1F7F27 − θγF7)/dF7F27)u3, u2 = −(F8/F7)u3, u3 = u3 > 0, u10 = (−κ/F27)u3 , and u4 = u5 = u6 = u7 = u8 = u9 = u11 = u12 = 0.

The left eigenvectors associated with the zero eigenvalue at β2 = β2∗ satisfying v.w = 1 given by v = (v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12) are v1 = v2 = v4 = v5 = v6 = v10 = v11 = v12=0, v3 = v3 > 0, v7 = ((−δ3δ4F17 + δ3F16F26 − F15F25F26)/F24F25F26)v3, v8 = (δ4F17 − F16F26)/F25F26v3, and v9 = −(F17/F26)v3.

After going through detailed computations and simplification, we have the following bifurcation coefficients a and b as

⟹a = 2β2∗v3u32[G1 − G2] where G1=θκ/d(d + θ) and G2 = ϵβ∗τ1V0 + β∗(d + τ1)Y10 + ϵβ∗dY20/d(d + τ1). Thus, the bifurcation coefficient a is positive whenever G1 > G2. Furthermore, b = v3u3(∂2f2(0, 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 ℛ3 = ℛ2 = 1 whenever G1 > G2.

Theorem 8 .

Model (4) will undergo backward bifurcation if a = G1 > G2 > 0

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 variable ℛ3 that depends differentiably on a parameter p is defined as SI(p) = (∂ℛ3/∂p)∗(p/ℛ3) [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 of ℛ1 and ℛ2 since ℛ3 = max{ℛ1, ℛ2}. 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.

Parameter	Value	Source
Λ	0.0413∗N0	Estimated
d	0.02	Estimated
δ 1	0.498	[7]
δ 2	0.08	[7]
δ 3	0.2885	[6]
δ 4	0.3105	[6]
ψ 1	1.1	Assumed
ψ 2	1.2	Assumed
ψ 3	1.4	Assumed
ν	1	Assumed
d P	0.1	[16]
θ	0.1	[18]
d A	0.333	[6]
π	0.2	[18]
τ 1	0.0025	[18]
ϵ	0.002	[18]
d AP	0.42	Assumed
κ	0.2	[18]
κ 1	0.2	[7]
κ 2	0.15	[7]
κ 3	0.13	Assumed
σ 1	0.498	[7]
σ 2	0.08	[7]
σ 3	0.230	Assumed
β 1	Variable	[6]
β 2	Variable	[6]
ρ 1, ρ2, ρ3, ω1, ω2, ω3	1.2,1,1,1,1,1	Assumed

Using the values of parameters in Table 3, the sensitivity indices are calculated in Tables 4 and 5.

Table 4

Sensitivity indices of ℛ3 = ℛ1.

Sensitivity index	Values
SI(β1)	+1
SI(ρ1)	+0.6134
SI(δ1)	- 0.0639
SI(d)	-0.3150
SI(κ1)	-0.1371
SI(κ2)	-0.0264
SI(δ2)	-0.0141

Table 5

Sensitivity indices of ℛ3 = ℛ2.

Sensitivity index	Values
SI(Λ)	+1
SI(β2)	+1
SI(d)	-0.4421
SI(κ)	-0.6559
SI(dP)	-0.3852
SI(π)	-0.3852
SI(ϵ)	-0.3852
SI(τ1)	-0.3852

In this paper, with parameter values in Table 3, we have got ℛ1 = 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 (ℛ1).

Similarly, with parameter values in Table 3, we have got ℛ2 = 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 ℛ2, 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 ℛ1 and ℛ2 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)

Figure 2 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

Figure 2

Local stability of endemic equilibrium point of the coepidemic model (4) whenever ℛ1 = 1.386 at β1 = 2 and ℛ2 = 9.69 at β2 = 0.2.

4.2.2. Effect of Parameters on the Threshold Parameter ℛ2

In this subsection, as we see in Figure 3, we have investigated the effect of pneumonia vaccination portion π on the pneumonia effective reproduction number  ℛ2. 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 ℛ2 < 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 on ℛ2.

In this subsection, as we see in Figure 4, we have investigated the effect of pneumonia spreading rate β2 on the pneumonia effective reproduction number ℛ2 by keeping the other rates as in Table 3. Figure 4 reflects that when the value of β2 increases, the pneumonia effective reproduction number ℛ2 increases, and whenever the value of β2 < 0.022 implies ℛ2 < 1. Therefore, public policymakers must concentrate on minimizing the values of pneumonia spreading rate β2 to minimize pneumonia effective reproduction number ℛ2.

Figure 4

Effect of pneumonia transmission on ℛ2.

4.2.3. Effect of Pneumonia Treatment Rate on Infectious Population

In this subsection, as we see in Figure 5, 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 Figures 6–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 Figure 9 the effect of the spreading rate of HIV/AIDS β1 on the acute HIV-pneumonia coepidemic population Y7. 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 Figures 10–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 Section 1, 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 is ℛ3 = max{ℛ1, ℛ2} = max{1.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.