Michael Byamukama^*, Mohammed Kizito
Department of Mathematics, Mabara University of Science and Technology, Mbarara, Uganda.
DOI: 10.4236/jamp.2025.135103   PDF    HTML   XML   99 Downloads   464 Views   Citations

Abstract

A mathematical model integrating both vaccination and treatment strategies for childhood pneumonia is developed to assess conditions for effective disease control. Special attention is given to individuals who do not respond to the vaccine, and their role in pneumonia transmission dynamics is analyzed. The effective reproduction number is derived, and conditions for the asymptotic stability of both the disease-free and endemic equilibria are established. An optimal control analysis is conducted to determine the most effective strategies for minimizing pneumonia cases. Three time-dependent control measures are considered: 1) pneumonia prevention, 2) identification and treatment of carriers, and 3) enhanced treatment of infected individuals. Numerical simulations and sensitivity analysis reveal that while vaccination and treatment significantly reduce pneumonia prevalence, the presence of vaccine-unresponsive serotypes poses a challenge to disease eradication. Additionally, results indicate that implementing multiple control strategies is more effective than relying on a single intervention. A combination of all three measures yields the best outcomes in reducing both infection and carriage rates.

Keywords

Pneumonia, Effective Reproduction Number, Optimal Protection, Treatment, Vaccination, Serotypes

Share and Cite:

1. Introduction

Childhood pneumonia remains a leading cause of mortality among infants and children under five years old, posing a significant disease burden worldwide [1]-[4]. While it affects children in both developed and developing countries, the impact is most severe in sub-Saharan Africa and South Asia due to high child populations and inadequate healthcare systems. Among the various forms of pneumonia, Streptococcus pneumoniae is the most common cause of severe and fatal bacterial pneumonia in children [1] [2]. The infection targets the lower respiratory tract, leading to fluid accumulation in the airspaces and symptoms such as tachypnea, dyspnea, hypoxia, and cough [5]. In addition to respiratory complications, pneumonia can cause life-threatening invasive infections such as sepsis and meningitis.

Bacterial pneumonia is primarily transmitted through airborne droplets, but reactivation due to incomplete treatment is also possible. Infants may contract the infection through blood, particularly during or shortly after birth. Fortunately, pneumonia is treatable, with amoxicillin being the first-line therapy, administered twice daily for 3 - 5 days [6]. Clotrimazole serves as an alternative treatment. However, treatment failure can occur, requiring hospitalization for children whose respiratory rates do not improve within 48 - 72 hours of therapy [6].

Effective preventive measures for childhood pneumonia include national immunization programs, improved nutrition, reduced exposure to tobacco smoke and indoor air pollution, and strengthened HIV prevention and treatment strategies. However, the effectiveness of these interventions varies across countries, highlighting the need for region-specific strategies to reduce pneumonia-related mortality [7].

Vaccination is a key strategy in pneumonia prevention, with most governments implementing routine childhood immunization campaigns. The Prevenar 13 (PCV 13) vaccine is recommended for infants, requiring four doses at 2, 4, 6, and 12 - 15 months of age [4] [8] [9]. However, challenges persist, as some children carry multiple serotypes of S. pneumoniae, making them unresponsive to vaccination. This complicates pneumonia dynamics by allowing vaccinated individuals to remain infectious. Additionally, studies have shown a high incidence of pneumonia within the first year of life despite good PCV coverage, influenced by HIV infection, malnutrition, and the need for full vaccination series completion for optimal protection [8] [10]. Several risk factors further increase childhood pneumonia susceptibility, including inadequate breastfeeding, incomplete immunization, indoor air pollution, low birth weight, and severe malnutrition [11]. Addressing these factors through comprehensive prevention and treatment strategies is essential for reducing the global burden of childhood pneumonia.

Without a doubt, numerous mathematical models have been developed to provide deeper insights into the dynamics of pneumonia, including [12]-[24], among others. Additionally, the opportunistic nature of pneumonia has been explored in several co-infection models [25]-[32].

However, a close evaluation of these studies reveals that individuals who do not respond to the PCV 13 vaccine have received little attention. In this model, optimal protection is considered through a combination of pneumonia vaccination, raising community awareness on risk-reducing behaviors, such as avoiding smoking and alcohol consumption, practicing proper hand hygiene with soap or alcohol-based sanitizers, and improving immune response through proper nutrition. On the treatment side, optimal efforts involve timely medical intervention for infected individuals, screening for carriers, and providing appropriate treatment. To achieve maximum protection and treatment, this model assumes that all children are screened, treated, vaccinated, and that parents or guardians take proactive measures to protect them from pneumonia.

2. Compartmental Model Formulation

The model is partitioned into five compartments, namely; $S\left(t\right)$ -individuals with no pneumonia but at risk of contracting it, $P_c\left(t\right)$ -individuals who carry pneumonia with no symptoms but are capable of spreading it, $P_I\left(t\right)$ -individuals with pneumonia symptoms who are infectious, $R\left(t\right)$ -individuals who have previously recovered from pneumonia and $V\left(t\right)$ -individuals immunized against pneumonia. Combining all compartments gives the total population as:

$N\left(t\right)=S\left(t\right)+P_c\left(t\right)+P_I\left(t\right)+R\left(t\right)+V\left(t\right)$

The population is assumed to mix homogeneously, giving everyone an equal chance of contracting the infection. Upon recovery, some individuals receive vaccination, while others become susceptible again. Additionally, some individuals carry pneumonia serotypes not covered by the vaccine, meaning that even those initially thought to be protected may still become infected and join either the carrier or infected class. Finally, vaccinated individuals receive timely boosters, ensuring the vaccine’s effectiveness does not wane.

2.1. Pneumonia Model Parameters

The parameters of the model and their descriptions are tabulated below.

Table 1. Parameter values and their description.

Parameter	Meaning	Value	Source
$\text{Λ}$	Constant increment rate into susceptible population	0.5000	Assumed
$\mu$	Non-pneumonia related death rate in all compartments	0.002	[13]
$\theta$	Upon infection, individuals join carriers at this rate	0.3380	[25]
$\sigma$	Pneumonia induced death rate	0.3300	[16]
$\beta$	Medication rate of pneumonia carriers	0.5150	Assumed
$\tau$	Medication rate of infected children	0.9000	[23]
$\pi$	Carriers develop symptoms at this rate	0.1000	[14]
$\gamma$	Susceptible individuals are vaccinated at this rate	0.3000	Assumed
$\psi$	Upon recovery, some individuals become susceptible at this rate	0.3000	Assumed
$\delta$	Proportion of serotypes that become carriers	0.2000	Assumed
$\epsilon$	Rate of vaccine non-responsiveness	0.0020	[19]
$\xi$	Infectious chance of pneumonia carriers	0.4102	[25] [31]
$\varphi$	Upon recovery, some individuals are vaccinated at this rate	0.7000	Assumed
$p$	Probability of transmission per contact	0.8900	[16]
$k$	Rate of contact with either a carrier or an infected person	1-10	[13]

