Optimal control strategies on HIV/AIDS and pneumonia co-infection with mathematical modelling approach

Abstract

In this paper, a compartmental model on the co-infection of pneumonia and HIV/AIDS with optimal control strategies was formulated using the system of ordinary differential equations. Using qualitative methods, we have analysed the mono-infection and HIV/AIDS and pneumonia co-infection models. We have computed effective reproduction numbers by applying the next-generation matrix method, applying Castillo Chavez criteria the models disease-free equilibrium points global stabilities were shown, while we have used the Centre manifold criteria to determine that the pneumonia infection and pneumonia and HIV/AIDS co-infection exhibit the phenomenon of backward bifurcation whenever the corresponding effective reproduction number is less than unity. We carried out the numerical simulations to investigate the behaviour of the co-infection model solutions. Furthermore, we have investigated various optimal control strategies to predict the best control strategy to minimize and possibly to eradicate the HIV/AIDS and pneumonia co-infection from the community.

KEYWORDS:

• HIV/AIDS
• pneumonia
• co-infection
• vaccination
• protection

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 34H05
• 92-10

1. Introduction

Infectious diseases caused by tiny microorganisms such as viruses, bacteria, fungi and parasites are clinically verified illnesses and have been reported as the leading cause of death in human beings throughout the world [1,2].

AIDS caused by human immunodeficiency virus (HIV) has been among the leading life-affecting disease which is expanding throughout nations in the world since 1981 [3–5]. It remains a significant world health issue that impacts almost 70 million people worldwide and has been a significant cause of morbidity and mortality [6]. HIV is transmitted through sexual intercourse, needle sharing, and direct contact with blood or other body fluids containing the virus and from mother to child during childbirth [3,7].

Pneumonia caused by various pathogenic microbial agents such as virus, bacteria, fungi and parasites is a major respiratory infectious disease identified as an inflammatory condition of the lungs [8–10]. Among the pathogenic microbial agents which have potential in causing pneumonia infection, bacteria especially Streptococcus pneumoniae has been reported as the leading cause [8–11]. The microbial agents enter the lungs, rapidly multiply its number, settle in the air passage called alveoli of the human being lung, the lung will be filled with fluid and pus and makes breathing difficult [8,10]. Pneumonia is commonly a highly transmitted disease and major causes of morbidity and mortality in both children and adults throughout nations in the world [1,11,12].

Many scholars in different nations have proposed and analysed mathematical model of the phenomenon of the transmission and control of infectious diseases. Mathematical models of infectious diseases are used for examining the spreading dynamics and also for making predictions of quantitative measures of different controlling strategies and their effectiveness as we reviewed in the literature [1–11,13–40]. Even though, mathematical epidemiologists did not give attention like the common HIV/AIDS and TB co-infection and other co-infections, pneumonia is the most opportunistic infection for HIV/AIDS infected individuals; the co-infection of HIV/AIDS and pneumonia in one host is a common phenomenon. From few essential mathematical epidemiological research studies that have been studied by scholars on the transmission dynamics of HIV/AIDS and pneumonia co-infection, the existence of pneumonia and HIV co-infection has been delineating in the literature [20,21].

In this study, we have examined research papers that have been done on the transmission dynamics of different infectious diseases. Huo and Chen [16] formulated a stage structure HIV/AIDS deterministic model to study the transmission of HIV/AIDS with treatment. The model analysis showed that ART at the asymptomatic initial stages of the HIV incidence or before-AIDS stage is the most effective to decrease its spreading rate. Omondi et al. [17] analysed a sex-structured community infection model and discussed male and female HIV infection trends with heterosexual activities. The analysis of the study deduced that ART treatment has a considerable impact on the HIV/AIDS spreading rate. Saha et al. [32] formulated and analysed a mathematical mode of HIV/AIDS with optimal control theory to prevent the transmission through pre-exposure prevention (PrEP) and limited treatment strategies. From the finding of the study, one can observe that the combined strategies effect is more powerful and effective than single strategies. Samanta et al. [33] formulated and analysed a non-autonomous HIV/AIDS epidemic model with distributed time delay. Teklu [27] constructed a mathematical model on COVID-19 in the presence of prevention and control strategies. The result emphasizes that the spreading of COVID-19 in the considered community will be minimized through applications of the intervention strategies presented on the study. Ahmed et al. [22] constructed a mathematical model of HIV and COVID-19 co-infection to globally assess the COVID-19 pandemic situation in different countries affected by both diseases, such as South Africa, Brazil and many other countries. Teklu and Mekonnen [20] proposed and analysed HIV/AIDS and pneumonia co-infection model with treatment at each infection stage and their model analysis verified that applying treatment mechanisms for both the single infections and co-infection models is the most effective strategy to minimize the co-infection transmission dynamics. Teklu and Rao [21] constructed and examined a deterministic mathematical model on the co-existence of HIV/AIDS and pneumonia with control measures such as pneumonia vaccination and treatments of pneumonia and HIV/AIDS infections. Even though researchers stated in the literature [20] and [21] invested much effort in studying HIV/AIDS and pneumonia co-infection, they did not considered pneumonia protection, pneumonia vaccination, pneumonia treatment, HIV protection by using condom, and HIV treatment as prevention and control strategies simultaneously in their proposed co-infection model formulation and analysis. Moreover, based on the findings and limitations in the studies stated in [20] and [21], we are motivated and inspired to undertake this study and to fulfil the gap.

The remaining part of this study has structured in the sequence: the model is formulated in Section 2 and is analysed in Section 3, optimal control analysis and numerical simulation, and discussion and conclusion of the study are carried out in Sections 4, 5 and 6 respectively.

2. Model descriptions and formulation

In this study, we partitioned the total human population at a given time t denoted by $N\left(t\right),$ into 11 distinct classes depending on their infection status: susceptible individuals to both pneumonia and HIV $D_1\left(t\right)$, pneumonia protected individuals (by good hygiene, quitting smoking, and getting regular physical activity and eating healthy etc.) $D_2\left(t\right)$, HIV protected individuals by using condom $D_3\left(t\right)$, pneumonia vaccinated individuals $D_4\left(t\right)$, pneumonia mono-infected individuals $D_5\left(t\right)$, HIV unaware mono-infected individuals $D_6\left(t\right),$ HIV aware mono-infected individuals $D_7\left(t\right)$, HIV unaware and pneumonia co-infected individuals $D_8\left(t\right)$, HIV aware and pneumonia co-infected individuals $D_9\left(t\right)$, individuals recovered from pneumonia $D_{10}\left(t\right)$ and HIV aware treated individuals $D_{11}\left(t\right)$ so that $N\left(t\right)=D_1\left(t\right)+D_2\left(t\right)+D_3\left(t\right)+D_4\left(t\right)+D_5\left(t\right)+D_6\left(t\right)+D_7\left(t\right)+D_8\left(t\right)+D_9\left(t\right)+D\left(t\right)+D_{11}\left(t\right)$. Since HIV/AIDS is a chronic infectious disease the susceptible individuals acquires HIV infection at the standard incidence rate given by (1) ${\lambda }_H\left(t\right)=\frac{{\text{β}}_1}{N}\left(D_6\left(t\right)+{\rho }_1{D}_7\left(t\right)+{\rho }_2{D}_8\left(t\right)+{\rho }_3{D}_9\left(t\right)\right).$(1) where ${\rho }_3\ge {\rho }_2\ge {\rho }_1\ge 1$ are the modification parameters that increase infectivity and ${\text{β}}_1$ is the HIV transmission rate.

Since pneumonia is an acute infectious disease, the susceptible individuals acquire pneumonia infection at the mass action incidence rate given by (2) ${\lambda }_{\text{P}}\left(\text{t}\right)={\text{β}}_2\left(D_5\left(\text{t}\right)+{\omega }_1{D}_8\left(\text{t}\right)+{\omega }_2{D}_9\left(\text{t}\right)\right).$(2) where ${\omega }_2\ge {\omega }_1\ge 1$ are the modification parameters that increase infectivity and ${\text{β}}_2$ is the pneumonia transmission rate.

Assume: $k_1$, $k_2$, $k_3,$ and $k_4$ where $k_1+k_2+k_3+k_4=1$ are portions of the recruited individuals who are entering to susceptible class, pneumonia protected class, HIV protected class and pneumonia vaccinated class respectively. The susceptible class is increased by individuals from pneumonia vaccinated class in which those individuals who are vaccinated against pneumonia but did not respond to vaccination with waning rate of $\rho$ and from pneumonia recovery class who lose their temporary immunity by the rate $\eta$, pneumonia vaccine may not be 100% efficient, so vaccinated individuals also have a chance of being infected with portion $ϵ$ of pneumonia serotype not covered by pneumonia vaccine where $0\le ϵ<1$, the human population distribution is homogeneous in each class, individuals who exactly apply protection measures such as good hygiene, quitting smoking, and getting regular physical activity and eating healthy and so on have fully protected from pneumonia infection, HIV-treated individuals do not transmit infection to others, individuals in each class are subject to natural death rate $\mu$, population of human being is variable, there is no dual-infection transmission simultaneously and no vertical HIV transmission.

In this section using parameters described in Table 1, model variables definitions given in Table 2, and the given assumptions the schematic diagram of the co-infection is given by Figure 1.

Figure 1. Schematic diagram of the pneumonia and HIV/AIDS co-infection transmission dynamics where ${\lambda }_H\left(t\right)$ and ${\lambda }_P\left(t\right)$ are given in (1) and (2) respectively.

Table 1. Model parameter descriptions.

Parameters	Descriptions
$\mu$	Human natural mortality rate
$\mathrm{\Delta }$	Number of recruited individuals
${\alpha }_1$	Pneumonia protection lose rate
${\alpha }_2$	HIV protection lose rate
$ϵ$	The proportion not covered by the pneumonia vaccine
$\theta$	Progression rate
${\varphi }_1$	Modification parameter
${\varphi }_2$	Modification parameter
$d_1$	Pneumonia death rate
$d_2$	HIV/AIDS death rate for unaware
$d_3$	HIV/AIDS death rate for aware
$\kappa$	Pneumonia recovery rate
$\gamma$	HIV aware infection treatment rate
$\rho$	Pneumonia vaccination waning rate
$\upsilon$	Modification parameter
${\text{β}}_1$	HIV/AIDS transmission rate
${\text{β}}_2$	Pneumonia transmission rate
${\text{k}}_1$	Portion of recruitment entered to susceptible
${\text{k}}_2$	Pneumonia protection portion
${\text{k}}_3$	HIV protection portion
${\text{k}}_4$	Pneumonia vaccination portion
$\delta$	Co-infection progression rate
${\theta }_1$	Pneumonia treatment rate for HIV unaware co-infected
${\theta }_2$	Pneumonia treatment rate for HIV aware co-infected
$\eta$	The temporary immunity development rate
$d_4$	HIV unaware and pneumonia co-infected death rate
$d_5$	HIV aware and pneumonia co-infected death rate

Table 2. Biological definitions of model variables.

Variables	Biological definitions
$D_1$	Susceptible class to both HIV and pneumonia
$D_2$	Individuals who are protected against pneumonia
$D_3$	Individuals who are protected against HIV
$D_4$	Pneumonia vaccinated class
$D_5$	Pneumonia infected class
$D_6$	Individuals mono-infected with HIV and unaware
$D_7$	Individuals mono-infected with HIV and aware
$D_8$	Co-infected individuals unaware of HIV infection
$D_9$	Co-infected individuals aware of HIV infection
$D_{10}$	Pneumonia recovery class
$D_{11}$	HIV/AIDS treatment class

Using Figure 1, the system of differential equations of the HIV/AIDS and pneumonia co-infection is given by (3) $\begin{array}{rl}\dot{D_1}& =k_1\mathrm{\Delta }+{\alpha }_1{D}_2+{\alpha }_2{D}_3+\rho D_4+\eta D_{10}-\left({\lambda }_H+{\lambda }_P+\mu \right)D_1,\\ \dot{D_2}& =k_2\mathrm{\Delta }-\left({\lambda }_H+{\alpha }_1+\mu \right)D_2,\\ \dot{D_3}& =k_3\mathrm{\Delta }-\left({\alpha }_2+\mu +{\lambda }_P\right)D_3,\\ \dot{D_4}& =k_4\mathrm{\Delta }-\left(\rho +\mu +{\lambda }_H+ϵ{\lambda }_P\right)D_4,\\ \dot{D_5}& ={\lambda }_P{D}_1+{\lambda }_P{D}_2+ϵ{\lambda }_P{D}_4-\left(\mu +d_1+\kappa +\upsilon {\lambda }_H\right)D_5,\\ \dot{D_6}& ={\lambda }_H{D}_1+{\lambda }_H{D}_2+{\lambda }_H{D}_4+{\theta }_1{D}_8-\left(\theta +\mu +d_2+{\varphi }_1{\lambda }_P\right)D_6,\\ \dot{D_7}& =\theta D_6+{\theta }_2{D}_9-\left(\gamma +d_3+\mu +{\varphi }_2{\lambda }_P\right)D_7,\\ \dot{D_8}& ={\varphi }_1{\lambda }_P{D}_6+\upsilon {\lambda }_H{D}_5-\left(\mu +d_4+\delta +{\theta }_1\right)D_8,\\ \dot{D_9}& =\delta D_8+{\varphi }_2{\lambda }_P{D}_7-\left(\mu +d_5+{\theta }_2\right)D_9,\\ \dot{D_{10}}& =\kappa D_5-\left(\mu +\eta \right)D_{10},\\ \dot{D_{11}}& =\gamma D_7-\mu D_{11},\end{array}$(3) with the corresponding initial conditions (4) $\begin{array}{rl}& D_1\left(0\right)>0,D_2\left(0\right)\ge 0,D_3\left(0\right)\ge 0,\\ & D_4\left(0\right)\ge 0,D_5\left(0\right)\ge 0,D_6\left(0\right)\ge 0,D_7\left(0\right)\ge 0,D_8\left(0\right)\ge 0,D_9\left(0\right)\ge 0,\\ & D_{10}>0,\text{and}{D}_{11}>0.\end{array}$(4) The sum of all the differential equations in (3) is (5) $\dot{N}=\mathrm{\Delta }-\mu N-(d_1{D}_5+d_2{D}_6+d_3{D}_7+d_4{D}_8+d_5{D}_9).$(5)

2.1. Two qualitative properties of the co-infection model

The model (3) is both biologically and mathematically meaningful if and only if all the model solutions are non-negative and bounded in the invariant region (6) $\mathrm{\Omega }=\left\{\begin{array}{c}\left(D_1,D_2,D_3,D_4,D_5,D_6,D_7,D_8,D_9,D_{10},D_{11}\right)\in {\mathbb{R}}_{+}^{11},N\le \frac{\mathbf{\Lambda }}{\mathit{d}}\end{array}\right\}$(6)

Theorem 2.1 (Non-negativity of the model solutions)

Let us given the initial data in Equation (4), then the solutions $D_1\left(t\right)$, $D_2\left(t\right),D_3\left(t\right),D_4\left(t\right),D_5\left(t\right),D_6\left(t\right)$, $D_7\left(t\right),D_8\left(t\right)$, $D_9\left(t\right),$ $D_{10}\left(t\right),$ and $D_{11}\left(t\right)$ of the COVID-19 and HIV/AIDS co-infection model (3) are nonnegative for all time $t>0$.

Proof:

Appendix A

Theorem 2.2 (The Invariant region)

All the feasible non-negative solutions of the co-infection model (3) are bounded in the region (6).

Proof:

Appendix B

Note: Since the model (3) solutions are both non-negative and bounded in the region (6) the pneumonia and HIV/AIDS co-infection model (3) is both mathematically and biologically meaningful then we can consider the two mono-infection models, namely, HIV mono-infection and pneumonia mono-infection models. This is fundamental for the analysis of the full pneumonia and HIV/AIDS co-infection model.

3. Qualitative analysis of the models

Before analysing the pneumonia and HIV/AIDS co-infection model given in Equation (3), it is very crucial to gain some basic backgrounds about the pneumonia and HIV/AIDS mono-infection models.

3.1. Analysis of HIV/AIDS mono-infection model

