Abstract

This paper proposes and analyses a mathematical model for the transmission dynamics of malaria with four-time dependent control measures in Kenya: insecticide treated bed nets (ITNs), treatment, indoor residual spray (IRS), and intermittent preventive treatment of malaria in pregnancy (IPTp). We first considered constant control parameters and calculate the basic reproduction number and investigate existence and stability of equilibria as well as stability analysis. We proved that if , the disease-free equilibrium is globally asymptotically stable in . If , the unique endemic equilibrium exists and is globally asymptotically stable. The model also exhibits backward bifurcation at . If , the model admits a unique endemic equilibrium which is globally asymptotically stable in the interior of feasible region . The sensitivity results showed that the most sensitive parameters are mosquito death rate and mosquito biting rates. We then consider the time-dependent control case and use Pontryagin’s Maximum Principle to derive the necessary conditions for the optimal control of the disease using the proposed model. The existence of optimal control problem is proved. Numerical simulations of the optimal control problem using a set of reasonable parameter values suggest that the optimal control strategy for malaria control in endemic areas is the combined use of treatment and IRS; for epidemic prone areas is the use of treatment and IRS; for seasonal areas is the use of treatment; and for low risk areas is the use of ITNs and treatment. Control programs that follow these strategies can effectively reduce the spread of malaria disease in different malaria transmission settings in Kenya.

1. Introduction

Malaria is a leading cause of mortality and morbidity among the under-five group and the pregnant women in Sub-Saharan Africa [1]. These groups are at high risk due to weakened and immature immunity, respectively. With the recent conversion of the Millennium Development Goals (MDGs) to Sustainable Development Goals (SDGs) as part of Global Malaria Action Plan for a malaria-free world by 2030, reducing malaria is critical to post-2015 malaria strategies and for achieving the SDGs such as ensuring healthy lives and promoting well-being for all at all ages. Most Kenyans are vulnerable to malaria because of poverty, inadequate health care infrastructures, and low income of the country. Prompt access to effective treatment for malaria is unacceptably low in Kenya due to the socioeconomic barriers to accessing health care. These challenges call for urgent development of effective and optimal strategies for preventing and controlling spread of malaria.

Malaria transmission is highly variable across Kenya because of the different transmission intensities driven by climate and temperature. Kenya has four malaria epidemiological zones: the endemic areas, the seasonal malaria transmission, the malaria epidemic prone areas, and the low risk malaria areas [2]. Malaria is caused by Plasmodium parasites and it is transmitted from one individual to another by the bite of infected female anopheline mosquitoes [1]. Malaria is an entirely preventable and treatable disease, provided the currently recommended interventions are properly implemented. Controlling malaria transmission involves interrupting the malaria transmission for specific transmission settings and for the most at-risk groups of malaria. World Health organization (WHO) recommended malaria intervention strategies include the use of long-lasting insecticide-treated bed nets (LLINs), indoor residual spray (IRS), chemoprevention for the most vulnerable such as intermittent preventive treatment for pregnant women (IPTp), confirmation of malaria diagnostics through rapid diagnostics tests (RDTs) and microscopy for every suspected case, and timely treatment with artemisinin-based combination therapies (ACTs) [1, 3].

Mathematical modelling has become an important tool in understanding the dynamics of disease transmission and in decision making processes regarding intervention programs for disease control. Mathematical models provide a framework for understanding the transmission dynamics for malaria and can be used for the optimal allocation of different interventions against malaria [4, 5]. Optimal control is a branch of mathematics developed to find optimal ways of controlling a dynamic system [6]. Application of optimal control theory can be an important tool for estimating the efficacy of various policies and control measures and the cost of implementing them. Optimal control approaches have been successfully previously used in decision making for the infectious diseases such as tuberculosis [7, 8], HIV [9], and influenza [10]. Optimal control theory has also been applied in malaria control to assess the impact of antimalarial control measures by formulating the model as an optimal control problem. Most malaria models for analyzing effect of interventions in optimal control used the standard Susceptible-Exposed-Infectious-Recovered (SEIR) model for humans and Susceptible-Exposed-Infectious (SEI) model for mosquitoes. Okosun et al. [11] considered three control variables in assessing the optimal control and cost effectiveness of the interventions but not for different settings. Mwamtobe et al. [12] used three control variables and for only one region in Malawi. Kim et al. [13] used two control efforts in the optimal control model for malaria transmission in South Korea. Agusto et al. [14] used three system control variables. Silva and Torres [15] used one control variable.