2.2. Compartmental Diagram of the Model

Figure 1. Pneumonia model compartmental diagram.

2.3. System of Equations Describing Pneumonia Model

The variables, parameters in Table 1, compartmental diagram in Figure 1 and the assumptions above are used to get the following system of equations.

$\frac{\text{d}S\left(t\right)}{\text{d}t}=\text{Λ}+\psi R\left(t\right)-\left(\alpha +\gamma +\mu \right)S\left(t\right),$

$\frac{\text{d}{P}_C\left(t\right)}{\text{d}t}=\theta \alpha S\left(t\right)+\delta \epsilon \alpha V\left(t\right)-k_1{P}_C\left(t\right),$

$\frac{\text{d}{P}_I\left(t\right)}{\text{d}t}=\left(1-\theta \right)\alpha S\left(t\right)+\left(1-\delta \right)\epsilon \alpha V\left(t\right)+\pi P_C\left(t\right)-k_2{P}_I\left(t\right),$ (1)

$\frac{\text{d}R\left(t\right)}{\text{d}t}=\beta P_C\left(t\right)+\tau P_I\left(t\right)-k_{3}R\left(t\right),$

$\frac{\text{d}V\left(t\right)}{\text{d}t}=\gamma S\left(t\right)+\varphi R\left(t\right)-\left(\epsilon \alpha +\mu \right)V\left(t\right),$

where $k_1=\mu +\beta +\pi ,k_2=\tau +\mu +\sigma$ and $k_3=\mu +\varphi +\psi$ , together with initial conditions $S\left(0\right)=S_0$ , $P_C\left(0\right)=P_{C0}$ , $P_I\left(0\right)=P_{I0}$ , $R\left(0\right)=R_0$ and $V\left(0\right)=V_0$ .

The entire population of system (1) $N\left(t\right)$ , alters at a rate

$\frac{\text{d}N}{\text{d}t}=\text{Λ}-\mu N-\sigma P_I.$ (2)

Children contract pneumonia at the force of infection given as

$\alpha =\frac{pk\left(P_I+\xi P_C\right)}{N},$ (3)

where $\xi$ is a scaling factor, and $pk$ refers to the pneumonia transmission potential.

3. Analytical Solutions of the Model

3.1. Basic Properties of the Model

Under this subsection, the focal properties of solutions of the pneumonia model which are significant in the proofs of stability are investigated. The state variables of the model are always non-negative and the solutions stay positive with respect to initial conditions in bounded region

$\Omega =\left\{\left(S,P_C,P_I,R,V\right)\in {ℝ}_{+}^5,N\le \frac{\text{Λ}}{\mu }\right\}.$ (4)

We seek to show that each state variable defined is non-negative for all $t>0$ in the bounded region (4) and the model is mathematically and epidemiologically relevant.

Lemma 1. The pneumonia compartment solutions $S\left(t\right),P_C\left(t\right),P_I\left(t\right),R\left(t\right)$ and $V\left(t\right)$ are non-negative for $t>0.$

Proof. Let the initial values be positive, that is, $S\left(0\right)>0$ , $P_C\left(0\right)>0$ , $P_I\left(0\right)>0$ , $R\left(0\right)>0$ and $V\left(0\right)>0$ , then for all $t>0$ we prove that $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ and $V\left(t\right)>0$ . In contrast, suppose that there is a time $t_a$ to the extent that $S\left(t_a\right)=0$ , $S^{\prime }\left(t_a\right)<0$ and $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ , $V\left(t\right)>0$ for $0<t_a<t$ or there is a time $t_b$ to the extent that $P_C\left(t_b\right)=0$ , ${P^{\prime }}_C\left(t_b\right)<0$ and $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ , $V\left(t\right)>0$ for $0<t_b<t$ or there is a time $t_c$ to the extent that $P_I\left(t_c\right)=0$ , ${P^{\prime }}_I\left(t_c\right)<0$ and $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ , $V\left(t\right)>0$ for $0<t_c<t$ or there is a time $t_d$ to the extent that $R\left(t_d\right)=0$ , $R^{\prime }\left(t_d\right)<0$ and $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ , $V\left(t\right)>0$ for $0<t_d<t$ or there is a time $t_e$ to the extent that $V\left(t_e\right)=0$ , $S^{\prime }\left(t_e\right)<0$ and $S\left(t\right)>0$ , $P_C\left(t\right)>0$ , $P_I\left(t\right)>0$ , $R\left(t\right)>0$ , $V\left(t\right)>0$ for $0<t_d<t$ . In reference to system (1), we obtain

$\frac{\text{d}S\left(t_a\right)}{\text{d}t}=\text{Λ}+\psi R\left(t_a\right)>0,$

$\frac{\text{d}{P}_C\left(t_b\right)}{\text{d}t}=\frac{kp{P}_I\left(t_b\right)}{N}\left(\theta S\left(t_b\right)+\delta \epsilon V\left(t_b\right)\right)>0,$

$\frac{\text{d}{P}_I\left(t_c\right)}{\text{d}t}=\frac{kp{P}_c\left(t_c\right)}{N}\left(\left(1-\theta \right)S\left(t_c\right)+\left(1-\delta \right)\epsilon V\left(t_c\right)\right)+\pi P_C\left(t_c\right)>0,$

$\frac{\text{d}R\left(t_d\right)}{\text{d}t}=\beta P_C\left(t_d\right)+\tau P_I\left(t_d\right)>0,$

$\frac{\text{d}V\left(t_e\right)}{\text{d}t}=\gamma S\left(t_e\right)+\varphi R\left(t_e\right)>0.$

This gives a contradiction and thus, the compartmental solutions are non-negative for $t>0$ . □

Lemma 2. The region $\Omega$ described in (4) is bounded in ${ℝ}_{+}^5$ .

Proof. If there are no fatalities due to pneumonia, Equation (2) becomes

$\frac{\text{d}N\left(t\right)}{\text{d}t}+\mu N\left(t\right)\le \text{Λ}.$ (5)

Given that $N\left(0\right)=N_0$ , integrate (5) to achieve

$N\left(t\right)\le \frac{\text{Λ}}{\mu }+\left(N_0-\frac{\text{Λ}}{\mu }\right){\text{e}}^{-\mu t}.$ (6)

After a very long time, it can be shown that

$N\left(t\right)\le \frac{\text{Λ}}{\mu }.$ (7)

Inequality (7) gives the population limit and the solution of system (1) is bounded in the region (4).

Therefore, by Lemmas 1 and 2, it is sufficient to study the dynamics of pneumonia model (1).

3.2. Pneumonia-Free Equilibrium Point

The pneumonia free equilibrium is possible when there are no carriers and infected individuals in the compartment and hence no recovered individuals, that is, $P_C=P_I=R=0.$ Equating the right hand side of system (1) to zero and solving

