Modelling and optimal control of pneumonia disease with cost-effective strategies

ABSTRACT We propose and analyse a nonlinear mathematical model for the transmission dynamics of pneumonia disease in a population of varying size. The deterministic compartmental model is studied using stability theory of differential equations. The effective reproduction number is obtained and also the asymptotic stability conditions for the disease free and as well as for the endemic equilibria are established. The possibility of bifurcation of the model and the sensitivity indices of the basic reproduction number to the key parameters are determined. Using Pontryagin's maximum principle, the optimal control problem is formulated with three control strategies: namely disease prevention through education, treatment and screening. The cost-effectiveness analysis of the adopted control strategies revealed that the combination of prevention and treatment is the most cost-effective intervention strategies to combat the pneumonia pandemic. Numerical simulation is performed and pertinent results are displayed graphically.


Introduction
In the report of WHO [18], 'infectious diseases are the leading cause of death in Human beings'. According to the fact sheet of WHO [18] sixteen percent of all deaths each year are from infectious diseases that means over 9.5 million deaths annually attribute to infectious diseases, with most of them in developing countries. From 9.5 million annual death, 'Pneumonia and other respiratory infections cause about 2 million child deaths yearly in developing countries' [19]. If we compare infectious diseases such as Malaria, HIV/AIDS, Measles and Pneumonia for under five-year children in Africa, pneumonia is the leading cause of deaths [19]. According to Institute for Health Metrics and Evaluation [4], every 35 seconds a child dies from pneumonia.
In Ethiopia, pneumonia is one of the leading cause of death. The reported cases shows that, it has been increasing aggregatively in the past 7 years (see Figure 1). A lot of scholars proposed models for understanding of infectious disease dynamics and also for making quantitative predictions of different intervention strategies and their effectiveness [10][11][12][13]. Very few essential research have been done on the dynamics of pneumonia have been done in the last decade. Some of them are, Melegaro et al. [8], Joseph Emaline [5], Ssebuliba, [16] and Okaka et al. [9], proposed a model on pneumonia dynamics. Additionally, Ong'ala et al. [14] studied and estimated the basic reproductive number as a random variable by first developing and analysing a deterministic model for transmission patterns of pneumonia.
All the above studies have developed a deterministic as well as stochastic mathematical model of pneumonia dynamics by subdividing the population into sub-classes of susceptible, infectious, vaccinated, treated, carrier and recovered. But none of them considered optimal control and cost-effectiveness strategies and also no study have been undertaken by applying optimal control. This, therefore motivated us to undertake this study to fulfil this gap. To estimate some parameters demographic data were collected from Health Minster of federal democratic republic of Ethiopia.

Model description and formulation
The model divides the total population into five sub-classes according to their disease status. Susceptible (S), vaccinated (V ), carrier (C), infected (I) and recovered (R). The model assumes that a fraction of the population has been vaccinated before the disease out break at the rate of (p) and (1 − p) fraction of population susceptible. (We consider this model due to the reason that, in African particularly in Ethiopian context all new born infants are not taking Pneumococcal conjugate vaccine (PCV). Only those mothers who are aware or who stay around town or city will go to their nearby health centre to vaccinate their infants but there are a lot of newborn left without vaccination).The Susceptible class is increased from vaccinated class in which those individuals who are vaccinated but did not respond to vaccination with waning rate of φ and from recovered class in which those individuals who loss their temporary immunity by δ rate. However, individuals from susceptible class move to vaccinated class with vaccination rate of ϑ. The susceptible class is infected either by carrier or symptomatically infected individuals with a force of infection λ = ξ((I(t) + ϒC(t))/N) where, ξ = kτ , k is contact rate, τ is the probability that a contact is effective to cause infection and ϒ is transmission coefficient for the carrier. If ϒ > 1 then, the carries infect susceptible more likely than infective. If ϒ = 1, then both carriers and infective have equal chance to infect the susceptible, but if ϒ < 1 then the infective have good chance to infect susceptible than carriers. The model assumes vaccination is not 100% effective, so vaccinated classes (V ) also have a chance of being infectious or carrier with small proportion and the force of infection for the vaccinated class is λ v = λ, where 0 ≤ < 1 and is the proportion of the serotype not covered by the vaccine. Newly infected individuals by the force of infection become either carrier with a probability of ρ to join the carrier class C or move to the infected class I with probability of 1 − ρ. The carrier class can develop disease symptom or can screen themselves and join the infected class with a rate of χ or recover by gaining natural immunity at β rate. Individuals in the infected class move to recovered compartment at a per capita rate of η by treatment, with treatment efficacy of q proportion of individuals join the recovered class or join the carrier class with (1 − q) proportion by adapting the treatment, or die from the disease at the rate α. In all compartments μ is the natural mortality rate of individuals and also all the parameters are positive. The above model description can be represented diagrammatically in Figure 2.
The above flow diagram can be written in to a system of five differential equations as follows: With initial condition