In this section, we assume there is no pneumonia infection in the community, i.e. $D_2=$ $D_4=$ $D_5=D_8=D_9=D_{10}=0$ in Equation (3), then the HIV/AIDS mono-infection model is given by (7) $\begin{array}{rl}\dot{D_1}& =k_1\mathrm{\Delta }+{\alpha }_2{D}_3-\left({\lambda }_H+\mu \right)D_1,\\ \dot{D_3}& =k_3\mathrm{\Delta }-\left({\alpha }_2+\mu \right)D_3,\\ \dot{D_6}& ={\lambda }_H{D}_1-\left(\theta +\mu +d_2\right)D_6,\\ \dot{D_7}& =\theta D_6-\left(\gamma +d_3+\mu \right)D_7,\\ \dot{D_{11}}& =\gamma D_7-\mu D_{11},\end{array}$(7) where the total population is $N_1\left(t\right)=D_1\left(t\right)+D_3\left(t\right)+D_6\left(t\right)+D_7\left(t\right)+D_{11}\left(t\right),$ the HIV mono-infection model force of infection is given by ${\lambda }_H=\frac{{\text{β}}_1}{N_1}\left(D_6+{\rho }_1{D}_7\right)$, and initial conditions $D_1\left(0\right)>0$, $D_3\left(0\right)\ge 0$, $D_6\left(0\right)\ge 0$, $D_7\left(0\right)\ge 0$ and $D_{11}\left(0\right)\ge 0$. In a similar manner of the full co-infection model (3) in the region ${\mathrm{\Omega }}_1=\left\{\left(D_1,D_3,D_6,D_7,D_{11}\right)\in {{\mathbb{R}}^5}_{+},N_1\le \frac{\mathrm{\Delta }}{\mu }\right\}$, it is sufficient to consider the dynamics of the model (7) in ${\mathrm{\Omega }}_1$ as biologically and mathematically well-posed.

3.1.1. Local stability of the model disease-free equilibrium point (DFE)

The disease-free equilibrium point of the HIV mono-infection system (7) is obtained by making its right-hand side as zero and setting the infected classes and treatment class to zero as $D_6=D_7=D_{11}=0$ we do have $D_1^0=\frac{k_1\mathrm{\Delta }\left({\alpha }_2+\mu \right)+{\alpha }_2{k}_3\mathrm{\Delta }}{\mu \left({\alpha }_2+\mu \right)}$, $D_3^0=\frac{k_3\mathrm{\Delta }}{{\alpha }_2+\mu }$. Hence the disease-free equilibrium point is given by $E_{HM}^0=\left(D_1^0,D_3^0,0,0,0\right)=\left(\frac{k_1\mathrm{\Delta }\left({\alpha }_2+\mu \right)+{\alpha }_2{k}_3\mathrm{\Delta }}{\mu \left({\alpha }_2+\mu \right)},\frac{k_3\mathrm{\Delta }}{{\alpha }_2+\mu },0,0,0\right)$.

The local stability of the HIV mono-infection model (7) DFE is examined by its basic reproduction number denoted by ${\mathcal{R}}_{HM}$, which is calculated by using the next-generation operator method stated by Van den Driesch and Warmouth [24]. Applying the method given in [24], the transmission matrix $F$ and the transition matrix $V$, i.e. for the new infection and the remaining transfer respectively, are given by $F=\left[\begin{array}{ccc}\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}& \begin{array}{c}\frac{{\rho }_1{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}\end{array}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],$and $V=\left[\begin{array}{ccc}\theta +\mu +d_2& 0& 0\\ -\theta & \gamma +d_3+\mu & 0\\ 0& -\gamma & \mu \end{array}\right].$Then we have calculated that $V^{-1}=\left[\begin{array}{ccc}\frac{1}{\left(\theta +\mu +d_2\right)}& 0& 0\\ \frac{\theta }{\left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}& \frac{1}{\left(\gamma +\mu +d_3\right)}& 0\\ \frac{\gamma \theta }{\mu \left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}& \frac{\gamma }{\mu \left(\gamma +\mu +d_3\right)}& \frac{1}{\mu }\end{array}\right],$and $F{V}^{-1}=\left[\begin{array}{ccc}\frac{{\text{β}}_1{D}_1^0}{(D_1^0+D_3^0)\left(\theta +\mu +d_2\right)}+\frac{{\text{β}}_1{\rho }_1\theta D_1^0}{(D_1^0+D_3^0)\left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}& \begin{array}{c}\frac{{\text{β}}_1{\rho }_1{D}_1^0}{(D_1^0+D_3^0)\left(\gamma +\mu +d_3\right)}\end{array}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right].$Then the effective reproduction number of the HIV mono-infection model (7) is defined as the largest eigenvalue in magnitude of the next generation matrix, $F{V}^{-1}$ given by ${\mathcal{R}}_{HM}=\frac{{\text{β}}_1\left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\alpha }_2{k}_3}{\left({\alpha }_2+\mu \right)\left(\theta +\mu +d_2\right)}+\frac{{\text{β}}_1{\rho }_1\theta \left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\rho }_1\theta {\alpha }_2{k}_3}{\left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}.$Here ${\mathcal{R}}_{HM}$ is defined as the total average number of secondary HIV unaware and HIV aware infection cases acquired from a typical HIV unaware or HIV aware individual during his/her effective infectious period in a susceptible population. The threshold result ${\mathcal{R}}_{HM}$ is the effective reproduction number for HIV mono-infection.

Theorem 3.1:

The DFE point of the HIV mono-infection model given in Equation (7) is locally asymptotically stable (LAS) if ${\mathcal{R}}_{HM}<1$, and it is unstable if ${\mathcal{R}}_{HM}>1$.

Proof:

The local stability of the disease-free equilibrium point of HIV mono-infection model (7) is studied by applying the Routh–Hurwitz stability criteria.

The Jacobian matrix of the HIV mono-infection model given in Equation (7) at the disease-free equilibrium point $E_{HM}^0$ is given by $J\left(E_{HM}^0\right)=\left[\begin{array}{ccccc}-\mu & {\alpha }_2& -\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}& -\frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}& 0\\ 0& -\left({\alpha }_2+\mu \right)& 0& 0& 0\\ 0& 0& \frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}-\left(\theta +\mu +d_2\right)& \frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}& 0\\ 0& 0& \theta & -\left(\gamma +d_3+\mu \right)& 0\\ 0& 0& 0& \gamma & -\mu \end{array}\right],$Then the corresponding characteristic equation of the Jacobian matrix $J\left(E_{HM}^0\right)$ is given by $\left|\begin{array}{ccccc}-\mu -\lambda & {\alpha }_2& -\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}& -\frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}& 0\\ 0& -\left({\alpha }_2+\mu \right)-\lambda & 0& 0& 0\\ 0& 0& \frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}-\left(\theta +\mu +d_2\right)-\lambda & \frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}& 0\\ 0& 0& \theta & -\left(\gamma +d_3+\mu \right)-\lambda & 0\\ 0& 0& 0& \gamma & -\mu -\lambda \end{array}\right|=0.$ $\begin{array}{rl}& \Rightarrow \left(-\mu -\lambda \right)\left(-\left({\alpha }_2+\mu \right)-\lambda \right)\left(-\mu -\lambda \right)\\ & \times \left[\left(\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}-\left(\theta +\mu +d_2\right)-\lambda \right)\right(-\left(\gamma +d_3+\mu \right)-\lambda )\\ & -\frac{{\text{β}}_1{\rho }_1{D}_1^0\theta }{D_1^0+D_3^0}]=0.\end{array}$Then we have determined $\left(-\mu -\lambda \right)\left(-\left({\alpha }_2+\mu \right)-\lambda \right)\left(-\mu -\lambda \right)\left({\lambda }^2+a\lambda +b\right)=0$, where $a=\left(\gamma +d_3+\mu \right)+\left(\theta +\mu +d_2\right)-\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0},$and $b=-\left(\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}-\left(\theta +\mu +d_2\right)\right)\left(\gamma +d_3+\mu \right)-\frac{{\text{β}}_1{\rho }_1{D}_1^0\theta }{D_1^0+D_3^0}.$Then we have got ${\lambda }_1=-\mu <0$ or ${\lambda }_2=-\left({\alpha }_2+\mu \right)<0$ or ${\lambda }_3=-\mu <0$ or (8) ${\lambda }^2+a\lambda +b=0.$(8) For Equation (8), we can apply Routh–Hurwitz stability criteria stated in [29] to determine that both eigenvalues are negative if ${\mathcal{R}}_{HM}<1$ and hence we can conclude that since all the eigenvalues are negative whenever ${\mathcal{R}}_{HM}<1$ the disease-free equilibrium point of the model (7) is locally asymptotically stable whenever ${\mathcal{R}}_{HM}<1.$ From Theorem 3.1, biologically one can conclude that HIV infection dies out from the population (whenever ${\mathcal{R}}_{HM}<1$) if the initial size of the sub-populations of the HIV mono-infection model given in Equation (7) is in the basin of attraction of the disease-free equilibrium point $E_{HM}^0$.

3.1.2. Existence of endemic equilibrium point(s)

Let $E_{HM}^{\ast }=\left(D_1^{\ast },D_3^{\ast },D_6^{\ast },D_7^{\ast },D_{11}^{\ast }\right)$ be an arbitrary endemic equilibrium point of the HIV mono-infection model (7) which can be determined by making the right-hand side of Equation (7) as zero. Then after a number of steps of computations, we have determined that (9) $\begin{array}{rl}{D}_1^{\ast }& =\frac{k_1\mathrm{\Delta }{m}_1+{\alpha }_2{k}_3\mathrm{\Delta }}{m_1\left({\lambda }_H^{\ast }+\mu \right)},D_3^{\ast }=\frac{k_3\mathrm{\Delta }}{m_1},D_6^{\ast }=\frac{k_1\mathrm{\Delta }{m}_1{\lambda }_H^{\ast }+{\alpha }_2{k}_3\mathrm{\Delta }{\lambda }_H^{\ast }}{m_1{m}_2\left({\lambda }_H^{\ast }+\mu \right)},\\ D_7^{\ast }& =\frac{k_1\mathrm{\Delta }\theta m_1{\lambda }_H^{\ast }+{\alpha }_2{k}_3\mathrm{\Delta }\mathrm{\Theta }{\lambda }_H^{\ast }}{m_1{m}_2{m}_3\left({\lambda }_H^{\ast }+\mu \right)},D_{11}^{\ast }=\frac{k_1\mathrm{\Delta }\theta \gamma m_1{\lambda }_H^{\ast }+{\alpha }_2{k}_3\mathrm{\Delta }\theta \gamma {\lambda }_H^{\ast }}{\mu m_1{m}_2{m}_3\left({\lambda }_H^{\ast }+\mu \right)},\end{array}$(9) where $m_1=\left({\alpha }_2+\mu \right)$, $m_2=\left(\theta +\mu +d_2\right)$, and $m_3=\left(\gamma +d_3+\mu \right)$.

Now substitute $D_6^{\ast }$ and $D_7^{\ast }$ which are given in Equation (9) into the HIV/AIDS force of infection ${\lambda }_H^{\ast }=\frac{{\text{β}}_1{D}_6^{\ast }+{\text{β}}_1{\rho }_1{D}_7^{\ast }}{D_1^{\ast }+D_3^{\ast }+D_6^{\ast }+D_7^{\ast }+D_{11}^{\ast }}$, then we have determined the result ${\lambda }_H^{\ast }=\frac{m_4{\lambda }_H^{\ast }}{m_5+m_6{\lambda }_H^{\ast }}$, (10) $\left(m_5+m_6{\lambda }_H^{\ast }-m_4\right){\lambda }_H^{\ast }=0.$(10) Then the non-zero solution of (10) is ${\lambda }_H^{\ast }=\frac{m_4-m_5}{m_6}$, where $m_4={\text{β}}_1{k}_1\mathrm{\Delta }{m}_1{m}_3{m}_3\mu$ + ${\text{β}}_1{\alpha }_2{k}_3\mathrm{\Delta }{m}_3{m}_3\mu$ +${\text{β}}_1{\rho }_1{k}_1\mathrm{\Delta }\theta m_1{m}_3\mu +$ ${\text{β}}_1{\rho }_1{\alpha }_2{k}_3\mathrm{\Delta }\theta m_3\mu$, $m_5=k_1\mathrm{\Delta }{m}_1{m}_2{m}_3\mu +{\alpha }_2{k}_3\mathrm{\Delta }{m}_2{m}_3\mu +k_3\mathrm{\Delta }{m}_2{m}_3\mu \mu$, $m_6=k_3\mathrm{\Delta }{m}_2{m}_3\mu$+$k_1\mathrm{\Delta }{m}_1{m}_3\mu +{\alpha }_2{k}_3\mathrm{\Delta }{m}_3\mu +k_1\mathrm{\Delta }\theta m_1\mu +{\alpha }_2{k}_3\mathrm{\Delta }\theta \mu +k_1\mathrm{\Delta }\theta \gamma m_1+{\alpha }_2{k}_3\mathrm{\Delta }\theta \gamma$.

Therefore, the required non-zero solution (force of infection) is obtained as ${\lambda }_H^{\ast }=\frac{\left[k_1\mathrm{\Delta }{m}_1{m}_2{m}_3\mu +k_3\mathrm{\Delta }{m}_2{m}_3\mu \left({\alpha }_2+\mu \right)\right]\left({\mathcal{R}}_{HM}-1\right)}{\begin{array}{c}{k}_3\mathrm{\Delta }{m}_2{m}_3\mu +k_1\mathrm{\Delta }{m}_1{m}_3\mu +{\alpha }_2{k}_3\mathrm{\Delta }{m}_3\mu +k_1\mathrm{\Delta }\theta m_1\mu +{\alpha }_2{k}_3\mathrm{\Delta }\mathrm{\Theta }\mu \\ +k_1\mathrm{\Delta }\theta \gamma m_1+{\alpha }_2{k}_3\mathrm{\Delta }\theta \gamma \end{array}}.$Then we do have ${\lambda }_H^{\ast }>0$ whenever ${\mathcal{R}}_{HM}>1$. Thus the HIV/AIDS mono-infection model (7) has a unique positive endemic equilibrium point if and only if ${\mathcal{R}}_{HM}>1$.

Theorem 3.2:

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

3.1.3. Global stability of disease-free equilibrium point

We now list two conditions which are sufficient to guarantee the global stability of the HIV only disease-free equilibrium point. Applying Castillo-Chavez et al. [25], let us rewrite the dynamical system (7) as $\dot{X}=F\left(X,Z\right)$$\dot{Z}=G\left(X,Z\right)$, $G\left(X,0\right)=0$, where $X=\left(D_1,D_3\right)$ and $Z=\left(D_6,D_7,D_{11}\right)$, with $X={\mathrm{\Re }}_{+}^2$ denoting the total number of uninfected individuals and $Z={\mathrm{\Re }}_{+}^3$ denoting the total number of infected individuals. The conditions ($H_1$) and $\left(H_2\right)$ below must met to guarantee global asymptotically stability.

$\left(H_1\right)$ for $\dot{X}=F\left(X,0\right)$, $X_0$ is globally asymptotically stable.

$\left(H_2\right)$ $\text{G}\left(X,Z\right)=\text{A}Z-\hat{\mathrm{G}}\left(X,Z\right),\hat{\mathrm{G}}\left(X,Z\right)\ge 0$ for $\left(X,Z\right)\in {\mathrm{\Omega }}_1$, where $A=D_{Z}G\left(X_0,0\right)$ is a Metzler matrix (the off diagonal elements of $A$ are nonnegative) and ${\mathrm{\Omega }}_1$ is the region where the model makes biological sense.

Theorem 3.3:

The fixed point $E_{HM}^0=\left(X_0,0\right)$ is a globally asymptotically stable equilibrium point of (7) provided ${\mathcal{R}}_{HM}<1$ and the assumptions ($H_1$) and $\left(H_2\right)$ are satisfied otherwise unstable.

Proof:

Consider $\dot{X}=F\left(X,Z\right),$ $\dot{Z}=G\left(X,Z\right),$where $F_1\left(X,Z\right)=k_1\mathrm{\Delta }+{\alpha }_2{D}_3-\left({\lambda }_H+\mu \right)D_1$ and $F_1\left(\mathcal{X},0\right)=\left[k_1\mathrm{\Delta }+{\alpha }_2{D}_3-\mu D_1\right]$,