for compartmental variables $S$ and $V$ , leads to; $S=\frac{\text{Λ}}{\gamma +\mu }$ and $V=\frac{\text{Λ}\gamma }{\mu \left(\gamma +\mu \right)}$ . Thus, the pneumonia-free equilibrium point is given as

$E_0=\left(\frac{\text{Λ}}{\gamma +\mu },0,0,0,\frac{\text{Λ}\gamma }{\mu \left(\gamma +\mu \right)}\right),$ (8)

implying that at equilibrium, the population consists of susceptible and vaccinated individuals only.

3.3. Effective Reproduction Number

The effective reproduction number is the mean number of secondary pneumonia cases caused by one infectious pneumonia carrier or infected individual in a community subjected to vaccination. We use the next generation matrix approach as applied in [16] [23] [25]. The infected classes are given a priority in this strategy such that system (1) is re-arranged to become:

$\frac{\text{d}{P}_C\left(t\right)}{\text{d}t}=\theta \alpha S\left(t\right)+\delta \epsilon \alpha V\left(t\right)-k_1{P}_C\left(t\right),$

$\frac{\text{d}{P}_I\left(t\right)}{\text{d}t}=\left(1-\theta \right)\alpha S\left(t\right)+\left(1-\delta \right)\epsilon \alpha V\left(t\right)+\pi P_C\left(t\right)-k_2{P}_I\left(t\right),$

$\frac{\text{d}S\left(t\right)}{\text{d}t}=\text{Λ}+\psi R\left(t\right)-\left(\alpha +\gamma +\mu \right)S\left(t\right),$ (9)

$\frac{\text{d}R\left(t\right)}{\text{d}t}=\beta P_C\left(t\right)+\tau P_I\left(t\right)-k_{3}R\left(t\right),$

$\frac{\text{d}V\left(t\right)}{\text{d}t}=\gamma S\left(t\right)+\varphi R\left(t\right)-\left(\epsilon \alpha +\mu \right)V\left(t\right).$

From system (9), apply the symbol $ℐ$ for the matrix of pneumonia infection terms and $\mathcal{T}$ for the matrix of other transfer terms in the infectious compartments such that

$ℐ=\left[\begin{array}{c}\theta \alpha S+\delta \epsilon \alpha V\\ \left(1-\theta \right)\alpha S+\left(1-\delta \right)\epsilon \alpha V\end{array}\right]\text{and}\mathcal{T}=\left[\begin{array}{c}{k}_1{P}_C\\ -\pi P_C+k_2{P}_I\end{array}\right].$ (10)

Evaluating the Jacobian matrices of $ℐ$ and $\mathcal{T}$ in (10) at pneumonia-free equilibrium point (8), we respectively obtain the following matrices

$I=\left[\begin{array}{cc}\frac{pk\xi }{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)& \frac{pk}{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)\\ \frac{pk\xi }{\gamma +\mu }\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)& \frac{pk}{\gamma +\mu }\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)\end{array}\right],$ (11)

$T=\left[\begin{array}{cc}{k}_1& 0\\ -\pi & k_2\end{array}\right].$ (12)

Matrices (11) and (12) are used to obtain the next generation matrix $I{T}^{-1}$ as

$I{T}^{-1}=\frac{pk}{k_1{k}_2}\left[\begin{array}{cc}\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \right)& k_1\left(\theta \mu +\delta \epsilon \gamma \right)\\ \left(\xi k_2+\pi \right)\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \right)& k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)\end{array}\right].$ (13)

The eigenvalues of matrix (13) are $\lambda =0$ and $\lambda =\frac{pk}{k_1{k}_2\left(\gamma +\mu \right)}\left(k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)+\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \gamma \right)\right)$ . Mathematically,

the effective reproduction number is the dominant eigenvalue of the next generation matrix (13). Thus, the effective reproduction number of the pneumonia model is given as

$R_e=\frac{pk}{k_1{k}_2\left(\gamma +\mu \right)}\left(k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)+\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \gamma \right)\right).$ (14)

The basic reproduction number $R_0$ is obtained on condition that there is no vaccination strategy, that is, $\gamma =0$ . This leads to

$R_0=\frac{pk}{k_1{k}_2}\left(k_1\left(1-\theta \right)+\theta \left(\xi k_2+\pi \right)\right).$ (15)

The relationship between $R_e$ and $R_0$ in (14) and (15) respectively is given as

$R_e=\frac{\mu R_0}{\gamma +\mu }+\frac{\epsilon \gamma pk}{k_1{k}_2\left(\gamma +\mu \right)}\left(k_1\left(1-\delta \right)+\left(\xi k_2+\pi \right)\delta \right).$ (16)

Furthermore, if there are no children that do not react to the vaccine, that is, $\epsilon =0$ , Equation (16) gives

$R_{ev}=\frac{\mu R_0}{\gamma +\mu }.$ (17)

Proposition 3. For all positive parameters in the effective reproduction number, then, $R_e<R_{ev}<R_0$ .

3.4. Local and Global Stability of Pneumonia-Free Equilibrium Point

Theorem 1. The pneumonia-free equilibrium point (8) is locally asymptotically stable as long as $R_e<1$ .

Proof. The Jacobian matrix of system (1) evaluated at the pneumonia-free equilibrium point is given by

$J_{E_0}=\left[\begin{array}{ccccc}-k_0& -pk\theta \xi \frac{\mu }{\gamma +\mu }& -pk\left(1-\theta \right)\frac{\mu }{\gamma +\mu }& \psi & 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)-k_1& \frac{pk}{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)& 0& 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(z\right)-\pi & \frac{pk}{\gamma +\mu }\left(z\right)-k_2& 0& 0\\ 0& \beta & \tau & -k_3& 0\\ \gamma & \frac{pk\epsilon \gamma \xi }{\gamma +\mu }& \frac{pk\epsilon \gamma }{\gamma +\mu }& \varphi & -\mu \end{array}\right],$ (18)

where $k_0=-\left(\mu +\gamma \right)$ , $k_1=\left(\mu +\beta +\pi \right)$ , $k_2=\left(\tau +\mu +\sigma \right)$ , $z=\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma$ and $k_3=\left(\mu +\varphi +\psi \right).$ The characteristic polynomial of the Jacobian matrix (18) is given by

$P\left(\lambda \right)=\left(\lambda +k_0\right)\left(\lambda +k_3\right)\left(\lambda +\mu \right)\left({\lambda }^2+d_0\lambda +d_1\right),$ (19)

where $d_0=k_1+k_2+\frac{kp}{\gamma +\mu }\left(\mu \left(\left(1-\xi \right)\theta -1\right)+ϵ\gamma \left(\delta \left(1-\xi \right)-1\right)\right)$ and $d_1=k_1{k}_2\left(1-R_e\right)>0$ , if $R_e<1$ . From the characteristic polynomial (19), we see that $\lambda =-k_0$ , $\lambda =-k_3$ , and $\lambda =-\mu$ are clearly negative. By the Routh-Hurwitz criterion, the sub-characteristic equation

