Mathematical Methods and Models in the Natural to the Life SciencesView this Special Issue
Optimal (Control of) Intervention Strategies for Malaria Epidemic in Karonga District, Malawi
Malaria is a public health problem for more than 2 billion people globally. About 219 million cases of malaria occur worldwide and 660,000 people die, mostly (91%) in the African Region despite decades of efforts to control the disease. Although the disease is preventable, it is life-threatening and parasitically transmitted by the bite of the female Anopheles mosquito. A deterministic mathematical model with intervention strategies is developed in order to investigate the effectiveness and optimal control strategies of indoor residual spraying (IRS), insecticide treated nets (ITNs) and treatment on the transmission dynamics of malaria in Karonga District, Malawi. The effective reproduction number is analytically computed, and the existence and stability conditions of the equilibria are explored. The model does not exhibit backward bifurcation. Pontryagin’s Maximum Principle which uses both the Lagrangian and Hamiltonian principles with respect to a time dependent constant is used to derive the necessary conditions for the optimal control of the disease. Numerical simulations indicate that the prevention strategies lead to the reduction of both the mosquito population and infected human individuals. Effective treatment consolidates the prevention strategies. Thus, malaria can be eradicated in Karonga District by concurrently applying vector control via ITNs and IRS complemented with timely treatment of infected people.
Malaria is a vector-borne infectious disease found mainly in tropical regions (Sub-Saharan Africa, Central and South America, the Indian subcontinent, Southeast Asia, and the Pacific islands) . It is a life-threatening disease transmitted through the bites of infected mosquitoes . There are four different types of Plasmodium parasites: Plasmodium falciparum (the only parasite which causes malignant malaria), Plasmodium vivax (causes benign malaria with less severe symptoms; the vector can remain in the liver for up to three years and can lead to a relapse), Plasmodium malariae (also causes benign malaria and is relatively rare), and Plasmodium ovale (causes benign malaria and can remain in the blood and liver for many years without causing symptoms). This study focuses mainly on malignant malaria. Severe malaria can affect the patient’s brain and central nervous system and can be fatal [2, 3]. Furthermore, Medicinenet  has reported another relatively new species Plasmodium knowlesi which has been causing malaria in Malaysia and areas of Southeast Asia. It is also a dangerous species that is typically found only in long-tailed and pigtail macaque monkeys. Like P. falciparum, P. knowlesi may be deadly to anyone infected .
Beyond the human toll, malaria wreaks significant economic havoc in endemic regions, decreasing gross domestic product (GDP) by as much as 1.3% in countries with high levels of transmission—the disease accounts for up to 40% of public health expenditures, 30–50% of in-patient hospital admissions, and up to 60% of out-patient health clinic visits . The Government of Malawi has put in place several control strategies through the National Malaria Control Programme to reduce and possibly (ultimately) eliminate malaria. The main strategic areas that have been identified for scaling-up of malaria control activities include malaria case management, intermittent preventive treatment (IPT) of pregnant women using sulfadoxine-pyrimethamine (SP), and malaria prevention with special emphasis on the use of ITNs as well as IRS .
Various compartmental models for the spread of malaria have been proposed [10–12]. Extensive review of mathematical models of malaria dynamics using -type models and their variants can be found in Chiyaka et al. [13, 14]. For information on the unexpected stability of malaria elimination and eradication and additional literature on malaria dynamics see Smith et al. , who also showed that if mosquito birth rate is increased by 25%, then the minimal effective spraying period is reduced by half, and if doubled, the period is reduced by three-quarters. Our goal is to assess the role of optimal control of preventive measures, namely, ITNs and IRS as well as treatment of the transmission dynamics of malaria in Karonga District, Malawi, without actually targeting a certain high-risk group. Despite a plethora of studies on the dynamics of malaria and its control (see Chiyaka et al.  and Smith et. al. ), to the best of our knowledge, the proposed model with exposed immigrants is seemingly new: the use of optimal values of a combination of three intervention strategies, namely, indoor residual spraying (IRS), insecticide treated nets (ITNs), and treatment to investigate the effectiveness and optimal control strategies in the transmission dynamics of malaria in a locality in a resource constrained setting.
2. Model Formulation and Analysis
We formulate an optimal control model for malaria with the population under study being subdivided into compartments according to individuals’ disease status. We consider the total population sizes denoted by and for the human hosts and Anopheles female mosquitoes, respectively. We employ the framework to describe a disease with temporary immunity on recovery from infection. The model indicates that the passage of individuals is from the susceptible class, , to the exposed class, , then to the infectious class, , and finally to the recovery class, . represents the number of individuals not yet infected with the malaria parasite at time . The latent or exposed class represents individuals who are infected but not yet infectious. Individuals in the class are infected with malaria and are capable of transmitting the disease to susceptible mosquitoes. represents the class of individuals who have temporarily recovered from the disease. The susceptible human population is increased by recruitment (birth) at a constant rate, , while others are generated through migration by , where is the proportion of infected immigrants into the exposed class and is the rate at which people migrate into Karonga District. All the recruited individuals are assumed to be naive when joining the community.
There is some finite probability, , that the parasites (in the form of sporozoites) will be passed onto the humans when an infectious female Anopheles mosquito bites a susceptible human. The parasite then moves to the liver where it develops into its next life stage, merozoites. Using the approach adopted in , susceptible individuals acquire malaria through contact with infectious mosquitoes at the rate . The infected person moves to the exposed class at the rate . The preventive variable represents the use of ITNs as a means of minimizing or eliminating mosquito-human contacts. After a certain period of time, the parasite (in the form of merozoites) enters the blood stream, usually signaling the clinical onset of malaria. Then the exposed individuals become infectious and progress to the infected state at a constant rate . The individuals who have experienced infection may recover with temporary immunity at a constant rate and move to the recovery class, while some infectious humans after recovery without immunity become immediately susceptible again at the rate . Infectious individuals recover due to treatment at a rate with representing the control effort on treatment and being the proportion of individuals who recover spontaneously. Recovered individuals lose immunity at a rate . The natural and disease induced death rates are and , respectively. The disease induced death rate is very small in comparison with the recovery rate.
The mosquito population is divided into three compartments: susceptible ; exposed ; and infectious . Female Anopheles mosquitoes enter the susceptible class through birth at a rate . The parasites in the form of gametocytes enter the mosquito population with probability . This happens when the mosquito bites an infectious human and the mosquito moves from the susceptible to the exposed class. Mosquitoes are assumed to suffer death due to natural causes at a rate . The exposed mosquitoes progress to the class of symptomatic mosquitoes at a rate . It is assumed that the disease does not induce death to the mosquito population. Finally, the mortality rate of the mosquito population increases at a rate proportional to , where is a rate constant when using IRS and is a constant rate of IRS. Figure 1 represents a schematic flow diagram of the proposed model.
The model state variables are represented in Table 1. Table 2 represents prevention and control strategies practised in the district. Table 3 shows parameters of the model.
The state variables in Table 1, the prevention and control parameters in Table 2, and the model parameters in Table 3 for the malaria model satisfy (1). It is assumed that all state variables and parameters of the model which monitors human and mosquito populations are positive for all . We will, therefore, analyse the model in a suitable region.
The assumptions lead to the following deterministic system of nonlinear ordinary differential equations which describe the evolutionary dynamics of a malaria model with a combination of interventions: where , .
The term denotes the rate at which the human hosts become infected by infectious mosquitoes and refers to the rate at which the susceptible mosquitoes are infected by the infectious human hosts .
Adding the first six equations of model (1) and assuming that there is no disease induced death, that is, , gives , so that as . Thus, is an upper bound of provided that . Further, if , then will decrease to this level, . Similar calculation for the vector equations shows that as .
The feasible region is positive-invariant and attracting.
Lemma 1. The region is positively invariant for the model system (1) with initial conditions in .
Proof. Let give . The first equation of model (1) gives which can be rewritten as Therefore, so that Similarly, it can be shown that , , , , , and , for all . This completes the proof.
2.1. Existence and Stability of Equilibrium Points
We analyse system (1) to obtain the equilibrium points of the system and their stability. Let be the equilibrium points of system (1). At an equilibrium point, we have
2.1.1. Disease-Free Equilibrium,
In the absence of malaria, that is, , model system (1) has an equilibrium point called the disease-free equilibrium, , and is given by To establish the linear stability of , we employ van den Driessche and Watmough’s next generation matrix approach . A reproduction number obtained this way determines the local stability of the disease-free equilibrium point for and instability for . Following van den Driessche and Watmough , the associated next generation matrices and of system (1) can be determined from and , respectively, where The Jacobian matrix of and is given by where Algebraic manipulation of the matrices leads to the effective reproduction number where is the contribution of the mosquito population when it infects the humans, and is the human contribution when they infect the mosquitoes.
The expression for the effective reproduction number, , has a biological meaning that is readily interpreted from terms under the square root sign. Consider the following terms.(i) represents the number of secondary human infections caused by one infected mosquito vector.(ii) represents the number of secondary mosquito infections caused by one infected human host.
The square root represents the geometric mean of the average number of secondary host infections produced by one vector and the average number of secondary vector infections produced by one host. This effective reproduction number serves as an invasion threshold both for predicting outbreaks and evaluating control strategies that would reduce the spread of the disease. The threshold quantity, , measures the average number of secondary cases generated by a single infected individual in a susceptible human population , where a fraction of the susceptible human population is under prevention and the infected class is under treatment. In the absence of any protective measure, the effective reproduction number with treatment is Also if ITNs are the only intervention strategy, then Similarly, if IRS is the only means of protection, then The value of the basic reproduction number can be obtained from the value of the effective reproduction number when there are no control measures : In general, it is easy to prove that for , , due to reduction of likelihood of infection by protection. This implies that ITNs and IRS have a positive impact on the malaria dynamics as they contribute to the reduction of secondary infections. Therefore, from van den Driessche and Watmough  (Theorem 2), the following result holds.
Lemma 2. The disease-free equilibrium, , of the malaria model with intervention strategies (1), given by (8), is locally asymptotically stable if and unstable if .
2.1.2. Endemic Equilibrium Point,
When malaria is present, the model system (1) has a steady state, , called the endemic equilibrium. In order to establish the stability of , we express system (1) in dimensionless variables, namely, , , , , , , and , where , . Then differentiating with respect to time, , respectively results in the following system: Let , , and .
Then The system can now be reduced to a nine-dimensional system by eliminating and since and . The following feasible region can be shown to be positively invariant. denotes the nonnegative cone of including its lower dimensional faces. Thus we have the following system of equations: We seek to establish whether a unique endemic equilibrium exists. This is done by making more realistic assumptions, namely, that the protective control measures may not be totally effective.
Existence and Uniqueness of Endemic Equilibrium, . To compute the steady states of the system (23), we set the derivative with respect to time in (23) equal to zero, and after simplification the following algebraic equations are obtained: To calculate the dimensionless proportions in terms of , consider the first equation in the system (24): where .
From the second equation we have Equating these two equations for we obtain Therefore, with and becomes Now solve for from the third equation Let , , , , , , and . Then substituting in the equation we obtain Using equation five of the system (24), we obtain Solving for gives Let , , , and . Substituting in the equation gives Substituting in the equation gives The existence of the endemic equilibrium in , can be determined when . After some algebraic manipulations we obtain where For , the existence of endemic equilibria is determined by the presence of positive real solutions of the quadratic expression (35). Consider Since and , then . Thus, there exists exactly one positive endemic equilibrium for whenever . This gives the threshold for the endemic persistence. Therefore, we have proved the existence and uniqueness of the endemic equilibrium, , for the system (1). This result is summarized in the following theorem.
Theorem 3. If , then the model system (23) has a unique endemic equilibrium .
The result in Theorem 3 indicates the impossibility of backward bifurcation in the malaria model, since it has no endemic equilibrium when . Thus, the model (1) has a globally asymptotically stable disease-free equilibrium whenever .
Next we need to investigate the global stability property of the endemic equilibrium for the case when the disease does not induce death.
2.1.3. Global Stability of the Endemic Equilibrium for
The model system (1), when , has a unique endemic equilibrium. By letting then the following can be claimed.
Theorem 4. The endemic equilibrium point of the malaria model (1) with is globally asymptotically stable whenever in and .
Proof. As for the case of Theorem 3, it can be shown that the unique endemic equilibrium for this special case exists only if . Additionally, as . Letting and and substituting into (1) give the limiting system
Dulac’s multiplier (see ) gives
since and in .
Thus, by Dulac’s criterion, there are no periodic orbits in and . Since is positively invariant and the endemic equilibrium exists whenever , then from the Poincare-Bendixson Theorem , it follows that all solutions of the limiting system originating in remain in for all . Furthermore, the absence of periodic orbits in implies that the unique endemic equilibrium of the special case of the malaria model is globally asymptotically stable whenever .
The malaria model has a locally asymptotically stable disease-free equilibrium whenever and a unique endemic equilibrium whenever . In addition, the unique endemic equilibrium is globally asymptotically stable for the case if .
In the next section, we apply the optimal control method using Pontryagin’s Maximum Principle to determine the necessary conditions for the combined optimal control of ITNs, IRS, and treatment effort which are being practised in Karonga District, Malawi.
3. Analysis of Optimal Control of the Malaria Model
The force of infection in the human population is reduced by a factor of where represents the use of ITNs as a means of minimizing or eliminating mosquito-human contact. represents the control effort on treatment of infectious individuals. This indeed represents the situation when individuals in the community seek treatment after visiting the hospitals or dispensary in their areas. For the mosquito population, we have a third control variable, . IRS affects the whole mosquito population by increasing its mortality rate by . We will use an approach similar to that in Lashari and Zaman  which consists of applying Pontryagin’s Maximum Principal to determine the conditions under which eradication of the disease can be achieved in finite time. Following the dynamics of the model system (1) with appropriate initial conditions, the bounded Lebesgue measurable control is used with the objective functional defined as subject to the differential equations in (1), where , , , , , are positive weights. We choose a quadratic cost on the controls in line with what is known in the literature on epidemic controls [21–23].
The purpose of this section is to minimize the cost functional (41). This functional includes the exposed and infectious human population and the total mosquito population. In addition, it has the cost of implementing personal protection using ITNs, , treatment of infected individuals, , and spraying of houses, . A linear function has been chosen for the cost incurred by exposed individuals, , infected individuals, , and the mosquito population, . A quadratic form is used for the cost on the controls , , and , such that the terms , , and describe the cost associated with the ITNs, treatment and IRS, respectively. We select to model the control efforts via a linear combination of quadratic terms , and the constants and , representing a measure of the relative cost of the interventions over the time horizon . We seek an optimal control , , and such that where is the control set, subject to the system (1) and appropriate initial conditions. We develop the optimal system for which the necessary conditions it must satisfy come from Pontryagin’s Maximum Principle .
3.1. Existence of an Optimal Control Problem
Pontryagin’s Maximum Principal converts (1) and (41) into a problem of minimizing pointwise the Lagrangian, , and Hamiltonian, , with respect to , , and . The Lagrangian of the control problem is given by We search for the minimal value of the Lagrangian. This can be done by defining the Hamiltonian, , for the control problem as
Next, we prove the existence of an optimal control for the model system (1).
Theorem 5. The model system (1) with the initial conditions at has control strategies and there exists an optimal control such that
Proof. The state and the control variables of the system (1) are nonnegative values. The control set is closed and convex. The integrand of the objective cost function expressed by (1) is a convex function of on the control set . The Lipschitz property of the state system with respect to the state variables is satisfied since the state solutions are bounded. It can easily be shown that there exist positive numbers and a constant such that This concludes the existence of an optimal control because the state variables are bounded.
3.2. Classification of the Optimal Control Problem
We use Pontryagin’s Maximum Principle to develop the necessary conditions for this optimal control since there exists an optimal control for minimizing the functional (41) subject to the system of equations in (1). From Lashari and Zaman , if is an optimal solution of an optimal control problem, then there exists a nontrivial vector function satisfying the following equations: Hence the necessary conditions of the Hamiltonian, , can be applied in (45).
Theorem 6. For the optimal control triple with their optimal state solutions that minimizes over , then there exist adjoint variables satisfying with transversality conditions Additionally, the optimal control triple that minimizes over satisfies the optimality condition
Proof. The adjoint equations can be determined by using the differential equations governing the adjoint variables. The Hamiltonian function, , is differentiated with respect to , , , , , , and . The adjoint equation is given by with the transversality conditions Solving , , and , evaluating at the optimal control on the interior of the control set, where , for , and letting , , , , , , and yields for which Then We achieve the uniqueness of the optimal control for small due to the prior boundedness of the state and adjoint functions and the resulting Lipschitz structure of the ordinary differential equations. The uniqueness of the optimal control triple trails from the uniqueness of the optimal system, which consists of (1), (49), and (50) with characterization of the optimal control (51).
The optimality system is comprised of the state system (1), the adjoint system (49), initial conditions at , boundary conditions (50), and the characterization of the optimal control (51). Hence the state and optimal control can be calculated using the optimality system. Hence using the fact that the second derivatives of the Lagrangian with respect to , , and , respectively, are positive indicates that the optimal problem is a minimum at controls , , and . Substituting , , and in the system (1), we obtain with at :