$F_2\left(X,Z\right)=k_3\mathrm{\Delta }-\left({\alpha }_2+\mu \right)D_3$ and $F_1\left(\mathcal{X},0\right)=\left[k_3\mathrm{\Delta }-\left({\alpha }_2+\mu \right)D_3\right]$.

Then, $G\left(\mathcal{X},\mathcal{Z}\right)=A\mathcal{Z}-\hat{G}\left(\mathcal{X},\mathcal{Z}\right)=\left(\begin{array}{c}{\lambda }_H{D}_1-\left(\theta +\mu +d_2\right)D_6\\ \theta D_6-\left(\gamma +d_3+\mu \right)D_7\\ \gamma D_7-\mu D_{11}\end{array}\right),$where $A=\left[\begin{array}{ccc}-\left(\theta +\mu +d_2\right)& \frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}& 0\\ \theta & -\left(\gamma +d_3+\mu \right)& 0\\ 0& \gamma & -\mu \end{array}\right],\mathcal{Z}=\left(\begin{array}{c}{D}_6\\ D_7\\ D_{11}\end{array}\right),$and $A\mathcal{Z}=\left(\begin{array}{c}\frac{{\text{β}}_1{D}_1^0}{D_1^0+D_3^0}{D}_6-\left(\theta +\mu +d_2\right)D_6+\frac{{\text{β}}_1{\rho }_1{D}_1^0}{D_1^0+D_3^0}{D}_7\\ \theta D_6-\left(\gamma +d_3+\mu \right)D_7\\ \gamma D_7-\mu D_{11}\end{array}\right).$We get $\hat{G}\left(\mathcal{X},\mathcal{Z}\right)=\left(\begin{array}{c}\widehat{G_1}\left(\mathcal{X},\mathcal{Z}\right)\\ \widehat{G_2}\left(\mathcal{X},\mathcal{Z}\right)\\ \widehat{G_3}\left(\mathcal{X},\mathcal{Z}\right)\end{array}\right)=\left(\begin{array}{c}\frac{{\text{β}}_1\left(D_6+{\rho }_1{D}_7\right)D_1^0}{D_1^0+D_3^0}+\frac{{\text{β}}_1\left(D_6+{\rho }_1{D}_7\right)D_1}{D_1+D_3}\\ 0\\ 0\end{array}\right)$implies that $\widehat{{\text{G}}_1}\left(\mathcal{X},\mathcal{Z}\right)\ge 0$ in the region ${\mathrm{\Omega }}_1$ whenever ${\mathcal{R}}_{HM}<1$. This implies that $E_{HM}^0$ is globally asymptotically stable whenever ${\mathcal{R}}_{HM}<1$.

3.2. Pneumonia mono-infection model analysis

The corresponding pneumonia mono-infection model of the system (3) is determined by making $D_3=$ $D_6=D_7=D_8=D_9$ = $D_{11}=D_{12}=0$, and it is given by (11) $\begin{array}{rl}\stackrel{}{{\dot{D}}_1}& =k_1\mathrm{\Delta }+{\alpha }_1{D}_2+\rho D_4+\eta D_{10}-\left({\lambda }_P+\mu \right)D_1,\\ \stackrel{}{{\dot{D}}_2}& =k_2\mathrm{\Delta }-\left({\alpha }_1+\mu \right)D_2,\\ \stackrel{}{{\dot{D}}_4}& =k_4\mathrm{\Delta }-\left(\rho +\mu +ϵ{\lambda }_P\right)D_4,\\ \stackrel{}{{\dot{D}}_5}& ={\lambda }_P{D}_1+ϵ{\lambda }_P{D}_4-\left(\mu +d_1+\kappa \right)D_5,\\ \stackrel{}{{\dot{D}}_{10}}& =\kappa D_5-\left(\mu +\eta \right)D_{10},\end{array}$(11) with pneumonia infection initial conditions $D_1\left(0\right)>0$, $D_2\left(0\right)\ge 0,D_4\left(0\right)\ge 0,$ $D_5\left(0\right)\ge 0$, $D_{10}\left(0\right)\ge 0$, the total population is given by $N_2\left(t\right)=D\left(t\right)+D_2\left(t\right)$+$D_4\left(t\right)+D_5\left(t\right)+D_{10}\left(t\right),$ and pneumonia force of infection given by ${\lambda }_P={\text{β}}_2{D}_5\left(t\right)$.

Here like the full model (3) and the HIV/AIDS sub-model (7) in the region ${\mathrm{\Omega }}_2=\left\{\left(D_1,D_2,D_4,D_5,D_{10}\right)\in {{\mathbb{R}}^5}_{+},N_2\le \frac{\mathrm{\Delta }}{\mu }\right\}$, it is sufficient to consider the dynamics of model (11) in ${\mathrm{\Omega }}_2$ be both biologically and mathematically meaningful.

3.2.1. Local stability of disease-free equilibrium point

Disease-free equilibrium point of the pneumonia mono-infection model (11) is obtained by making its right-hand side as zero and setting the infected class and isolation with treatment class to zero as $D_5=D_{10}=0$ and after some simple steps of calculations we have got $D_1^0=\frac{k_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+k_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)}$, $D_2^0=\frac{k_2\mathrm{\Delta }}{{\alpha }_1+\mu }$, and $D_4^0=\frac{k_4\mathrm{\Delta }}{\rho +\mu }$. Hence the pneumonia mono-infection model (11), disease-free equilibrium point is given by $\begin{array}{rl}{E}_{PM}^0& =(D_1^0,D_2^0,D_3^0,D_4^0,D_5^0,D_{10}^0)\\ & =(\frac{k_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+k_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)},\\ & \frac{k_2\mathrm{\Delta }}{{\alpha }_1+\mu },\frac{k_4\mathrm{\Delta }}{\rho +\mu },0,0).\end{array}$Here we apply the van den Driesch and Warmouth next-generation matrix approach [24] to determine the pneumonia mono-infection model (11) effective reproduction number ${\mathcal{R}}_{PM}$. After long computations, we have got the transmission matrix given by $F=\left[\begin{array}{cc}{\text{β}}_2{D}_1^0+ϵ{\text{β}}_2{D}_4^0& 0\\ 0& 0\end{array}\right],$and the transition matrix given by $V=\left[\begin{array}{c}\mu +d_1+\kappa \\ -\kappa \end{array}\begin{array}{c}0\\ \mu +\eta \end{array}\right].$Then we have determined that $V^{-1}=\left[\begin{array}{cc}\frac{1}{\mu +d_1+\kappa }& 0\\ \frac{\kappa }{\left(\mu +d_1+\kappa \right)\left(\mu +\eta \right)}& \frac{1}{\mu +\eta }\end{array}\right],$and $F{V}^{-1}=\left[\begin{array}{cc}\frac{{\text{β}}_2{D}_2^0+ϵ{\text{β}}_2{D}_4^0}{\mu +d_1+\kappa }& 0\\ 0& 0\end{array}\right].$Then the spectral radius (effective reproduction number ${\mathcal{R}}_{PM}$) of $F{V}^{-1}$ of the pneumonia mono-infection model (11) is $\begin{array}{rl}{\mathcal{R}}_{PM}& =\frac{{\text{β}}_2{D}_1^0+ϵ{\text{β}}_2{D}_4^0}{\mu +d_1+\kappa }\\ & =\frac{\begin{array}{c}{\text{β}}_2{k}_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\text{β}}_2{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)\\ +{\text{β}}_2{k}_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)+{\text{β}}_2ϵk_4\mathrm{\Delta }\mu \left({\alpha }_1+\mu \right)\end{array}}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right)}.\end{array}$

Theorem 3.4:

The disease-free equilibrium point (DFE) $E_{PM}^0$ of the pneumonia mono-infection model (11) is locally asymptotically stable if ${\mathcal{R}}_{PM}<1$ otherwise unstable.

Proof:

The local stability of the disease-free equilibrium of the system (11) at point $\begin{array}{rl}{E}_{PM}^0& =(\frac{k_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+k_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)},\\ & \frac{k_2\mathrm{\Delta }}{{\alpha }_1+\mu },\frac{k_4\mathrm{\Delta }}{\rho +\mu },0,0)\end{array}$can be studied from its Jacobian matrix and Routh–Hurwitz stability criteria. The Jacobian matrix of the dynamical system at the disease-free equilibrium point is given by $\begin{array}{rl}& J\left(E_{PM}^0\right)\\ & =\left[\begin{array}{ccc}-\mu & {\alpha }_1& \rho \\ 0& -\left({\alpha }_1+\mu \right)& 0\\ 0& 0& -\left(\rho +\mu \right)\\ 0& 0& 0\\ 0& 0& 0\end{array}\begin{array}{cc}-{\text{β}}_2{D}_1^0& \eta \\ 0& 0\\ -{\text{β}}_2ϵD_4^0& 0\\ {\text{β}}_2{D}_1^0+{\text{β}}_2ϵD_4^0-\left(\mu +d_1+\kappa \right)& 0\\ \kappa & -\left(\mu +\eta \right)\end{array}\right].\end{array}$Then the characteristic equation of the above Jacobian matrix is given by $\left|\begin{array}{ccc}-\mu & {\alpha }_1& \rho \\ 0& -\left({\alpha }_1+\mu \right)& 0\\ 0& 0& -\left(\rho +\mu \right)\\ 0& 0& 0\\ 0& 0& 0\end{array}\begin{array}{cc}-{\text{β}}_2{D}_1^0& \eta \\ 0& 0\\ -{\text{β}}_2ϵD_4^0& 0\\ M& 0\\ \kappa & -\left(\mu +\eta \right)\end{array}\right|=0,$where $M={\text{β}}_2{D}_1^0+{\text{β}}_2ϵD_4^0-\left(\mu +d_1+\kappa \right)$ and after some steps of computations, we have got ${\lambda }_1=-\mu <0$ or ${\lambda }_2=-\left({\alpha }_1+\mu \right)<0$ or ${\lambda }_3=-\left(\rho +\mu \right)<0$ or ${\lambda }_4={\text{β}}_2{D}_1^0+{\text{β}}_2ϵD_4^0-\left(\mu +d_1+\kappa \right)=\left(\mu +d_1+\kappa \right)\left[\frac{{\text{β}}_2{D}_1^0+{\text{β}}_2ϵD_4^0}{\mu +d_1+\kappa }-1\right]=\mu +d_1+\kappa \left)\right[{\mathcal{R}}_{PM}-1]<0$ if ${\mathcal{R}}_{PM}<1$ or ${\lambda }_5=-\left(\mu +\eta \right)<0$. Therefore, since all the eigenvalues of the characteristics polynomials of the system (11) are negative if ${\mathcal{R}}_{PM}<1$ then disease-free equilibrium point of the pneumonia mono-infection model (11) is locally asymptotically stable.

3.2.2. Existence of endemic equilibrium point (s)

Before checking the global stability of the disease-free equilibrium point of the pneumonia mono-infection model (11), we shall find the possible number of endemic equilibrium point(s) of the model (11). Let $E_P^{\ast }=\left(D_1^{\ast },D_2^{\ast },D_4^{\ast },D_5^{\ast },D_{10}^{\ast }\right)$ be the endemic equilibrium point of pneumonia mono-infection and ${\lambda }_P^{\ast }={\text{β}}_2{D}_5^{\ast }$ be the pneumonia mono-infection mass action incidence rate (‘force of infection’) at the equilibrium point. To find equilibrium point(s) for which pneumonia mono-infection is endemic in the population, the equations are solved in terms of ${\lambda }_P^{\ast }={\text{β}}_2{D}_5^{\ast }$ at an endemic equilibrium point. Now setting the right-hand side of the equations of the model to zero (at steady state) gives $\begin{array}{rl}{D}_1^{\ast }& =\frac{b_5{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2+b_6{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2+b_7\left(b_2+ϵ{\lambda }_P^{\ast }\right)+b_8{\lambda }_P^{\ast }}{b_1{b}_3{b}_4{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2\left({\lambda }_P^{\ast }+\mu \right)-b_1\eta \kappa {\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }},\\ D_2^{\ast }& =\frac{k_2\mathrm{\Delta }}{b_1},D_4^{\ast }=\frac{k_4\mathrm{\Delta }}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)},\\ D_5^{\ast }& =\frac{b_5{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }+b_6{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }+\left(b_2{b}_7{\lambda }_P^{\ast }+b_7ϵ\lambda {{}_P^{\ast }}^2\right)}{b_{12}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2\left({\lambda }_P^{\ast }+\mu \right)-b_{13}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }}\\ & +\frac{b_8\lambda {{}_P^{\ast }}^2+\left(b_2{b}_9+b_9ϵ{\lambda }_P^{\ast }\right)\left(\lambda {{}_P^{\ast }}^2+\mu {\lambda }_P^{\ast }\right)-b_{10}\lambda {{}_P^{\ast }}^2-b_{11}\lambda {{}_P^{\ast }}^3}{b_{12}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2\left({\lambda }_P^{\ast }+\mu \right)-b_{13}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }},\end{array}$And $D_{10}^{\ast }=\frac{\kappa D_5^{\ast }}{b_4},$where $b_1={\alpha }_1+\mu$, $b_2=\rho +\mu$, $b_3=\mu +d_1+\kappa$, $b_4=\mu +\eta$, $b_5=k_1\mathrm{\Delta }{b}_1{b}_3{b}_4$, $b_6={\alpha }_1{k}_2\mathrm{\Delta }{b}_3{b}_4$, $b_7=\rho k_4\mathrm{\Delta }{b}_1{b}_3{b}_4$, $b_8=k_4\mathrm{\Delta }{b}_1\eta \kappa ϵ$, $b_9=b_1{b}_3{b}_4{k}_4\mathrm{\Delta }ϵ$, $b_{10}=b_2{b}_1\eta \kappa k_4\mathrm{\Delta }ϵ$, $b_{11}=b_1{k}_4\mathrm{\Delta }ϵ\eta \kappa ϵ$, $b_{12}=b_1{b}_3{b}_3{b}_4$, $b_{13}=b_1{b}_3\eta \kappa$.

Then we have substituted $\begin{array}{rl}{D}_5^{\ast }& =\frac{b_5{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }+b_6{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }+\left(b_2{b}_7{\lambda }_P^{\ast }+b_7ϵ\lambda {{}_P^{\ast }}^2\right)}{b_{12}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2\left({\lambda }_P^{\ast }+\mu \right)-b_{13}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }}\\ & +\frac{b_8\lambda {{}_P^{\ast }}^2+\left(b_2{b}_9+b_9ϵ{\lambda }_P^{\ast }\right)\left(\lambda {{}_P^{\ast }}^2+\mu {\lambda }_P^{\ast }\right)-b_{10}\lambda {{}_P^{\ast }}^2-b_{11}\lambda {{}_P^{\ast }}^3}{b_{12}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2\left({\lambda }_P^{\ast }+\mu \right)-b_{13}{\left(b_2+ϵ{\lambda }_P^{\ast }\right)}^2{\lambda }_P^{\ast }}\end{array}$in the pneumonia force of infection given by ${\lambda }_P^{\ast }={\text{β}}_2{D}_5^{\ast }$ and we have computed the non-zero solution of ${\lambda }_P^{\ast }$ from the cubic equation (12) $f_2\lambda {}_P^{\ast }2+f_1{\lambda }_P^{\ast }+f_0=0,$(12) where (13) $\begin{array}{l}{f}_2=b_{12}{ϵ}^2-b_{13}{ϵ}^2>0,\\ f_1=2b_2{b}_{12}ϵ+b_{12}\mu {ϵ}^2-2b_2{b}_{13}ϵ-b_5{ϵ}^2-b_6{ϵ}^2-b_9ϵ+b_{11}\end{array}$(13) $f_0=b_1{b}_2{b}_3{b}_4\left[1-{\mathcal{R}}_{PM}\right]>0$ if ${\mathcal{R}}_{PM}<1$.

It can be seen from (13) that $f_2>0$ (since the entire model parameters are non-negative). Furthermore, $f_0>0$ whenever ${\mathcal{R}}_{PM}<1$. Thus the number of possible positive real roots the polynomial (12) can have depends on the sign of f₁. This can be analysed using the Descartes’ rule of signs on the cubic $g\left(x\right)$ = $f_2{x}^2+f_{1}x+f_0$ (with x = ${\lambda }_p^{\ast }$). Hence, the following results are established.