${\lambda }^2+d_0\lambda +d_1=0$

has negative eigenvalues if $d_0>0,d_1>0$ , which is only possible when $R_e<1$ . Therefore, we conclude that the pneumonia-free equilibrium point is locally asymptotically stable if $R_e<1$ . □

Theorem 2. The pneumonia-free equilibrium point $E_0$ is globally stable if $R_e<1$ .

Proof. Consider the Lyapunov function candidate:

$W={\vartheta }_1{P}_C+{\vartheta }_2{P}_I,$ (20)

where ${\vartheta }_1$ and ${\vartheta }_2$ are positive constants. The time derivative of $W$ in (20) is given by

$\frac{\text{d}W}{\text{d}t}={\vartheta }_1\frac{\text{d}{P}_C}{\text{d}t}+{\vartheta }_2\frac{\text{d}{P}_I}{\text{d}t}.$ (21)

Substituting the equations for $\frac{\text{d}{P}_C}{\text{d}t}$ and $\frac{\text{d}{P}_I}{\text{d}t}$ into Equation (21) gives

$\frac{\text{d}W}{\text{d}t}={\vartheta }_1\left(\alpha \theta S+\delta \epsilon \alpha V-k_1{P}_C\right)+{\vartheta }_2\left(\alpha \left(1-\theta \right)S+\pi P_C+\alpha \epsilon \left(1-\delta \right)V-k_2{P}_I\right).$ (22)

Substituting for $\alpha =pk\left(\frac{P_I+\xi P_C}{N}\right)$ in (22) and noting that $\frac{S}{N}\approx \frac{\mu }{\gamma +\mu }$ , $\frac{V}{N}\approx \frac{\gamma }{\gamma +\mu }$ at the pneumonia-free equilibrium point (8), we obtain

$\begin{array}{c}\frac{\text{d}W}{\text{d}t}\le \left(\left(pk\xi \left(\mu \theta +\gamma \delta \epsilon \right)-k_1\right){\vartheta }_1+\left(pk\xi \left(\mu \left(1-\theta \right)+\gamma \left(1-\delta \right)ϵ\right)+\pi \right){\vartheta }_2\right)P_C,\\ +\left(\left(pk\left(\mu \theta +\gamma \delta \epsilon \right)\right){\vartheta }_1+\left(pk\left(\left(1-\theta \right)\mu +\left(1-\delta \right)ϵ\gamma \right)-k_2\right){\vartheta }_2\right)P_I.\end{array}$ (23)

Let ${\vartheta }_1>0$ and ${\vartheta }_2=\frac{\left(pk\xi \left(\mu \theta +\gamma \delta \epsilon \right)-k_1\right){\vartheta }_1}{\left(pk\xi \left(\mu \left(1-\theta \right)+\gamma \left(1-\delta \right)\epsilon \right)+\pi \right)}$ .

Thus, inequality (23) simplifies to

$\frac{\text{d}W}{\text{d}t}\le \left(k_1{k}_2{\vartheta }_1{P}_C+\frac{k_1{\vartheta }_1{P}_I}{\gamma +\mu }\right)\left(R_e-1\right)$

which confirms the global stability of the pneumonia-free equilibrium point when $R_e<1$ . □

3.5. Pneumonia Endemic State

The pneumonia endemic equilibrium point $E_e$ exists when pneumonia regularly occurs in the community. It is determined from system (1) by equating its right-hand side to zero:

$\{\begin{array}{l}\alpha \theta S+\delta \epsilon \alpha V-k_1{P}_C=0,\hfill \\ \alpha \left(1-\theta \right)S+\pi P_C+\alpha \epsilon \left(1-\delta \right)V-k_2{P}_I=0,\hfill \\ \text{Λ}+\varphi R-\left(\alpha +\mu +\gamma \right)S=0,\hfill \\ \beta P_C+\tau P_I-k_{3}R=0,\hfill \\ \gamma S+\phi R-\left(\mu +ϵ\alpha \right)V=0.\hfill \end{array}$ (24)

Solving system (24) together with $\text{Λ}-\mu N-\sigma P_I=0$ and $\alpha =pk\left(\frac{P_I+\xi P_C}{N}\right)$ leads to the endemic steady-state. The endemic steady-state $E_e=\left(S_e,P_{Ce},P_{Ie},R_e,V_e\right)$ is given by:

$S_e=\frac{k_1{k}_2{N}^{\text{*}}\left(\gamma +\mu \right)}{pk\left(k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)+\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \gamma \right)\right)},$

$P_{Ce}=\frac{-B+\sqrt{B^2-4AC}}{2A},$

$P_{Ie}=\frac{\text{Λ}-\mu N^{\text{*}}}{\sigma },$

$R_e=\frac{\sigma \beta P_{Ce}+\tau \left(\text{Λ}-\mu N^{\text{*}}\right)}{\sigma k_3},$

$V_e=\frac{\gamma N^{\text{*}2}+\beta \phi N^{\text{*}}{P}_{Ce}+\tau \phi \left(\text{Λ}-\mu N^{\text{*}}\right)}{\sigma k_3\left(\mu N^{\text{*}}+m\epsilon \left(\xi P_{Ce}+\text{Λ}-\mu N^{\text{*}}\right)\right)},$

where $A=D_{1}D{k}_{3}m\epsilon {\xi }^2+d_1\beta \phi N^{\text{*}}-k_1{k}_3{D}_1{N}^{\text{*}}m\xi \epsilon$ , $B=\left(k_3{D}_{1}D\left(m\xi \epsilon \text{Λ}+\xi \mu N^{\text{*}}\right)+d_1\beta \phi N^{\text{*}}\text{Λ}+d_1\gamma \xi N^{\text{*}}-k_1{k}_3{D}_{1}m\epsilon N^{\text{*}}\text{Λ}\right)$ , $C=\left(k_3{D}_{1}D\mu +d_1\gamma \right)\left(\text{Λ}-\mu N^{\text{*}}\right)N^{\text{*}}$ , $D=m\theta D_1$ , $d_1=m\delta \epsilon$ , $m=pk$ and $D_1=\frac{k_1{k}_2{N}^{\text{*}}\left(\gamma +\mu \right)}{pk\left(k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)+\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \gamma \right)\right)}$ .

3.6. Bifurcation Analysis of the Pneumonia Model

Let $S=r_1,P_C=r_2,P_I=r_3,R=r_4,V=r_5$ and $N=r_1+r_2+r_3+r_4+r_5$ . Re-write model (1) as

$\frac{\text{d}{r}_1}{\text{d}t}=\text{Λ}+\psi r_4-\left(\alpha +\gamma +\mu \right)r_1,$

$\frac{\text{d}{r}_2}{\text{d}t}=\theta \alpha r_1+\delta \epsilon \alpha r_5-k_1{r}_2,$

$\frac{\text{d}{r}_3}{\text{d}t}=\left(1-\theta \right)\alpha r_1+\left(1-\delta \right)\epsilon \alpha r_5+\pi r_2-k_2{r}_3,$ (25)