IPTp is one of the WHO recommended prevention therapies for the pregnant women. IPTp has been shown to be effective in reducing maternal and infant mortality that are related to malaria for the most at-risk group for malaria [1619]. Very few studies have been carried out in applying optimal control theory to malaria transmission models for different transmission settings. The combined effect of ITNs, IRS, and natural death on reducing the mosquito population has not been demonstrated in optimal control theory in malaria control. The effect of IPTp which is WHO recommended preventive therapy for the most at-risk group for malaria (pregnant women) has not been studied in optimal control theory. No optimal control model for malaria interventions for different transmission settings exits for Kenya. No optimal control model for four control variables incorporating the IPTp malaria intervention studies exits for Kenya. No optimal control model has been stratified by the age group (under five) and specific categories (pregnant women).

In this paper, a model for malaria transmission dynamics with four control strategies is formulated and analyzed. We then formulate an optimal control problem and derive expressions for the optimal control for the malaria transmission model with four control variables and then use the optimal control theory to study the effectiveness of all possible combinations of four malaria preventive measures among the pregnant women and children under five years of age.

2. Model Formulation

The model is formulated by considering the human and mosquito subgroups. The considered model consists of population of susceptible , Exposed humans , infected humans , recovered humans , susceptible mosquitoes , exposed mosquitoes, and infected mosquitoes . The total population sizes at time for humans and mosquitoes are denoted by and , respectively. We employ the SEIRS type model for humans to describe a disease with temporary immunity on recovery from infection. Mosquitoes are assumed not to recover from the parasites so the mosquito population can be described by the SEI model. In the model we incorporate four time-dependent control measures simultaneously: (i) the use of treated bed nets , (ii) treatment of infective humans , (iii) spray of insecticides , and (iv) treatment to protect pregnant women and their newborn children: intermittent preventive treatment for pregnant women (IPTp) . represents the number of individuals not yet infected with the malaria parasite at time , represents individuals who are infected but not yet infectious, is the class representing those who are infected with malaria and are capable of transmitting the disease to susceptible mosquitoes, and represents the class of individuals who have temporarily recovered from the disease.

Figure 1 describes the dynamics of malaria in human and mosquito populations together with interventions.

The susceptible humans (pregnant women and children under the age of five) () are recruited at the rate . They either die from natural causes (at a rate ) or move to the exposed class () by acquiring malaria through contact with infectious mosquitoes at a rate or , where is the transmission probability per bite, is the per capita biting rate of mosquitoes, is the contact rate of vector per human per unit time, is the preventive measure using ITNs, is the preventive measure using IPTp, is the infectious mosquitoes at time , is the total number of individuals (pregnant women and children under the age of five), and is the total number of pregnant women. Susceptible class is divided into whole population (children under the age of five and pregnant women) being exposed and the population for the pregnant women being exposed. Exposed individuals move to the infectious class after the development of clinical symptoms at the rate . Infectious individuals are assumed to recover at a rate , where is the rate of spontaneous recovery, is the control on treatment of infected individuals, and is the efficacy of treatment. Infectious individuals who do not recover die at a rate . Individuals infected with malaria suffer a disease induced death at rate of and natural death . Infected individuals then progress to partially immune group where upon recovery the partially immune individual loses immunity at the rate and becomes susceptible again.

Susceptible mosquitoes () are recruited at the rate and acquire malaria infection (following contact with humans infected with malaria) at the rate . They either die from natural causes (at a rate ) or move to the exposed class by acquiring malaria through contacts with infected humans at a rate , where is the probability for a vector to get infected after biting an infectious human and are individuals infected with malaria at time . The mosquito population is reduced, due to the use of insecticides spray, at a rate , where represents the control due to IRS and represents the efficacy of IRS. Mosquito population is also reduced as a result of natural death () and at the rate , where represents the control due to ITNs and is the efficacy due to ITNs. Newly infected mosquitoes are moved into the exposed class () at a rate and progresses to the class of symptomatic mosquitoes ().

The state variables of the model are represented and described in Table 1. Table 2 describes and shows parameters of the model. Table 3 describes and represents malaria prevention and control strategies practiced in Kenya.

The following assumptions have been used in the formulation of the model:(i)Population for human and mosquito being constant (no immigrants).(ii)No recovery for infected mosquitoes.(iii)Mosquitoes not dying due to disease infection.(iv)All parameters in the model being nonnegative.Putting the above formulations and assumptions together gives the following system of nonlinear differential equations describing the dynamics of malaria in human and mosquito populations together with interventions:With initial conditions is the per capita incidence rate among mosquitoes (force of infection for susceptible mosquitoes), is the force of infection for susceptible humans (pregnant and under 5), and is the force of infection for susceptible pregnant humans.

The total population size for the human is and for mosquito is and their differential equations are given by and , respectively.

2.1. Mathematical Analysis of the Malaria Model with Intervention Strategies

We will assume that the control parameters are constant so as to determine the basic reproduction number, steady states, and their stability as well as the bifurcation analysis.