Theorem 3.5:

The pneumonia mono-infection model (11) could have

1. a unique endemic equilibrium point if ${\mathcal{R}}_{PM}>1$ either of the following holds:

   1. $f_1>0$.

   2. $f_1<0$.

2. two endemic equilibrium points if ${\mathcal{R}}_{PM}<1$, f₁<0.

Here, the statement in (b) shows the happening of the backward bifurcation in the model (11), i.e. the locally asymptotically stable disease-free equilibrium point co-exists with a locally asymptotically stable endemic equilibrium point if ${\mathcal{R}}_{PM}<1$; examples of the existence of backward bifurcation phenomenon in mathematical epidemiological models, and the causes, can be seen in [3,21,23,27,39]. The epidemiological consequence is that the classical epidemiological requirement of having the reproduction number ${\mathcal{R}}_{PM}$ to be less than 1, even though necessary, is not sufficient for the effective control of the disease. The existence of the backward bifurcation phenomenon in sub-model (11) is now explored.

Theorem 3.6:

The pneumonia mono-infection model (11) exhibits backward bifurcation at ${\mathcal{R}}_{PM}=1$ whenever the inequality $H_2>H_1$ holds, whenever H₁ = $\frac{-{\text{β}}_2{\text{β}}^{\ast }{x}_1^0\left(\rho +\mu \right)\left(\mu +\eta \right)-{\text{β}}_2{\text{β}}^{\ast }ϵx_3^0\rho \left(\mu +\eta \right)-{\text{β}}_2ϵ{\text{β}}^{\ast }ϵx_3^0\mu \left(\mu +\eta \right)}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)}$ and $H_2=\frac{{\text{β}}_2\kappa \eta \left(\rho +\mu \right)}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)}$.

Proof:

Appendix C.

3.2.3. Global stability of the model disease-free equilibrium point (DFE)

We now list two conditions which are sufficient to guarantee the global stability of the pneumonia only disease-free equilibrium point. Following Castillo-Chavez et al. [25], let us rewrite the dynamical system (11) as $\dot{X}=F\left(X,Z\right)$ $\dot{Z}=G\left(X,Z\right),G\left(X,0\right)=0,$where $X=\left(D_1,D_2,D_4\right)$ and $Z=\left(D_5,D_{10}\right)$, with $X={\mathrm{\Re }}_{+}^3$ denoting the total number of uninfected individuals and $Z={\mathrm{\Re }}_{+}^3$ denoting the total number of infected individuals. The condition ($H_1$) and $\left(H_2\right)$ below must met to guarantee global asymptotically stability.

$\left(H_1\right)$ for $\dot{X}=F\left(X,0\right)$, $X_0$ is globally asymptotically stable.

$\left(H_2\right)$ $\text{G}\left(X,Z\right)=\text{A}Z-\hat{\mathrm{G}}\left(X,Z\right),\hat{\mathrm{G}}\left(X,Z\right)\ge 0$ for $\left(X,Z\right)\in {\mathrm{\Omega }}_2$, where $A=D_{Z}G\left(X_0,0\right)$ is a Metzler matrix (the off diagonal elements of $A$ are nonnegative) and ${\mathrm{\Omega }}_1$ is the region where the model makes biological sense.

Theorem 3.7:

The fixed point $E_{PM}^0=\left(X_0,0\right)$ is a globally asymptotically stable equilibrium point of (11) provided ${\mathcal{R}}_{PM}<{\mathcal{R}}_C$ and the assumptions (H₁) and (H₂) are satisfied otherwise unstable.

Proof:

Consider $\dot{X}=F\left(X,Z\right),$ $\dot{Z}=G\left(X,Z\right),$where $F_1\left(X,Z\right)=k_1\mathrm{\Delta }+{\alpha }_1{D}_2+\rho D_4+\eta D_{10}-\left({\lambda }_P+\mu \right)D_1$ and $F_1\left(\mathcal{X},0\right)=\left[k_1\mathrm{\Delta }+{\alpha }_1{D}_2+\rho D_4+\eta D_{10}-\mu D_1\right]$, $F_2\left(X,Z\right)=k_2\mathrm{\Delta }-\left({\alpha }_1+\mu \right)D_2$ and $F_2\left(\mathcal{X},0\right)=\left[k_2\mathrm{\Delta }-\left({\alpha }_1+\mu \right)D_2\right]$, $F_3\left(X,Z\right)=k_4\mathrm{\Delta }-\left(\rho +\mu \right)D_4$ and $F_3\left(\mathcal{X},0\right)=\left[k_4\mathrm{\Delta }-\left(\rho +\mu \right)D_4\right]$.

Then, $G\left(\mathcal{X},\mathcal{Z}\right)=A\mathcal{Z}-\hat{G}\left(\mathcal{X},\mathcal{Z}\right)=\left(\begin{array}{c}{\lambda }_P{D}_1+ϵ{\lambda }_P{D}_4-\left(\mu +d_1+\kappa \right)D_5\\ \kappa D_5-\left(\mu +\eta \right)D_{10}\end{array}\right)$, where, $A=\left[\begin{array}{cc}{\text{β}}_2{D}_1^0+{\text{β}}_2ϵD_4^0-\left(\mu +d_1+\kappa \right)& 0\\ \kappa & -\left(\mu +\eta \right)\end{array}\right]$, $\mathcal{Z}=\left(\begin{array}{c}{D}_5\\ D_{10}\end{array}\right)$, and $A\mathcal{Z}=\left(\begin{array}{c}{\text{β}}_2{D}_1^0{D}_5+{\text{β}}_2ϵD_4^0{D}_5-\left(\mu +d_1+\kappa \right)D_5\\ \kappa D_5-\left(\mu +\eta \right)D_{10}\end{array}\right)$.

We get, $\hat{G}\left(\mathcal{X},\mathcal{Z}\right)=\left(\begin{array}{c}\widehat{G_1}\left(\mathcal{X},\mathcal{Z}\right)\\ \widehat{G_2}\left(\mathcal{X},\mathcal{Z}\right)\end{array}\right)=\left(\begin{array}{c}{\text{β}}_2{D}_1^0{D}_5+{\text{β}}_2ϵD_4^0{D}_5-{\lambda }_P{D}_1+ϵ{\lambda }_P{D}_4\\ 0\end{array}\right)$.

Thus $\widehat{{\text{G}}_1}\left(\mathcal{X},\mathcal{Z}\right)\ge 0$ in the region ${\mathrm{\Omega }}_1$, whenever $\begin{array}{rl}& D_1+ϵD_4\\ & \le \frac{k_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+k_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)+k_4\mathrm{\Delta }\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)},\end{array}$that is whenever ${\mathcal{R}}_{PM}<{\mathcal{R}}_C$, where ${\mathcal{R}}_C=\frac{\mu \left({\alpha }_1+\mu \right)\left(\mu +d_1+\kappa \right)+2{\text{β}}_2{\alpha }_1{k}_2\mathrm{\Delta }}{\mu \left({\alpha }_1+\mu \right)\left(\mu +d_1+\kappa \right)}$. This implies that $E_{HM}^0$ is globally asymptotically stable whenever ${\mathcal{R}}_{PM}<{\mathcal{R}}_C$.

3.3. Pneumonia and HIV/AIDS co-infection model analysis

3.3.1. Disease-free equilibrium point

The disease-free equilibrium point of the dynamical system (3) whenever the state variable $D_5=D_6=D_7=D_8=D_9=0$ is given by $\begin{array}{rl}{E}_0& =\left(D_1^0,D_2^0,D_3^0,D_4^0,D_5^0,D_6^0,D_7^0,D_8^0,D_9^0,D_{12}^0,D_{12}^0\right)\\ & =\left(\frac{k_1\mathrm{\Delta }}{\mu }+\frac{{\alpha }_1{k}_2\mathrm{\Delta }}{{\alpha }_1+\mu }+\frac{{\alpha }_2{k}_3\mathrm{\Delta }}{{\alpha }_2+\mu }+\frac{\rho k_4\mathrm{\Delta }}{\rho +\mu },\frac{k_2\mathrm{\Delta }}{{\alpha }_1+\mu },\frac{k_3\mathrm{\Delta }}{{\alpha }_2+\mu },\frac{k_4\mathrm{\Delta }}{\rho +\mu },0,0,0,0,0,0,0\right).\end{array}$

3.3.2. Effective reproduction number

Using similar approach to Sections 3.1.1 and 3.2.1, the largest (dominant) eigenvalue (spectral radius) of the matrix: $F{V}^{-1}=\left[\frac{\mathrm{\partial }{\mathcal{F}}_i\left(E_0\right)}{\mathrm{\partial }{X}_j}\right]\left[\frac{\mathrm{\partial }{\nu }_i\left(E_0\right)}{\mathrm{\partial }{X}_j}\right]$, where ${\mathcal{F}}_i$ is the rate of appearance of new infection in compartment i, ${\nu }_i$ is the transfer of infections from one compartment $i$ to another and $E_0$ is the disease-free equilibrium point. After some steps of calculations we have got $F=\left[\begin{array}{ccccccc}\mathrm{\aleph }& 0& {\omega }_1\mathrm{\aleph }& {\omega }_2\mathrm{\aleph }& 0& 0& 0\\ 0& {\mathrm{\aleph }}_2& {\rho }_2{\mathrm{\aleph }}_2& {\rho }_3{\mathrm{\aleph }}_2& {\rho }_1{\mathrm{\aleph }}_2& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right],$and $\begin{array}{rl}V& =[\begin{array}{ccc}\left(\mu +d_1+\kappa \right)& 0& 0\\ 0& \left(\theta +\mu +d_2\right)& -{\theta }_1\\ 0& 0& \left(\mu +d_4+\delta +{\theta }_1\right)\\ 0& 0& -\delta \\ 0& -\theta & 0\\ -\kappa & 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \left(\mu +d_5+{\theta }_2\right)& 0& 0& 0\\ -{\theta }_2& \left(\gamma +d_3+\mu \right)& 0& 0\\ 0& 0& \left(\mu +\eta \right)& 0\\ 0& -\gamma & 0& \mu \end{array}],\end{array}$where ${\mathrm{\aleph }}_1={\text{β}}_2{D}_1+{\text{β}}_2{D}_3+ϵ{\text{β}}_2{D}_4,{\mathrm{\aleph }}_2=\frac{{\text{β}}_1}{N}\left(D_1+D_2+D_4\right)$ and we have computed as $F{V}^{-1}=\left[\begin{array}{ccccccc}\frac{{\mathrm{\aleph }}_1}{\left(\mu +d_1+\kappa \right)}& a_{21}& 0& 0& 0& 0& 0\\ 0& \frac{{\mathrm{\aleph }}_2}{\left(\theta +\mu +d_2\right)}+\frac{{\rho }_1{\mathrm{\aleph }}_2\theta }{\left(\theta +\mu +d_2\right)\left(\gamma +d_3+\mu \right)}& 0& 0& 0& 0& 0\\ 0& a_{23}& 0& 0& 0& 0& 0\\ 0& a_{24}& 0& 0& 0& 0& 0\\ 0& a_{25}& 0& 0& 0& 0& 0\\ 0& a_{26}& 0& 0& 0& 0& 0\\ 0& a_{27}& 0& 0& 0& 0& 0\end{array}\right].$After some computations and simplifications, the HIV/AIDS and pneumonia co-infection model effective reproduction number is given by $\begin{array}{rl}{\mathcal{R}}_{HP}& =\text{max}\left\{{\mathcal{R}}_{CM},{\mathcal{R}}_{HM}\right\}\\ & =\text{max}\{\frac{\begin{array}{c}{\text{β}}_2{k}_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\text{β}}_2{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)\\ +{\text{β}}_2{k}_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)+{\text{β}}_2ϵk_4\mathrm{\Delta }\mu \left({\alpha }_1+\mu \right)\end{array}}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right)},\end{array}$$\frac{{\text{β}}_1\left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\alpha }_2{k}_3}{\left({\alpha }_2+\mu \right)\left(\theta +\mu +d_2\right)}+\frac{{\text{β}}_1{\rho }_1\theta \left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\rho }_1\theta {\alpha }_2{k}_3}{\left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}\}$ where ${\mathcal{R}}_{PM}=\frac{{\text{β}}_2{D}_1^0+ϵ{\text{β}}_2{D}_4^0}{\mu +d_1+\kappa }=\frac{{\text{β}}_2{k}_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\text{β}}_2{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+{\text{β}}_2{k}_4\mathrm{\Delta }\rho \left({\alpha }_1+\mu \right)+{\text{β}}_2ϵk_4\mathrm{\Delta }\mu \left({\alpha }_1+\mu \right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right)}$ is the pneumonia effective reproduction number and ${\mathcal{R}}_{HM}=\frac{{\text{β}}_1\left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\alpha }_2{k}_3}{\left({\alpha }_2+\mu \right)\left(\theta +\mu +d_2\right)}+\frac{{\text{β}}_1{\rho }_1\theta \left(1-k_3\right)\left({\alpha }_2+\mu \right)+{\text{β}}_1{\rho }_1\theta {\alpha }_2{k}_3}{\left(\theta +\mu +d_2\right)\left(\gamma +\mu +d_3\right)}$ is the HIV/AIDS effective reproduction number.

3.3.3. Locally stability of the disease-free equilibrium (DFE)

The Jacobian matrix of the system (3) at the DFE given by $E_0$ is computed as $\begin{array}{rl}J\left(E_0\right)& =(\begin{array}{cccccc}-\mu & {\alpha }_1& {\alpha }_2& \rho & -{\text{β}}_2{D}_1^0& -\frac{{\text{β}}_1}{N}{D}_1^0\\ 0& -\left({\alpha }_1+\mu \right)& 0& 0& 0& -\frac{{\text{β}}_1}{N}{D}_2^0\\ 0& 0& -\left({\alpha }_2+\mu \right)& 0& -{\text{β}}_2{D}_3^0& 0\\ 0& 0& 0& -\left(\rho +\mu \right)& -ϵ{\text{β}}_2{D}_4^0& -\frac{{\text{β}}_1}{N}{D}_4^0\\ 0& 0& 0& 0& {\mathrm{\aleph }}_3& 0\\ 0& 0& 0& 0& 0& {\mathrm{\aleph }}_8\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\\ & \begin{array}{ccccc}-\frac{{\text{β}}_1}{N}{D}_1^0{\rho }_1& -{\mathrm{\aleph }}_4& -{\mathrm{\aleph }}_6& \eta & 0\\ -\frac{{\text{β}}_1}{N}{D}_2^0{\rho }_1& -\frac{{\text{β}}_1}{N}{D}_2^0{\rho }_2& -\frac{{\text{β}}_1}{N}{D}_2^0{\rho }_3& 0& 0\\ 0& -{\text{β}}_2{D}_3^0{\omega }_1& -{\text{β}}_2{D}_3^0{\omega }_2& 0& 0\\ -\frac{{\text{β}}_1}{N}{D}_4^0{\rho }_1& -{\mathrm{\aleph }}_5& -{\mathrm{\aleph }}_7& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ {\mathrm{\aleph }}_9& 0& 0& 0& 0\\ 0& {\mathrm{\aleph }}_{10}& 0& 0& 0\\ 0& 0& {\mathrm{\aleph }}_{11}& 0& 0\\ 0& 0& 0& -\left(\mu +\eta \right)& 0\\ 0& 0& 0& 0& -\mu \end{array}),\end{array}$where ${\mathrm{\aleph }}_3={\text{β}}_2\left(D_1^0+D_3^0+ϵD_4^0\right)-\left(\mu +d_1+\kappa \right)$, ${\mathrm{\aleph }}_4=\left(\frac{{\text{β}}_1}{N}{\rho }_2+{\text{β}}_2{\omega }_1\right)D_1^0$, ${\mathrm{\aleph }}_5=\left(\frac{{\text{β}}_1}{N}{\rho }_2+ϵ{\text{β}}_2{\omega }_1\right)D_4^0$, ${\mathrm{\aleph }}_6=\left(\frac{{\text{β}}_1}{N}{\rho }_3+{\text{β}}_2{\omega }_2\right)D_1^0$ and ${\mathrm{\aleph }}_7=\left(\frac{{\text{β}}_1}{N}{\rho }_3+ϵ{\text{β}}_2{\omega }_2\right)D_4^0$, ${\mathrm{\aleph }}_8=\frac{{\text{β}}_1}{N}\left(D_1^0+D_2^0+D_4^0\right)-\left(\theta +\mu +d_2\right)$, ${\mathrm{\aleph }}_9=-\left(\gamma +d_3+\mu \right),{\mathrm{\aleph }}_{10}=-\left(\mu +d_4+\delta +{\theta }_1\right),{\mathrm{\aleph }}_{11}=-\left(\mu +d_5+{\theta }_2\right)$.