$\frac{\text{d}{r}_4}{\text{d}t}=\beta r_2+\tau r_3-k_3{r}_4,$

$\frac{\text{d}{r}_5}{\text{d}t}=\gamma r_1+\varphi r_4-\left(\epsilon \alpha +\mu \right)r_5,$

where $\alpha =\frac{pk\left(r_3+\xi r_2\right)}{r_1+r_2+r_3+r_4+r_5}$ is the pneumonia force of infection. Evaluating the Jacobian matrix of system (25) at pneumonia-free equilibrium, leads to

$J_{E_0}=\left[\begin{array}{ccccc}-k_0& -pk\theta \xi \frac{\mu }{\gamma +\mu }& -pk\left(1-\theta \right)\frac{\mu }{\gamma +\mu }& \psi & 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)-k_1& \frac{pk}{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)& 0& 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(z\right)-\pi & \frac{pk}{\gamma +\mu }\left(z\right)-k_2& 0& 0\\ 0& \beta & \tau & -k_3& 0\\ \gamma & \frac{pk\epsilon \gamma \xi }{\gamma +\mu }& \frac{pk\epsilon \gamma }{\gamma +\mu }& \varphi & -\mu \end{array}\right].$ (26)

Make $k=k^{\text{*}}$ the bifurcation parameter and assess the model when $R_e=1$ such that $k^{\text{*}}=\frac{k_1{k}_2\left(\gamma +\mu \right)}{p\left(k_1\left(1-\delta \right)\epsilon \gamma +\left(1-\theta \right)\mu +\left(\pi +\xi k_2\right)\left(\delta \epsilon {\gamma }_{\mu }\theta \right)\right)}$ . The charactereistic polynomial of the Jacobian (26) when $R_e=1$ is obtained as

$P\left(\lambda \right)=\lambda \left(\lambda +k_0\right)\left(\lambda +\mu \right)\left(\lambda +k_3\right)\left(\lambda +k_1+k_2+\frac{kp}{\gamma +\mu }\left(\mu \left(\left(1-\xi \right)\theta -1\right)+\epsilon \gamma \left(\delta \left(1-\xi \right)-1\right)\right)\right).$ (27)

Considering the characteristic polynomial (27), there is a simple zero eigenvalue and all other eigenvalues are less than zero. It is now sufficient to apply the Centre-Manifold theory, as applied in [25] [30] [31] in order to deduce whether there is possibility of forward or backward bifurcation at $R_e=1$ .

The right eigenvector $x=\left(x_1,x_2,x_3,x_4,x_5\right)$ that matches up with the zero eigenvalue is got from

$J_{E_0}{x}^{\text{T}}=\left[\begin{array}{ccccc}-k_0& -pk\theta \xi \frac{\mu }{\gamma +\mu }& -pk\left(1-\theta \right)\frac{\mu }{\gamma +\mu }& \psi & 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)-k_1& \frac{pk}{\gamma +\mu }\left(\theta \mu +\delta \epsilon \gamma \right)& 0& 0\\ 0& \frac{pk\xi }{\gamma +\mu }\left(z\right)-\pi & \frac{pk}{\gamma +\mu }\left(z\right)-k_2& 0& 0\\ 0& \beta & \tau & -k_3& 0\\ \gamma & \frac{pk\epsilon \gamma \xi }{\gamma +\mu }& \frac{pk\epsilon \gamma }{\gamma +\mu }& \varphi & -\mu \end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right]=\left[\begin{array}{l}0\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \end{array}\right].$ (28)

Evaluating system (28) leads to $x_1=\left(\frac{M\left(\theta \mu +\delta \epsilon \gamma \right)-m{k}_1\left(\left(1-\theta \right)\mu +\gamma \tau \xi \right)}{k_0\left(m\xi \left(\theta \mu +\delta \epsilon \gamma \right)-k_1\right)}\right)x_3$ ,$x_3>0$ , $x_2=\left(\frac{m\left(\theta \mu +\delta \epsilon \gamma \right)}{m\xi \left(\theta \mu +\delta \epsilon \gamma \right)-k_1}\right)x_3$ $x_4=\left(\frac{\beta m\left(\theta \mu +\delta \epsilon \gamma \right)}{k_3\left(m\xi \left(\theta \mu +\delta \epsilon \gamma \right)-k_1\right)}+\frac{\tau }{k_3}\right)x_3$ and $x_5=\left(\frac{A\left(\theta \mu +\delta \epsilon \gamma \right)+k_1{k}_{3}B}{\mu k_0{k}_3\left(m\xi \left(\theta \mu +\delta \epsilon \right)-k_1\right)}\right)x_3$ , where $M=m^2\mu \xi \left(2\theta -1\right)+\psi m\left(\tau \xi +\beta \right)$ , $A=k_3\gamma m\left(m\mu \xi \left(2\theta -1\right)+\psi \left(\tau \xi +\beta \right)\right)+m^2\epsilon \gamma \xi k_0{k}_3+\beta m{k}_0+\mu k_3\left(m\epsilon \gamma k_3+\tau \right)$ , $B=\gamma m\left(\left(1-\theta \right)\mu +\psi \tau \xi \right)-\mu \left(m\epsilon \gamma k_3+\tau \right)$ and $m=\frac{p{k}^{\text{*}}}{\gamma +\mu }.$

In the same way, the left eigenvector $y=\left(y_1,y_2,y_3,y_4,y_5\right)$ is obtained from $y{J}_{E_0}=0$ , which on evaluation, leads to $y_1=y_4=y_5=0$ , $y_2=\left(\frac{m\left(\theta \mu +\delta \epsilon \gamma \right)}{n\xi \left(\theta \mu +\delta \epsilon \gamma \right)-k_1}\right)y_3$ and $y_3>0$ . The bifurcation quantities $a_1$ and $a_2$ are acquired as:

$\begin{array}{c}{a}_1={\displaystyle \sum }_{k,i,j=1}^5{y}_k{x}_i{x}_j\frac{{\partial }^2{m}_k}{\partial x_i\partial x_j}\\ =\frac{-2p{k}^{\text{*}}\mu x_3{y}_3}{\left(\gamma +\mu \right){\left(k_1-m\xi M\right)}^3}\left(k_1^{3}A+k_1^{2}MB+k_1{M}^{2}C+M^{3}D\right),\end{array}$ (29)

$\begin{array}{c}{a}_2={\displaystyle \sum }_{k,i=1}^5{y}_k{x}_i\frac{{\partial }^2{m}_k}{\partial x_i\partial k^{\text{*}}}\\ =\frac{p{x}_3{y}_3\left(m\left(1+\xi \right)M-k_1\right)}{\left(\gamma +\mu \right){\left(m\xi M-k_1\right)}^2}\times \left(\left(m\xi M-k_1\right)\left(\mu \left(1-\theta \right)+\epsilon \gamma \left(1-\delta \right)+mM\right)\right),\end{array}$ (30)

