Cost effectiveness and optimal control analysis for bimodal pneumonia dynamics with the effect of children's breastfeeding

The global impact of exclusive versus inclusive nursing on particular baby mortalities and morbidities from conception to 6 months is examined in this study. Exclusive breastfeeding practices are more crucial and effective in preventing illness outbreaks when there is no access to appropriate medications or vaccinations. Additionally, this study takes optimal control theory into account, applying it to a system of differential equations that uses Pontryagin's Maximum Principle to describe a bimodal pneumonia transmission behavior in a vulnerable compartment. The proposed pneumonia transmission model was then updated to include two control variables. These include preventing illness exposure in susceptible children through various preventative measures and treating infected children through antibiotics, hospital care, and other treatments. If the threshold number (ℜ₀) is less than one, then treatment and prevention rates are increased, and the disease will be wiped out of the population. However, when (ℜ₀) is greater than one, then the disease persists in the population, which indicates that prevention and treatment rates are low. To evaluate the cost-effectiveness of all potential control techniques and their combinations, the incremental cost-effectiveness ratio (ICER) was determined. The simulation results of the identified model show that the interventions of prevention and treatment scenarios were the most successful in eradicating the dynamics of the pneumonia disease's propagation during the epidemic, but they were ineffective from a cost-saving perspective. Therefore, limiting pneumonia transmission to prevention alone during an outbreak is the most economical course of action.

1. Introduction

Infant (child) disability and death are the primary continuing public health issues worldwide. However, newborn (child) mortality and morbidity rates are greatly impacted by deaths caused by infectious illnesses. Infectious diseases can take 7 from 10 childhood deaths throughout the world. Pneumonia is one of the most common causes of death worldwide among acute respiratory infections, accounting for 30% of all child fatalities. Ninety-five percent of cases of pneumonia occur in developing countries. As a result, infectious illnesses are more likely to kill newborn babies in these countries [1, 2].

Among acute respiratory infection (ARI) diseases, pneumonia is the one that affects children's lungs. Approximately 740,180 children aged 05 years died because of pneumonia in 2019, which accounts for 14 and 22% of all deaths of children below 5 years and 1–5 year(s) old, respectively, and deaths are higher in Asia and Africa [3]. Hence, of infectious diseases, pneumonia causes the most children's deaths worldwide [4]. Pneumonia can be caused either by viruses, bacteria, or fungi; among these, bacterial pneumonia is the leading cause of death for children under 6 years of age. By immunizing against the disease, providing appropriate nutrition (EBF), and decreasing environmental variables, pneumonia can be avoided [5]. Additionally, using many control measures, such as prevention, treatment, and reducing indoor air pollution, can halt the spread of pneumonia. The following research has been carried out to address non-exclusive EBF or a lack of EBF, one of the main risk factors for infectious illnesses.

The first natural diet for infants is their mother's milk, which contains all the nutrients and energy required for a baby throughout the first 6 months of life [6]. According to WHO recommendations, newborns should receive only breast milk for the first 6 months of their lives. Thereafter, additional (complementary) foods are allowed for 18 months or more, followed by breastfeeding. Hence, infants (children) can achieve good growth and development [7]. Therefore, for children in the first months of life up to 6 months, any additional food or liquid (even water) is not permitted except vitamins, mineral supplements, and medicine[7, 8].

An intervention of double control has been offered to help eliminate the mortality and disability rates among children because of infectious diseases. Breastfeeding is one of the most popular and cost-effective strategies (interventions) for preventing pediatric pneumonia and all other causes of death [9–11]. Furthermore, the WHO, UNICEF, AAP, AAFP, and NNPE advocate starting breastfeeding promptly within the first hour after birth and continuing to exclusively breastfeed with human milk for the following 6 months to reduce the baby (child) death and disability rate. Continual breastfeeding with other appropriate foods will follow for the first 2 years of life to ensure that the children have healthy optimal growth and development [2].

Most of the studies assure that over two-thirds of the deaths occurring globally in the first year of life of children are often associated with a loss of exclusive breastfeeding or inappropriate feeding exercises [10]. Sub-optimal breastfeeding contributes to 18% of acute respiratory disease deaths among children under 5 years old in low-income countries [6].

Evidence suggests that if the EBF length is properly maintained, it can significantly increase immunity and lower the risk of death and disability from communicable and non-communicable diseases in both the early and advanced phases [12, 13]. EBF throughout the first 6 months of a baby's (or child's) life can typically lower the likelihood of developing any infectious diseases [14]. For the first 6 months of their lives, infants (children) who were nursed exclusively had a higher risk of contracting infectious diseases than those who were not [9, 15].

According to [10, 16], 1.24 million or 96% of child deaths occur during the first 6 months of life due to inappropriate EBF practices, and the mortality rate is higher in Africa and Asia. Additionally, poor breastfeeding results in more than 236,000 child deaths annually in a select few nations, including Nigeria, China, Mexico, Indonesia, and India [17]. Furthermore, in low- and middle-income nations, inadequate breastfeeding was found to be responsible for 18 and 30% of acute respiratory and diarrheal mortalities, respectively [18]. To reduce child mortality among children under the age of five, the WHO advises that an EBF of 90% is needed globally. Furthermore, the Sustainable Development Goals (SDGs) plan envisaged an increase in EBF of 50% by 2025 [19, 20]. According to the study by [12, 20], raising the EBF rate in middle-income and developing nations to an ideal level can reduce infant mortality among children under the age of five by 13 to 15%.

Mathematical models are frequently used to (i) analyse the dynamics of the spread of infectious diseases like cholera, bronchiolitis, pneumonia, and others; (ii) employ a variety of control methods to reduce or stop the spread of infectious diseases; and (iii) predict the effects of these diseases on people's lives, socio-economic systems, and national health programmes and policies. However, none of the aforementioned studies take into account a mathematical model method to illustrate the transmission behavior of infectious diseases, particularly pneumonia.

Several mathematical modeling studies have been conducted to estimate the potential burden of the endemic and the various control approaches for the endemic disease of pneumonia in children. Tilahun et al. [21] considered a non-linear deterministic model for the transmission of the pneumonia disease in a population of variable size, together with optimal control and cost-effectiveness measures. Agusto et al. [22] studied the advantage of isolation strategies and quarantine effectiveness measures against outbreaks of disease in the absence of appropriate medicines or vaccines.

Swai et al. [23] formulated an optimal control of pneumonia transmission in two strains by incorporating drug resistance. Additionally, how measures such as vaccination, public awareness campaigns, and therapy can reduce pneumonia transmission patterns should be considered. Tessema et al. [24] also developed a deterministic mathematical model of drug-resistant pneumonia with ideal preventive measures and cost-effectiveness evaluations. Based on the simulation values of optimal controls for the proposed model, they concluded that the combination of prevention, treatment, and screening of infectious persons is the most efficient and cost-effective way to remove pneumonia infections from the community. The diagnostic problem of distinguishing between bacterial and non-bacterial pneumonia is the main reason antibiotics are used to treat pneumonia in children. Consequently, Wu et al. [25] present causal Bayesian networks (BNs) in their model as useful tools for resolving this problem because they provide succinct maps of the probabilistic relationships between variables and produce results in a way that is understandable and justified by incorporating domain expert knowledge and numerical data.

Kotola and Mekonnen [26] created a deterministic model to demonstrate the efficacy of interventions for pneumonia and meningitis co-infection and provide a reasoned recommendation to public health officials, decision-makers in government policy, and programme implementers. Owing to their shared clinical characteristics and significant effects on human morbidity and mortality, pneumonia and tuberculosis are two of the most frequent airborne infections. Therefore, in a community of populations with both diseases, co-infection of the two diseases becomes inevitable. Owing to a lack of resources, the significant illness burden that these endemics together impose necessitates an efficient intervention to mitigate the impact. Thus, the authors in Gweryina et al. [27] use a pragmatic approach to create an SEIR model for the co-dynamics of tuberculosis and pneumonia. Using a variety of parameters, Naveed et al. [28] investigated the dynamics of delayed pneumonia-like infectious illnesses. Kassa et al. [29] and Rafiq et al. [30] offer a mathematical model of COVID-19 that includes bimodal virus transmission in a susceptible compartment.