Invariant region
In this section, a region in which solutions of the model system are uniformly bounded is the proper subset ⊂ 5 + . The total population at any time t is given by N = S+V +C+I+R. and dN/dt = π − μN − αI(t).
In the absence of mortality due to pneumonia equation, it becomes By solving Equation (2), we obtain 0 ≤ N ≤ π/μ. Therefore, the feasible solution set of the system equation (1) of the model remain in the region:

Positivity of the solutions
In this section, to obtain the solution of the model is non-negative we stated and proved the following theorem. Proof: From the system of differential equation (1), let us take the first equation; Similarly, we obtained

Disease-free equilibrium
In this section we obtain the equilibrium point at which the epidemic is eradicated from the population. Letting the right-hand side of Equation (1) to zero and letting C = I = 0, leads to

The effective reproductive number (R eff )
In this section we obtained the threshold parameter that governs the spread of a disease which is called the effective reproduction number is determined. To obtain the effective reproduction number, we used the next-generation matrix method so that it is the spectral radius of the next-generation matrix [17] and we obtain

Local stability of disease-free equilibrium
Theorem 3.2: The disease-free equilibrium point is locally asymptotically stable if R eff < 1 and unstable if R eff > 1.

Proof:
To prove local stability of disease-free equilibrium, we obtained the Jacobian matrix of the system (1) at the disease-free equilibrium E 0 : To obtain the eigenvalue of Equation (5), when we expand Equation (6) Then by Routh-Hurwitz criteria equation (8) have strictly negative root. The determinant of Equation (7) can be obtained when we rearrange Equation (9), it becomes where By Routh-Hurwitz criteria, and also Thus, the disease-free equilibrium is locally asymptotically stable if R eff < 1.

The endemic equilibrium
The endemic equilibrium is denoted by E * and defined as a steady-state solutions for the model (1). This can occur when there is a persistence of the disease. It can be obtained by equating the system of Equation (1) to zero. Then we obtained . where is the endemic equilibrium of the model 1.

Lemma 3.3:
For R eff > 1 a unique endemic equilibrium point E * exist and no endemic equilibrium otherwise.
Proof: For the disease to endemic, dC/dt > 0 and dI/dt > 0, that is: From the second inequality of Equation (10), From the first inequality of Equation (10), From the fact (S + V)/N ≤ 1, By substituting Equation (12) into Equation (11) we can get, Then, by rearranging and cancelling of I in both sides, we can get Thus a unique endemic equilibrium exist when R eff > 1.
Using expression for I * in the endemic equilibrium, we plot Figure 10, that shows there is a backward bifurcation. We used a set of estimated parameters in Table 2. The figure reflects the co-existence of the disease free with endemic equilibrium.

Local stability of the endemic equilibrium
Theorem 3.4: If R eff > 1 then the endemic equilibrium E * of system (1) is locally asymptotically stable in .
Proof: To Prove the local stability of endemic equilibrium, we obtain the Jacobian matrix at endemic equilibrium in Equation (14): The characteristic polynomial of Equation (14) in terms of force of infection is obtained as where , For, λ is the force of infection evaluated at endemic equilibrium.
According to the Routh-Hurwitz criterion, for R eff > 1, the endemic equilibrium (E * ) is locally asymptotically stable if,

The global stability of the endemic equilibrium
Theorem 3.5: If R eff > 1, the endemic equilibrium E * of the model (1) is globally asymptotically stable.
Proof: To prove the global asymptotic stability of the endemic equilibrium, we use the method of Lyapunov functions. Define.
By direct calculating the derivative of L along the solution of Equation (1), we have Thus collecting positive terms together and negative terms together from Equation (15) leads to Thus if Q < K, then dL/dt ≤ 0; Noting that dL/dt = 0 if and only if S = S * , V = V * , C = C * , I = I * , R = R * Therefore, the largest compact invariant set in {(S * , V * , C * , I * , R * ) ∈ : dL/dt = 0} is the singleton E * , where E * is the endemic equilibrium of the system (1).
By LaSalle's invariant principle ( [6], it implies that E * is globally asymptotically stable in if Q < K.

Sensitivity analysis of the model parameters
We carried out the sensitivity analysis to determine the model robustness to parameter values. The normalized forward sensitivity index of a variable to a parameter is a ratio of the relative change in the variable to the relative change in the parameter. If a variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives.