where

$M=\theta \mu +\delta \epsilon \gamma ,$

$\begin{array}{l}A=\left(1-\theta \right)\left(2\gamma +5\mu \right)+\epsilon \left(1-\delta \right)\left(\mu \left(1-\xi \right)+\gamma \right)+\mu \xi \left(1-\epsilon \right)\\ +{\xi }^{2}m\mu M^2\left(\left(1-\xi \right)+{\xi }^2\left(\delta \epsilon +\theta \right)\right),\end{array}$

$\begin{array}{l}B=m(\xi \left(2m\mu {\xi }^2+3\mu \xi \left(1-\epsilon \right)+4\gamma \left(1-\theta \right)+\gamma \mu \left(1-\delta \right)\right)\\ +\epsilon \left(\gamma \left(1-\xi \right)+\delta \left(\mu +m\gamma \xi \right)\right)),\end{array}$

$\begin{array}{l}C=m^2(\left(\mu \xi +\gamma \left(1+{\xi }^2\right)\right)\left(1-\theta \right)+\left(\mu {\xi }^2\left(1+3\theta \right)+\delta \epsilon \gamma \right)\left(1-\xi \right)\\ +\epsilon \gamma \xi \left(1-\delta \right)+{\xi }^2\delta \epsilon \left(\gamma +\mu \right)),\end{array}$

and

$D=m^3{\xi }^2\left(2\mu +\xi \left(2\gamma +\mu \left(2+\xi \epsilon \left(\theta +\delta \right)\right)\right)\right).$

From the expressions of (29) and (30), we realize that $a_1<0$ , $a_2>0$ . This implies that model (25) possesses a forward bifurcation at the pneumonia-free equilibrium point and thus, there is at least one positive endemic equilibrium point at $R_e=1$ . By this result, the following theorems hold.

Theorem 3. If the effective reproduction number, $R_e$ is equal to one, then the system undergoes a forward bifurcation at the pneumonia-free equilibrium, guaranteeing the existence of at least one positive endemic equilibrium point.

Theorem 4. The endemic state $E_e$ is locally asymptotically stable if $R_e>1$ .

3.7. Sensitivity Investigation of the Effective Reproduction Number

In the struggle to reduce pneumonia burden, the first step is to ensure that the effective reproduction number is smaller than one, that is, $R_e<1$ . To achieve this, we need to determine the most vital parameters that can be reduced or increased by analyzing the sensitivity indices of the effective reproduction number with respect to the parameters. The methodology in [16] [23] [25] is applied, where the relation for a normalized sensitivity index of the effective reproduction number, $R_e$ with respect to a parameter $X$ is defined as:

$D_X^{\left(R_e\right)}=\frac{\partial R_e}{\partial X}\times \frac{X}{R_e}.$

For example, consider the rate of contact parameter $p$ with

$R_e=\frac{pk}{k_1{k}_2\left(\gamma +\mu \right)}\left(k_1\left(\left(1-\theta \right)\mu +\left(1-\delta \right)\epsilon \gamma \right)+\left(\xi k_2+\pi \right)\left(\theta \mu +\delta \epsilon \gamma \right)\right),$

then the sensitivity index of $p$ is given as:

$D_p^{\left(R_e\right)}=\frac{\partial R_e}{\partial p}\times \frac{p}{R_e}=1.$

Applying parameter values in Table 1, the sensitivity indices are generated as shown in Table 2.

Highlights from Sensitivity Analysis

From Table 2 and Figure 2, it is evident that parameters with positive sensitivity

Table 2. Sensitivity indices of the effective reproduction number.

Parameter	Sensitivity index	Parameter	Sensitivity index
$p$	+1.0000	$\tau$	−0.6949
$k$	+1.0000	$\sigma$	−0.0645
$\theta$	+0.0574	$\beta$	−0.2248
$\delta$	+0.0287	$\pi$	−0.0140
$\epsilon$	+0.3333	$\gamma$	−0.6627
$\xi$	+0.2393	$\tau$	−0.6949

Figure 2. Visualization of model sensitivity indices.

indices, namely $p,k,\epsilon ,\xi ,\delta$ and $\theta$ contribute to an increased prevalence of pneumonia, as their increase leads to a higher effective reproduction number. Conversely, parameters $\beta ,\sigma ,\gamma ,\pi$ and $\tau$ have a negative impact on the effective reproduction number, meaning that increasing their values directly reduces $R_e$ , thereby lowering the prevalence of pneumonia. By focusing on reducing parameters that increase $R_e$ and enhancing those that decrease it, policymakers and healthcare practitioners can develop targeted interventions to lower the prevalence of pneumonia and improve public health outcomes.

4. Optimal Protection and Treatment in the Dynamics of Pneumonia

4.1. Introduction

Optimal control theory has demonstrated to be an efficient strategy in apprehending mechanisms through which infectious diseases can be restricted by concocting the optimal intervention plans. The priority is always to minimize the cost of infection or the cost of implementing the control strategy, or both. Various infectious studies have embraced optimal control theory, for example check [15] [19] [21] [22] [28] [31], where the results have been magnificent in proposing the best plans for disease reduction if not elimination. Besides vaccination and treatment in the current pneumonia model (1), we further apply optimal control theory, where we seek the following three controls in the formulation of the desired optimal control model.

Control 1. $u_1$ : Pneumonia protection effort, that prevents susceptible and the unreactive vaccinated individuals from getting infected with pneumonia. We seek to achieve this strategy through community awareness programs based on proper sanitation, use of protective gears like masks when in polluted environments, taking the required diet and educating the sick to self-isolate.

Control 2. $u_2$ : Effort aimed at reducing pneumonia infected individuals through testing and providing timely and proper treatment. All children with pneumonia symptoms should be screened and treated promptly.

Control 3. $u_3$ : Pneumonia carriage reduction effort. This involves identifying carriers and subjecting them to treatment. This can be achieved through contact tracing, where all contacts of a previously infected individuals are screened and treated if found positive.

These controls will be instructive in identification of the most suitable intervention plan that will probably lead to significant reduction of pneumonia cases.

4.2. Model Formulation

Let $\left\{S\left(t\right),P_C\left(t\right),P_I\left(t\right),R\left(t\right),V\left(t\right)\right\}$ be the state variables of the system at time $t$ , and $\left\{u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\right\}$ be the control inputs to model (1) at time $t$ . Then, the optimal system is described by the following five first-order differential equations:

$\frac{\text{d}S}{\text{d}t}=\text{Λ}+\psi R-\left(\left(1-u_1\right)\alpha +\mu +\gamma \right)S,$

$\frac{\text{d}{P}_C}{\text{d}t}=\left(1-u_1\right)\alpha \theta S+\delta \epsilon \left(1-u_1\right)\alpha V-\left(\mu +\pi +\beta +u_3\right)P_C,$