We describe the basic properties and analysis of the formulated malaria model with control strategies through mathematical analysis of the formulated model.

2.1.1. Basic Properties of the Model: Positivity and Invariant Regions

All the state variables and parameters for model (1) are nonnegative for all .

The feasible solutions set for model (1) given by is positively invariant and hence model (1) is biologically, epidemiologically meaningful and mathematically well posed in the domain .

System (1) has always a disease-free equilibrium given by

2.1.2. Basic Reproduction Number

The matrices and for the new infection terms and the remaining transfer terms at disease-free equilibrium [26], respectively, are given byIt follows that the basic reproduction of model (1), denoted by , is given by

2.1.3. Stability Analysis of Disease-Free Equilibrium Point

(1) Local Stability of Disease-Free Equilibrium Point

Theorem 1. The disease-free equilibrium point for system (1) is locally asymptotically stable if .

Proof. The Jacobian matrix of the malaria model (1) at the disease-free equilibrium point is given byThe eigenvalues of the Jacobian matrix are the solutions of the characteristic equationExpanding the determinant into a characteristic equation we haveHence we havewhereThus, applying the Routh-Hurwitz criteria [27] to polynomial (10), we have that all determinants of the Hurwitz matrices are positive. Hence all the eigenvalues of the Jacobian have negative real part, implying that the DFE point is (at least) locally asymptotically stable .

Next, we study the global behavior of the disease-free equilibrium for system (1).

(2) Global Stability of Disease-Free Equilibrium Point

Theorem 2. The DFE, , of system of (1) is globally asymptotically stable if .

Proof. We consider the following Lyapunov function:whereComputing the derivative of along the solution of the system of differential equation (1),Thus we have established that if and the equality holds if and only if and . If then when and are sufficiently close to and , respectively, except when
On the boundary when , and and , as .
Therefore the largest compact invariant when is the singleton in . By LaSalle’s invariant principle [28], is globally asymptotically stable for in .

Next, we investigate the endemic equilibrium and its stability of system (1).

2.1.4. Stability Analysis of Endemic Equilibrium Point

First we determine the existence of the endemic equilibrium points.

The endemic equilibrium () of the model is given byTo derive , we solve model (1) by equating it to zero.

Substituting and solving for (as an expression of parameters only) through some algebraic manipulation giveorwhereWe use the quadratic formula to find the roots of (17); that is,HenceFrom the quadratic equation (17) we analyze the possibility of multiple endemic equilibria. It is important to note that the coefficient is always positive with and having different signs. All the negative terms in are in the form of where is the control that is bounded by 1.

It follows that(i)there is a unique endemic equilibrium if and or ,(ii)there is a unique endemic equilibrium if (i.e., if ),(iii)there are two endemic equilibria if , , and ,(iv)there are no endemic equilibria otherwise.Note that the hypothesis is equivalent to .

Thus the results of this section can be summarized in the following theorem.

Theorem 3. If , is an equilibrium of system (1) and it is locally asymptotically stable. Furthermore, there exists an endemic equilibrium if conditions in (i) are satisfied or two endemic equilibria if conditions in (iii) are satisfied. If , then is unstable and there exists a unique endemic equilibrium.

Item (iii) indicates the possibility of backward bifurcation in model (1) when . In the next section we will prove the occurrence of multiple equilibria for comes from the backward bifurcation and this will give information on the local stability of the endemic equilibrium.

(1) Local Stability Analysis of Endemic Equilibrium Point. We use centre manifold theory [29] to investigate the stability of the endemic equilibria for model (1) where we carry out bifurcation analysis of system (1) at . We consider a transmission rate as bifurcation parameter so that .

To apply the theory, we introduce dimensionless state variables into system (1).

The system of (1) can be written as We let , , , , , , and

Therefore system (1) in vector form can be written aswhere and are transposed matrices and with , Let be the bifurcation parameter; the system is linearized at disease-free equilibrium point when with . That is, The Jacobian matrix of (1) calculated at is given by Centre manifold approach [29] is then applied.

A right eigenvector associated with the eigenvalue zero is . Solving the system we have the following right eigenvector:The left eigenvectors satisfying are . Solving the system, the left eigenvector will be given byComputing for the sign of and as indicated in the theorem givesso that is always positive.

Therefore the following result is established.

Theorem 4. Model (1) exhibits backward bifurcation at whenever and and . Whenever and , then model (1) exhibits a forward bifurcation at .

Finally, we will investigate the global stability of the endemic equilibrium in the feasible region.

(2) Global Stability Analysis of Endemic Equilibrium Point. Global stability results for the endemic equilibrium can be obtained when it is unique and whenever it exists. We have established in Theorem 4 that if this implies the existence and uniqueness of the endemic equilibrium.