Definition:
The normalized forward sensitivity index of a variable, u, which depends differentiability on index of a parameter, p is defined as u p = (∂u/∂p)(p/u) From an explicit formula for (R eff ) in Equation (4), we derive an analytical expression for the sensitivity of R eff as R eff p = (∂R eff /∂p)(p/R eff ) to each of the parameter involved in (R eff ). For example, the sensitivity index of R eff with respect to k is α where obtained and evaluated at, p = 0.6, φ = 0.001, ϑ = 0.9, = 0.4, χ = 0.00274, q = 0.5, η = 0.0238, β = 0.0115, k = 6, τ = 0.89, ρ = 0.338, μ = 0.002, α = 0.33 to obtain the following results. Table 1 shows the sensitivity indices of R eff to the parameters for the pneumonia model, evaluated at the baseline parameter values given in Table 2. The parameters are ordered from most sensitive to least. The most sensitive parameter is the contact rate, and the least sensitive parameter is the progression proportion of the disease. This result implies that, when the parameters k, , τ , φ and χ are increased keeping other parameters constant, they increase the value of R eff thus, they increase the endemicity of the disease as they have positive indices. While the parameters χ , p, ϑ, μ, α, ρ, β, η and q decrease the value of R eff

Extension of the model into an optimal control
In this section, we apply optimal control strategies on the model (1). This helped us to identify the best intervention strategies that helps to eradicate the disease in the specified time. The optimal control model is an extension pneumonia model by including the following three controls defined as (i) u 1 a prevention effort, that protect susceptible from contacting the disease.
(ii) u 2 a treatment effort, to minimize infection by treating infectious. (iii) u 3 a screening effort,to help carriers to screen themselves.
After incorporating, u 1 , u 2 and u 3 in pneumonia model (1), we obtain the following optimal control model of pneumonia: To study the optimal levels of the controls, the control set U is Lebesgue measurable and it is defined as Our aim is to obtain a control u and S,V,C,I and R that minimize the proposed objective function J and the form of the objective functional is taken in line with literature on epidemic models [1], given by where b 1 , b 2 and w i are positive. The expression 1 2 w i u 2 i represents cost which is associated with the controls u i .The form is quadratic because we assume that costs are nonlinear in its nature. Our aim is to minimize the number of carriers, infectives and costs. Thus, we seek to find an optimal triple controls (u * 1 , u * 2 , u * 3 ) such that

The Hamiltonian and optimality system
By using the principle of Pontryagin et al. [15], 'Pontryagin's Maximum Principle Pontryagin', we got the necessary conditions which is satisfied by optimal pair. Therefore, by this principle, we obtained a Hamiltonian (H) defined as where L(C, I, u 1 , 2, 3, 4, 5 are the adjoint variable functions to be determined suitably by applying Pontryagin's maximal principle [15] and also using [3] for existence of the optimal control pairs. Theorem 5.1: For an optimal control set u 1 , u 2 , u 3 that minimizes J over U, there is an adjoint variables, λ 1 , . . . , λ 5 such that: where,

Proof:
The form of the adjoint equation and transversality conditions are standard results from Pontryagin's maximum principle [15]. We differentiate Hamiltonian 5.1 with respect to states S, V, C, I and R, respectively, and then the adjoint system can be written as Similarly by following the approach of Pontryagin et al. [15], to get the controls, we solved the equation, ∂H/∂u i = 0 at u * i , for i = 1, 2, 3 and obtained: When we write by using standard control arguments involving the bounds on the controls, we conclude In compact notation The optimality system is formed from the optimal control system (the state system) and the adjoint variable system by incorporating the characterized control set and initial and transversal condition.

Uniqueness of the optimality system
Due to the priori boundedness of the state, adjoint functions and the resulting Lipschitz structure of the ODEs, we can obtain the uniqueness of solutions of the optimality system for the small time interval. Hence the following theorem. For the proof of the theorem, see [2].

Numerical simulations
In this section, we perform some numerical experimentation on the basic model (1) and the resulting optimality system consisting of the state equations (16) and the adjoint system (18). We make use of the parameter values given in Table 2 for the simulation. An iterative scheme is used to find the optimal solution of the optimality system. Since the state system (1) has initial conditions and the adjoint systems (18) have final conditions, we solve the state system using a forward fourth-order Runge-Kutta method and solve the adjoint system using a backward fourth-order Runge-Kutta method. The solution iterative scheme involves making a guess of the controls and using that guess to solve the state system. The initial guess of the control together with the solution of the state systems is used to solve the adjoint systems. The controls are then updated using a convex combination of the previous controls and the values obtained using the characterizations. The updated controls are then used to repeat the solution of the state and adjoint systems. This process is repeated until the values in the current iteration are close enough to the previous iteration values [7].
Using different combinations of the controls, such as one control only at a time, two controls at a time and also all controls at a time, that we analyse and compare numerical results from simulations with the following scenarios.
(vii) Using all the three controls, prevention u 1 treatment of infective (u 2 ) and screening of carriers u 3 .