$\frac{\text{d}{P}_I}{\text{d}t}=\left(1-u_1\right)\alpha \left(1-\theta \right)S+\pi P_C+\alpha \left(1-u_1\right)\epsilon \left(1-\delta \right)V-\left(\mu +\sigma +\tau +u_2\right)P_I,$ (31)

$\frac{\text{d}R}{\text{d}t}=\left(\beta +u_3\right)P_C+\left(\tau +u_2\right)P_I-\left(\varphi +\psi +\mu \right)R,$

$\frac{\text{d}V}{\text{d}t}=\gamma S+\varphi R-\left(\mu +\epsilon \left(1-u_1\right)\alpha \right)V.$

The goal is to minimize the number of infected individuals $P_I$ and carriers $P_C$ in the community while keeping the cost associated with controls $u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)$ at a minimum value. Therefore, we seek to minimize the objective function $J$ , defined over a feasible set of controls $u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)$ applied over the pre-defined finite time interval $\left[t_0,t_m\right]$ , given by:

$J={\displaystyle {\int }_{t_0}^{t_m}\left[{\mathcal{C}}_1{P}_C+{\mathcal{C}}_2{P}_I+\frac{1}{2}\left(W_1{u}_1^2+W_2{u}_2^2+W_3{u}_3^2\right)\right]\text{d}t},$ (32)

where $C_1,C_2,W_1,W_2$ and $W_3$ are positive coefficients. In the objective function (32), $C_1{P}_C$ represents the cost associated with pneumonia carriers, while $C_2{P}_I$ accounts for the cost related to individuals infected with pneumonia. The coefficients $W_1,W_2$ and $W_3$ represent the additional costs incurred for implementing each respective control strategy $u_1\left(t\right),u_2\left(t\right)$ and $u_3\left(t\right)$ . Thus, we look for the optimal plans $\left(u_1^{\text{*}},u_2^{\text{*}},u_3^{\text{*}}\right)$ such that:

$J\left(u_1^{\text{*}},u_2^{\text{*}},u_3^{\text{*}}\right)=\underset{\left(u_1,u_2,u_3\right)\in U}{\text{min}}J\left(u_1,u_2,u_3\right),$

where the admissible control set is defined as:

$U=\left\{\left(u_1,u_2,u_3\right):0\le u_1<1,0\le u_2<1,0\le u_3<1\right\}.$

The Pontryagin’s maximum principle is used to determine the necessary conditions for the satisfaction of an optimal control problem. By this principle, we define the Hamiltonian, $h$ as:

$\begin{array}{l}h=C_1{P}_C+C_2{P}_I+\frac{1}{2}\left(W_1{u}_1^2+W_2{u}_2^2+W_3{u}_3^2\right)\\ +L_1\frac{\text{d}S}{\text{d}t}+L_2\frac{\text{d}{P}_C}{\text{d}t}+L_3\frac{\text{d}{P}_I}{\text{d}t}+L_4\frac{\text{d}R}{\text{d}t}+L_5\frac{\text{d}V}{\text{d}t},\end{array}$

where $L_1,L_2,L_3,L_4$ and $L_5$ are the costate variables corresponding to the state variables. The characterization of the optimal control problem is presented in the theorem below.

Theorem 5. Let $u^{*}=\left(u_1^{*},u_2^{*},u_3^{*}\right)\in U$ be a set of optimal plans and $X=\left(S,P_C,P_I,R,V\right)$ the corresponding set of solutions that minimizes $J$ over $U$ . Then, there exist costate variables $L_X$ such that:

$\frac{\text{d}{L}_X}{\text{d}t}=-\frac{\partial h}{\partial X},$ (33)

$L_X\left(t_m\right)=0,$ (34)

$\frac{\partial k}{\partial u_i}\left(u_i^{\text{*}}\right)=0,i=1,2,3.$ (35)

Equations (33), (34), and (35) represent the adjoint conditions, transversality conditions, and optimality conditions, respectively.

Proof. We apply Pontryagin’s maximum principle for bounded controls as applied in [19] [22] [31]. The adjoint system is obtained from Equation (33) as follows:

$\frac{\text{d}{L}_1}{\text{d}t}=-\frac{\partial h}{\partial S}=\left(L_1-L_5\right)\gamma +\left(L_1-\theta L_2-\left(1-\theta \right)L_3\right)\left(1-u_1\right)\alpha +L_1\mu ,$

$\begin{array}{l}\frac{\text{d}{L}_2}{\text{d}t}=-\frac{\partial h}{\partial P_C}=\left(\left(L_1-\theta L_2-\left(1-\theta \right)L_3\right)\frac{S}{N}+\left(L_1-\delta L_2-\left(1-\delta \right)L_3\right)\frac{V}{N}\right)\\ \times \left(1-u_1\right)\xi pk+\left(L_2-L_4\right)\left(\beta +u_3\right)+\left(L_2-L_3\right)\pi +L_2\mu -C_1,\end{array}$

$\begin{array}{l}\frac{\text{d}{L}_3}{\text{d}t}=-\frac{\partial h}{\partial P_I}=\left(\left(L_1-\theta L_2-\left(1-\theta \right)L_3\right)\frac{S}{N}+\left(L_1-\delta L_2-\left(1-\delta \right)L_3\right)\frac{V}{N}\right)\\ \times \left(1-u_1\right)pk+\left(L_3-L_4\right)\left(\tau +u_2\right)+L_3\left(\mu +\sigma \right)-C_2,\end{array}$ (36)

$\frac{\text{d}{L}_4}{\text{d}t}=-\frac{\partial h}{\partial R}=\left(L_4-L_1\right)\psi +\left(L_4-L_5\right)\varphi +L_4\mu ,$

$\frac{\text{d}{L}_5}{\text{d}t}=-\frac{\partial h}{\partial V}=\left(L_5-\delta L_2-\left(1-\delta \right)L_3\right)\epsilon \alpha +L_5\mu .$

The transversality conditions are:

$L_1\left(t_m\right)=L_2\left(t_m\right)=L_3\left(t_m\right)=L_4\left(t_m\right)=L_5\left(t_m\right)=0.$

The optimal conditions are obtained from Equation (35) by differentiating the Hamiltonian equation with respect to the control plans. This leads to:

$\frac{\partial h}{\partial u_1}=W_1{u}_1+\left(L_1-\theta L_2-\left(1-\theta \right)L_3\right)\alpha S+\left(L_5-\delta L_2-\left(1-\delta \right)L_3\right)\alpha \epsilon V,$

$\frac{\partial h}{\partial u_2}=W_2{u}_2-\left(L_4-L_3\right)P_I,$ (37)

$\frac{\partial h}{\partial u_3}=W_3{u}_3-\left(L_4-L_2\right)P_C.$

From system (37), we have the following boundary conditions:

$u_1^{\text{*}}=\text{max}\left\{0,\text{min}\left(\frac{\left(\theta L_2+\left(1-\theta \right)L_3-L_1\right)\alpha S+\left(\delta L_2+\left(1-\delta \right)L_3-L_5\right)\epsilon \alpha V}{W_1}\right)\right\},$