The global behavior of the endemic equilibrium of model (1) when it exists is explored by proving that such an equilibrium is globally asymptotic stable in the interior of the feasible region . We will use the geometric approach to global stability as described by Li and Muldowney [30]. The following conditions are required for the global stability of the endemic equilibrium, : (i) the uniqueness of in the interior of (condition ); (ii) the existence of an absorbing compact set in the interior of (condition ); and (iii) the fulfillment of a Bendixson criterion (i.e., inequality (A.6)).

Theorem 5. If and and , then the endemic equilibrium of the malaria model (1) is globally asymptotically stable in the interior of .

Proof. Following Li and Muldowney [30], for system (1), under the assumption of , satisfies conditions , the existence of the endemic equilibrium has also been shown, and the instability of DFE implies the uniform persistence [31]; that is, there exists a constant such that any solutions with in the interior of satisfyThe uniform persistence together with boundedness of is equivalent to the existence of a compact set in the interior of which is absorbing for (4) as described by Hutson and Schmitt [32]. Thus, is verified. Moreover, is the only equilibrium in the interior of , so that is also verified.
What remains is to find conditions for which the Bendixson criterion given by (A.6) is verified. To this aim, let us begin by observing that, from the Jacobian matrix associated with a general solution of reduced system (1), we get the second additive compound matrix :where where
Choose now matrix . Then , and the matrix can be written in block form aswhereThe vector norm in is here chosen to be Let denote the Lozinskii measure with respect to this norm. Using the method of estimating in [29], we havewhere and are matrix norms with respect to the vector norm and denotes the Lozinskii measure with respect to norm. Since is a scalar, its Lozinskii measure with respect to any norm in is equal to .
ThereforeThereforeWe rewrite the last two equations of system (1) for and as Substituting (39) into (37) and (40) into (38) we haveFor the uniform persistence constant , there exists a time independent of , the compact absorbing set, such thatfor we haveRelations (41) imply whereAlong each solution to (1) with where is the compact absorbing set, we have, for ,Which implies . This proves that the unique endemic equilibrium is globally asymptotically stable whenever it exist, thus completing the proof.

2.1.5. Sensitivity Analysis

Sensitivity analysis is to assess the relative impact of each of the parameters of the basic reproductive number. The normalized forward sensitivity index of the reproduction number with respect to these parameters given in Table 4 is computed. The index measures the relative change in a variable with respect to relative changes in parameters.

Definition 6. Following Chitnis et al. [33], the normalized forward sensitivity index of a variable, , that depends on a parameter, , is defined as .
The sensitivity index of the model parameters is given bySensitivity indices for the control parameters are given byThrough sensitivity analysis, it is observed that the most sensitive parameters to across all the settings were the mosquito’s natural death rate, , and mosquito biting rate, ; this was followed by the by the mosquito contact rate with humans, , probability of mosquito getting infected, , the probability of humans getting infected, , and the recruitment rate of mosquitoes and humans (see Table 4).

Sensitivity Indices of . Sensitivity indices are calculated using parameters in Table 5.

3. Analysis of Optimal Control

3.1. Analysis of Optimal Control of the Malaria Model with Intervention Strategies

We consider the case of time-dependent control variables. The malaria dynamics model is extended and an optimal control problem is formulated. We formulate an optimal control model for malaria disease in order to determine optimal malaria control strategies (ITNs, IRS, IPTp, and treatment) with minimal implementation cost.

For the optimal control problem of the given system, we consider the control variables relative to the state variables , , , , , , where control variables are bounded and measured withWe define the objective function assubject toIn the objective function is the final time and quantities , , and are weights constants of the total mosquito population, infected individuals, and exposed individuals, respectively, while , , , and are weight constants for use with ITNs, treatment effort, IRS, and IPTp efforts, respectively. The total mosquito population is part of the objective function because it is affected by the use of IRS and ITNs. In addition, and are included in the objective function because individuals in these classes are affected by the use of ITNs, IPTps, and treatment, respectively. A quadratic cost on the controls was chosen in line with what is known in the literature on epidemic optimal controls for malaria [11, 12]. The cost of implementing personal protection using ITNs is , treatment of infected individuals is , spraying of houses with IRS is , and preventive method of IPTp is . 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, mosquito control (IRS), and chemoprevention (IPTp), respectively.

Our aim with the given objective function is to minimize total number of mosquito population , the number of exposed humans , and infected humans while minimizing the cost of control , , , and . We select to model the control efforts via a linear combination of quadratic terms and the constants which represents a measure of the relative cost of the interventions over .