Until now, only Legesse et al. [31] formulated a S₁S₂CIR deterministic mathematical model by grouping susceptible children as inclusively and exclusively breastfeed children and verify that inclusive breastfeeding children are more exposed to pneumonia than those children breastfeed exclusively. However, they did not take into account optimal control analysis in their research. Furthermore, no research has been carried out so far to assess the impact of EBF practice on child mortality rates and the efficacy of EBF practice in lowering pediatric mortality due to infectious disease (pneumonia). With this as a backdrop, the study's objective is to apply mathematical models with optimal control and accessible methods to treat pneumonia in infants between the ages of 0 and 6 months who do not participate in EBF. By increasing the prevalence of EBF and stepping up efforts to reduce non-exclusive breastfeeding, the findings of this study will help in making decisions that will reduce child mortality and impairment from pneumonia.

The article is organized as follows. The proposed model is formulated in the Construction of a Bimodal Pneumonia Model section and its analysis is presented in the Analyzing the Model Qualitatively section. Stability analysis of the equilibria is then discussed in the Equilibrium Point Stability section. Extension of the proposed model into optimal control is presented in The Proposed Model Under Optimal Control section. Numerical simulations are performed to support the analytical results discussed in the Analyzing the Model Qualitatively section and are presented in the Results and Discussion. Cost-effective analysis is performed in the subsequent section followed, finally, by the Conclusion.

2. Construction of the bimodal pneumonia model

In this model, the overall population size N(t) is divided into five mutually exclusive compartments based on the disease condition of the population as a whole. Furthermore, the total population size N(t) at any given time t is given by:

$\begin{array}{ll}N(t)=S_I(t)+S_E(t)+E(t)+I(t)+R(t)& (1)\end{array}$

At any time instant t ∈ [0, ∞), the real valued differentiable state variables S_I(t),S_E(t),E(t),I(t), and R(t) represent the number of susceptible children that are not exclusively breastfed, susceptible children that are exclusively breastfed, children exposed to the disease, children that are seriously infected, and children who have obtained temporary immunity from pneumonia, respectively. This research assumes that the two susceptible classes S_I(t) and S_E(t) are enlisted into the population at rates of Λ₁ and Λ₂, respectively. They acquire pneumonia infection through effective contact with the infected humans I(t) or via inhalation of contaminated air droplets at a force of infection given by

$f_i=\frac{{\beta }_{i}I}{N}$, where i = 1, 2 and β_j = kP_j, where j = 1, 2.

Here β_j = kP_j for j = 1, 2 denotes the transmission rates. However, k stands for the number of contacts, and P_j is the probability of close contact rates between two susceptible humans with the infected individuals causing infection.

Humans exposed to pneumonia advance at a γ rate to the infected compartment I(t). The sub-populations are all reduced at the same time because a consistent natural mortality rate of μ is taken into account for each compartment. The parameters σ and α at the infected stage indicate the mortality rate from pneumonia disease, which only falls in the infected class, and the percentage of children who recover due to therapy or innate immunity, respectively. Those individuals that have recovered from pneumonia are assumed to have partial immunity and again become susceptible at a rate of δ. This study also assumes that a child who has obtained partial immunity does not again join exclusively breastfed children because as one individual is infected with infectious diseases they cannot regain their original immunity [31]. Using the parameter values, basic model assumptions, and state variables described above, we have generated a systematic diagram (Figure 1), and the corresponding model equation is given by Equation (2).