Control with prevention only
We simulate the model by preventive intervention only. From Figure 3, we see that the decrease of infectious and carrier population due to implementation of prevention. This can be attribute the fact that prevention minimizes the rate of joining of individuals in to infective as well as carrier compartments. This implies that optimized prevention reduces the burden of the infection of pneumonia. Figure 4 shows a decrease of infectious population up to 4 month, then after start to go up. Those individuals, who were previously with the disease are being treated and that is why the number of infective population goes down for the first four months. Then, due to lack of prevention newly infected individuals start to join the infective as well as the carrier classes. That is why the number of infective start to goes up after four months of going down and the number of carrier also starts to go up after five months.

Control with screening only
Screening helps carriers to move into the infective classes and start to get treatment. Figure 5 shows a decrease in carrier population up to five months and then start to increase because due to lack of prevention. Susceptible start to be infected and joins carrier as well as infective classes. As a result of this screening only might not be sufficient to eradicate the burden of the infection of pneumonia.

Control with prevention and treatment
We used prevention and treatment as intervention strategy, and Figure 7 show that, the number of infective and also carriers goes down in the specified time. Therefore, this  strategies is effective in eradicating the disease from the community in a specified period of time ( Figure 6).

Control with prevention and screening
In this strategy we used prevention and screening. The first Figure 7 shows that the curve for optimal control is above the curve of with out control. Due to the reason that there is no treatment but individuals from carrier groups are joining infective compartment by screening and also there are a number of infected people in the compartment before prevention with out getting treatment so this situation make the curve to goes up for a time being. After some time the number of infectious goes down because due to prevention strategies, new infection is no more coming and also since there is no treatment the number of infective population start to goes down by disease causing death and natural death rates.

Control with treatment and screening
We used treatment and screening controls as intervention. From Figures 8, we observe that optimal control of the combination of treatment and screening helps to bring down the infectious as well as the carrier population which helps to eradicate the disease in the community.

Control with prevention, treatment and screening
We implement all control the three controls interventions that helps to minimize the objective function. From Figure 9, we observe that the number of the infectious and carrier populations decrease at the specified time due to the intervention strategies. Therefore, applying this strategy helps to eradicate pneumonia disease in specified period of time ( Figure 10).

Cost-effectiveness analysis
Cost-effectiveness analysis used to rank the implemented strategies interims of their cost. Applying one intervention only might to be effective to eradicate the disease from the community. Therefore, we analysed strategies that used more than one intervention method. To achieve, this we used incremental cost-effectiveness ratio (ICER), stated by Baba and Makinde [1]; "ICER = Difference in costs between strategies Difference in health effects between strategies .  In Table 3 we obtain the total number of infectious averted and total cost for the implemented strategies. The difference between the total infectious individuals without control and the total infectious individuals with control is used to obtain the total number of infectious averted. And also to find the total cost for the implemented strategies we used the cost function, which is 1 2 w 1 u 2 1 , 1 2 w 2 u 2 2 and 1 2 w 3 u 2 3 over time. We used the parameter values in Table 2 and to apply ICER technique first we ordered the intervention strategies for pairwise comparison as in Table 2 from A to D with increasing order of effectiveness.  Similarly, from ICER (B) and ICER (C) we can see that strategy C saves 0.0015 than strategy B. Therefore, we exclude strategy B, because it is a bit expensive and finally, we compared strategy C and D. From ICER (C) and ICER (D) we can see that strategy C saves 0.71 than strategy D. Therefore, we exclude strategy D, because it is a bit expensive. Therefore, we conclude that strategy C the cheapest of all compared strategies, that meant it is the most cost-effective for pneumonia disease control interventions strategy's. For further elaboration, Figure 11 shows that applying only one intervention costs the least interims of price but we did not consider this, due to the reason that a single intervention is not effective to eradicate the disease. And additionally we observe from the figure, applying all the three intervention at once is the most expensive of all the applied intervention strategies.

Discussions and conclusions
In Section 2 we described and proposed a pneumonia model, which is deterministic in its nature and also the population is assumed to be variable in size. In Section 3 we analysed the model by obtaining the feasible region, positivity of the solution set, effective reproductive number, equilibria points and their stability. Moreover, we described the possibility of backward bifurcation in Section 3. In Section 4, sensitivity analysis and interpretation of the sensitivity index is done. In Section 5 we extended the model by applying optimal control interventions and we obtained the Hamiltonian, the adjoint variables, the characterization of the controls and the optimality system. In Section 6 we the optimality system by considering different strategies as follows: • By applying a single control, prevention, treatment or screening.
• By applying two control, prevention and treatment.
• By applying two control, prevention and screening.
• By applying two control, treatment and screening.
• By applying all control, prevention, treatment and screening.
In Section 7 numerically we investigated cost effectiveness to determine, the least and the most expensive strategies by using ICER. From the pairwise result in this study, we conclude that, the combination of prevention and treatment is the best cost-effective strategies interims of cost as well as health benefits.