We seek an optimal control , , , and such thatPontryagin’s Maximum Principle is used to solve this optimal control problem and the derivation of necessary conditions that an optimal control must satisfy [6]. Pontryagin’s Maximum Principle converts the state system (1) and objective function (51) into a problem of minimizing pointwise the Lagrangian, , and Hamiltonian, , with respect to , , , and .

3.2. Existence of Optimal Control Problem

The existence of an optimal control can be proved by using the theorem given in Fleming and Rishel [34]. It can be clearly seen that the system of equations given by (1) is bounded from above by a linear system. According to the result in Fleming and Rishel [34], the solution exists if the following hypotheses are met:(H1):the set of controls and corresponding state variables is nonempty.(H2):the control set is convex and closed.(H3):right hand side of each equation of the state system in (1) is continuous, bounded by a linear functional in the state and control, and can be written as a linear function of with coefficients depending on time and state.(H4):there exist constants and such that the integrand of the objective functional is convex and satisfiesIn order to confirm the above hypotheses, a result given by Lukes [35] is used to establish the existence of solutions of state system (1). Since the coefficients are bounded, is satisfied. The solutions are bounded and hence the control set satisfies as well. The system is bilinear in and hence the right hand side of state system (1) satisfies condition (since the solutions are bounded). Note that the integrand of the objective functional is convex. The last condition is also satisfied.

The state and the control variables of system (1) are nonnegative values and nonempty. The control set is closed and convex. The integrand of the objective cost function expressed by (51) 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 where , , , , , , , , and

This concludes existence of an optimal control since the state variables are bounded. Hence we have the following theorem.

Theorem 7. Given the objective functional , where subject to (1) with initial conditions, then there exists an optimal control such that .

Lagrangian for a problem discusses how the techniques come and Hamiltonian helps in solving for the adjoint variable. In order to find an optimal solution, first we find the Lagrangian and Hamiltonian for the optimal control problem (51). The Lagrangian of the optimal problem is given byWe seek to find the minimal value of the Lagrangian. To do this, we define the Hamiltonian for the control problem which consists of the integrand of the objective functional (Lagrangian, ) and the inner product of the right hand sides of the state equations and the costate variables or adjoint variables asTaking , , and we obtain the Hamiltonian given by

3.3. The Optimality System

In order to find the necessary conditions for this optimal control, we apply Pontryagin’s Maximum Principle [6] as described by Lenhart and Workman [36] as follows.

If are optimal solution of an optimal control problem, then there exists a nontrivial vector function satisfying the following conditions.

The state equation isThe optimality condition isand the adjoint equation isNow, we apply the necessary conditions to the Hamiltonian.

Theorem 8. Given the optimal controls , , , and solutions , , , , , , of the corresponding state system (1), there exist adjoint variables , , , , , , satisfyingwith transversality conditionsFurthermore , , , are represented by

Proof. To determine the adjoint equations and the transversality conditions we use the Hamiltonian . The Hamiltonian function, , is differentiated with respect to , , , , , , and . The adjoint/costate equation is given bywith transversality conditions

In order to minimize Hamiltonian, , with respect to the controls at the optimal controls, is differentiated with respect to , , , and on the set , and the solution for the optimal control point is obtained after equating to zero. This is the optimality condition.

Solving , , , and , evaluating at the optimal control on the interior of the control set, where , for , and letting , , , , , , and yield for whichBy applying the boundary condition of each control, the solution of (67) becomes The optimality system is comprised of state system (1), adjoint system (61), initial conditions at , boundary conditions (62), and the characterization of the optimal control (63). 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 .

The optimality system is solved using the forward-backward fourth-order Runge-Kutta scheme in statistical computing platform [37]. The optimal strategy is obtained by solving the state and adjoint systems and the transversality conditions. First we start to solve state (1) using the Runge-Kutta fourth-order forward in time with a guess for the controls , , , and over the simulated time. Then, using the current iteration of the state equations with the initial guess for the controls, the adjoint equations in system (57) are solved by a backward method with the transversality conditions (62). Then the controls are updated by using a convex combination of the previous controls and the value from characterizations (63). This process is repeated and iterations stopped if the values of the unknowns at the previous iterations are very close to the ones at the present iterations [36].

4. Numerical Results on Optimal Control Analysis and Discussion

The parameters in model (1) were estimated using clinical malaria data and demographics statistics of Kenya. Those that were not available were obtained from literature published by researchers in malaria endemic countries which have similar environmental conditions compared to Kenya. Table 5 provides a summary of the estimated values of all parameters of the malaria model. Data was collected from the literature, Division of Malaria Control (DOMC), Kenya National Bureau of Statistics, Malaria Indicator Survey for Kenya, Demographic Health Survey (DHS) for Kenya, WHO websites, and hospital records (from Kisumu, Kisii, Chuka (Tharake-Nithi), and Nyeri counties representing the four different transmission settings/epidemiological zones in Kenya). The rates are per day.