Then the eigenvalues of the matrix $J\left(E_0\right)$ are ${\lambda }_1=-\mu <0$ or ${\lambda }_2=-\left({\alpha }_1+\mu \right)<0$ or ${\lambda }_3=-\left({\alpha }_2+\mu \right)<0$ or ${\lambda }_4=-\left(\rho +\mu \right)<0$ or ${\lambda }_5=-\mu <0$ or ${\lambda }_6=-\left(\mu +\eta \right)<0$ or ${\lambda }_7=\frac{{\text{β}}_2ϵk_4\mathrm{\Delta }}{\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right)}\left({\mathcal{R}}_{CM}-1\right)<0$ or ${\lambda }_8=-\left(\mu +d_4+\delta +{\theta }_1\right)<0$ or ${\lambda }_9=-\left(\mu +d_5+{\mathrm{\Theta }}_2\right)<0$ or ${\lambda }^2+\left[\left(\gamma +d_3+\mu \right)+\left(\theta +\mu +d_2\right)-{\mathrm{\aleph }}_8\right]\lambda -\left[\left({\mathrm{\aleph }}_8-\left(\theta +\mu +{\text{d}}_2\right)\right)\left(\gamma +{\text{d}}_3+\mu \right)+\theta {\rho }_1{\mathrm{\aleph }}_8\right]=0$. Then after some calculations, we have determined that the last two eigenvalues of the quadratic equation as ${\lambda }_{10}<0$ and ${\lambda }_{11}<0$ whenever ${\mathcal{R}}_{\text{HP}}=max\left\{{\mathcal{R}}_{PM},{\mathcal{R}}_{HM}\right\}<1$. Thus, since all the 10 eigenvalues are negative, the disease-free equilibrium point of the full model (3) is locally asymptotically stable whenever ${\mathcal{R}}_{\text{HP}}=max\left\{{\mathcal{R}}_{PM},{\mathcal{R}}_{HM}\right\}<1$.

4. Analysis of the optimal control strategies

In this section, we analyse the optimal control strategies on the co-infection model (3) to identify the best control strategy that reduces the number of HIV/AIDS infected individuals, the number of pneumonia infected individuals and number of pneumonia and HIVAIDS Co-infected infected individuals. The objective is to find the optimal values ${\mathfrak{u}}^{\ast }=\left({\mathfrak{u}}_1^{\ast },{\mathfrak{u}}_2^{\ast },{\mathfrak{u}}_3^{\ast },{\mathfrak{u}}_4^{\ast }\right)$ of the controls $\mathfrak{u}=\left({\mathfrak{u}}_1,{\mathfrak{u}}_2,{\mathfrak{u}}_3,{\mathfrak{u}}_4\right)$ such that the associated state trajectories $D_1^{\ast },D_2^{\ast },D_3^{\ast },D_4^{\ast },D_5^{\ast },D_6^{\ast },D_7^{\ast },D_8^{\ast },D_9^{\ast },D_{10}^{\ast },D_{11}^{\ast }$ are solution of the system (3) in the intervention time interval $\left[0,T_f\right]$ with initial conditions as given in (4) and minimize the objective function. The control ${\mathfrak{u}}_1\left(t\right)$ represents the efforts on preventing HIV/AIDS infection using condom that helps to reduce contact rate HIV and ${\mathfrak{u}}_2\left(t\right)$ represents the efforts on preventing pneumonia infections that helps to reduce contact rate of pneumonia. Control to increase the recovery of pneumonia infected individuals which means to increase pneumonia recovered individuals ${\mathfrak{u}}_3\left(t\right)$ satisfies $0\le {\mathfrak{u}}_3\left(t\right)\le 1$, due to the efficacy of anti-pneumonia drug used for the treatment of pneumonia infected individuals and control on the treatment of HIV/AIDS to increase treated number of HIV/AIDS infected individuals ${\mathfrak{u}}_4\left(t\right)$ satisfies $0\le {\mathfrak{u}}_4\left(t\right)\le 1$, due to the efficacy of anti-HIV drugs used for the treatment of HIV/AIDS infected individuals.

Then the model (3) is now rewritten as (17) $\begin{array}{rl}\dot{D_1}& =k_1\mathrm{\Delta }+{\alpha }_1{D}_2+{\alpha }_2{D}_3+\rho D_4+\eta D_{10}-\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_1-\left(1-{\mathfrak{u}}_2\right){\lambda }_P{D}_1-\mu D_1,\\ \dot{D_2}& =k_2\mathrm{\Delta }-\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_2-\left({\alpha }_1+\mu \right)D_2,\\ \dot{D_3}& =k_3\mathrm{\Delta }-\left(1-{\mathfrak{u}}_2\right){\lambda }_P{D}_3-\left({\alpha }_2+\mu \right)D_3,\\ \dot{D_4}& =k_4\mathrm{\Delta }-\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_4-\left(1-{\mathfrak{u}}_2\right)ϵ{\lambda }_P{D}_4-\left(\rho +\mu \right)D_4,\\ \dot{D_5}& =\left(1-{\mathfrak{u}}_2\right){\lambda }_P{D}_1+\left(1-{\mathfrak{u}}_2\right){\lambda }_P{D}_3+\left(1-{\mathfrak{u}}_2\right)ϵ{\lambda }_P{D}_4-\left(1-{\mathfrak{u}}_1\right)\upsilon {\lambda }_H{D}_5\\ & -\left(\mu +d_1+{\mathfrak{u}}_3\kappa \right)D_5,\\ \dot{D_6}& =\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_1+\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_2+\left(1-{\mathfrak{u}}_1\right){\lambda }_H{D}_4+{\mathfrak{u}}_3{\theta }_1{D}_8\\ & -\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\lambda }_P{D}_6-\left(\theta +\mu +d_2\right)D_6,\\ \dot{D_7}& =\theta D_6+{\mathfrak{u}}_3{\theta }_2{D}_9-\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\lambda }_P{D}_7-\left({\mathfrak{u}}_4\gamma +d_3+\mu \right)D_7,\\ \dot{D_8}& =\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\lambda }_P{D}_6+\left(1-{\mathfrak{u}}_1\right)\upsilon {\lambda }_H{D}_5-\left(\mu +d_4+\delta +{\mathfrak{u}}_3{\theta }_1\right)D_8,\\ \dot{D_9}& =\delta D_8+\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\lambda }_P{D}_7-\left(\mu +d_5+{\mathfrak{u}}_3{\theta }_2\right)D_9,\\ \dot{D_{1}0}& ={\mathfrak{u}}_3\kappa D_5-\left(\mu +\eta \right)D_{10},\\ \dot{D_{1}1}& ={\mathfrak{u}}_4\gamma D_7-\mu D_{11},\end{array}$(17) with the corresponding initial conditions (18) $\begin{array}{rl}{D}_1\left(0\right)& >0,D_2\left(0\right)\ge 0,D_3\left(0\right)\ge 0,D_4\left(0\right)\ge 0,D_5\left(0\right)\ge 0,D_6\left(0\right)\ge 0,D_7\left(0\right)\ge 0,\\ D_8\left(0\right)& \ge 0,D_9\left(0\right)\ge 0,D_{10}>0,\text{and}{D}_{11}>0.\end{array}$(18) The optimal control problem is to minimize the objective functional (19) $\begin{array}{rl}& J\left({\mathfrak{u}}_1,{\mathfrak{u}}_2,{\mathfrak{u}}_3,{\mathfrak{u}}_4\right)\\ & =\underset{0}{\overset{T_f}{\int }}\left({\mathfrak{w}}_1{D}_5+{\mathfrak{w}}_2{D}_7+{\mathfrak{w}}_3{D}_8+{\mathfrak{w}}_4{D}_9+\frac{{\mathfrak{B}}_1}{2}{\mathfrak{u}}_1^2+\frac{{\mathfrak{B}}_2}{2}{\mathfrak{u}}_2^2+\frac{{\mathfrak{B}}_3}{2}{\mathfrak{u}}_3^2+\frac{{\mathfrak{B}}_4}{2}{\mathfrak{u}}_4^2\right)dt.\end{array}$(19) The coefficients ${\mathfrak{w}}_1,{\mathfrak{w}}_2,{\mathfrak{w}}_3,$ and ${\mathfrak{w}}_4$ are positive weight constants and $\frac{{\mathfrak{B}}_1}{2},\frac{{\mathfrak{B}}_2}{2},\frac{{\mathfrak{B}}_3}{2}$ and $\frac{{\mathfrak{B}}_4}{2}$ are the measure of relative costs of interventions associated with the controls ${\mathfrak{u}}_1,{\mathfrak{u}}_2,{\mathfrak{u}}_3$ and ${\mathfrak{u}}_4$, respectively, and also balances the units of integrand. In the cost functional, the term ${\mathfrak{w}}_1{D}_5$ refers to the cost related to pneumonia infected class, the term ${\mathfrak{w}}_2{D}_7$ refers to the cost related to individuals mono-infected with HIV and aware, the term ${\mathfrak{w}}_3{D}_8$ refers to the cost related to co-infected individuals unaware of HIV infection and the term ${\mathfrak{w}}_4{A}_9$ refers to the cost related to co-infected individuals aware of HIV infection.

$I\left(D_1,D_2,D_3,D_4,D_5,D_6,D_7,D_8,D_9,D_{10},D_{11},\mathfrak{u}\right)={\mathfrak{w}}_1{D}_5+{\mathfrak{w}}_2{D}_7+{\mathfrak{w}}_3{D}_8+{\mathfrak{w}}_4{D}_9+\frac{{\mathfrak{B}}_1}{2}{\mathfrak{u}}_1^2+\frac{{\mathfrak{B}}_2}{2}{\mathfrak{u}}_2^2+\frac{{\mathfrak{B}}_3}{2}{\mathfrak{u}}_3^2+\frac{{\mathfrak{B}}_4}{2}{\mathfrak{u}}_4^2$, measures the current cost at time t. The set of admissible control functions is defined by (20) ${\mathrm{\Omega }}_{\mathfrak{u}}=\left\{\left({\mathfrak{u}}_1\left(t\right),{\mathfrak{u}}_2\left(t\right),{\mathfrak{u}}_3\left(t\right),{\mathfrak{u}}_4\left(t\right)\right)\in L^4:0\le {\mathfrak{u}}_1\left(t\right),{\mathfrak{u}}_2\left(t\right),{\mathfrak{u}}_3\left(t\right),{\mathfrak{u}}_4\left(t\right)\le 1,t\in \left[0,T_f\right]\right\}.$(20) More precisely, we seek an optimal control pair (21) $J\left({\mathfrak{u}}_1^{\ast },{\mathfrak{u}}_2^{\ast },{\mathfrak{u}}_3^{\ast },{\mathfrak{u}}_4^{\ast }\right)=\underset{\mathrm{\Omega }}{min}J\left({\mathfrak{u}}_1,{\mathfrak{u}}_2,{\mathfrak{u}}_3,{\mathfrak{u}}_4\right).$(21)

4.1. Characterization of optimal control

In this section, we present optimality conditions for the optimal control problem defined above and detail its properties. According to the Pontryagin’s Maximum Principle as stated in [29], if ${\mathfrak{u}}^{\ast }\left(t\right)\in {\mathrm{\Omega }}_{\mathfrak{u}}$ is optimal for dynamical system (17) with initial value (18) and (21) with fixed final time $T_f$, then there exists a non-trivial absolutely continuous mapping

$\lambda :\left[0,T_f\right]\to {\mathbb{R}}^{14},\lambda =\left({\lambda }_1\right(t),{\lambda }_2(t),{\lambda }_3(t),{\lambda }_4(t),{\lambda }_5(t),{\lambda }_6(t),{\lambda }_7(t),{\lambda }_8(t),{\lambda }_9(t),{\lambda }_{10}(t),{\lambda }_{11}(t\left)\right)$ called the adjoint vector, such that

1. The Hamiltonian function is defined as (22) $\begin{array}{rl}\mathcal{H}& ={\mathfrak{w}}_1{D}_5+{\mathfrak{w}}_2{D}_7+{\mathfrak{w}}_3{D}_8+{\mathfrak{w}}_4{D}_9+\frac{{\mathfrak{B}}_1}{2}{\mathfrak{u}}_1^2+\frac{{\mathfrak{B}}_2}{2}{\mathfrak{u}}_2^2+\frac{{\mathfrak{B}}_3}{2}{\mathfrak{u}}_3^2+\frac{{\mathfrak{B}}_4}{2}{\mathfrak{u}}_4^2\\ & +\sum _{i=1}^{11}{\lambda }_i{\mathcal{G}}_i.\end{array}$(22) where ${\mathcal{G}}_i$ stands for the right-hand side of model (16) which is the ith-state variable equation.

2. The control system (23) $\begin{array}{rl}& \frac{d{D}_1}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_1},\frac{d{D}_2}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_2},\frac{d{D}_3}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_3},\frac{d{D}_4}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_4},\frac{d{D}_5}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_5},\\ & \frac{d{D}_6}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_6},\frac{d{D}_7}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_7},\frac{d{D}_8}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_8},\\ & \frac{d{D}_9}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_9},\frac{d{D}_{10}}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_{10}},\frac{d{D}_{11}}{dt}=\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{\lambda }_{11}}.\end{array}$(23)

3. The adjoint system (24) $\begin{array}{rl}& \frac{d{\lambda }_1}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_1},\frac{d{\lambda }_2}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_2},\frac{d{\lambda }_3}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_3},\\ & \frac{d{\lambda }_4}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_4},\frac{d{\lambda }_5}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_5},\frac{d{\lambda }_6}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_6},\frac{d{\lambda }_7}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_7},\\ & \frac{d{\lambda }_8}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_8},\frac{d{\lambda }_9}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_9},\frac{d{\lambda }_{10}}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_{10}},\frac{d{\lambda }_{11}}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_{11}}.\end{array}$(24)

4. And the optimality condition (25) $\mathcal{H}\left({\mathbb{E}}^{\ast },\mathfrak{u},{\lambda }^{\ast }\right)=\underset{\mathfrak{u}\in {\mathrm{\Omega }}_{\mathfrak{u}}}{min}\left({\mathbb{E}}^{\ast },{\mathfrak{u}}^{\ast },{\lambda }^{\ast }\right),$(25) holds for almost all $t\in \left[0,T_f\right]$.

5. Moreover, the transversality condition (26) ${\lambda }_i\left(T_f\right)=0,i=1,2,3,\dots ,14$(26) also holds true.

   Furthermore, discussions and characterization of optimal controls and adjoint variables are carried out in Appendix D.

5. Numerical simulations

In this section, using values of the parameters given in Table 3 and some initial population values we carried out, some numerical simulations regarding on the complete HIV/AIDS and pneumonia co-infection model (3) and the resulting optimal control based dynamical system consisting of the state equations (17) and the adjoint system (27). Applying ordinary differential equation solvers coded in MATLAB programming language, we have verified some basic qualitative results of the co-infection model (3) by performed numerical simulation of the models using some initial population and parameter values collected from the ministry of health of Ethiopia (EMOH) and others from the literature and some of them by reasonable assumptions given in Table 3.

Table 3. Parameter values used for numerical simulation.

Symbol	Baseline parameter value per year	Source
$\mathrm{\Delta }$	369 individuals	Calculated
$\mu$	0.02	[21]
${\alpha }_1$	0.31	Estimated from [6]
${\alpha }_2$	0.59	Estimated from [6]
$d_1$	0.23	Calculated
$d_2$	0.34	Calculated
$\theta$	0.27	Calculated
${\theta }_1$	0.29	Assumed
${\theta }_2$	0.30	Assumed
$\upsilon$	1.1	Assumed
${\varphi }_1$	1.21	Assumed
${\varphi }_2$	1.25	Assumed
$\gamma$	0.31	Calculated
$d_3$	0.33	[21]
$d_4$	0.42	[21]
$d_5$	0.51	Assumed
$\eta$	0.42	Estimated from [1]
$\delta$	0.53	Assumed
$ϵ$	0.043	Calculated
${\text{β}}_1$	0.38	Calculated
${\text{β}}_2$	2.35	Calculated
$k_1$	0.30	Calculated
$k_2$	0.30	Assumed
$k_3$	0.20	Estimated from [6]
$k_4$	0.20	Assumed
$\rho$	0.30	Estimated from [1]
$\kappa$	0.62	Estimated from [1]