$u_2^{\text{*}}=\text{max}\left\{0,\text{min}\left(\frac{\left(L_2-L_4\right)P_C}{W_3}\right)\right\},$

$u_3^{\text{*}}=\text{max}\left\{0,\text{min}\left(\frac{\left(L_2-L_4\right)P_C}{W_3}\right)\right\}.$

Thus, the optimality system is a combination of systems (31) and (36). □

5. Numerical Investigation

In this section, we numerically investigate system (1) using MATLAB’s built-in ode45 solver. The baseline parameters from Table 1 are applied, along with initial state values of $S\left(0\right)=1000$ , $P_C\left(0\right)=20$ , $P_I\left(0\right)=50$ , $R\left(0\right)=10$ , and $V\left(0\right)=10$ . Additionally, a fourth-order Runge-Kutta scheme is employed to numerically solve the optimality system. The iterative process is initialized with weight factors $W_1=2$ , $W_2=2$ , and $W_3=6$ , as well as the costs $C_1=10$ and $C_2=15$ , and is continued until the state, adjoint, and control values converge. We propose the following intervention plans:

• Plan 1: Using all control strategies ($u_1\ne 0,u_2\ne 0,u_3\ne 0$ ).

• Plan 2: Pneumonia protection and reduction of infected individuals ($u_1\ne 0,u_2\ne 0,u_3=0$ ).

• Plan 3: Pneumonia protection and reduction of carriers ($u_1\ne 0,u_2=0,u_3\ne 0$ ).

• Plan 4: Reduction of carriers and infected individuals ($u_1=0,u_2\ne 0,u_3\ne 0$ ).

• Plan 5: Pneumonia protection only ($u_1\ne 0,u_2=0,u_3=0$ ).

• Plan 6: Carrier reduction only ($u_1=0,u_2=0,u_3\ne 0$ ).

• Plan 7: Reduction of infected individuals only ($u_1=0,u_2\ne 0,u_3=0$ ).

5.1. Role of Vaccination and Treatment in the Dynamics of Pneumonia

Figure 3. Role of vaccination and treatment on the dynamics of pneumonia.

In Figure 3, we observe the effect of treatment only ($R_0$ ), the effect of vaccination and treatment in presence of serotypes uncovered by the vaccine ($R_e$ ) and the effect of vaccination and treatment in absence of serotypes uncovered by the vaccine ($R_{ev}$ ). The result in Figure 3, tells that a combination of vaccination and treatment would be enough to significantly reduce the burden pneumonia puts on children, however, the presence of serotypes uncovered by the vaccine undermine the effort.

5.2. Effects of Serotypes Uncovered by the Vaccine on the Dynamics of Pneumonia

Figure 4. Effect of serotypes uncovered by the vaccine on pneumonia dynamics.

As shown in Figure 4, a higher prevalence of serotypes not covered by the vaccine leads to an increase in both pneumonia carriers and infected children. This highlights the impact of serotyping in pneumonia, as the presence of uncovered serotypes can reduce vaccine efficacy and contribute to a higher disease burden.

5.3. Optimal Protection and Treatment on the Dynamics of Pneumonia

Figure 5. Intervention with all control strategies.

In the presence of pneumonia protection, screening, and treatment, Figure 5 demonstrates the potential eradication of both pneumonia carriers and infected individuals in the proposed optimal control model.

Implementing Plan 2, as shown in Figure 6, leads to a significant reduction in the number of pneumonia-infected individuals. Similarly, Figure 7 illustrates a notable decline in pneumonia carriers from the onset of the epidemic, suggesting that over time, the number of infectious individuals can be drastically reduced.

Figure 6. Intervention with pneumonia protection and reduction of infected individuals.

Figure 7. Intervention with pneumonia protection and reduction of carriers.

Advocating for the identification and treatment of carriers and infected children, Figure 8 indicates the possibility of eventually eliminating the epidemic. While there is an initial increase in the number of infected individuals at the onset of the disease, this declines over time.

Figure 8. Intervention with reduction of carriers and infected individuals.

At the beginning of the pneumonia outbreak, implementing Plan 5 results in only a slight decrease in the number of infectious individuals, as seen in Figure 9. This suggests that protection alone may not be sufficient to significantly reduce pneumonia cases.

The control strategy targeting only carriers, illustrated in Figure 10, shows an immediate reduction in carriers from the start of the outbreak. Although the number of infected individuals is initially high, it gradually declines over time.

Figure 9. Intervention aimed at pneumonia protection only.

Figure 10. Intervention with carrier identification and treatment.

Figure 11. Intervention with treating infected individuals.

Figure 11 reveals that the pneumonia carriage trajectory in Plan 7 closely resembles that of the no-control scenario at the beginning of the outbreak. However, over time, the number of infected individuals decreases with the implementation of Plan 7.

Overall, all proposed control plans demonstrate effectiveness in reducing the number of infectious individuals. However, Plan 1 emerges as the most effective strategy for eradicating childhood pneumonia. Additionally, the results indicate that single intervention strategies are less effective than a combination of two or more control measures.

6. Discussion and Conclusion

In this study, we developed and analyzed a deterministic model for pneumonia transmission, which was further extended to incorporate optimal control strategies for pneumonia prevention, identification, and treatment of both carriers and infected individuals. The primary objective was to determine the most effective approach to minimizing the number of infected and carrier individuals while keeping intervention costs minimal.

Sensitivity analysis revealed that certain parameters have a direct influence on the effective reproduction number ($R_e$ ). Specifically, parameters such as the contact rate ($k$ ), the rate of vaccine non-responsiveness ($\epsilon$ ), the probability of transmission per contact ($p$ ), and additional factors ($\xi$ , $\delta$ , and $\theta$ ) contribute to an increase in pneumonia prevalence. Reducing these rates could therefore play a crucial role in controlling childhood pneumonia. Conversely, parameters such as the treatment rate of carriers ($\beta$ ), the medication rate of infected children ($\tau$ ), and the vaccination rate ($\gamma$ ) are negatively correlated with pneumonia prevalence, indicating that increasing these rates would lower $R_e$ and help curb disease spread.

Numerical simulations further demonstrated that a combination of treatment and vaccination is more effective than treatment alone. However, this strategy is challenged by the presence of children who do not respond to the vaccine, reducing its overall impact. Additionally, numerical results provided valuable insights into optimal control plans, showing that implementing a combination of all control measures would be the most effective approach to eliminating the pneumonia burden.

In conclusion, the developed model offers valuable insights into the dynamics of pneumonia transmission and the effects of various control interventions. By implementing a comprehensive control plan, the model demonstrates that it is possible to significantly reduce the prevalence of pneumonia in the population. This work underscores the importance of carefully selecting key parameters in the control of pneumonia and provides a framework for optimizing intervention strategies. Future studies could expand upon this model by incorporating factors such as cost-effectiveness, drug resistance, multi-strain pneumonia models, and more realistic stochastic population dynamics.