In addition the effect of the different intervention strategies is estimated as , , , and and the cost of intervention is for , , , and (White et al. [38] and Hansen et al. [39]). The initial state variables are constant across all the epidemiological zones and are chosen as , , , , , , and . The following weight factors were also fixed for the different epidemiological scenarios as , , , , , , and .

4.1. Dynamics of Human Population with Intervention Strategies

We simulated malaria model with intervention strategies to find the dynamics of human and mosquito populations. It is observed that the control strategy leads to decrease in the number of infected human () as shown in Figure 2.

The figure shows steady decrease in susceptible human population at the initial period as the exposure of humans to disease increases. Thereafter graph of susceptible humans increases as the exposed and infected human population decreases due to positive effect of the intervention strategies being implemented.

4.2. Numerical Results on Optimal Control Analysis

In this section, we investigate numerically the effect of the several optimal control strategies on the spread of malaria. We compare the numerical results from the simulation using one control and various combinations of two, three, and four control strategies. This was done by comparing when there were not any intervention strategies and when there were intervention strategies. There are 15 different control strategies for each of the four different epidemiological zones in Kenya that are explored. We use the case of endemic zone with the case of one control variable, two control variables, three control variables, and all the four control variables being in use for the illustration purpose.

The results in Figures 36 show a significant difference in and with the control strategy compared to and without the control strategy. It is observed that the control strategy leads to decrease in the number of infected humans (). The uncontrollable case leads to a decrease in the number of infected mosquitoes (), while the control strategy leads to decrease in the infected number. The control profiles show the upper bound time for each strategy for each setting before dropping to the lower bound.

There are several combinations for the different settings as described below.

Results of only one intervention strategy for the endemic epidemiological zones is as follows.

(a) Optimal Protection Using ITN. Only the control () on ITNs is used to optimize the objective function , while the control on treatment (), the control on IRS (), and control on IPTp () are set to zero.

(b) Optimal Treatment. Only the control () on treatment is used to optimize the objective function , while the control on ITNs (), the control on IRS (), and control on IPTp () are set to zero.

(c) Optimal IRS. Only the control () on IRS is used to optimize the objective function , while the control on treatment (), the control on ITNs (), and control on IPTp () are set to zero.

(d) Optimal IPTp. Only the control () on IPTp is used to optimize the objective function , while the control on treatment (), the control on IRS (), and control on ITNs () are set to zero.

Results of combining 2 intervention strategies for the endemic epidemiological zones are as follows.

(a) Optimal ITNs and Treatment. With this strategy, the control on ITNs () and the treatment () are used to optimize the objective function while setting the control on IRS () and control on IPTp () to zero.

(b) Optimal ITN and IRS. With this strategy, the control on ITNs () and the IRS () are used to optimize the objective function while setting the control on treatment () and control on IPTp () to zero.

(c) Optimal ITN and IPTp. With this strategy, the control on ITNs () and IPTp () are used to optimize the objective function while setting the control on treatment () and control on IRS () to zero.

(d) Optimal Treatment and IRS. With this strategy, the control on treatment () and the IRS () are used to optimize the objective function while setting the control on ITNs () and control on IPTp () to zero.

(e) Optimal Treatment and IPTp. With this strategy, the control on treatment () and the IPTp () are used to optimize the objective function while setting the control on IRS () and control on ITNs () to zero.

(f) Optimal IRS and IPTp. With this strategy, the control on IRS () and the IPTp () are used to optimize the objective function while setting the control on treatment () and control on ITNs () to zero.

Results of combining three intervention strategies for the endemic epidemiological zones are as follows.

(a) Optimal ITN, Treatment, and IRS. In this case ITN control (), treatment control (), and IRS control () are used to optimize the objective function , while IPTp control () is set to zero.

(b) Optimal ITN, Treatment, and IPTp. In this case ITNs control (), treatment control (), and IPTp control () are used to optimize the objective function , while IRS control () is set to zero.

(c) Optimal ITN, IRS, and IPTp. In this case ITNs control (), IRS control (), and IPTp control () are used to optimize the objective function , while treatment control () is set to zero.

(d) Optimal Treatment, IRS, and IPTp. In this case treatment control (), IRS control (), and IPTp control () are used to optimize the objective function , while ITNs control () is set to zero.

Results of combining the four intervention strategies for the endemic epidemiological zones are as follows.

Optimal ITN, Treatment, IRS, and IPTp. In this case all the control function ITNs control (), treatment control (), IRS control (), and IPTp control () are used to optimize the objective function .