5.1. Numerical simulations for the co-infection model (3)

In this section, using some initial population values, parameter values given in Table 3 and applying the forward and backward Range–Kutta scheme we investigate the direction of the co-infection model bifurcation diagram and the behaviours of the co-infection model solutions whenever the corresponding effective reproduction number is greater than unity. We implement the numerical simulations in the maximum time level to be 5 years.

5.1.1. Backward bifurcation diagram

The diagram representation of the HIV/AIDS and pneumonia co-infection model bifurcation diagram given in Figure 2 shows the phenomenon of the HIV/AIDS and pneumonia co-infection model backward bifurcation, which results in the co-existence of positive HIV/AIDS and pneumonia co-infection endemic and disease-free equilibrium points whenever ${\mathcal{R}}_{PM}<1$. In such a case, the common conditions of HIV/AIDS and pneumonia co-infection minimization such as making ${\mathcal{R}}_{PM}<1$ will not be necessarily effective, and the initial number of HIV/AIDS and pneumonia co-infected individuals also plays a crucial role.

Figure 2. Backward bifurcation diagram of the co-infection model.

5.1.2. Trajectories of the model (3) solutions whenever ${\mathcal{R}}_{\mathbf{\text{H}}P}>1$

In this section, we have carried out numerical simulation of the complete co-infection model (3) by applying some initial population values and parameter values given in Table 3 and the simulation illustrated by Figure 3 shows the HIV/AIDS and pneumonia co-infection dynamical system (3) solutions are approaching to the complete model endemic equilibrium point whenever ${\mathcal{R}}_{\text{HP}}=\text{max}\left\{{\mathcal{R}}_{HM},{\mathcal{R}}_{PM}\right\}=max\left\{1.23,2.91\right\}=2.91>1$. Biologically it means that the outbreak of HIV/AIDS and pneumonia co-infection throughout the community is consistently present in the study area.

Figure 3. HIV/AIDS and pneumonia co-infection model (3) solutions behaviour whenever ${\mathcal{R}}_{\text{HP}}$ = 2.91 > 1.

5.2. Optimal control problem numerical simulations

In this section, we investigate the effect of intervention strategies on the transmission of HIV/AIDs and pneumonia co-infection in a population of the study area. The optimal control problem stated in (17)–(21) is solved numerically using the forward and backward Range–Kutta scheme. To verify the qualitative analysis and to investigate the most effective controlling strategy to minimize the number of HIV/AIDs and pneumonia co-infection individuals in a population we implemented the numerical simulation of the co-infection model using some initial population values and parameter values illustrated in Table 3. We implement the numerical simulations in the maximum time level to be 5 years by considering the following possible illustrated optimal control strategies:

1. Single strategy at a time

   • A: use of HIV protection by condom only

   • B: use of pneumonia protection only

   • C: use of HIV/AIDS treatment only

   • D: use of pneumonia treatment only.

2. Double strategies simultaneously

   • E: use HIV protection + pneumonia protection only

   • F: use HIV/AIDS treatment + pneumonia treatment only

   • G: use HIV protection + HIV/AIDS treatment only

   • H: use pneumonia protection + pneumonia treatment only

   • I: use HIV protection + pneumonia treatment only

   • J: use pneumonia protection + HIV/AIDS treatment only.

3. Triple strategies simultaneously

   • K: use HIV protection + pneumonia protection + HIV/AIDS treatment only

   • L: use HIV protection + pneumonia protection + pneumonia treatment only

   • M: use HIV protection + pneumonia treatment + HIV/AIDS treatment only

   • N: use pneumonia protection + pneumonia treatment + HIV/AIDS treatment only.

4. All the possible mentioned strategies simultaneously

   • O: use all the four strategies simultaneously.

5.2.1. Impacts of single strategies on the co-infection

In this section, simulation is done when there is no control strategy in place and considering the following controlling strategies: Strategy A: use of HIV protection by condom, we present the simulation of optimal control system (17) with condom use $\left(u_1\right)$ as a protection against HIV infection by Figure 4(A). Strategy B: Use of pneumonia protection strategy, we present the simulation of optimal control system (17) with protection mechanism $\left(u_2\right)$ as a protection against pneumonia infection by Figure 4(B). Strategy C: Use of HIV/AIDS treatment strategy, we present the simulation of optimal control system (17) with HIV/AIDS treatment mechanism $\left(u_4\right)$ as a treatment against HIV infection by Figure 4(C). Strategy D: Use of pneumonia treatment strategy, we present the simulation of optimal control system (17) with pneumonia treatment mechanism $\left(u_3\right)$ as a treatment against pneumonia infection by Figure 4(D).

Figure 4. Impact of single strategies on the HIV/AIDS and pneumonia co-infection.

From Figure 4, we observed that the protective strategies illustrated in B and C are more effective strategies compared to the treatment strategies illustrated in A and D. But we recommend that strategy B is the most effective strategy to tackle the co-infection problem in the community.

5.2.2. Impacts of double strategies on the co-infection