$\begin{array}{ll}\{\begin{array}{l}\frac{d{S}_I}{dt}={\Lambda }_1+\delta R-f_1{S}_I-\mu S_I\\ \frac{d{S}_E}{dt}={\Lambda }_2-f_2{S}_E-\mu S_E\\ \frac{dE}{dt}=f_1{S}_I+f_2{S}_E-(\gamma +\mu )E\\ \frac{dI}{dt}=\gamma E-(\sigma +\alpha +\mu )I\\ \frac{dR}{dt}=\sigma I-(\mu +\delta )R\end{array}& (2)\end{array}$

With the following initial conditions:

$\begin{array}{ll}{S}_I(0)\ge 0,S_E(0)\ge 0,E(0)\ge 0,I(0)\ge 0,R(0)\ge 0& (3)\end{array}$

FIGURE 1

Figure 1. Model flow chart.

3. Analyzing the model qualitatively

This subsection explains the qualitative behavior of the model being considered for the long run.

3.1. Positivity and boundedness of solution

To ensure that the generated dynamical system's (2) positivity of solution is both epidemiologically meaningful and theoretically well-posed, we must show that all the state variables of the dynamical systems are non-negative.

Theorem 3.1. All the solutions of Equation (2) with the positive initial condition given on Equation (3) are non-negative.

Proof. From Equation (3), all the state variables are positive or zero at the initial time, then T > 0. To show the positivity of all the state variables select any equation of Equation (2), randomly let it be

$\begin{array}{l}\frac{d{S}_I}{dt}={\Lambda }_1+\delta R-f_1{S}_I-\mu S_I\end{array}$

$\begin{array}{l}\frac{d{S}_I}{dt}+(f_1+\mu )S_I={\Lambda }_1+\delta R\end{array}$

$\begin{array}{ll}\frac{d}{dt}[e^{{\int }_0^T(f_1+\mu )d{t}^{\prime }}{S}_I]=(e^{{\int }_0^T(f_1+\mu )d{t}^{\prime }})\left[{\Lambda }_1+\delta R\right]& (4)\end{array}$

where t′ ∈ [0, T] and each state variable are non-negative at t′.

Equation (4) is integrated with regard to time to produce

$\begin{array}{ll}{S}_I(t)=k_1{S}_I(0)+k_1\left[{\displaystyle {\int }_0^T}(e^{{\int }_0^T(f_1+\mu )d{t}^{\prime }})\left[{\Lambda }_1+\delta R\right]dt\right]\ge 0& (5)\end{array}$

where $k_1=e^{-{\int }_0^T(f_1+\mu )d{t}^{\prime }}$. From Equation (4), we observe that S_I(t) is non-negative for all t > 0. In a similar fashion, one can show S_E(t) ≥ 0, E(t) ≥ 0, I(t) ≥ 0 and R(t) ≥ 0.

Theorem 3.2. The closed positive invariant set Ω is a biologically and mathematically well-posed region of the initial value problems defined on Equations (2), (3), where

$\begin{array}{l}\Omega =\{(S_I,S_E,C,E,R)\in R_{+}^5:\\ 0<N(S_I,S_E,E,I,R)\le \frac{{\Lambda }_1+{\Lambda }_2}{\mu }\}& (6)\end{array}$

Proof. For convenience, we let S₁ = S_I, S₂ = S_E, r₁ = γ + μ, r₂ = σ + α + μ, r₃ = μ + δ throughout this study. Differentiating Equation (1) with respect to t gives

$\begin{array}{l}\frac{dN}{dt}=\frac{d{S}_I}{dt}+\frac{d{S}_E}{dt}+\frac{dE}{dt}+\frac{dI}{dt}+\frac{dR}{dt}\end{array}$

$\begin{array}{l}\frac{dN}{dt}={\Lambda }_1+{\Lambda }_2-\mu (S_I+S_E+E+I+R)-\alpha I\end{array}$

$\begin{array}{ll}\frac{dN}{dt}={\Lambda }_1+{\Lambda }_2-\mu N-\alpha I& (7)\end{array}$

In the absence of infectious rate Equation (7) reduced to

$\begin{array}{ll}\frac{dN}{dt}\le {\Lambda }_1+{\Lambda }_2-\mu N.& (8)\end{array}$

Integrating both sides of Equation (8) with regard to t and taking the limit of Equation (8) as t → ∞, we obtain

$\begin{array}{ll}N(t)\le \frac{{\Lambda }_1+{\Lambda }_2}{\mu }-\frac{{\Lambda }_1+{\Lambda }_2-\mu N_0}{\mu }{e}^{-\mu t}& (9)\end{array}$

$\begin{array}{ll}N(t)\le \frac{{\Lambda }_1+{\Lambda }_2}{\mu }.& (10)\end{array}$

Therefore, each solution of the initial value problems on Equations (2) and (3) remains in Equation (6) for all t > 0. This result can be summarized as lemma below.

Lemma 3.1. Ω is a positively invariant region for the Equation (2) with initial condition Equation (3) in $R_{+}^5$.

3.2. Threshold parameter

Before calculating the expression for threshold quantity (ℜ₀), determine pneumonia free of Equation (2). For this aim, equate the right hand side of Equation (2) to zero. After that, substitute in S₁(t) = S_I(0) > 0, S_E(t) = S_E(0) > 0, E(t) = E₀ = I(t) = I₀ and R(t) = R₀ = 0. Thus,

$\begin{array}{ll}{E}_0=\left(\frac{{\Lambda }_1}{\mu },\frac{{\Lambda }_2}{\mu },0,0,0\right).& (11)\end{array}$

Hence, E₀ is the pneumonia free-equilibrium of Equation (2). By using DFE we can account the threshold number (ℜ₀) following the work in Agusto [22], and we used the method of next generation matrix to obtain the required threshold number and from the transmission matrix

$\begin{array}{l}F=DF(E_0)=\left[\begin{array}{cc}\frac{\partial F_1(E_0)}{\partial E}& \frac{\partial F_1(E_0)}{\partial I}\\ \frac{\partial F_2(E_0)}{\partial E}& \frac{\partial F_2(E_0)}{\partial I}\end{array}\right].\end{array}$

Where F₁(t) = f₁S_I + f₂S_E and F₂(t) = γE

$\begin{array}{ll}F=\left[\begin{array}{cc}0& \frac{{\Lambda }_1{\beta }_1+{\Lambda }_2{\beta }_2}{{\Lambda }_1+{\Lambda }_2}\\ \gamma & 0\end{array}\right]& (12)\end{array}$

and the transition matrix V is given by

$\begin{array}{l}V=DF(E_0)=\left[\begin{array}{cc}\frac{\partial V_1(E_0)}{\partial E}& \frac{\partial V_1(E_0)}{\partial I}\\ \frac{\partial V_2(E_0)}{\partial E}& \frac{\partial V_2(E_0)}{\partial I}\end{array}\right].\end{array}$

where V₁ = r₁E and V₂ = r₂I.

$\begin{array}{ll}\Rightarrow V=\left[\begin{array}{cc}{r}_1& 0\\ 0& r_2\end{array}\right]& (13)\end{array}$

Hence, using the next generation matrix calculated from Equations (12), (13) we get

$\begin{array}{ll}F{V}^{-1}=\frac{1}{r_1{r}_2}\left[\begin{array}{cc}0& r_1\frac{{\beta }_1{\Lambda }_1+{\beta }_2{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}\\ \gamma r_2& 0\end{array}\right]& (14)\end{array}$

Now, the governing eigenvalue of Equation (14) represents ℜ₀ of Equation (2), which is given by

$\begin{array}{ll}{\Re }_0=\sqrt[]{\frac{\gamma ({\beta }_1{\Lambda }_1+{\beta }_2{\Lambda }_2)}{r_1{r}_2({\Lambda }_1+{\Lambda }_2)}}& (15)\end{array}$

The threshold number (ℜ₀) is a quantity that determines how pneumonia spreads within the population or fades out of the society. If ℜ₀ < 1 then the disease will fade out of the community. This shows that more exclusivity in breastfeeding children is added to the susceptible class. Because exclusively breastfed individuals have high natural immunity, they are less exposed to the diseases. ℜ₀ > 1 shows that there is a continuation of disease spread within the population.

3.3. Existence of the model's endemic equilibrium point

In this part, we examine the condition known as EE of Equation (2). The fundamental motivation for this equilibrium is that it is utilized to estimate how long pneumonia will continue to affect the population. To identify the prerequisites for an equilibrium in which community pneumonia is endemic (that is, at least one of E^* ≠ 0 or I^* ≠ 0), denoted by $E_e=(S_I^{*},S_E^{*},E^{*},I^{*},R^{*})$. To find E_e, equate each equation in Equation (2) to zero and express each state variable in terms of the force of infection at the steady state ($f_i^{*}$ where i=1,2), given by

$\begin{array}{l}{f}_1^{*}=\frac{{\beta }_{1}I}{N^{*}},f_2^{*}=\frac{{\beta }_{2}I}{N^{*}}\\ \begin{array}{l}{S}_I^{*}=\frac{{\Lambda }_1{r}_3+\delta \sigma I^{*}}{r_3(f_1^{*}-\mu )},\\ S_E^{*}=\frac{{\Lambda }_2}{f_1^{*}-\mu }=\frac{{\Lambda }_2({\Lambda }_1+{\Lambda }_2)}{\mu ({\beta }_2{I}^{*}-({\Lambda }_1+{\Lambda }_2))},\\ E^{*}=\frac{r_2{I}^{*}}{\gamma }{\Re }_0^2,\\ I^{*}=\frac{\gamma E^{*}}{r_2},\\ R^{*}=\frac{\sigma }{r_3}{\Re }_0^2\end{array}& (16)\end{array}$

Therefore, the existence of E_e in Equation (16) depends on ℜ₀, meaning that E_e from Equation (2) exists if ℜ₀ > 1.

4. Equilibrium point stability analysis

The two equilibria of Equation (2) are shown in this subsection to have both local and global asymptotic stability. We employ the Jacobian matrices of system on Equation (2) at DFE and EE for local stability and the Lyapunov function for the global stability of both equilibria to confirm this stability.

4.1. Local stability analyses

Theorem 4.1. The disease free-equilibrium point (E₀), of Equation (2) corresponding to the considered model is locally asymptotically stable if ℜ₀ < 1 and not stable otherwise.

Proof. To prove this, first determine the Jacobian matrix evaluated at E₀ becomes

$\begin{array}{ll}J(E_o)=\left[\begin{array}{ccccc}-\mu & 0& 0& -\frac{{\beta }_1{\Lambda }_1}{{\Lambda }_1+{\Lambda }_2}& \delta \\ 0& -\mu & 0& -\frac{{\beta }_2{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}& 0\\ 0& 0& -r_1& \frac{{\beta }_1{\Lambda }_1+{\beta }_2{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}& 0\\ 0& 0& \gamma & -r_2& 0\\ 0& 0& 0& \sigma & -r_3\end{array}\right]& (17)\end{array}$

The characteristic polynomial of Equation (17) becomes

$\begin{array}{ll}\text{Ψ}(\lambda )={(\lambda +\mu )}^2(\lambda +r_3)({\lambda }^2+D_1\lambda +D_2)& (18)\end{array}$

The first three eigenvalues of Equation (18) are λ = −μ a double root, λ = −r₃. All are negative, and we use the RouthHurwitz criterion to confirm the presence of the remaining eigenvalues in the manner described below:

$\begin{array}{l}{D}_1=r_1+r_2>0,\end{array}$

$\begin{array}{l}{D}_2=r_1{r}_2-\gamma \frac{{\beta }_1{\Lambda }_2+{\beta }_2\alpha }{{\Lambda }_1+{\Lambda }_2}=r_1{r}_2(1-{\Re }_0^2)>1\end{array}$

As a result, the RouthHurwitz criteria's required condition is confirmed whenever ℜ₀ < 1. Therefore, the DFE (E₀) of Equation (2) is locally asymptotically stable (LAS) when ℜ₀ < 1.

Theorem 4.2. The disease endemic equilibrium point (E_e), of Equation (2) is LAS in Ω if ℜ₀ > 1 and unstable otherwise.

Proof. To prove the local stability of E_e, first determine the desired Jacobean matrix J(E_e) of system (2) at the endemic equilibrium, which is given as Equation (19)

$\begin{array}{ll}J(E_e)=\left[\begin{array}{ccccc}-(f_1^{*}+\mu )& 0& 0& -\frac{{\beta }_1{S}_I^{*}}{N^{*}}& \delta \\ 0& -(f_2^{*}+\mu )& 0& -\frac{{\beta }_2{S}_E^{*}}{N^{*}}& 0\\ f_1^{*}& f_2^{*}& -r_1& \frac{{\beta }_1{S}_I^{*}+{\beta }_2{S}_E^{*}}{N^{*}}& 0\\ 0& 0& \gamma & -r_2& 0\\ 0& 0& 0& \sigma & -r_3\end{array}\right]& (19)\end{array}$

The characteristics polynomial corresponding to Equation (19) is

$\begin{array}{ll}(\lambda +(f_1^{*}+\mu ))(\lambda +(f_2^{*}+\mu ))(\lambda +r_1)(\lambda +r_2)(\lambda +r_3)=0& (20)\end{array}$

The first three roots of Equation (20) are λ = −r₁ < 0, λ = −r₂ < 0, and λ = −r₃ < 0 and the remaining roots can be calculated from

$\begin{array}{l}{\lambda }^2+a_1\lambda +a_2\end{array}$

where

$\begin{array}{l}{a}_1=f_1^{*}+f_2^{*}+2\mu and{a}_2=f_1^{*}{f}_2^{*}+2\mu (f_1^{*}+f_2^{*})+{\mu }^2\end{array}$

and $f_1^{*},f_2^{*}$ are defined as the force of infection at the endemic equilibrium.

As ${\lambda }^2+a_1\lambda +a_2$ has both roots with a negative real part (and the system with characteristic equation $P(\lambda )={\lambda }^2+a_1\lambda +a_2=0$ is stable) if and only if a₁, a₂ > 0, clearly a₁, a₂ > 0. Hence by RouthHurwitz criteria, for ℜ₀ > 1, the endemic equilibrium (E_e) is LAS.

4.2. Global stability analysis

In this section, we use LaSalle's invariant principle to analyse the global stability of both equilibria of Equation (2) by creating suitable Lyapunov functions.

Theorem 4.3. If ℜ₀ < 1, then the disease free-equilibrium (E₀) of Equation (2) is GAS in Ω and unstable otherwise.

Proof. We first create a suitable Lyapunov function of the type

$\begin{array}{ll}L(t)=k_{1}E(t)+k_{2}I(t)& (21)\end{array}$

where k_i, i = 1, 2 are positive real numbers to be chosen later. Upon differentiating Equation (21) along its trajectories with respect to t and simplifying, the result yields

$\begin{array}{l}\frac{dL}{dt}=k_1\frac{dE}{dt}+k_2\frac{dI}{dt},\end{array}$

$\begin{array}{ll}\frac{dL}{dt}=k_1(f_1{S}_I+f_1{S}_E-r_{1}E)+k_2(\gamma E-r_{2}I)& (22)\end{array}$

Now, we choose k₁ = γ and k₂ = r₁, and simplification of Equation (22) yields

$\begin{array}{l}\frac{dL}{dt}=\gamma (f_1{S}_I+f_1{S}_E)-r_1{r}_{2}I,\\ \frac{dL}{dt}=\left[\gamma \left(\frac{{\beta }_1{\Lambda }_1+{\beta }_2{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}\right)-r_1{r}_1\right]I& (23)\end{array}$

Simplification and some rearrangement of Equation (23) will give:

$\begin{array}{ll}\frac{dL}{dt}=-r_1{r}_2(1-{\Re }_0^2)I& (24)\end{array}$

Thus, $\frac{dL}{dt}<0$ whenever ℜ₀ < 1. Additionally, $\frac{dL}{dt}=0$ if and only if E(t) = 0 and I(t) = 0. Hence, the largest compact invariant set $\{(S_I,S_E,E,I,R)\in \Omega :\frac{dL}{dt}=0\}$ is the singleton E₀, which is the disease-free equilibrium. Therefore, using LaSalle's invariant principle [32], we conclude that the point E₀ is globally asymptotically stable in Ω if ℜ₀ < 1.

Theorem 4.4. The disease endemic equilibrium point (E_e) of Equation (2) is GAS in the invariant region stated in Theorem 3.2 as Ω if ℜ₀ > 1.

Proof. To prove the global behavior of E_e, we systematically construct a Lyapunov function V of the form as Legesse et al. [31]

$\begin{array}{ll}V(x_i)={\displaystyle \sum _1^n}\left(x_i-x_i^{*}-x_i^{*}ln\left(\frac{x_i}{x_i^{*}}\right)\right)& (25)\end{array}$

where x_i represents the compartments in the model and i = 1,...5 and $x_i^{*}$ is the endemic equilibrium point. This is defined as

$\begin{array}{l}\begin{array}{l}V\left(S_I^{*},S_E^{*},C^{*},I^{*},R^{*}\right)=\left(S_I-S_I^{*}-S_I^{*}ln\left(\frac{S_I}{S_I^{*}}\right)\right)+\\ \left(S_E-S_E^{*}-S_E^{*}ln\left(\frac{S_E}{S_E^{*}}\right)\right)+\\ +\left(E-E^{*}-E^{*}ln\left(\frac{E}{E^{*}}\right)\right)+\\ \left(I-I^{*}-I^{*}ln\left(\frac{I}{I^{*}}\right)\right)+\left(R-R^{*}-R^{*}ln\left(\frac{R}{R^{*}}\right)\right)\end{array}\end{array}$

Then, after differentiating V with regard to time t, the following is obtained.

$\begin{array}{l}\frac{dV}{dt}=\left(1-\frac{S_I^{*}}{S_I}\right)\frac{d{S}_I}{dt}+\left(1-\frac{S_E^{*}}{S_2}\right)\frac{d{S}_E}{dt}+\left(1-\frac{E^{*}}{E}\right)\frac{dE}{dt}+\\ \left(1-\frac{I^{*}}{I}\right)\frac{dI}{dt}+\left(1-\frac{R^{*}}{R}\right)\frac{dR}{dt}& (26)\end{array}$

Next, substituting $\frac{d{S}_I}{dt},\frac{d{S}_E}{dt},\frac{dE}{dt},\frac{dI}{dt},\frac{dR}{dt}$ in Equation (26) using Equation (2) gives

$\begin{array}{l}\begin{array}{l}\frac{dV}{dt}=\left(1-\frac{S_I^{*}}{S_I}\right)\left({\Lambda }_1+\delta R-f_1{S}_I-\mu S_I\right)+\left(1-\frac{S_E^{*}}{S_E}\right)\\ \left({\Lambda }_2-f_2{S}_E-\mu S_E\right)+\left(1-\frac{E^{*}}{E}\right)\left(f_1{S}_I+f_2{S}_E-\left(\gamma +\mu \right)E\right)+\\ \left(1-\frac{I^{*}}{I}\right)\left(\gamma E-\left(\sigma +\alpha +\mu \right)I\right)+\left(1-\frac{R^{*}}{R}\right)\left(\sigma I-\left(\mu +\delta \right)R\right)\end{array}\\ \begin{array}{l}=\left(\frac{S_I-S_I^{*}}{S_I}\right)\left({\Lambda }_1+\delta R-\left(f_1+\mu \right)\left(S_I-S_I^{*}\right)-\left(f_1+\mu \right)S_I^{*}\right)\\ +\left(\frac{S_E-S_E^{*}}{S_E}\right)\left({\Lambda }_2-\left(f_2+\mu \right)\left(S_E-S_E^{*}\right)-\left(f_2+\mu \right)S_E^{*}\right)+\\ \left(\frac{E-E^{*}}{E}\right)\left(f_1{S}_I+f_2{S}_E-\left(\gamma +\mu \right)\left(E-E^{*}\right)-\left(\gamma +\mu \right)E^{*}\right)\\ +\left(\frac{I-I^{*}}{I}\right)\left(\gamma E-\left(\sigma +\alpha +\mu \right)\left(I-I^{*}\right)-\left(\sigma +\alpha +\mu \right)I^{*}\right)\\ +\left(\frac{R-R^{*}}{R}\right)\left(\sigma I-\left(\mu +\delta \right)\left(R-R^{*}\right)-\left(\mu +\delta \right)R^{*}\right).\end{array}& (27)\end{array}$

We can put $\frac{dV}{dt}$ as $\frac{dV}{dt}={\text{Ψ}}_1-{\text{Ψ}}_2$ where

$\begin{array}{l}\begin{array}{l}{\text{Ψ}}_1={\Lambda }_1+{\Lambda }_2+\delta R+f_1{S}_I+\gamma E+\sigma I+(f_1+\mu )\frac{{S_I^{*}}^2}{S_I}+\\ (f_2+\mu )\frac{{S_E^{*}}^2}{S_E}+(\gamma +\mu )\frac{{E{}^{*}}^2}{E}+(\mu +\delta )\frac{{R{}^{*}}^2}{R}.\end{array}\end{array}$

$\begin{array}{l}{\text{Ψ}}_2=\frac{{(S_I-S_I^{*})}^2}{S_I}(f_1+\mu )+{\Lambda }_1\frac{S_I^{*}}{S_I}+\delta R\frac{S_I^{*}}{S_I}+(f_1+\mu )S_I^{*}+\\ \frac{{(S_E-S_E^{*})}^2}{S_E}(f_2+\mu )+\frac{S_E^{*}}{S_E}{\Lambda }_2-(f_1+\mu )S_E^{*}+\frac{f_1{S}_I{E}^{*}}{E}+f_2{S}_E+\\ \frac{(f_2{S}_E{E}^{*}}{E}-(\gamma +\mu )\frac{{(E-E^{*})}^2}{E}-(\gamma +\mu )E^{*}-\gamma E\frac{I^{*}}{I}-\\ (\sigma +\alpha +\mu )\frac{{(I-I^{*})}^2}{I}+(\sigma +\alpha +\mu )I^{*}+(\sigma +\alpha +\mu )\frac{{I{}^{*}}^2}{I}+\\ \sigma I+\sigma I\frac{R^{*}}{R}-(\mu +\delta )\frac{{(R-R^{*})}^2}{R}+(\mu +\delta )R^{*}.\end{array}$

Thus, if P < N, then $\frac{dV}{dt}\le 0$. Hence, $\frac{dV}{dt}\le 0$ when ℜ₀ > 1. Clearly, $\frac{dV}{dt}=0$ if and only if $S_I=S_I^{*},S_E=S_E^{*},E=E^{*},I=I^{*}$, and R = R^*. Therefore, the largest compact positive invariant in set $\{(S_I,S_E,E,I,R)\in \Omega :\frac{dV}{dt}=0\}$ is the singleton E_e, which is a disease endemic equilibrium of Equation (2). Generally, by LaSalle's invariant principle, E_e is GAS in the biologically feasible region when ℜ₀ > 1.

5. The proposed model under optimal control

This section focuses on using optimum control techniques with the model under consideration from Equation (2). In a short amount of time, we were able to manage or reduce the diseases in the community with the use of these strategies. The pneumonia model is expanded to include the following two control variables, each of which is defined as follows:

u₁: a campaign to prevent the spread of the disease among people who are vulnerable.

u₂: by treating infectious diseases, a treatment effort is made to minimize infection or maximize recovery.

After incorporating u₁ and u₂ in Equation (2), we obtain the following optimal control model Equation (28).

$\begin{array}{ll}\{\begin{array}{l}\frac{d{S}_I}{dt}={\Lambda }_1+\delta R-(1-u_1)f_1{S}_I-\mu S_I\\ \frac{d{S}_E}{dt}={\Lambda }_2-(1-u_1)f_2{S}_E-\mu S_E\\ \frac{dE}{dt}=(1-u_1)(f_1{S}_I+f_2{S}_E)-(\gamma +\mu )E\\ \frac{dI}{dt}=\gamma E-(\sigma +u_2)I-(\alpha +\mu )I\\ \frac{dR}{dt}=(\sigma +u_2)I-(\mu +\delta )R\end{array}& (28)\end{array}$

The control set U is Lebesgue measurable and has the following definition in order to explore the optimal levels of the controls: U = {(u₁(t), u₂(t))}: {0 ≤ u₁ < 1, 0 ≤ u₂ < 1, 0 ≤ t ≤ T} where {0 ≤ u₁ < 1, 0 ≤ u₂ < 1, 0 ≤ t ≤ T} is the set of admissible controls. Our goal is to find a control u and S_I, S_E, E, I, and R that minimize the proposed objective function J given below, while maintaining the lowest cost of control implementation in Equation (2). The proposed objective functional J should follow the epidemic Equation (2), which is given by

$\begin{array}{ll}J(u_1,u_2)={\displaystyle \underset{u_1,u_2}{min}}{\displaystyle {\int }_0^{t_f}}(b_{1}E+b_{2}I+\frac{1}{2}{\displaystyle \sum _{i=1}^2}{w}_i{u}_i^2)dt& (29)\end{array}$

subject to Equation (3), where b₁ and b₂ are the weight positive constants associated with the number of exposed children and infected children, respectively, while w₁ and w₂ are positive constants, present the relative cost weight, which is associated with control measures u₁ and u₂, respectively. We assume costs are non-linear in nature; hence, the control variables in J are in second degree polynomial form [21, 23]. The major thing that is required of us is to reduce the number of exposed and affected children while maintaining a low cost. Thus, we are going to find optimal controls $(u_1^{*},u_2^{*})$, such that

$\begin{array}{l}J(u_1^{*},u_2^{*})=min\left\{J(u_1,u_2)/u_i\in U\right\},\end{array}$

where U = (u₁, u₂): each u_i is measurable with 0 ≤ u_i < 1 ,i = 1, 2 for t ∈ [0, t_f].

5.1. The Hamiltonian and optimality system

Here, applying the principle of Pontryagin [34], Maximum Principle, we can drive the necessary conditions that the optimal control solution must satisfy [35]. Therefore, this principle converts the model Equations (28), (29) into a problem of minimizing a Hamiltonian, H, point-wise with respect to u₁ and u₂, and we obtained a Hamiltonian (H) defined as:

$\begin{array}{l}H(t,x(t),u(t),\lambda (t))=f(t,x(t),u(t))+\lambda g(t,x(t),u(t))\end{array}$

where

$\begin{array}{l}f(t,x(t),u(t))=b_{l}E+b_{2}I+\frac{1}{2}{w}_1{u}_1^2+\frac{1}{2}{w}_2{u}_2^2,\end{array}$

$\begin{array}{l}g(t,x(t),u(t))={(g_1,g_2,g_3,g_4,g_5)}^T,\end{array}$

where

$\begin{array}{l}\begin{array}{l}{g}_1={\Lambda }_1+\delta R-(1-u_1)f_1{S}_I-\mu S_I,\\ g_2={\Lambda }_2-(1-u_1)f_2{S}_E-\mu S_E,\\ g_3=(1-u_1)(f_1{S}_I+f_2{S}_E)-(\gamma +\mu )E,\\ g_4=\gamma E-(r_2+u_2)I,\\ g_5=(\sigma +u_2)I-(\mu +\delta )R.\end{array}.\end{array}$

Hence the Hamiltonian becomes

$\begin{array}{c}H(S_I,S_E,E,I,R,t)=f(E,I,u_1,u_2,t)+{\lambda }_1\frac{d{S}_I}{dt}+{\lambda }_2\frac{d{S}_E}{dt}+{\lambda }_3\frac{dE}{dt}+\\ {\lambda }_4\frac{dI}{dt}+{\lambda }_5\frac{dR}{dt}\end{array}$

$\begin{array}{l}H(S_I,S_E,E,I,R,t)=f(E,I,u_1,u_2,t)+{\lambda }_1{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3+{\lambda }_4{g}_4+{\lambda }_5{g}_5\end{array}$

$\begin{array}{ll}H=b_{l}E+b_{2}I+\frac{1}{2}{w}_1{u}_1^2+\frac{1}{2}{w}_2{u}_2^2+{\lambda }_1{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3+{\lambda }_4{g}_4+{\lambda }_5{g}_5& (30)\end{array}$

where $f(E,I,u_1,u_2,t)=b_{1}E+b_{2}I+\frac{1}{2}{\sum }_{i=1}^2{w}_i{u}_i^2,{\lambda }_i,i=1,2$ are the adjoint variable functions which are determined by using Pontryagin's maximal principle [34] and use Swai et al. [23] for verification of existence of the optimal control pairs.

Theorem 5.1. There exists adjoint variable λ_i, where i = 1, ..., 5 with transversality conditions λ_i(t_f) = 0, i = 1, ..., 5 for an optimal control $(u_1^{*},u_2^{*})$ that minimizes J(u₁, u₂) such that:

$\begin{array}{l}\frac{d\lambda }{dt}=-\frac{\partial H}{\partial X},\end{array}$

where $X={(S_I,S_E,E,I,R)}^T$ and $\lambda ={({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5)}^T$ λ(T) = 0 transidentality condition.

Now,

$\begin{array}{c}\frac{d{\lambda }_1}{dt}=-\frac{\partial H}{\partial S_I}=-({\lambda }_1(0-(1-u_1)f_1-\mu )+\\ {\lambda }_2(0)+{\lambda }_3((1-u_1)(f_1)+\\ {\lambda }_4(0)+{\lambda }_5(0))\\ ={\lambda }_1((1-u_1)f_1+\mu )-{\lambda }_3(1-u_1)f_1\end{array}$

$\begin{array}{c}\frac{d{\lambda }_2}{dt}=-\frac{\partial H}{\partial S_E}=-({\lambda }_2({\lambda }_1(0)-(1-u_1)f_2-\mu )+{\lambda }_3((1-u_1)(f_2)+\\ {\lambda }_4(0)+{\lambda }_5(0))\\ ={\lambda }_2((1-u_1)f_2+\mu )-{\lambda }_3(1-u)f_2\end{array}$

$\begin{array}{c}\frac{d{\lambda }_3}{dt}=-\frac{\partial H}{\partial E}=-(b_1+{\lambda }_1(0)+{\lambda }_2(0)+{\lambda }_3(-(\gamma +\mu ))+{\lambda }_4(\gamma )+{\lambda }_5(0))\\ =-b_1+{\lambda }_3(\gamma +\mu )-{\lambda }_4\gamma =-b_1+{\lambda }_3{\gamma }_1-{\lambda }_4\gamma \end{array}$

$\begin{array}{c}\begin{array}{c}\frac{d{\lambda }_4}{dt}=-\frac{\partial H}{\partial I}=-(b_2-\frac{{\lambda }_1{\beta }_1(1-u_1)S_I}{N}-\frac{{\lambda }_2{\beta }_2(1-u_1)S_E}{N}+\\ \frac{{\lambda }_3(1-u_1)({\beta }_1{S}_I+{\beta }_2{S}_E)}{N}-{\lambda }_4(\sigma +u_2+\alpha +\mu )-{\lambda }_5(\sigma +u_2)\end{array}\\ \begin{array}{c}=-b_2+{\lambda }_1{\beta }_1(1-u_1)\frac{{\Lambda }_1}{{\Lambda }_1+{\Lambda }_2}+{\lambda }_2{\beta }_2(1-u_1)\frac{{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}-\\ {\lambda }_3(1-u_1)\left(\frac{{\beta }_1{\Lambda }_1+{\beta }_2{\Lambda }_2}{{\Lambda }_1+{\Lambda }_2}\right)+{\lambda }_4({\gamma }_2+u_2)-{\lambda }_5(\sigma +u_2)\end{array}\end{array}$

$\begin{array}{c}\begin{array}{c}\frac{d{\lambda }_5}{dt}=-\frac{\partial H}{\partial R}=-({\lambda }_1\sigma +{\lambda }_2(0)+{\lambda }_3(0)+{\lambda }_4(0)-{\lambda }_5(\mu +\delta ))\end{array}\\ =-({\lambda }_1\delta -{\lambda }_5(\mu +\delta ))=-{\lambda }_1\delta +{\lambda }_5{\gamma }_3\end{array}$

In a similar manner, we obtained the controls by solving the equation $\frac{\partial H}{\partial u_i}=0$ at $u_i^{*}$, for i = 1, 2 in accordance with Pontryagin [34]'s methodology and obtained:

$\begin{array}{l}{w}_1{u}_1+{\lambda }_1{f}_1{S}_I+{\lambda }_2{f}_2{S}_E-{\lambda }_3(f_1{S}_I+f_2{S}_E)=0\end{array}$

$\begin{array}{l}{u}_1=\frac{({\lambda }_3-{\lambda }_1){\beta }_1{\Lambda }_1+({\lambda }_3-{\lambda }_2){\beta }_2{\Lambda }_2}{w_1({\Lambda }_1+{\Lambda }_2)}I\end{array}$

$\begin{array}{l}{u}_1^{*}=max\left\{0,min\left\{1,\frac{({\lambda }_3-{\lambda }_1){\beta }_1{S}_I+({\lambda }_3-{\lambda }_2){\beta }_2{S}_E}{w_{1}N}I\right\}\right\}\end{array}$

Similarly $\frac{\partial H}{\partial u_2}=0$

$\begin{array}{l}{w}_2{u}_2-{\lambda }_{4}I+{\lambda }_{5}I=0\end{array}$

From this

$\begin{array}{l}{u}_2=\frac{({\lambda }_4-{\lambda }_5)I}{w_2}\end{array}$

This implies that

$\begin{array}{l}{u}_2^{*}=\{\begin{array}{l}{u}_2,if0<u_2<1\\ 0,if{u}_2<0\\ 1,if{u}_2>1\end{array}\end{array}$

The above equation in compact notation is

$\begin{array}{l}{u}_2^{*}=max\left\{0,min\left\{1,\frac{({\lambda }_4-{\lambda }_5)I}{w_2}\right\}\right\}\end{array}$

Considering the bounds of the control quintuple, we have

$\begin{array}{ll}\begin{array}{l}{u}_1^{*}=max\left\{0,min\left\{1,\frac{({\lambda }_3-{\lambda }_1){\beta }_1{\Lambda }_1+({\lambda }_3-{\lambda }_2){\beta }_2{\Lambda }_2}{w_1({\Lambda }_1+{\Lambda }_2)}I\right\}\right\}\\ u_2^{*}=max\left\{0,min\left\{1,\frac{({\lambda }_4-{\lambda }_5)I}{w_2}\right\}\right\}\end{array}& (31)\end{array}$

The optimality system is obtained from the state Equation (29) together with adjoint variables and the transversality condition in Theorem 5.1 by including the characterized control set and initial condition.

$\begin{array}{ll}\begin{array}{l}\frac{d{S}_I}{dt}={\Lambda }_1+\delta R-(1-u_1)f_1{S}_I-\mu S_I\\ \frac{d{S}_E}{dt}={\Lambda }_2-(1-u_1)f_2{S}_E-\mu S_E\\ \frac{dE}{dt}=(1-u_1)(f_1{S}_I+f_2{S}_E)-(\gamma +\mu )E\\ \frac{dI}{dt}=\gamma E-(\sigma +u_2)I-(\alpha +\mu )I\\ \frac{dR}{dt}=(\sigma +u_2)I-(\mu +\delta )R\\ \frac{d{\lambda }_1}{dt}=(1-u_1)f_1({\lambda }_1-{\lambda }_3)+{\lambda }_1\mu \\ \frac{d{\lambda }_2}{dt}=(1-u_1)f_2({\lambda }_2-{\lambda }_3)+{\lambda }_2\mu \\ \frac{d{\lambda }_3}{dt}=-b_1+{\lambda }_3{\gamma }_1-{\lambda }_4\gamma \\ \frac{d{\lambda }_4}{dt}=-b_2+(1-u_1)({\lambda }_1-{\lambda }_3)\frac{{\beta }_1{S}_I}{N}+(1-u_1)\\ ({\lambda }_2-{\lambda }_3)\frac{{\beta }_2{S}_E}{N}+{\lambda }_4(r_2+u_2)-{\lambda }_5(\sigma +u_2)\\ \frac{d{\lambda }_5}{dt}=-{\lambda }_1\delta +{\lambda }_5(\mu +\delta )\end{array}& (32)\end{array}$

$\begin{array}{r}{\lambda }_i(t_f)=0,i=1,...,5\\ S_I(0)=S_{I0},S_E(0)=S_{E0},E(0)=E_0,I(0)=I_0,R(0)=R_0\end{array}$

Therefore, using the optimality system 32, it is possible to calculate the optimal control. Consequently, the optimal problem is minimal at control $u_1^{*}$ and $u_2^{*}$, as shown by the fact that the second derivatives of the Lagrangian with regard to u₁ and u₂, respectively, are positive.

6. Results and discussion

To analyse the dynamics of pneumonia disease with or without control measures, numerical simulations are performed on the suggested model and optimality system using the parameter values indicated in Table 1. In addition, we assumed the initial population size to be S_I(0) = 40; S_E(0) = 100; E(0) = 50; I(0) = 15; and R(0) = 1 for the purpose of numerical simulation. The weight constant values are chosen as b₁ = 3; b₂ = 3; w₁ = 0.05 and w₂ = 0.03. First, we simulate the pneumonia model for the case R₀ = 0.8513 < 1, which indicates that the pneumonia disease dies out from the society. As a result, the pneumonia model's solution trajectory moves toward a disease-free equilibrium point. The disease-free equilibrium point is demonstrated to be locally asymptotically stable as all the trajectories of the model converge to DFE, see Figure 2A. Next, we plotted the graphics for the case R₀ = 1.4232 > 1, which implies that the disease is endemic. In this case, the solution curves are converging to the endemic equilibrium point, which verifies the linear stability of the EE point (see Figure 2B).

TABLE 1

Table 1. Pneumonia model parameter values with their source.

FIGURE 2

Figure 2. The pneumonia model's trajectory of solutions converges to (A) DFEP and (B) EEP.

Now, to extend the proposed model to optimal control, we focus on the parameter values and initial population, which give R₀ = 1.4232 > 1 to analyze the model. In light of the fact that diseases are still prevalent in society, adding control factors to the mode is appropriate. Figures 3–6 demonstrate the impact of prevention and treatment on the dynamics of pneumonia.

FIGURE 3

Figure 3. (A) Dynamics of sub-populations for the DFE point. (B) The phase portrait for S_I(t) and S_E(t) vs. E(t).

The plot in Figure 3A illustrates that subpopulations converge to the DFE point, which indicates that pneumonia has been eliminated from the community. Moreover, it can be observed that the two susceptible populations decrease while the exposed and infected children increase for a few years and decrease rapidly afterward to the DFE point. Figure 3B reveals that even if controls are applied, non-exclusively breastfed children are more exposed to pneumonia than exclusively breastfed children. In general, from Figures 2A, 3A, we can easily see the impact of control variables on the transmission dynamics of pneumonia.

6.1. Contingency plans

We utilized the following scenarios to assess how each regulation would affect the dynamics of pneumonia spread:

(i) Optimal use of prevention (u₁ only).

(ii) Optimal use of treatment (u₂ only).

(iii) Optimal use of prevention (u₁) and treatment (u₂) intervention.

6.1.1. Scenario A: control of pneumonia with prevention only

This scenario shows the use of only one control measure, prevention (u₁), and the other controls were set to zero. As clearly observed from Figures 4A, B, with the optimal use of a prevention strategy, the two susceptible individuals increase due to the prevention strategy, and when we compare it with the case free of prevention, the number of susceptibilities of individuals to the diseases is less. Moreover, the number of total exposed humans decreases more with control than when there is no control, as depicted in Figure 4C. Since the number of infection averted human from pneumonia disease due to this strategy is less in number, hence additional intervention is required.

FIGURE 4

Figure 4. Simulation of the optimal model showing the effect of prevention on (A) not exclusively breastfeed individuals, (B) exclusively breastfeed individuals, and (C) exposed individuals.

6.1.2. Scenario B: control of pneumonia with treatment only

Scenario B is shown in Figures 5A, B, which illustrate that treatment has a significant impact in reducing the number of children infected with pneumonia after 14 years. It can be noted that the number of infected individuals slightly decreases and becomes effective after some time; hence, more interventions are needed to eliminate the disease from the community.

FIGURE 5

Figure 5. Simulations showing (A) the optimal use of treatment only (u2) and (B) its control profile.

6.1.3. Scenario C: optimal use of the two controls

This strategy demonstrates the effect of the optimal use of prevention for the exposed humans and treatment for the infectious humans to decrease the number of exposed and infected individuals in the society. Additionally, this intervention reduces the spread of pneumonia dynamics governed by model (2) in the population. The numbers of exposed individuals and infectious individuals decrease more rapidly when the two control scenarios are in use compared with when controls are not used or one control is used, as depicted in Figures 6A, B. Figure 6F reveals that the optimal use of prevention. u₁(t) is maximum at 100% throughout the proposed days until reaching the final time that maximum prevention is applied to control pneumonia. Optimal use of treatment u₂(t) is kept at the maximum level for 48 days before arriving at the minimum at the final intervention time. Figures 6C–E reveal, respectively, the size of S_I, S_E, and R increases compared with non-control and one control intervention. This confirms that a maximum number of children's pneumonia diseases are averted due to the intervention of the two controls.

FIGURE 6

Figure 6. Simulations demonstrating optimal use of prevention (u₁) and treatment (u₂) on (A) E (B) I (C) S_E (D) S_I (E) R and (F) its control profile.

7. Cost-effectiveness analysis

In this section, we present cost-effectiveness analysis, which is used to evaluate the benefits related to a health intervention(s) or strategy (strategies) (for instance treatment and prevention), to elaborate the strategy's costs [22]. The number of infections averted is given as the difference between total infectious individuals without control and total infectious individuals with control. Using the parameter values in Table 1 and initial conditions of state variables with the weight constant values chosen, the ICER is determined for each intervention labeled as prevention, treatment, and a combination of both. The prevention strategy includes vaccination (immunization), personal hygiene, avoiding exposure to people who are ill, covering a cough, and adequate nutrition (scenario A), while the treatment intervention involves antibiotics that stop the infection from progressing (these medicines are used to treat bacterial pneumonia), hospital treatment (allowed for more severe cases), rest, etc. (scenario B). The combination of prevention and treatment scenario C. This is obtained by balancing the change between the costs and health outcomes of these intervention strategies; usually obtained by using the incremental cost-effectiveness ratio (ICER), which is described as:

$\begin{array}{ll}ICER=\frac{\text{change in total costs between strategies}}{\text{change in health benefits between strategies}}& (33)\end{array}$

where the numerator of the ICER represents the difference in cost-benefit and the denominator measures the change in health benefit.

According to the simulation outcomes of the optimality system, the control scenarios are then ranked in ascending order of total number of infections averted, i.e., prevention of infections in susceptible children using vaccines, personal hygiene and others (strategy A), treatment of infected individuals with antibiotics (strategy B), and a combination of prevention and treatment (strategy C), as shown in Table 2.

TABLE 2

Table 2. Incremental cost-effectiveness ratio in increasing order of total infections averted.

The ICER is obtained through the following computation:

$\begin{array}{r}\begin{array}{r}ICER(B)=\frac{1102.5}{701.9053}=1.5707\\ ICER(A)=\frac{346.4642-1102.5}{1665.791-701.9053}=-0.7827\\ ICER(C)=\frac{1753.1540-346.4642}{2178.746-1665.791}=2.741\end{array}\end{array}$

Now, comparing strategy A and B incrementally, the ICER for the two competing strategies is calculated as above and it shows that ICER (B) > ICER (A). From this, we can see that strategy A saves 0.7827 more than strategy B, and strategy B is a bit more expensive. Hence, we excluded strategy B from the set of competing strategies, and finally, we compared strategies A and C as depicted in Table 3. From ICER (A) and ICER (C) in Table 3 we can see that strategy C saves 2.741 than strategy A. Hence, we exclude strategy C, because it is a bit expensive. Therefore, we conclude that strategy A the cheapest of all compared strategies, that meant it is the most cost-effective for pneumonia disease control intervention strategies.

TABLE 3

Table 3. Comparison between intervention strategies A and C.

8. Conclusion

This study is concerned with the mathematical analysis of a pneumonia transmission model with naturally acquired immunity in the presence of effective exclusively breastfed infants and a lack of naturally acquired immunity due to the loss of exclusively breastfed infants. This work also shows that if the threshold number is smaller than unity, then the pneumonia-free equilibrium point is both locally and globally asymptotically stable, which means pneumonia is wiped out of the community. If the threshold number is greater than unity, then an endemic equilibrium of the model occurs, which shows the persistence of the diseases in the population.

To control pneumonia spread dynamics in a population, multiple time-dependent control variables, including prevention using vaccines, personal hygiene, etc., treatment of infectious humans using antibiotics, hospital treatment, and rest are considered. An analysis of the optimal control model is carried out theoretically, and the model is simulated to determine the effects of combining the two control intervention strategies on the spread dynamics of pneumonia in the community. It is shown that the number of infected children is minimized through prevention and treatment intervention strategies. Throughout this work, based on the results in Table 3, we recommend the prevention of susceptible children from being exposed to the diseases using vaccination, public health education, etc., to reduce new exposed cases and the number of infected children due to pneumonia in our society with the least cost.

In general, we considered cleanliness as a method of preventing pneumonia in children under the age of five in the earlier studies on the dynamics of bimodal pneumonia [31]. However, in the present study, we considered the extension of the bimodal pneumonia model to optimal control using two time-dependent control measures, namely prevention and treatment. In addition, we analyzed the cost-effectiveness of intervention strategies. The result of the analysis reveals that prevention strategies are the most cost-effective way of eradicating pneumonia. Therefore, the present study is more effective and cost-effective in preventing pneumonia transmission than the previous study.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

FL conceptualized, planned, and prepared the article as well as the figures. The final manuscript's review, modification, and literature search were all aided by the efforts of all writers. The article's submission was reviewed and approved by all authors.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Elyas L, Mekasha A, Admasie A, Assefa E. Exclusive breastfeeding practice and associated factors among mothers attending private pediatric and child clinics, Addis Ababa, Ethiopia: a cross-sectional study. Int J Pediatr. (2017) 2017:8546192. doi: 10.1155/2017/8546192

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Abdulla F, Hossain MM, Karimuzzaman M, Ali M, Rahman A. Likelihood of infectious diseases due to lack of exclusive breastfeeding among infants in Bangladesh. PLoS ONE. (2022) 17:e0263890. doi: 10.1371/journal.pone.0263890

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Morty RE. World health day observances in November 2021: advocating for adult and pediatric pneumonia, preterm birth, and chronic obstructive pulmonary disease. Am J Physiol Lung Cell Mol Physiol. (2021) 321:L954L957. doi: 10.1152/ajplung.00423.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Otoo D, Opoku P, Charles S, Kingsley AP. Deterministic epidemic model for (SVCSyCAsyIR) pneumonia dynamics, with vaccination and temporal immunity. Infect Dis Model. (2020) 5:42–60. doi: 10.1016/j.idm.2019.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kizito M, Tumwiine J, A. mathematical model of treatment and vaccination interventions of pneumococcal pneumonia infection dynamics. J Appl Mathem. (2018) 2018:2539465. doi: 10.1155/2018/2539465

CrossRef Full Text | Google Scholar

6. Alebel A, Tesma C, Temesgen B, Ferede A, Kibret GD. Exclusive breastfeeding practice in Ethiopia and its association with antenatal care and institutional delivery: a systematic review and meta-analysis. Int Breastfeed J. (2018) 13:1–12. doi: 10.1186/s13006-018-0173-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Rajeshwari K, Bang A, Chaturvedi P, Kumar V, Yadav B, Bharadva K, et al. Infant and young child feeding guidelines: 2010. Indian Pediatr. (2010) 47:995–1004.

PubMed Abstract | Google Scholar

8. Arage G, Gedamu H. Exclusive breastfeeding practice and its associated factors among mothers of infants less than six months of age in Debre Tabor town, Northwest Ethiopia: a cross-sectional study. Adv Public Health. (2016) 2016:3426249. doi: 10.1155/2016/3426249

CrossRef Full Text | Google Scholar

9. Turin CG, Ochoa TJ. The role of maternal breast milk in preventing infantile diarrhea in the developing world. Curr Trop Med Rep. (2014) 1:97–105. doi: 10.1007/s40475-014-0015-x

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tewabe T, Mandesh A, Gualu T, Alem G, Mekuria G, Zeleke H. Exclusive breastfeeding practice and associated factors among mothers in Motta town, East Gojjam zone, Amhara Regional State, Ethiopia, 2015: a cross-sectional study. Int Breastfeed J. (2016) 12:1–7. doi: 10.1186/s13006-017-0103-3

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Mgongo M, Mosha MV, Uriyo JG, Msuya SE, Stray-Pedersen B. Prevalence and predictors of exclusive breastfeeding among women in Kilimanjaro region, Northern Tanzania: a population based cross-sectional study. Int Breastfeed J. (2013) 8:1–8. doi: 10.1186/1746-4358-8-12

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Mogre V, Dery M, Gaa PK. Knowledge, attitudes and determinants of exclusive breastfeeding practice among Ghanaian rural lactating mothers. Int Breastfeed J. (2016) 11:1–8. doi: 10.1186/s13006-016-0071-z

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ajetunmobi OM, Whyte B, Chalmers J, Tappin DM, Wolfson L, Fleming M, et al. Breastfeeding is associated with reduced childhood hospitalization: evidence from a Scottish Birth Cohort (1997-2009). J Pediatr. (2015) 166:620–5. doi: 10.1016/j.jpeds.2014.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Martin CR, Ling PR, Blackburn GL. Review of infant feeding: key features of breast milk and infant formula. Nutrients. (2016) 8:279. doi: 10.3390/nu8050279

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Arifeen S, Black RE, Antelman G, Baqui A, Caulfield L, Becker S. Exclusive breastfeeding reduces acute respiratory infection and diarrhea deaths among infants in Dhaka slums. Pediatrics. (2001) 108:e67–e67. doi: 10.1542/peds.108.4.e67

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sefene A, Birhanu D, Awoke W, Taye T. Determinants of exclusive breastfeeding practice among mothers of children age less than 6 month in Bahir Dar city administration, Northwest Ethiopia; a community based cross-sectional survey. Sci J Clin Med. (2013) 2:153–9. doi: 10.11648/j.sjcm.20130206.12

CrossRef Full Text | Google Scholar

17. Collective GB UNICEF. Nurturing the health and wealth of nations: the investment case for breastfeeding. Technical Report, World Health Organization. (2017).

Google Scholar

18. Bhandari N, Chowdhury R. Infant and young child feeding. Proc Indian Nat Sci Acad. (2016) 82:1507–17. doi: 10.16943/ptinsa/2016/48883

CrossRef Full Text | Google Scholar

19. Organization WH. Tracking universal health coverage: first global monitoring report. World Health Organization. (2015).

Google Scholar

20. Sajjad S, Roshan R, Tanvir S. Impact of maternal education and source of knowledge on breast feeding practices in Rawalpindi city. MOJCRR. (2018) 1:212–42. doi: 10.15406/mojcrr.2018.01.00035

CrossRef Full Text | Google Scholar

21. Tilahun GT, Makinde OD, Malonza D. Modelling and optimal control of pneumonia disease with cost-effective strategies. J Biol Dyn. (2017) 11:400–26. doi: 10.1080/17513758.2017.1337245

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Agusto FB. Optimal isolation control strategies and cost-effectiveness analysis of a two-strain avian influenza model. Biosystems. (2013) 113:155–64. doi: 10.1016/j.biosystems.2013.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Swai MC, Shaban N, Marijani T. Optimal control in two strain pneumonia transmission dynamics. J Appl Mathem. (2021) 2021:8835918. doi: 10.1155/2021/8835918

CrossRef Full Text | Google Scholar

24. Tessema FS, Bole BK, Rao PK. Optimal control strategies and cost effectiveness analysis of Pneumonia disease with drug resistance. Int J Nonl Analy Appl. (2023) 14:903–17. doi: 10.22075/ijnaa.2022.26746.3402

CrossRef Full Text | Google Scholar

25. Wu Y, Mascaro S, Bhuiyan M, Fathima P, Mace AO, Nicol MP, et al. Predicting the causative pathogen among children with pneumonia using a causal Bayesian network. PLoS Comput Biol. (2023) 19:e1010967. doi: 10.1371/journal.pcbi.1010967

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kotola BS, Mekonnen TT. Mathematical model analysis and numerical simulation for codynamics of meningitis and pneumonia infection with intervention. Sci Rep. (2022) 12:1–22. doi: 10.1038/s41598-022-06253-0

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Gweryina RI, Madubueze CE, Bajiya VP, Esla FE. Modeling and analysis of tuberculosis and pneumonia co-infection dynamics with cost-effective strategies. Results Control Optimiz. (2023) 10:100210. doi: 10.1016/j.rico.2023.100210

CrossRef Full Text | Google Scholar

28. Naveed M, Baleanu D, Raza A, Rafiq M, Soori AH, Mohsin M. Modeling the transmission dynamics of delayed pneumonia-like diseases with a sensitivity of parameters. Adv Differ Equat. (2021) 2021:1–19. doi: 10.1186/s13662-021-03618-z

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kassa SM, Njagarah JB, Terefe YA. Analysis of the mitigation strategies for COVID-19: from mathematical modelling perspective. Chaos, Solitons Fractals. (2020) 138:109968. doi: 10.1016/j.chaos.2020.109968

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Rafiq M, Ali J, Riaz MB, Awrejcewicz J. Numerical analysis of a bi-modal COVID-19 sitr model. Alexandria Eng J. (2022) 61:227–35. doi: 10.1016/j.aej.2021.04.102

CrossRef Full Text | Google Scholar

31. Legesse FM, Rao KP, Keno TD. Mathematical Modeling of a Bimodal Pneumonia Epidemic with Non-breastfeeding Class. Appl Math. (2023) 17:95–107. doi: 10.18576/amis/170111

CrossRef Full Text | Google Scholar

32. Dano LB, Rao KP, Keno TD. Modeling the combined effect of hepatitis b infection and heavy alcohol consumption on the progression dynamics of liver cirrhosis. J Mathematics. (2022) 2022:6936396. doi: 10.1155/2022/6936396

CrossRef Full Text | Google Scholar

33. Otieno O, Joseph M, John O. “Mathematical Model for Pneumonia Dynamics among Children,” in The 2012 southern Africa mathematical sciences association conference (SAMSA 2012). (2012).

Google Scholar

34. Pontryagin LS. Mathematical Theory of Optimal Processes. London: CRC press. (1987).

Google Scholar

35. Lenhart S, Workman JT. Optimal Control Applied to Biological Models. New York: Chapman and Hall/CRC. (2007). doi: 10.1201/9781420011418

CrossRef Full Text | Google Scholar