The same procedure is repeated for other combination for strategies and for other epidemiological zones (epidemic, seasonal, and low). Based on the simulation findings for the highest number of infections being inverted at a lower cost, the combined use of treatment and IRS reduces the infected human and mosquito population faster at a lower cost for the endemic settings (105 infections at 368.258). For the epidemic prone settings the use of treatment and IRS (111.03 infections at 388.6051) has more impact in reducing the infected human and mosquito population. For seasonal areas much impact will be felt when treatment is used (115.6983 infections at 231.3967). For the low risk areas, just the use of ITNs and treatment (119.0659 infections at 595.32) will be sufficient to reduce infected human and mosquito population. This is deduced from the intervention which takes shorter time to start reducing the number of infected mosquitoes and humans.

5. Discussion

In this paper, we formulated a mathematical model for the transmission dynamics of malaria with four time-dependent control measures in Kenya. We first consider control parameters to be constant and perform mathematical analysis of the model. The analysis showed that there exists a domain where the model is epidemiologically and mathematically well posed. Stability analysis of the model showed that the disease-free equilibrium is globally asymptotically stable if in . If , the unique endemic equilibrium exists and is globally asymptotically stable in . The model exhibited backward bifurcation backward bifurcation at implying the existence of multiple endemic equilibria for . The existence of multiple endemic equilibria emphasizes the fact that is not sufficient to eradicate disease from the population and the need to lower much below one to make the disease-free equilibrium stable globally. This behavior has important public health implications because it means that bringing below 1 is not enough to eradicate malaria.

We then consider the case of time-dependent control variable from where we formulated an optimal control problem and derived expressions for the optimal control for the malaria model with four control variables with an aim of minimizing the number of malaria infections in humans (derive optimal prevention and treatment strategies) while keeping the cost low. In the optimal control problem, use of one control at a time and the different combination of interventions can be explored to investigate and compare the effects of the control strategies on malaria eradication for different transmission settings. The analysis of the model showed that the state and optimal control can be calculated using the optimality system. The optimality system is comprised of the state system, the adjoint/costate system, initial conditions at , boundary conditions (transversality), and the characterization of the optimal control [11, 12, 14].

The results of the optimal control problem will be able to show which intervention or combination of the different intervention strategies has the highest impact on the control of the disease especially for different transmission settings. Our results shows that an effective IRS use and treatment will be beneficial to the community for the control of malaria disease (infected human and mosquito population) faster at a lower cost for the endemic settings. This is slightly different from the findings of Agusto et al. [14] who found that the combination of the personal protection, treatment, and insecticides spray had the highest impact on the control of the disease. This could be in endemic settings where both preventive and treatment measures work better which implies that the effect of protection using IRS is better. Griffin et al. [40] found that use of treatment, long-lasting insecticides treated bed nets (LLITNs), and IRS with high levels coverage would result in reducing malaria transmission for high settings though the study did not consider the cost aspect.

The findings show that, for the epidemic prone areas, the optimal control strategy for reducing the infected human and mosquito population was the use of treatment and IRS. This is slightly different from Agusto et al. [14] findings on resource limited settings in which the study recommended the use of personal protection and insecticides. This was further different from the findings of Mwamtobe et al. [12] who noted that the prevention strategies (use of ITNs and IRS) lead to the reduction of both the mosquito population and infected human individuals. This is because in epidemic areas emphasis is usually more placed on preventive strategies.

The results show that for seasonal areas much impact will be felt when treatment is used which is different from Mwamtobe et al. [12] who recommended IRS and ITNs. This is also comparable to Kim et al. [13] findings that mosquito-reduction strategies are more effective than personal protection. This is because in seasonal areas malaria transmission is usually not so high and hence if the mosquito-reduction strategies can be implemented then malaria transmission can be reduced. Griffin et al. [40] found that for the high seasonal transmission settings the use of LLITNs, IRS, and treatment would help reduce the transmission of malaria.

The results show that, for the low risk areas, just the use of ITNs and treatment will be sufficient to reduce infected human and mosquito population. This is comparable to Silva and Torres [15] who found the optimal use of ITNS would prevent malaria transmission the same as Kim et al. [13]. The findings are comparable to those by Griffin et al. [40]. In low transmission areas prevention strategies seem to be better because the population is not infected.

These findings support the WHO concerns on the capability of only one intervention strategy in reducing malaria transmission. The findings are however applicable to the designing of intervention strategies for malaria especially when costs aspects are of concern. This modelling approach also addresses effectiveness of the recommended intervention for at-risk group of malaria (pregnant women) by the WHO. The modelling approach has also been implemented in the statistical computing platform which is free statistical software.

To the best of our knowledge, this is the first ever optimal control modelling and simulation of malaria intervention strategies in free statistical computing platform; future testing and refinement of the model together with simulation with data from other representative settings should be done to improve the results and the model. These findings were based on the use of secondary data; a more designed study may be needed to ascertain the findings of these studies. Regardless, it does not invalidate the findings.