In this section, simulation is done when there is no control strategy in place and considering the following controlling strategies: Strategy E: use HIV and pneumonia protection strategies simultaneously, we present the simulation of optimal control system (17) with HIV and pneumonia protection strategies $(u_1$ and $u_2)$ simultaneously as a protection against HIV and pneumonia infections respectively and is illustrated by Figure 5(E). Strategy F: Use HIV and pneumonia treatment strategies simultaneously, we present the simulation of optimal control system (17) with HIV and pneumonia treatment strategies $(u_3$ and $u_4)$ simultaneously as treatments against pneumonia and HIV infections respectively and is illustrated by Figure 5(F). Strategy G: Use HIV protection and HIV treatment strategies simultaneously, we present the simulation of optimal control system (17) with HIV protection and HIV treatment strategies $(u_1$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated by Figure 5(G). Strategy H: Use pneumonia protection and pneumonia treatment strategies simultaneously, we present the simulation of optimal control system (17) with pneumonia protection and pneumonia treatment strategies $(u_1$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated by Figure 5(H). Strategy I: Use HIV protection and pneumonia treatment strategies simultaneously, we present the simulation of optimal control system (17) with HIV protection and pneumonia treatment strategies $(u_1$ and $u_3)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated by Figure 5(I). Strategy J: Use pneumonia protection and HIV treatment strategies simultaneously, we present the simulation of optimal control system (17) with pneumonia protection and HIV treatment strategies $(u_2$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated by Figure 5(J).

Figure 5. Effects of double strategies on the HIV/AIDS and pneumonia co-infection.

From Figure 5, we observed that the protective and treatment strategies illustrated in G, H, I and J are more effective strategies compared to other strategies illustrated in E and F. But we recommend that strategy J is the most effective strategy to tackle the co-infection problem in the community.

5.2.3. Impacts of triple strategies on the co-infection

In this section, simulation is done when there is no control strategy in place and considering the following controlling strategies: Strategy K: use both HIV and pneumonia protections and HIV treatment strategies simultaneously, we present the simulation of optimal control system (17) with both HIV and pneumonia protections and HIV treatment strategies $(u_1$, $u_2$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated in Figure 6(K). Strategy L: Use both HIV and pneumonia protections and pneumonia treatment strategies simultaneously, we present the simulation of optimal control system (17) with both HIV and pneumonia protections and pneumonia treatment strategies $(u_1$, $u_2$ and $u_3)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated in Figure 6(L). Strategy M: Use both HIV and pneumonia treatments and HIV protection strategies simultaneously, we present the simulation of optimal control system (17) with both HIV and pneumonia treatments and HIV protection strategies $(u_1$, $u_3$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated in Figure 6(M). Strategy N: Use both HIV and pneumonia treatments and pneumonia protection strategies simultaneously, we present the simulation of optimal control system (17) with both HIV and pneumonia treatments and pneumonia protection strategies $(u_2$, $u_3$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated in Figure 6(N).

Figure 6. Effects of triple strategies on the HIV/AIDS and pneumonia co-infection.

From Figure 6, we observed that the protections and treatment strategies illustrated in K and L are more effective strategies compared to the strategies illustrated in M and N. But we recommend that strategy K is the most effective strategy to tackle the co-infection problem in the community.

5.2.4. Impacts of combination of all strategies on the co-infection

Strategy O: Use all the four strategies HIV protection, pneumonia protection, pneumonia treatment and HIV treatment strategies $(u_1$, $u_2,u_3$ and $u_4)$ simultaneously, we present the simulation of optimal control system (17) with all strategies $(u_1$, $u_2,u_3$ and $u_4)$ simultaneously as a control strategy against HIV and pneumonia co-infection and is illustrated by Figure 7.

Figure 7. Effects of triple strategies on the HIV/AIDS and pneumonia co-infection.

In this section, numerical simulation is carried out when there is no control strategy in place and when there are controls involving protection and treatment strategies for both pneumonia and HIV/AIDS single infections. Figure 7 shows the result that all the protection and treatment strategies efforts are implemented, and the number of individuals co-infected with HIV/AIDS and pneumonia decreases drastically to 0 after 3. Using the result given in Figure 7, we also compared strategy O to each of other strategies and found out that the strategy shows a significant decline in the number of HIV/AIDS and pneumonia co-infected individuals and hence strategy O is the most effective strategies to tackle the co-infection spreading in the community.

6. Discussion and conclusion

In this proposed study, a new 11 mathematical model for the transmission dynamics of pneumonia and HIV/AIDS co-infection with pneumonia protection, pneumonia vaccination, HIV/AIDS protection by using condom, pneumonia treatment and HIV/AIDS treatment measures in a considered community is constructed and assessed the effects of the prevention and control measures such as pneumonia protection, pneumonia vaccination, HIV/AIDS protection, pneumonia treatment, and HIV/AIDS treatment for the pneumonia and HIV/AIDS co-infected individuals in a region where the model is both mathematically and biologically well-posed. Both analytical and numerical analyses of the pneumonia and HIV/AIDS co-infection model (3) were performed to assess the effect of various prevention and controlling strategies for the dynamics of both single infections and co-infection in a considered community.

In this study, the theoretical results are verified by the numerical results and we can summarize results as: the disease-free equilibrium points of all the mono-infection models (7) and (11) and the complete co-infection model (3) are locally asymptotically stable whenever each corresponding effective reproduction numbers are less than unity respectively also the study discusses the global stability of disease-free equilibrium points of the two-sub-models. Using Centre manifold criteria, we have shown the backward bifurcation of the pneumonia infection model and also the HIV/AIDS and pneumonia co-infection model which means there is co-existence of the disease-free equilibrium point(s) together with positive endemic equilibrium point(s) in the region where the co-infection model effective reproduction number is less than unity.

From the numerical simulations, we have observed the following results: from Figure 2 of the numerical simulation, one can observe that the HIV/AIDS and pneumonia co-infection model exhibit the phenomenon of backward bifurcation means whenever the effective reproduction number is less than unity both positive endemic equilibrium and disease-free equilibrium points exist. Figure 3 shows that the solutions of the co-infection model (3) converging to its endemic equilibrium point whenever its effective reproduction number ${\mathcal{R}}_{HP}=max\left\{{\mathcal{R}}_{PM},{\mathcal{R}}_{HM}\right\}=${2.91, 1.23$\}=2.91>1.$

Figure 4 shows the impact of every single strategy separately to tackle the HIV/AIDS and pneumonia co-infection spreading dynamics. Figure 5 shows the impact of every double strategy separately to tackle the HIV/AIDS and pneumonia co-infection spreading dynamics. Figure 6 shows the impact of every triple strategy separately to tackle the HIV/AIDS and pneumonia co-infection spreading dynamics. Figure 7 shows the impact of the combination of all strategies simultaneously to tackle the HIV/AIDS and pneumonia co-infection spreading dynamics. Compared to the other control strategies investigated and discussed applying all the protections and treatments controlling strategies simultaneously as stated in the discussion of Figure 7 is the most effective strategy to decrease the HIV/AIDS and pneumonia co-infection transmission mechanisms in the community under the study.

The limitations of this study and the next potential researchers can incorporate and extend the model under this study will be: stochastic approach, fractional order derivative approach, environmental impacts, age and spatial structure, validate the model by applying appropriate real infection data.

Authors’ contributions

All authors have read and approved the final manuscript.

Math Classification

92-10 Mathematical modelling or simulation for problems pertaining to biology.

Acknowledgments

The authors thank Mr. Sitotaw Eshete for his valuable Wi-Fi contribution.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

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

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.

Appendices

Appendix A. Proof of Theorem 2.1

Proof:

Let us consider $D_1\left(0\right)>0$, $D_2\left(0\right)>0,D_3\left(0\right)>0$, $D_4\left(0\right)>0$, $D_5\left(0\right)>0$, $D_6\left(0\right)>0$, $D_7\left(0\right)>0$, $D_8\left(0\right)>0,D_9\left(0\right)>0$, $D_{10}\left(0\right)>0,$ and $D_{11}\left(0\right)>0$ then for all t > 0, we have to show that $D_1$ (t) > 0, $D_2\left(t\right)>0,D_3\left(t\right)>0$, $D_4\left(t\right)$> 0, $D_5\left(t\right)$>0, $D_6\left(t\right)$>0, $D_7\left(t\right)$>0, $D_8\left(t\right)$>0, $D_9\left(t\right)$>0, $D_{10}\left(t\right)>0,$ and $D_{11}\left(t\right)$ > 0.

Define: $\tau$= sup $\{t>0:D_1(\text{t})>0,D_2(t)>0,D_3(t)>0$, $D_4\left(t\right)>0,D_5\left(t\right)>0$, $D_6\left(t\right)$>0, $D_7\left(t\right)$>0, $D_8\left(t\right)$> 0, $D_9\left(t\right)$> 0$,D\left(t\right)>0,\text{and }{D}_{11}\left(t\right)>0\}$. Since all the state variables $D_1\left(t\right),D_2\left(t\right)$, $D_3\left(t\right)$, $D_4\left(t\right),D_5\left(t\right),D_6\left(t\right)$, $D_7\left(t\right),D_8\left(t\right)$, $D_9\left(t\right),D_{10}\left(t\right)$, and $D_{11}\left(t\right)$ are continuous and we can justify that $\tau >0$. If $\tau$ = +∞, then non-negativity holds. But, if 0 <$\tau$ < +∞, we will have $D_1\left(\tau \right)=0$ or $D_2\left(\tau \right)=0$ or $D_3\left(\tau \right)=0$ or $D_4\left(\tau \right)=0$ or $D_5\left(\tau \right)=0$ or $D_6\left(\tau \right)=$ or $D_7\left(\tau \right)=0$ or $D_8\left(\tau \right)=0$ or $D_9\left(\tau \right)=0$ or $D_{10}\left(\tau \right)=0$ or $D_{11}\left(\tau \right)=0$.

From the first equation of the model (3), we do have $\dot{D_1}+\left({\lambda }_H+{\lambda }_C+\mu \right)D_1=k_1\mathrm{\Delta }+{\alpha }_1{D}_2+{\alpha }_2{D}_3+\rho D_4+\eta D_{10}.$And integrate both sides using integrating factor, we have determined the constant value $D_1\left(\tau \right)=M_1{D}_1\left(0\right)+M_1{\int }_0^{\tau }{\exp }^{\int \left(\mu +{\lambda }_H\left(t\right)+{\lambda }_P\left(t\right)\right)dt}\left(k_1\mathrm{\Delta }+{\alpha }_1{D}_2+{\alpha }_2{D}_3+\rho D_4+\eta D_{10}\right)dt>0$ where $M_1=ex{p}^{-\left(\mu \tau +\underset{0}{\overset{\tau }{\int }}\left({\lambda }_H\right(w)+{\lambda }_P(w)\right)}>0,D_1\left(0\right)>0,$ and from the meaning of $\tau$_, the solutions $D_2\left(t\right)>0,$ $D_3\left(t\right)>0,D_4\left(t\right)>0$, $D_{10}\left(t\right)>0$, also the exponential function always is positive, then the solution $A_1\left(\tau \right)>0$ hence $A_1\left(\tau \right)\ne 0$.

Again from the second equation of the model (3), we do have ${\dot{D}}_2+\left({\lambda }_H+{\alpha }_1+\mu \right)D_2=k_2\mathrm{\Delta }$.

And also using t integrating factor after some calculations we obtained that

$D_2\left(\tau \right)=M_1{D}_2\left(0\right)+M_1{\int }_0^{\tau }{\exp }^{\int \left({\lambda }_H+{\alpha }_1+\mu )\right)dt}{k}_2\mathrm{\Delta }dt>0$ where $M_1=\text{ex}{\text{p}}^{-\left({\alpha }_1\tau +\mu \tau +\underset{0}{\overset{\tau }{\int }}\left({\lambda }_H\right(w)\right)}>0,D_2\left(0\right)>0$, and from the meaning of $\tau$, the solution $D_2\left(\tau \right)>0$ hence $D_2\left(\tau \right)\ne 0$.

Similarly, $D_3\left(\tau \right)>0$ hence $D_3\left(\tau \right)\ne 0$, $D_4\left(\tau \right)>0$ hence $D_4\left(\tau \right)\ne 0$, $D_5\left(\tau \right)>0$ hence $D_5\left(\tau \right)\ne 0$, $D_6\left(\tau \right)>0$ hence $D_6\left(\tau \right)\ne 0$, $D_7\left(\tau \right)>0$ hence $D_7\left(\tau \right)\ne 0$, $D_8\left(\tau \right)>0$ hence $D_8\left(\tau \right)\ne 0$, $D_9\left(\tau \right)>0$ hence $D_9\left(\tau \right)\ne 0$, $D_{10}\left(\tau \right)>0$ hence $D_{10}\left(\tau \right)\ne 0$, and $D_{11}\left(\tau \right)>0$ hence $D_{11}\left(\tau \right)\ne 0$.

Thus $\tau =+\mathrm{\infty }$, and hence all the solutions of the pneumonia and HIV/AIDS co-infection model (3) are non-negative.

Appendix B. Proof of theorem 2.2

Proof: Let $\left(D_1,D_2,D_3,D_4,D_5,D_6,D_7,D_8,D_9,D_{10},D_{11}\right)\in {\mathbb{R}}_{+}^{11}$ be an arbitrary non-negative solution of the system (3) with initial conditions given in Equation (4).

Now adding all the differential equations given in Equation (3), we do have the derivative of the total population $N$ which is given in Equation (5) as $\dot{N}=\mathrm{\Delta }-\mu N-(d_1{A}_5+d_2{A}_6+d_3{A}_7+d_4{A}_8+d_5{A}_9).$ Then by ignoring all the infection classes, we have determined that $\dot{N}\le \mathrm{\Delta }-\mu N,$ and using the separation of variables whenever $t\to \mathrm{\infty },$ we have obtained that $0\le N\le \frac{\mathrm{\Delta }}{\mu }$. Hence, all the positive feasible non-negative solutions of the co-infection model (3) entering into the region are given in Equation (6).

Appendix C. Proof of theorem 3.6

In this section, we apply the centre manifold theory [19] to ascertain the local stability of the endemic equilibrium due to the convolution of the first approach (eigenvalues of the Jacobian matrix). To make use of the centre manifold theory, the following change of variables is made by symbolizing $D_1=x_1$, $D_2=x_2,D_4=x_3$, $D_5=x_4$, and $D_{10}=x_5$ such that $N_2=x_1+x_2+x_3+x_4+x_5$. Furthermore, by using vector notation $X=(x_1,x_2,x_3,x_4,x_5{)}^T$, the COVID-19 mono-infection model (11) can be written in the form $\frac{dX}{dt}=F\left(X\right)$ with $F=(f_1,f_2,f_3,f_4,f_5{)}^T$, as follows: (14) $\begin{array}{rl}\frac{d{x}_1}{dt}& =f_1=k_1\mathrm{\Delta }+{\alpha }_1{x}_2+\rho x_3+\eta x_5-\mu x_1-{\lambda }_P{x}_1,\\ \frac{d{x}_2}{dt}& =f_2=k_2\mathrm{\Delta }-\left({\alpha }_1+\mu \right)x_2,\\ \frac{d{x}_3}{dt}& =f_4=k_4\mathrm{\Delta }-\left(\rho +\mu +ϵ{\lambda }_P\right)x_3,\\ \frac{d{x}_4}{dt}& ={\lambda }_P{x}_1+ϵ{\lambda }_P{x}_3-\left(\mu +d_1+\kappa \right)x_4,\\ \frac{d{x}_5}{dt}& =\kappa x_4-\left(\mu +\eta \right)x_5,\end{array}$(14) with ${\lambda }_P={\text{β}}_2{x}_5$ then the method entails evaluating the Jacobian of the system (14) at the DFE point $E_{PM}^0$, denoted by $J\left(E_{PM}^0\right)$ and this gives us $J\left(E_{PM}^0\right)=\left(\begin{array}{ccccc}-\mu & {\alpha }_2& \rho & -{\text{β}}_2{x}_1^0& \eta \\ 0& -\left({\alpha }_1+\mu \right)& & 0& 0\\ 0& 0& -\left(\rho +\mu \right)& -{\text{β}}_2ϵx_3^0& 0\\ 0& 0& 0& {\text{β}}_2{x}_1^0+{\text{β}}_2ϵx_3^0-\left(\mu +d_1+\kappa \right)& 0\\ 0& 0& 0& \kappa & -\left(\mu +\eta \right)\end{array}\right).$Consider, ${\mathcal{R}}_{PM}=1$ and suppose that ${\text{β}}_2$ $={\text{β}}^{\ast }$ is chosen as a bifurcation parameter. From ${\mathcal{R}}_{PM}=1$ as ${\mathcal{R}}_{PM}=\frac{{\text{β}}_2{x}_2^0+ϵ{\text{β}}_2{x}_4^0}{\mu +d_1+\kappa }=\frac{{\text{β}}_2{k}_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\text{β}}_2{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+{\text{β}}_2{k}_4\mathrm{\Delta }\left({\alpha }_1+\mu \right)\left(\rho +\mu ϵ\right)}{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right)}=1$.

Solving for β₂ we have got ${\text{β}}_2={\text{β}}^{\ast }=\frac{\mu \left({\alpha }_1+\mu \right)\left(\rho +\mu \right)\left(\mu +d_1+\kappa \right))}{k_1\mathrm{\Delta }({\alpha }_1+\mu )\left(\rho +\mu \right)+{\alpha }_1{k}_2\mathrm{\Delta }\left(\rho +\mu \right)+k_4\mathrm{\Delta }\left({\alpha }_1+\mu \right)\left(\rho +\mu ϵ\right)}$. $J_{{\text{β}}^{\ast }}=\left(\begin{array}{ccccc}-\mu & {\alpha }_2& \rho & -{\text{β}}^{\ast }{x}_1^0& \eta \\ 0& -\left({\alpha }_1+\mu \right)& & 0& 0\\ 0& 0& -\left(\rho +\mu \right)& -{\text{β}}^{\ast }ϵx_3^0& 0\\ 0& 0& 0& {\text{β}}^{\ast }{x}_1^0+{\text{β}}^{\ast }ϵx_3^0-\left(\mu +d_1+\kappa \right)& 0\\ 0& 0& 0& \kappa & -\left(\mu +\eta \right)\end{array}\right).$After some steps of the calculation, we have got the eigenvalues of $J_{{\text{β}}^{\ast }}$ as ${\lambda }_1=-\mu$, ${\lambda }_2=-\left({\alpha }_1+\mu \right)$ or ${\lambda }_3=-\left(\rho +\mu \right)$ or ${\lambda }_4=0$ or ${\lambda }_5=-\left(\mu +\eta \right)$. It follows that the Jacobian $J\left(E_{PM}^0\right)$ of Equation (14) at the disease-free equilibrium with ${\text{β}}_2={\text{β}}^{\ast }$, denoted by $J_{{\text{β}}^{\ast }}$, has a simple zero eigenvalue with all the remaining eigenvalues have negative real part. Hence, Theorem 2.2 of Castillo-Chavez and Song [19] can be used to analyse the dynamics of the model to show that the model (11) undergoes backward bifurcation at ${\mathcal{R}}_{PM}=1$.

Eigenvectors of $J_{{\text{β}}^{\ast }}$: For the case ${\mathcal{R}}_{PM}=1$, it can be shown that the Jacobian of the system (14) at ${\text{β}}_2={\text{β}}^{\ast }$ (denoted by $J_{{\text{β}}^{\ast }})$ has a right eigenvectors associated with the zero eigenvalue given by $u=(u_1,u_2,u_3,u_4,u_5{)}^T$ as (15) $\begin{array}{rl}& \left(\begin{array}{ccccc}-\mu & {\alpha }_2& \rho & -{\text{β}}^{\ast }{x}_1^0& \eta \\ 0& -\left({\alpha }_1+\mu \right)& & 0& 0\\ 0& 0& -\left(\rho +\mu \right)& -{\text{β}}^{\ast }ϵx_3^0& 0\\ 0& 0& 0& {\text{β}}^{\ast }{x}_1^0+{\text{β}}^{\ast }ϵx_3^0-\left(\mu +d_1+\kappa \right)& 0\\ 0& 0& 0& \kappa & -\left(\mu +\eta \right)\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\end{array}\right)\\ & =\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right).\end{array}$(15) Then solving equation (15) the right eigenvectors associated with the zero eigenvalue are given by $\begin{array}{rl}{u}_1& =\frac{-{\text{β}}^{\ast }{x}_1^0{u}_4\left(\rho +\mu \right)\left(\mu +\eta \right)-{\text{β}}^{\ast }ϵx_3^0\rho \left(\mu +\eta \right)u_4+\kappa \eta \left(\rho +\mu \right)u_4}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)},\\ u_2& =0,u_3=-\frac{{\text{β}}^{\ast }ϵx_3^0}{\left(\rho +\mu \right)}{u}_4,u_4=u_4>0,u_5=\frac{\kappa }{\mu +\eta }{u}_4.\end{array}$Similarly, the left eigenvector associated with the zero eigenvalues at ${\text{β}}_2={\text{β}}^{\ast }$ given by $v=(v_1,v_2,v_3,v_4,v_5{)}^T$ as (16) ${\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\\ v_4\\ v_5\end{array}\right)}^T\ast \left(\begin{array}{ccccc}-\mu & {\alpha }_2& \rho & -{\text{β}}^{\ast }{x}_1^0& \eta \\ 0& -\left({\alpha }_1+\mu \right)& & 0& 0\\ 0& 0& -\left(\rho +\mu \right)& -{\text{β}}^{\ast }ϵx_3^0& 0\\ 0& 0& 0& D& 0\\ 0& 0& 0& \kappa & -\left(\mu +\eta \right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),$(16) where $D={\text{β}}^{\ast }{x}_1^0+{\text{β}}^{\ast }ϵx_3^0-\left(\mu +d_1+\kappa \right)$.

Then solving Equation (16), the left eigenvectors associated with the zero eigenvalue are given by $v_1=v_2=v_3=v_4=0$ and $v_4=v_4>0$. After long steps of calculations, the bifurcation coefficients $a$ and $b$ are obtained as $\begin{array}{rl}a& =\sum _{i,j,k=1}^5{v}_4{u}_i{u}_j{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j=2v_4{u}_1{u}_4{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_1\mathrm{\partial }{x}_4+2v_4{u}_3{u}_4{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_3\mathrm{\partial }{x}_4\\ & =2v_4{u}_4\left[u_1{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_1\mathrm{\partial }{x}_4+u_3{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_3\mathrm{\partial }{x}_4\right]\\ & =2v_4{u}_4\left[{\text{β}}_2{u}_1+{\text{β}}_2ϵu_3\right]\\ =& 2v_4{u_4}^2\left[\frac{\begin{array}{c}-{\text{β}}_2{\text{β}}^{\ast }{x}_1^0\left(\rho +\mu \right)\left(\mu +\eta \right)-{\text{β}}_2{\text{β}}^{\ast }ϵx_3^0\rho \left(\mu +\eta \right)+{\text{β}}_2\kappa \eta \left(\rho +\mu \right)\\ -{\text{β}}_2ϵ{\text{β}}^{\ast }ϵx_3^0\mu \left(\mu +\eta \right)\end{array}}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)}\right],\\ & =2v_4{u}_4\left[D_2-D_1\right],\end{array}$where $H_1$ = $\frac{-{\text{β}}_2{\text{β}}^{\ast }{x}_1^0\left(\rho +\mu \right)\left(\mu +\eta \right)-{\text{β}}_2{\text{β}}^{\ast }ϵx_3^0\rho \left(\mu +\eta \right)-{\text{β}}_2ϵ{\text{β}}^{\ast }ϵx_3^0\mu \left(\mu +\eta \right)}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)}$ and $H_2=\frac{{\text{β}}_2\kappa \eta \left(\rho +\mu \right)}{\mu \left(\rho +\mu \right)\left(\mu +\eta \right)}$.

Thus the bifurcation coefficient $a$ is positive whenever $D_2>D_1$.

Moreover $\begin{array}{rl}b& =\sum _{i,k=1}^5{v}_k{u}_i{\mathrm{\partial }}^2{f}_k/\mathrm{\partial }{x}_i\mathrm{\partial }\text{β}\left(E_{CM}^0\right)\\ & =\sum _{i=1}^5{v}_4{u}_i{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_i\mathrm{\partial }\text{β}=v_4{u}_4{\mathrm{\partial }}^2{f}_4/\mathrm{\partial }{x}_4\mathrm{\partial }\text{β}=v_4{u}_4\left[x_1^0{u}_1+ϵx_3^0{u}_3\right]>0.\end{array}$Hence, from in Castillo-Chavez and Song [19] the pneumonia mono-infection model (11) exhibits a backward bifurcation at ${\mathcal{R}}_{PM}=1$ and whenever $H_2>H_1$.

Appendix D. Discussion and characterization of optimal controls and adjoint variables

Theorem 3.8:

Let ${\mathfrak{u}}^{\ast }=\left({\mathfrak{u}}_1^{\ast },{\mathfrak{u}}_2^{\ast },{\mathfrak{u}}_3^{\ast },{\mathfrak{u}}_4^{\ast }\right)$ be the optimal control and ${\dot{D}}_5^{\ast }=\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{\text{D}}_1^{\ast }+\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{\text{D}}_3^{\ast }+ϵ\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{\text{D}}_3^{\ast }-\upsilon \left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{\text{D}}_5^{\ast }-\left(\mu +d_1+{\mathfrak{u}}_3^{\ast }\kappa \right){\text{D}}_5^{\ast }$ be the associated unique optimal solutions of the optimal control problem (17) with initial condition (18) and objective functional (19) with fixed final time $T_f$ (20). Then there exists adjoint function ${\lambda }_i^{\ast }(\cdot ),i=1,\dots ,11$ satisfying the following canonical equations: $\begin{array}{rl}& \frac{d{\lambda }_1}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_1}=\left(1-{\mathfrak{u}}_1\right){\lambda }_H^{\ast }\left({\lambda }_1-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right){\lambda }_P^{\ast }\left({\lambda }_1-{\lambda }_5\right)+\mu {\lambda }_1,\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_2}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_2}=\left(1-{\mathfrak{u}}_1\right){\lambda }_H^{\ast }\left({\lambda }_2-{\lambda }_6\right)+{\alpha }_1\left({\lambda }_2-{\lambda }_1\right)+\mu {\lambda }_2,\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_3}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_3}=\left(1-{\mathfrak{u}}_2\right){\lambda }_P^{\ast }\left({\lambda }_3-{\lambda }_5\right)+{\alpha }_2\left({\lambda }_3-{\lambda }_1\right)+\mu {\lambda }_3,\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_4}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_4}=\left(1-{\mathfrak{u}}_1\right){\lambda }_H^{\ast }\left({\lambda }_4-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right)ϵ{\lambda }_P^{\ast }\left({\lambda }_4-{\lambda }_5\right)+\rho \left({\lambda }_4-{\lambda }_1\right)+\mu {\lambda }_4,\end{array}$ (27) $\begin{array}{rl}& \frac{d{\lambda }_6}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_6}=\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_1^{\ast }\left({\lambda }_1-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_2^{\ast }\left({\lambda }_2-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_4^{\ast }\left({\lambda }_4-{\lambda }_6\right)\\ & +\left({\lambda }_4-{\lambda }_5\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\text{β}}_2{D}_6^{\ast }\left({\lambda }_6-{\lambda }_8\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\text{β}}_2{D}_7^{\ast }\left({\lambda }_7-{\lambda }_9\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\upsilon {\lambda }_H\left({\lambda }_5-{\lambda }_8\right)+\left(\mu +d_1\right){\lambda }_5+{\mathfrak{u}}_3\kappa \left({\lambda }_5-{\lambda }_{10}\right),\end{array}$(27) $\begin{array}{rl}& \frac{d{\lambda }_6}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_6}=\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_1^{\ast }\left({\lambda }_1-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_2^{\ast }\left({\lambda }_2-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1}{N}{D}_4^{\ast }\left({\lambda }_4-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\upsilon \frac{{\text{β}}_1}{N}{D}_5^{\ast }\left({\lambda }_5-{\lambda }_8\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\lambda }_P^{\ast }\left({\lambda }_6-{\lambda }_8\right)+\left(\mu +d_2\right){\lambda }_6+\mathrm{\Theta }\left({\lambda }_6-{\lambda }_7\right)\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_7}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_7}=-{\mathfrak{w}}_2+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_1}{N}{D}_1^{\ast }\left({\lambda }_1-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_1}{N}{D}_2^{\ast }\left({\lambda }_2-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_1}{N}{D}_4^{\ast }\left({\lambda }_4-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\upsilon \frac{{\text{β}}_1{\rho }_1}{N}{D}_5^{\ast }\left({\lambda }_5-{\lambda }_8\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\lambda }_P^{\ast }\left({\lambda }_7-{\lambda }_9\right)+\left(d_3+\mu \right){\lambda }_7+{\mathfrak{u}}_4\gamma \left({\lambda }_7-{\lambda }_{11}\right),\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_8}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_8}\\ & =-{\mathfrak{w}}_3+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_2}{N}{D}_1^{\ast }\left({\lambda }_1-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right){\text{β}}_2{\omega }_1{D}_1^{\ast }\left({\lambda }_1-{\lambda }_5\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_2}{N}{D}_2^{\ast }\left({\lambda }_2-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_2\right){\text{β}}_2{\omega }_1{D}_3^{\ast }\left({\lambda }_3-{\lambda }_5\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_2}{N}{D}_4^{\ast }\left({\lambda }_4-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right)ϵ{\text{β}}_2{\omega }_1{D}_4^{\ast }\left({\lambda }_4-{\lambda }_5\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\upsilon \frac{{\text{β}}_1{\rho }_2}{N}{D}_5^{\ast }\left({\lambda }_5-{\lambda }_8\right)+{\mathfrak{u}}_3{\theta }_1\left({\lambda }_8-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\text{β}}_2{\omega }_1{D}_6^{\ast }\left({\lambda }_6-{\lambda }_8\right)\\ & +\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\text{β}}_2{\omega }_1{D}_7^{\ast }\left({\lambda }_7-{\lambda }_9\right)+\left(\mu +d_4\right){\lambda }_8+\delta \left({\lambda }_8-{\lambda }_9\right),\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_9}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_9}=-{\mathfrak{w}}_4+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_3}{N}{D}_1^{\ast }\left({\lambda }_1-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_2\right){\text{β}}_2{\omega }_2{D}_1^{\ast }\left({\lambda }_1-{\lambda }_5\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_3}{N}{D}_2^{\ast }\left({\lambda }_2-{\lambda }_6\right)\\ & +\left(1-{\mathfrak{u}}_2\right){\text{β}}_2{\omega }_2{D}_3^{\ast }\left({\lambda }_3-{\lambda }_5\right)+\left(1-{\mathfrak{u}}_1\right)\frac{{\text{β}}_1{\rho }_3}{N}{D}_4^{\ast }\left({\lambda }_4-{\lambda }_6\right)+\left(1-{\mathfrak{u}}_2\right)ϵ{\text{β}}_2{\omega }_2{D}_4^{\ast }\left({\lambda }_4-{\lambda }_5\right)\\ & +\left(1-{\mathfrak{u}}_1\right)\upsilon \frac{{\text{β}}_1{\rho }_3}{N}{D}_5^{\ast }\left({\lambda }_5-{\lambda }_8\right)+\left(1-{\mathfrak{u}}_2\right){\varphi }_1{\text{β}}_2{\omega }_2{D}_6^{\ast }\left({\lambda }_6-{\lambda }_8\right)+{\mathfrak{u}}_3{\theta }_2\left({\lambda }_9-{\lambda }_7\right)\\ & +\left(1-{\mathfrak{u}}_2\right){\varphi }_2{\text{β}}_2{\omega }_2{D}_7^{\ast }\left({\lambda }_7-{\lambda }_9\right)+\left(\mu +d_5\right){\lambda }_9,\end{array}$ $\begin{array}{rl}& \frac{d{\lambda }_{10}}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{A}_{10}}=-\eta {\lambda }_1+\left(\mu +\eta \right){\lambda }_{10},\end{array}$

$\frac{d{\lambda }_{11}}{dt}=-\frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }{D}_{11}}=\mu {\lambda }_{11}$ with transversality conditions (28) ${\lambda }_i^{\ast }\left(T_f\right)=0,i=1,2,\dots ,11.$(28) Moreover, the corresponding optimal controls ${\mathfrak{u}}_1^{\ast }\left(t\right),{\mathfrak{u}}_2^{\ast }\left(t\right),{\mathfrak{u}}_3^{\ast }\left(t\right),$ and ${\mathfrak{u}}_4^{\ast }\left(t\right)$ are given by (29) $\begin{array}{rl}{\mathfrak{u}}_1^{\ast }\left(t\right)& =max\left\{0,min\left\{\frac{\begin{array}{c}{\lambda }_H^{\ast }{D}_1^{\ast }\left({\lambda }_6-{\lambda }_1\right)+{\lambda }_H^{\ast }{D}_2^{\ast }\left({\lambda }_6-{\lambda }_2\right)+{\lambda }_H^{\ast }{D}_4^{\ast }\left({\lambda }_6-{\lambda }_4\right)\\ +\upsilon {\lambda }_H^{\ast }{D}_5^{\ast }\left({\lambda }_8-{\lambda }_5\right)\end{array}}{{\mathfrak{B}}_1},1\right\}\right\},\\ {\mathfrak{u}}_2^{\ast }\left(t\right)& =max\left\{0,min\left\{\frac{\begin{array}{c}{\lambda }_P^{\ast }{D}_1^{\ast }\left({\lambda }_5-{\lambda }_1\right)+{\lambda }_P^{\ast }{D}_3^{\ast }\left({\lambda }_5-{\lambda }_3\right)+ϵ{\lambda }_P^{\ast }{D}_4^{\ast }\left({\lambda }_5-{\lambda }_4\right)\\ +{\varphi }_1{\lambda }_P^{\ast }{D}_6^{\ast }\left({\lambda }_8-{\lambda }_6\right)+{\varphi }_2{\lambda }_P^{\ast }{D}_7^{\ast }\left({\lambda }_9-{\lambda }_7\right)\end{array}}{{\mathfrak{B}}_2},1\right\}\right\},\\ {\mathfrak{u}}_3^{\ast }\left(t\right)& =max\left\{0,min\left\{\frac{{\mathrm{\Theta }}_1{D}_8^{\ast }\left({\lambda }_8-{\lambda }_6\right)+{\mathrm{\Theta }}_2{D}_9^{\ast }\left({\lambda }_9-{\lambda }_7\right)+\kappa D_5^{\ast }\left({\lambda }_5-{\lambda }_{10}\right)}{{\mathfrak{B}}_3},1\right\}\right\},\\ {\mathfrak{u}}_4^{\ast }\left(t\right)& =max\left\{0,min\left\{\frac{\gamma D_7^{\ast }\left({\lambda }_7-{\lambda }_{11}\right)}{{\mathfrak{B}}_4},1\right\}\right\}.\end{array}$(29) From the previous analysis, to get the optimal point, we have to solve the system $\begin{array}{rl}{\dot{D}}_1^{\ast }& =k_1\mathrm{\Delta }+{\alpha }_1{D}_2^{\ast }+{\alpha }_2{D}_3^{\ast }+\rho D_4^{\ast }+\eta D_{10}^{\ast }-\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_1^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{D}_1^{\ast }-\mu D_1^{\ast },\\ {\dot{D}}_2^{\ast }& =k_2\mathrm{\Delta }-\left(\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }+{\alpha }_1+\mu \right)D_2^{\ast },\\ {\dot{D}}_3^{\ast }& =k_3\mathrm{\Delta }-\left({\alpha }_2+\mu +\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }\right)D_3^{\ast },\\ {\dot{D}}_4^{\ast }& =k_4\mathrm{\Delta }-\left(\rho +\mu +\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }+ϵ\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }\right)D_4^{\ast },\\ {\dot{D}}_5^{\ast }& =\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{D}_1^{\ast }+\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{D}_3^{\ast }+ϵ\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{D}_3^{\ast }\\ & -\upsilon \left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_5^{\ast }-\left(\mu +d_1+{\mathfrak{u}}_3^{\ast }\kappa \right)D_5^{\ast },\\ {\dot{D}}_6^{\ast }& =\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_1^{\ast }+\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_2^{\ast }+\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_4^{\ast }\\ & +{\mathfrak{u}}_3^{\ast }{\theta }_1{D}_8^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_1{\lambda }_P^{\ast }{D}_6^{\ast }-\left(\theta +\mu +d_2\right)D_6^{\ast },\\ {\dot{D}}_7^{\ast }& =\mathrm{\Theta }{D}_6^{\ast }+{\mathfrak{u}}_3^{\ast }{\mathrm{\Theta }}_2{D}_9^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_2{\lambda }_P^{\ast }{D}_7^{\ast }-\left({\mathfrak{u}}_4^{\ast }\gamma +d_3+\mu \right)D_7^{\ast },\\ {\dot{D}}_8^{\ast }& =\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_1{\lambda }_P^{\ast }{D}_6^{\ast }+\left(1-{\mathfrak{u}}_1^{\ast }\right)\upsilon {\lambda }_H^{\ast }{D}_5^{\ast }-\left(\mu +d_4+\delta +{\mathfrak{u}}_3^{\ast }{\mathrm{\Theta }}_1\right)D_8^{\ast },\\ {\dot{D}}_9^{\ast }& =\delta D_8^{\ast }+\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_2{\lambda }_P^{\ast }{D}_7^{\ast }-\left(\mu +d_5+{\mathfrak{u}}_3^{\ast }{\theta }_2\right)D_9^{\ast },\\ {\dot{D}}_{10}^{\ast }& ={\mathfrak{u}}_3^{\ast }\kappa D_5^{\ast }-\left(\mu +\eta \right)D_{10}^{\ast },\\ {\dot{D}}_{11}^{\ast }& ={\mathfrak{u}}_4^{\ast }\gamma D_7^{\ast }-\mu D_{11}^{\ast }.\end{array}$With the Hamiltonian $\begin{array}{rl}\mathcal{H}& ={\mathfrak{w}}_1{D}_5^{\ast }+{\mathfrak{w}}_2{D}_7^{\ast }+{\mathfrak{w}}_3{D}_8^{\ast }+{\mathfrak{w}}_4{D}_9^{\ast }+\frac{{\mathfrak{B}}_1}{2}({\mathfrak{u}}_1^{\ast }{)}^2+\frac{{\mathfrak{B}}_2}{2}({\mathfrak{u}}_2^{\ast }{)}^2+\frac{{\mathfrak{B}}_3}{2}({\mathfrak{u}}_3^{\ast }{)}^2+\frac{{\mathfrak{B}}_4}{2}({\mathfrak{u}}_4^{\ast }{)}^2\\ & +{\lambda }_1\left(k_1\mathrm{\Delta }+{\alpha }_1{D}_2^{\ast }+{\alpha }_2{D}_3^{\ast }+\rho D_4^{\ast }+\eta D_{10}^{\ast }-\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }{D}_1^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }{D}_1^{\ast }-\mu D_1^{\ast }\right)\\ & +{\lambda }_2\left(k_2\mathrm{\Delta }-\left(\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }+{\alpha }_1+\mu \right)D_2^{\ast }\right)\\ & +{\lambda }_3\left(k_3\mathrm{\Delta }-\left({\alpha }_2+\mu +\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_C^{\ast }\right)D_3^{\ast }\right)\\ & +{\lambda }_4\left(k_4\mathrm{\Delta }-\left(\rho +\mu +\left(1-{\mathfrak{u}}_1^{\ast }\right){\lambda }_H^{\ast }+ϵ\left(1-{\mathfrak{u}}_2^{\ast }\right){\lambda }_P^{\ast }\right)D_4^{\ast }\right)\\ & +{\lambda }_5\left(\right(1-{\mathfrak{u}}_2^{\ast }){\lambda }_P^{\ast }{D}_1^{\ast }+(1-{\mathfrak{u}}_2^{\ast }){\lambda }_P^{\ast }{D}_3^{\ast }+ϵ(1-{\mathfrak{u}}_2^{\ast }){\lambda }_P^{\ast }{D}_3^{\ast }-\upsilon (1-{\mathfrak{u}}_1^{\ast }){\lambda }_H^{\ast }{D}_5^{\ast }\\ & -\left(\mu +d_1+{\mathfrak{u}}_3^{\ast }\kappa \right)D_5^{\ast })\\ & +{\lambda }_6\left(\right(1-{\mathfrak{u}}_1^{\ast }){\lambda }_H^{\ast }{D}_1^{\ast }+(1-{\mathfrak{u}}_1^{\ast }){\lambda }_H^{\ast }{D}_2^{\ast }+(1-{\mathfrak{u}}_1^{\ast }){\lambda }_H^{\ast }{D}_4^{\ast }\\ & +{\mathfrak{u}}_3^{\ast }{\theta }_1{D}_8^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_1{\lambda }_P^{\ast }{D}_6^{\ast }-\left(\theta +\mu +d_2\right)D_6^{\ast })\\ & +{\lambda }_7\left(\theta D_6^{\ast }+{\mathfrak{u}}_3^{\ast }{\theta }_2{D}_9^{\ast }-\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_2{\lambda }_P^{\ast }{D}_7^{\ast }-\left({\mathfrak{u}}_4^{\ast }\gamma +d_3+\mu \right)D_7^{\ast }\right)\\ & +{\lambda }_8\left(\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_1{\lambda }_P^{\ast }{D}_6^{\ast }+\left(1-{\mathfrak{u}}_1^{\ast }\right)\upsilon {\lambda }_H^{\ast }{D}_5^{\ast }-\left(\mu +d_4+\delta +{\mathfrak{u}}_3^{\ast }{\theta }_1\right)D_8^{\ast }\right)\\ & +{\lambda }_9\left(\delta D_8^{\ast }+\left(1-{\mathfrak{u}}_2^{\ast }\right){\varphi }_2{\lambda }_P^{\ast }{D}_7^{\ast }-\left(\mu +d_5+{\mathfrak{u}}_3^{\ast }{\theta }_2\right)D_9^{\ast }\right)+{\lambda }_{10}\left({\mathfrak{u}}_3^{\ast }\kappa D_5^{\ast }-\left(\mu +\eta \right)D_{10}^{\ast }\right)\\ & +{\lambda }_{11}\left({\mathfrak{u}}_4^{\ast }\gamma D_7^{\ast }-\mu D_{11}^{\ast }\right),\end{array}$where

${\lambda }_H^{\ast }=\frac{{\text{β}}_1}{N}\left(D_6^{\ast }\left(t\right)+{\rho }_1{D}_7^{\ast }\left(t\right)+{\rho }_2{D}_8^{\ast }\left(t\right)+{\rho }_3{D}_9^{\ast }\left(t\right)\right)$ and ${\lambda }_P^{\ast }={\text{β}}_2\left(D_5^{\ast }\left(t\right)+{\omega }_1{D}_8^{\ast }\left(t\right)+{\omega }_2{D}_9^{\ast }\left(t\right)\right)$.