6. Conclusion

In this paper, we derived and analyzed a deterministic model for the transmission of malaria disease which incorporated the use of insecticide-treated bed nets (ITNs), treatment, indoor residual spray (IRS), and intermittent preventive treatment for pregnant women (IPTp) and performed optimal control analysis of the model. We first consider constant control parameters from where we investigate existence and stability of equilibria as well as stability analysis. We proved that if , the disease-free equilibrium is globally asymptotically stable in . If , the unique endemic equilibrium exists and is globally asymptotically stable. The model also exhibits backward bifurcation at .

We then consider the time-dependent control case from where we derived and analyzed the necessary conditions for the optimal control of the disease. There were 15 different control strategies for each of the four different epidemiological zones in Kenya that were explored using control plots (numerical simulations) which compared when there were not any intervention strategies and when there were intervention strategies. Using the optimal control approach we found that the combined use of treatment and IRS would reduce the highest number of infected human and mosquito population faster at a lower cost for the endemic settings. For the epidemic prone settings the use of treatment and IRS has more impact in reducing the infected human and mosquito population. For seasonal areas much impact will be felt when treatment is used. For the low risk areas, just the use of ITNs and treatment will be sufficient to reduce infected human and mosquito population. This was deduced from the intervention which takes shorter time to start reducing the number of infected mosquitoes and humans.

We conclude that, according to our model, the optimal control strategy for malaria control in endemic areas was the combined use of treatment and IRS; in epidemic prone areas it was the use of treatment and IRS; in seasonal areas it was the use of treatment; and in low risk areas was the use of ITNs and treatment. Control programs that follow these strategies can effectively reduce the spread of malaria disease in different malaria transmission settings in Kenya.

Appendix

The geometric approach to global stability is as described by Li and Muldowney for proving Theorem 4.

The general method considered is the one developed by Li & Muldowney [30]. Consider the autonomous dynamical system:where , is open set and simply connected and . Let be an equilibrium of (A.1); that is, . We recall that is said to be globally stable in if it is locally stable and all trajectories in converge to .

Assume that the following hypotheses hold: there exists a compact absorbing set ;: equation (A.1) has a unique equilibrium in .The basic idea of this method is that if equilibrium is (locally) stable, then the global stability is assured provided that hold and no nonconstant periodic solution of (A.1) exists.

Bendixson Criterion. Li and Muldowney showed that if hold and (A.1) satisfies a Bendixson criterion that is robust under local -perturbations of at all nonequilibrium nonwandering points for (A.1), then is globally stable in provided it is stable. Then, a new Bendixson criterion robust under local -perturbation and based on the use of the Lozinskiǐ measure is introduced.

Function is called local -perturbations of at , if there exists an open neighbourhood of in such that support and , where .

Point is said to be nonwondering for (A.1) for any neighborhood of in and there exists arbitrary large such that .

Let be a matrix-valued function that is on and considerwhere matrix is And is the second additive compound matrix of the Jacobian matrix, ; that is, . Generally speaking, for matrix, , is matrix (for a survey on compound matrices) and their relations to differential equations as described by Muldowney [41] and in the special case , one hasConsider Lozinskii measure of with respect to vector norm in , [42]:It is proved in [30] that if and hold, conditionguarantees that there are no orbits giving rise to a simple closed rectifiable curve in which is invariant for (A.1), that is, periodic orbits, homoclinic orbits, and heteroclinic cycles. In particular, condition (A.6) is proved to be a robust Bendixson criterion for (A.1). Besides, it is remarked that under assumptions , condition (A.6) also implies the local stability of .

As a consequence, the following theorem holds [30].

Theorem A.1. Assume that conditions hold. Then is globally asymptotically stable in provided that function and Lozinskiǐ measure exist such that condition (A.6) is satisfied.

Additional Points

Recommendation: in order to reduce malaria transmission in Kenya, the National Malaria Control Programme should tailor-make different intervention strategies for different epidemiological zones since some strategies work better than others in a resource limited setting. In addition the designed interventions that are preventive in nature should be implemented across all the different settings. There will be a need for the National Control Programme to create awareness on the malaria preventive measures.

Competing Interests

No competing interests exist.

Authors’ Contributions

All the authors have significant contribution to this paper and the final form of this paper is approved by all authors.

Acknowledgments

This work was (partially) funded by an African Doctoral Dissertation Research Fellowship (ADDRF) award offered by African Population and Health Research Centre (APHRC) in partnership with the International Development Research Centre (IDRC). Phase 4 of the ADDRF program is funded by IDRC (Grant no. 107508-001). This work was also supported by National Commission for Science and Technology (NACOSTI/ RCD/ST&I 6th Call PhD/225).