# Transmission Dynamics and Optimal Control of Malaria in Kenya

Authors: Gabriel Otieno; Joseph K. Koske; John M. Mutiso
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## 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 ifR
0
≤
1, the disease-free equilibrium is globally asymptotically stable in D. If R
0
>
1, the unique endemic equilibrium exists and is globally asymptotically stable. The model also exhibits backward bifurcation at R
0
=
1. If R
0
>
1, the model admits a unique endemic equilibrium which is globally asymptotically stable in the interior of feasible region D. 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.


---

## Body

## 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 byPlasmodium 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 [16–19]. 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 susceptibleS
h, Exposed humans E
h, infected humans I
h, recovered humans R
h, susceptible mosquitoes S
m, exposed mosquitoes, E
m and infected mosquitoes I
m. The total population sizes at time t for humans and mosquitoes are denoted by N
h
(
t
) and N
m
(
t
), 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 u
1
(
t
), (ii) treatment of infective humans u
2
(
t
), (iii) spray of insecticides u
3
(
t
), and (iv) treatment to protect pregnant women and their newborn children: intermittent preventive treatment for pregnant women (IPTp) u
4
(
t
). S
h
(
t
) represents the number of individuals not yet infected with the malaria parasite at time t, E
h
(
t
) represents individuals who are infected but not yet infectious, I
h
(
t
) is the class representing those who are infected with malaria and are capable of transmitting the disease to susceptible mosquitoes, and R
h
(
t
) represents the class of individuals who have temporarily recovered from the disease.Figure1 describes the dynamics of malaria in human and mosquito populations together with interventions.Figure 1
Malaria model with interventions.The susceptible humans (pregnant women and children under the age of five) (S
h) are recruited at the rate Λ
h. They either die from natural causes (at a rate μ
h) or move to the exposed class (E
h) by acquiring malaria through contact with infectious mosquitoes at a rate 1
-
u
1
(
β
ϵ
ϕ
I
m
/
N
h
)
S
h or 1
-
u
4
(
β
ϵ
ϕ
I
m
/
N
h
w
)
S
h, 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, u
1
(
t
)
∈
0,1 is the preventive measure using ITNs, u
4
(
t
)
∈
0,1 is the preventive measure using IPTp, I
m
(
t
) is the infectious mosquitoes at time t, N
h
(
t
) is the total number of individuals (pregnant women and children under the age of five), and N
h
w
(
t
) is the total number of pregnant women. Susceptible class S
h 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 α
h. Infectious individuals are assumed to recover at a rate b
+
τ
u
2, where b is the rate of spontaneous recovery, u
2
(
t
)
∈
0,1 is the control on treatment of infected individuals, and τ
∈
0,1 is the efficacy of treatment. Infectious individuals who do not recover die at a rate δ
h
+
μ
h. Individuals infected with malaria suffer a disease induced death at rate of δ
h and natural death μ
h. 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 (S
m) are recruited at the rate Λ
m and acquire malaria infection (following contact with humans infected with malaria) at the rate λ
m. They either die from natural causes (at a rate μ
m) or move to the exposed class by acquiring malaria through contacts with infected humans at a rate 1
-
u
1
(
λ
ϵ
ϕ
I
h
/
N
h
)
S
m, where λ is the probability for a vector to get infected after biting an infectious human and I
h
(
t
) are individuals infected with malaria at time t. The mosquito population is reduced, due to the use of insecticides spray, at a rate p
u
3, where u
3
(
t
)
∈
0,1 represents the control due to IRS and p represents the efficacy of IRS. Mosquito population is also reduced as a result of natural death (μ
m) and at the rate a
u
1, where u
1
(
t
) represents the control due to ITNs and a is the efficacy due to ITNs. Newly infected mosquitoes are moved into the exposed class (E
m) at a rate α
m and progresses to the class of symptomatic mosquitoes (I
m).The state variables of the model are represented and described in Table1. Table 2 describes and shows parameters of the model. Table 3 describes and represents malaria prevention and control strategies practiced in Kenya.Table 1
State variables of the malaria model.
Variable
Description
S
h
(
t
)
Number of susceptible individuals (pregnant and under 5) at timet
E
h
(
t
)
Number of exposed individuals (pregnant and under 5) at timet
I
h
(
t
)
Number of infectious humans (pregnant and under 5) at timet
R
h
(
t
)
Number of recovered humans (pregnant and under 5) at timet
S
m
(
t
)
Number of susceptible mosquitoes at timet
E
m
(
t
)
Number of exposed mosquitoes at timet
I
m
(
t
)
Number of infectious mosquitoes at timet
N
h
(
t
)
Total number of individuals (pregnant and under 5) at timet
N
h
w
(
t
)
Total number of pregnant individuals at timet
N
m
(
t
)
Total mosquito population at timetTable 2
Description of parameter variables of the malaria model.
Parameter
Description
ϕ
Mosquito contact rate with human
ϵ
Mosquito biting rate
β
Probability of human getting infected (the probability of transmission of infection from an infectious mosquito to a susceptible human provided there is a bite)
λ
Probability of a mosquito getting infected (the probability of transmission of infection from an infectious human to a susceptible mosquito provided there is a bite)
μ
h
Per capita natural death rate of humans
μ
m
Per capita natural death rate of mosquitoes
ψ
Per capita rate of loss of immunity by recovered individuals
α
h
Humans progression rate from exposed to infected
α
m
Mosquitoes progression rate from exposed to infected
Λ
h
Recruitment rate of human by birth and by getting pregnant
Λ
m
Recruitment of mosquitoes by birth
δ
h
Per capita disease induced death rate for humans (pregnant and under 5)
b
Proportion of spontaneous individual recovery
λ
h
Force of infection for susceptible humans (pregnant and under 5) to exposed individuals
λ
h
w
Force of infection for susceptible pregnant humans to exposed individuals
λ
m
Force of infection for susceptible mosquitoes to exposed mosquitoesTable 3
Prevention and control variables in the model.
Variable
Description
u
1
(
t
)
Preventive measure using insecticide-treated bed nets (ITNs)
u
2
(
t
)
The control effort on treatment of infectious individuals
u
3
(
t
)
Preventing measure using indoor residual spray (IRS)
u
4
(
t
)
Preventive measure using intermittent preventive treatment for pregnant women (IPTp)
p
Rate constant due to use of indoor residual spray
τ
Rate constant due to use of treatment effort
a
Rate constant due to use of insecticide-treated bed netsThe 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:(1)dShdt=Λh+ψRh-1-u1λhSh-1-u4λhwSh-μhSh,dEhdt=1-u1λhSh+1-u4λhwSh-αh+μhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-ψ+μhRh,dSmdt=Λm-1-u1λmSm-μm+au1+pu3Sm,dEmdt=1-u1λmSm-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im.With initial conditions (2)Sh0≥0,Eh0≥0,Ih0≥0,Rh0≥0,Sm0≥0,Em0≥0,Im0≥0,
λ
m
=
λ
ϵ
ϕ
I
h
/
N
h is the per capita incidence rate among mosquitoes (force of infection for susceptible mosquitoes), λ
h
=
β
ϵ
ϕ
I
m
/
N
h is the force of infection for susceptible humans (pregnant and under 5), and λ
h
w
=
β
ϵ
ϕ
I
m
/
N
h
w is the force of infection for susceptible pregnant humans.The total population size for the human isN
h
=
S
h
+
E
h
+
I
h
+
R
h and for mosquito is N
m
=
S
m
+
E
m
+
I
m and their differential equations are given by d
N
h
/
d
t
=
Λ
h
-
μ
h
N
h
-
δ
h
I
h and d
N
m
/
d
t
=
Λ
m
-
μ
m
+
a
u
1
+
p
u
3
N
m, 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 t
≥
0.The feasible solutions set for model (1) given by (3)
D
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
S
h
,
S
m
>
0
,
E
h
,
I
h
,
R
h
,
E
m
,
I
m
≥
0
;
N
h
≤
Λ
h
μ
h
;
N
m
≤
Λ
m
μ
m
+
a
u
1
+
p
u
3is positively invariant and hence model (1) is biologically, epidemiologically meaningful and mathematically well posed in the domain D.System (1) has always a disease-free equilibrium given by(4)E0=Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗=Λhμh,0,0,0,Λmμm+au1+pu3,0,0.


#### 2.1.2. Basic Reproduction Number

The matricesF and V for the new infection terms and the remaining transfer terms at disease-free equilibrium [26], respectively, are given by(5)F=0001-u1βϵϕ+1-u4βϵϕ000001-u1λϵϕμhΛmαhμm+au1+pu3000000,V=μh+αh000-αhδh+μh+b+τu20000αm+μm+au1+pu3000-αmμm+au1+pu3.It follows that the basic reproduction of model (1), denoted by R
0, is given by(6)R0=ρFV-1,R0=αh1-u1λϵϕΛmμh1-u1βϵϕαm+αh1-u1λϵϕΛmμh1-u4βϵϕαmμm+au1+pu3μh+αhδh+μh+b+τu2Λhαm+μm+au1+pu3μm+au1+pu3.


#### 2.1.3. Stability Analysis of Disease-Free Equilibrium Point

(
1) Local Stability of Disease-Free Equilibrium PointTheorem 1.
The disease-free equilibrium point for system (1) is locally asymptotically stable if R
0
<
1.Proof.
The Jacobian matrix(
J
) of the malaria model (1) at the disease-free equilibrium point is given by(7)-αh+μh0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu20000b+τu2-ψ+μh0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu30000αm-μm+au1+pu3.The eigenvalues of the Jacobian matrix are the solutions of the characteristic equation(8)J-λI=0,-αh+μh+λ0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu2+λ0000b+τu2-ψ+μh+λ0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu3+λ0000αm-μm+au1+pu3+λ=0.Expanding the determinant into a characteristic equation we have(9)-δh+μhαh+μh+λδh+μh+b+τu2+λαm+μm+au1+pu3+λμm+au1+pu3+λ-1-u12λϵ2ϕ2Λmμhαmβαh+1-u41-u1λϵ2ϕ2Λmμhαmβαhμm+au1+pu3Λm=0.Hence we have(10)
A
1
=
μ
m
+
a
u
1
+
p
u
3
,
A
2
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
,
A
3
=
δ
h
+
μ
h
+
b
+
τ
u
2
,
A
4
=
α
h
+
μ
h
,
Q
=
1
-
u
1
2
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
+
1
-
u
4
1
-
u
1
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
·
μ
m
+
a
u
1
+
p
u
3
Λ
m
-
1
,
λ
+
A
1
λ
+
A
2
λ
+
A
3
λ
+
A
4
-
Q
=
0
,
λ
4
+
B
1
λ
3
+
B
2
λ
2
+
B
3
λ
+
B
4
=
0
,where(11)B1=A4+A3+A2+A1,B2=A4A3+A2+A1+A3A2+A1+A1A2,B3=A4A3A2+A4A3A1+A4A2A1+A3A2A1,B4=A4A3A2A1-Q.Thus, 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 R
0
<
1.Next, we study the global behavior of the disease-free equilibrium for system (1).(
2) Global Stability of Disease-Free Equilibrium PointTheorem 2.
The DFE,E
0, of system of (1) is globally asymptotically stable if R
0
<
1.Proof.
We consider the following Lyapunov function:(12)L=c1Eh+c2Ih+c3Em+c4Im,where(13)c1=αhμm+au1+pu3αh+μhb+τu2+δh+μh,c2=1μm+au1+pu3b+τu2+δh+μh,c3=11-u1ϵϕλΛm,c4=αm+μm+au1+pu31-u1ϵϕλΛmαm.Computing the derivative of L along the solution of the system of differential equation (1),(14)
L
˙
=
α
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
E
h
+
1
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
I
h
+
1
1
-
u
1
ϵ
ϕ
λ
Λ
m
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
I
m
,
L
˙
=
α
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
E
h
+
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
I
h
+
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
I
m
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
1
β
ϵ
ϕ
α
m
+
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
4
β
ϵ
ϕ
α
m
μ
m
+
a
u
1
+
p
u
3
μ
h
+
α
h
δ
h
+
μ
h
+
b
+
τ
u
2
Λ
h
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
-
1
I
h
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
R
0
-
1
I
h
≤
0
iff
R
0
≤
1
.Thus we have established that L
˙
≤
0 if R
0
<
1 and the equality L
˙
=
0 holds if and only if R
0
=
1 and E
h
=
I
h
=
E
m
=
I
m
=
0. If R
0
>
1 then L
˙
>
0 when S
h
(
t
) and S
m
(
t
) are sufficiently close to Λ
h
/
μ
h and Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
), respectively, except when E
h
=
I
h
=
E
m
=
I
m
=
0
.
On the boundary whenE
h
=
I
h
=
E
m
=
I
m
=
0, N
h
˙
t
=
Λ
h
-
μ
h
N
h and N
m
˙
t
=
Λ
m
-
(
μ
m
+
a
u
1
+
p
u
3
) and N
h
(
t
)
→
Λ
h
/
μ
h, N
m
(
t
)
→
Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
) as t
→
∞.
Therefore the largest compact invariantD
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
L
˙
=
0 when R
0
<
1 is the singleton E
0 in D. By LaSalle’s invariant principle [28], E
0 is globally asymptotically stable for R
0
<
1 in D.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 (E
1) of the model is given by(15)E1=Sh∗∗,Eh∗∗,Ih∗∗,Rh∗∗,Sm∗∗,Em∗∗,Im∗∗.To derive E
1, we solve model (1) by equating it to zero.Substituting and solving forI
h
∗
∗ (as an expression of parameters only) through some algebraic manipulation give(16)Λh+ψb+τu2Ih∗∗ψ+μh+1-u1βϵϕNh1-u1μm+au1+pu3R0mIh∗∗1-u1λϵϕIh∗∗+Nhμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh+1-u4βϵϕNhw1-u4μm+au1+pu3R0mIh∗∗1-u4λϵϕIh∗∗+Nhwμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh-μhNhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh=0or(17)AIh∗∗2+BIh∗∗+C=0,where(18)
A
=
ψ
+
μ
h
Λ
h
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
+
ψ
b
+
τ
u
2
I
h
∗
∗
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
·
N
h
w
N
h
,
B
=
ψ
+
μ
h
1
-
u
1
β
ϵ
ϕ
1
-
u
1
·
μ
m
+
a
u
1
+
p
u
3
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
·
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
w
+
ψ
+
μ
h
1
-
u
4
β
ϵ
ϕ
1
-
u
4
μ
m
+
a
u
1
+
p
u
3
·
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
,
C
=
ψ
+
μ
h
μ
h
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
.We use the quadratic formula to find the roots of (17); that is,(19)Ih∗=-B±B2-4AC2A=-B+B2-4AC2Aor-B-B2-4AC2AIh∗=-B+B2-4AC2A=-B-B2-4AC2A=Φ.Hence(20)Sh∗∗=Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Eh∗∗=1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦNhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Ih∗∗=Φ=αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh≥0,Rh∗∗=b+τu2αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦμh+ψδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Sm∗∗=ΛmNh1-u1λϵϕΦ+μm+au1+pu3Nh,Em∗∗=1-u1αhλϵϕ1-u1μm+au1+pu3R0mΦΛmNhNhαm+μm+au1+pu31-u1λϵϕΦ+Nhμm+au1+pu31-u1λϵϕΦ+μm+au1+pu3Nh,Im∗∗=1-u1μm+au1+pu3R0mΦ1-u1λϵϕΦ+Nhμm+au1+pu3.From the quadratic equation (17) we analyze the possibility of multiple endemic equilibria. It is important to note that the coefficient A is always positive with B and C having different signs. All the negative terms in C are in the form of (
1
-
u
) where u is the control that is bounded by 1.It follows that(i)
there is a unique endemic equilibrium ifB
<
0 and C
=
0 or B
2
-
4
A
C
=
0,
(ii)
there is a unique endemic equilibrium ifC
<
0 (i.e., if R
0
>
1),
(iii)
there are two endemic equilibria ifC
>
0, B
<
0, and B
2
-
4
A
C
>
0,
(iv)
there are no endemic equilibria otherwise.Note that the hypothesis C
>
0 is equivalent to R
0
<
1.Thus the results of this section can be summarized in the following theorem.Theorem 3.
IfR
0
<
1, E
0 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 R
0
>
1, then E
0 is unstable and there exists a unique endemic equilibrium.Item (iii) indicates the possibility of backward bifurcation in model (1) when R
0
<
1. In the next section we will prove the occurrence of multiple equilibria for R
0
<
1 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 outbifurcation analysis of system (1) at R
0
=
1. We consider a transmission rate β as bifurcation parameter Ψ so that R
0
=
1.To apply the theory, we introduce dimensionless state variables into system (1).The system of (1) can be written as (21)dShdt=Λh+ψRh-μhSh-1-u1βϵϕShImNh-1-u4βϵϕShImNhw,dEhdt=1-u1βϵϕShImNh+1-u4βϵϕShImNhw-μhEh-αhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-μhRh-ψRh,dSmdt=Λm-1-u1λϵϕIhSmNh-μm+au1+pu3Sm,dEmdt=1-u1λϵϕIhSmNh-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im.We let x
1
=
S
h, x
2
=
E
h, x
3
=
I
h, x
4
=
R
h, x
5
=
S
m, x
6
=
E
m, and x
7
=
I
m
.Therefore system (1) in vector form can be written as(22)dXidt=HXi,where X
i
=
x
1
,
x
2
,
…
,
x
7
T and H
=
h
1
,
h
2
,
…
,
h
7
T are transposed matrices and N
h
=
Λ
h
/
μ
h with Ψ
∗
=
β, (23)dx1dt=Λh+ψx4-μhSh-1-u1Ψ∗ϵϕx7x1μhΛh-1-u4Ψ∗ϵϕx7x1μhΛh=h1,dx2dt=1-u1Ψ∗ϵϕx7x1μhΛh+1-u4Ψ∗ϵϕx7x1μhΛh-μh+αhx2=h2,dx3dt=αhx2-δh+μh+b+τu2x3=h3,dx4dt=b+τu2x3-μh+ψx4=h4,dx5dt=Λm-1-u1λϵϕx3x5μhΛh-μm+au1+pu3x5=h5,dx6dt=1-u1λϵϕx3x5μhΛh-αm+μm+au1+pu3x6=h6,dx7dt=αmx6-μm+au1+pu3x7=h7.Let Ψ
∗ be the bifurcation parameter; the system is linearized at disease-free equilibrium point when β
=
Ψ
∗ with R
0
=
1. That is, (24)R0=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,12=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,Ψ∗=Λhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2αhαmΛmμh1-u12ϕ2ϵλ+αhαmΛmμh1-u11-u4ϕ2ϵ2.The Jacobian matrix of (1) calculated at Ψ
∗ is given by (25)
-
μ
h
0
0
ψ
0
0
-
Ψ
∗
ϵ
ϕ
0
-
α
h
-
μ
h
0
0
0
0
Ψ
∗
ϵ
ϕ
0
α
h
-
δ
h
-
μ
h
-
b
-
τ
u
2
0
0
0
0
0
0
b
+
τ
u
2
-
μ
h
-
ψ
0
0
0
0
0
-
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
-
μ
m
+
a
u
1
+
p
u
3
0
0
0
0
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
0
0
0
0
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
.Centre manifold approach [29] is then applied.A right eigenvector associated with the eigenvalue zero isw
=
w
1
,
w
2
,
…
,
w
7. Solving the system we have the following right eigenvector:(26)w1=ψw4-Ψ∗ϵϕw7μh,w2=Ψ∗ϵϕw7αh+μh,w3=αhw2b+τu2+μh+δh,w4=b+τu2w3μh+ψ,w5=-1-u1λϵϕΛmμhw3Λhμm+au1+pu32,w6=1-u1λϵϕΛmμhw3Λhμm+au1+pu3αm+μm+au1+pu3,w7=αmw6μm+au1+pu3>0.The left eigenvectors satisfying v
·
w
=
1 are v
=
v
1
,
v
2
,
…
,
v
7. Solving the system, the left eigenvector will be given by(27)v1=0,v2=αhv3αh+μh,v3=v61-u1λϵϕΛmμhΛhμm+au1+pu3-δh-uh-b-τu2,v4=0,v5=0,v6=αmv7-μm-αm-au1-pu3,v7=Ψ∗ϵϕv2μm+au1+pu3.Computing for the sign of a and b as indicated in the theorem gives(28)
a
=
v
2
w
7
-
1
-
u
1
Ψ
∗
ϵ
ϕ
μ
h
-
1
-
u
4
Ψ
∗
ϵ
ϕ
μ
h
Λ
h
·
Ψ
∗
ϵ
w
7
α
h
+
μ
h
+
α
h
w
2
b
+
τ
u
2
+
μ
h
+
δ
h
+
v
6
w
3
-
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
h
·
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
m
w
3
Λ
h
μ
m
+
a
u
1
+
p
u
3
α
m
+
μ
m
+
a
u
1
+
p
u
3
+
α
m
w
3
μ
m
+
a
u
1
+
p
u
3
b
=
v
2
w
7
ϵ
ϕ
>
0so that b is always positive.Therefore the following result is established.Theorem 4.
Model (1) exhibits backward bifurcation at R
0
=
1 whenever a
>
0 and b
>
0 and R
0
<
1. Whenever a
<
0 and b
>
0, then model (1) exhibits a forward bifurcation at R
0
=
1.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 R
0
>
1 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 D. 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, E
1: (i) the uniqueness of E
1 in the interior of D (condition H
1); (ii) the existence of an absorbing compact set in the interior of D (condition H
2); and (iii) the fulfillment of a Bendixson criterion (i.e., inequality (A.6)).Theorem 5.
Ifa
<
0 and b
>
0 and R
0
>
1, then the endemic equilibrium of the malaria model (1) is globally asymptotically stable in the interior of D.Proof.
Following Li and Muldowney [30], for system (1), under the assumption of R
0
>
1, satisfies conditions H
1
-
(
H
2
), 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 c
>
0 such that any solutions (
S
h
t
,
I
h
t
,
I
m
t
) with (
S
h
0
,
I
h
0
,
I
m
0
) in the interior of D satisfy(29)min⁡limt→∞⁡inf⁡Sht,limt→∞⁡inf⁡Iht,limt→∞⁡inf⁡Imt.The uniform persistence together with boundedness of D is equivalent to the existence of a compact set in the interior of D which is absorbing for (4) as described by Hutson and Schmitt [32]. Thus, (
H
1
) is verified. Moreover, E
1 is the only equilibrium in the interior of D, so that (
H
2
) 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 S
h
,
I
h
,
I
m of reduced system (1), we get the second additive compound matrix J
2:(30)J2Sh,Ih,Im=-a11a121-u1ϵϕβSh+1-u4ϵϕβSh1-u1ϵϕλNm-Em-Im-a22-δh01-u1ϵϕβIm+1-u4ϵϕβIm-a33,where (31)a11=μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+δh+b+τu2,a12=1-u1ϵϕβSh+1-u4ϵϕβSh,a22=μh+δh+1-u1ϵϕβIm+1-u4ϵϕβIm+1-u1ϵϕλIh+μm+au1+pu3,a33=m+μm+au1+pu3+1-u1ϵϕλIh+1-u1ϵϕβIm+1-u4ϵϕβIm,where m
=
Λ
m
/
μ
m
+
a
u
1
+
p
u
3
/
Λ
h
/
μ
h
.
Choose now matrixP
(
x
)
=
P
S
h
,
I
h
,
I
m
=
diag
⁡
(
1
,
I
h
/
I
m
,
I
h
/
I
m
). Then P
f
P
-
1
=
diag
⁡
(
0
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
), and the matrix B
=
P
f
P
-
1
+
P
J
2
P
-
1 can be written in block form as(32)B=B11B12B21B22,where(33)B11=-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+b+τu2+δh,B12=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh,1-u1ϵϕβSh·ImIh+1-u4ϵϕβSh·ImIh,B21=1-u1ϵλϕImNm-Em+ImIhImT,B22=Ih˙Ih-Im˙Im-a2201-u1ϵϕβIm-1-u4ϵϕβIm-1-u1ϵϕλIhIh˙Ih-Im˙Im-a33.The vector norm · in R
+
3 is here chosen to be (34)x,y,z=max⁡x,y,z.Let σ
(
·
) denote the Lozinskii measure with respect to this norm. Using the method of estimating σ
(
·
) in [29], we have(35)σB≤sup⁡g1,g2=sup⁡σ1B11+B12,σ1B22+B21,where B
12 and B
21 are matrix norms with respect to the L vector norm and σ
1 denotes the Lozinskii measure with respect to L norm. Since B
11 is a scalar, its Lozinskii measure with respect to any norm in R
1 is equal to B
11.
Therefore(36)
σ
1
B
11
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
·
ϵ
ϕ
β
I
m
+
μ
h
+
b
+
τ
u
2
+
δ
h
,
σ
1
B
22
=
max
⁡
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
I
h
˙
I
h
-
I
m
˙
I
m
-
1
-
u
1
ϵ
ϕ
β
I
m
-
1
-
u
4
ϵ
ϕ
β
I
m
-
1
-
u
1
ϵ
ϕ
λ
I
h
-
μ
h
-
δ
h
-
b
-
τ
u
2
-
μ
m
-
a
u
1
-
p
u
3
,
σ
1
B
22
=
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
B
12
=
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
B
21
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
.Therefore(37)
g
1
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
μ
h
+
δ
h
+
b
+
τ
u
2
+
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
(38)
g
2
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
+
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
μ
m
+
a
u
1
+
p
u
3
+
1
-
u
1
ϵ
ϕ
λ
I
h
.We rewrite the last two equations of system (1) for I
h
˙ and I
m
˙ as (39)Ih˙Ih=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh-μh+δh+b+τu2,(40)Im˙Im=1-u1ϵϕλNm-Em-Im·IhIm-μm-au1-pu3.Substituting (39) into (37) and (40) into (38) we have(41)g1t=Ih˙Ih-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm.g2t=Ih˙Ih-μh+αh+1-u1ϵϕλIh.For the uniform persistence constant ε
>
0, there exists a time T
0
>
0 independent of x
(
0
)
∈
K, the compact absorbing set, such that(42)Iht>ε,Imt>ε;for t
>
T
0 we have(43)g1t≤Ih˙Ih-μh+αh+1-u1ϵϕβε+1-u4ϵϕβε,g2t≤Ih˙Ih-μh+αh+1-u1ϵϕλε.Relations (41) imply (44)σB≤Ih˙Ih-μfort>T0,where(45)
μ
=
min
⁡
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
ε
+
1
-
u
4
ϵ
ϕ
β
ε
,
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
λ
ε
.Along each solution S
h
t
,
I
h
t
,
I
m
(
t
) to (1) with S
h
0
,
I
h
0
,
I
m
(
0
)
∈
K where K is the compact absorbing set, we have, for t
>
T
0,(46)1t∫0tσBds≤1t∫0T0σBds+1tln⁡IhtIhT0-μt-T0t.Which implies q
-
2
<
μ
/
2
<
0. 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 Table4 is computed. The index measures the relative change in a variable with respect to relative changes in parameters.Table 4
Sensitivity indices (SI) ofR
0 to parameters for the malaria model.
Parameter
Sensitivity indices
Endemic
Seasonal
Epidemic
Low risk
μ
_
h
−0.0402531
−0.0402531
−0.0402531
−0.0402531
μ
_
m
−1.07211
−1.07211
−1.07211
−1.07211
α
_
h
0.00038817
0.00038817
0.00038817
0.00038817
α
_
m
0.22445
0.22445
0.22445
0.22445
λ
0.5
0.5
0.5
0.5
β
0.5
0.5
0.5
0.5
ϵ
1
1
1
1
Λ
_
h
−0.4987
−0.4987
−0.4987
−0.4987
Λ
_
m
0.5
0.5
0.5
0.5
b
−0.01818
−0.02048
−0.01639
−0.02563
τ
−0.322497
−0.322497
−0.322497
−0.322497
δ
_
h
−0.13695
−0.10336
−0.11321
−0.03508
ϕ
1
1
1
1Definition 6.
Following Chitnis et al. [33], the normalized forward sensitivity index of a variable, h, that depends on a parameter, l, is defined as ξ
l
=
(
l
/
h
)
×
(
∂
h
/
∂
l
).
The sensitivity index of the model parameters is given by(47)ξαhR0=μh2αh+μh,ξαmR0=pu3+μm+au12αm+pu3+au1+μm,ξμmR0=-μm2αm+3μm+au1+pu32αm+μm+au1+pu3μm+au1+pu3,ξδhR0=-δh2b+τu2+δh+μh,ξbR0=-b2b+τu2+δh+μh,ξμhR0=-μh2+αhδh+αhb+αhτu22μh+δh+b+τu2μh+αh.Sensitivity indices for the control parameters are given by(48)ξu1R0=-μ11-μ1,ξu2R0=-τμ22μh+ϵ+b+τμ2,ξu3R0=-p3pu3+3μm+au1+2αmu32pu3+au1+μmpu3+μm+αm,ξu4R0=-u41-u4.Through sensitivity analysis, it is observed that the most sensitive parameters to R
0 across all the settings were the mosquito’s natural death rate, μ
m, 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
R
0. Sensitivity indices are calculated using parameters in Table 5.Table 5
Parameter values for the full malaria model.
Parameter
Estimated value
Source
Endemic
Epidemic
Seasonal
Low risk
μ
h
0.00005447
0.00004644
0.00004281
0.00004566
KNBS (2009 Census estimates) [20]
μ
m
0.04
Estimated
α
h
0.07143
Estimated
α
m
0.0909
Chitnis et al. [21]
λ
0.42
Estimated
β
0.0655
Estimated
ϵ
0.2
Blayneh et al. [22]
ψ
0.01095
Estimated
Λ
h
0.00000575
0.00000575
0.00000548
0.00000438
KNBS (estimates based on 2009 Census) [20]
Λ
m
0.071
Niger and Gumel [23]
b
0.005
Chiyaka et al. [24]
τ
0.5
Assumed
δ
h
0.05
KNBS and ICF Macro [25]
p
0.25
Assumed
a
0.25
Assumed
ϕ
0.502
Blayneh et al. [22]
λ
h
0.00000149
0.00000123
0.00000445
0.00000226
Estimated
λ
h
w
0.00000247
0.00000203
0.00000693
0.00000328
Estimated
λ
m
0.00000048
0.00000394
0.00000143
0.00000073
Estimated
N
h
440169
535093
147832
290703
KNBS (2009 Census estimate) [20]
N
h
w
266343
324658
94857
200216
KNBS (2009 Census estimate) [20]
N
m
4401690
5350930
1478320
2907030
Estimated


## 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 t
≥
0.The feasible solutions set for model (1) given by (3)
D
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
S
h
,
S
m
>
0
,
E
h
,
I
h
,
R
h
,
E
m
,
I
m
≥
0
;
N
h
≤
Λ
h
μ
h
;
N
m
≤
Λ
m
μ
m
+
a
u
1
+
p
u
3is positively invariant and hence model (1) is biologically, epidemiologically meaningful and mathematically well posed in the domain D.System (1) has always a disease-free equilibrium given by(4)E0=Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗=Λhμh,0,0,0,Λmμm+au1+pu3,0,0.


### 2.1.2. Basic Reproduction Number

The matricesF and V for the new infection terms and the remaining transfer terms at disease-free equilibrium [26], respectively, are given by(5)F=0001-u1βϵϕ+1-u4βϵϕ000001-u1λϵϕμhΛmαhμm+au1+pu3000000,V=μh+αh000-αhδh+μh+b+τu20000αm+μm+au1+pu3000-αmμm+au1+pu3.It follows that the basic reproduction of model (1), denoted by R
0, is given by(6)R0=ρFV-1,R0=αh1-u1λϵϕΛmμh1-u1βϵϕαm+αh1-u1λϵϕΛmμh1-u4βϵϕαmμm+au1+pu3μh+αhδh+μh+b+τu2Λhαm+μm+au1+pu3μm+au1+pu3.


### 2.1.3. Stability Analysis of Disease-Free Equilibrium Point

(
1) Local Stability of Disease-Free Equilibrium PointTheorem 1.
The disease-free equilibrium point for system (1) is locally asymptotically stable if R
0
<
1.Proof.
The Jacobian matrix(
J
) of the malaria model (1) at the disease-free equilibrium point is given by(7)-αh+μh0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu20000b+τu2-ψ+μh0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu30000αm-μm+au1+pu3.The eigenvalues of the Jacobian matrix are the solutions of the characteristic equation(8)J-λI=0,-αh+μh+λ0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu2+λ0000b+τu2-ψ+μh+λ0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu3+λ0000αm-μm+au1+pu3+λ=0.Expanding the determinant into a characteristic equation we have(9)-δh+μhαh+μh+λδh+μh+b+τu2+λαm+μm+au1+pu3+λμm+au1+pu3+λ-1-u12λϵ2ϕ2Λmμhαmβαh+1-u41-u1λϵ2ϕ2Λmμhαmβαhμm+au1+pu3Λm=0.Hence we have(10)
A
1
=
μ
m
+
a
u
1
+
p
u
3
,
A
2
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
,
A
3
=
δ
h
+
μ
h
+
b
+
τ
u
2
,
A
4
=
α
h
+
μ
h
,
Q
=
1
-
u
1
2
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
+
1
-
u
4
1
-
u
1
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
·
μ
m
+
a
u
1
+
p
u
3
Λ
m
-
1
,
λ
+
A
1
λ
+
A
2
λ
+
A
3
λ
+
A
4
-
Q
=
0
,
λ
4
+
B
1
λ
3
+
B
2
λ
2
+
B
3
λ
+
B
4
=
0
,where(11)B1=A4+A3+A2+A1,B2=A4A3+A2+A1+A3A2+A1+A1A2,B3=A4A3A2+A4A3A1+A4A2A1+A3A2A1,B4=A4A3A2A1-Q.Thus, 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 R
0
<
1.Next, we study the global behavior of the disease-free equilibrium for system (1).(
2) Global Stability of Disease-Free Equilibrium PointTheorem 2.
The DFE,E
0, of system of (1) is globally asymptotically stable if R
0
<
1.Proof.
We consider the following Lyapunov function:(12)L=c1Eh+c2Ih+c3Em+c4Im,where(13)c1=αhμm+au1+pu3αh+μhb+τu2+δh+μh,c2=1μm+au1+pu3b+τu2+δh+μh,c3=11-u1ϵϕλΛm,c4=αm+μm+au1+pu31-u1ϵϕλΛmαm.Computing the derivative of L along the solution of the system of differential equation (1),(14)
L
˙
=
α
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
E
h
+
1
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
I
h
+
1
1
-
u
1
ϵ
ϕ
λ
Λ
m
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
I
m
,
L
˙
=
α
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
E
h
+
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
I
h
+
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
I
m
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
1
β
ϵ
ϕ
α
m
+
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
4
β
ϵ
ϕ
α
m
μ
m
+
a
u
1
+
p
u
3
μ
h
+
α
h
δ
h
+
μ
h
+
b
+
τ
u
2
Λ
h
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
-
1
I
h
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
R
0
-
1
I
h
≤
0
iff
R
0
≤
1
.Thus we have established that L
˙
≤
0 if R
0
<
1 and the equality L
˙
=
0 holds if and only if R
0
=
1 and E
h
=
I
h
=
E
m
=
I
m
=
0. If R
0
>
1 then L
˙
>
0 when S
h
(
t
) and S
m
(
t
) are sufficiently close to Λ
h
/
μ
h and Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
), respectively, except when E
h
=
I
h
=
E
m
=
I
m
=
0
.
On the boundary whenE
h
=
I
h
=
E
m
=
I
m
=
0, N
h
˙
t
=
Λ
h
-
μ
h
N
h and N
m
˙
t
=
Λ
m
-
(
μ
m
+
a
u
1
+
p
u
3
) and N
h
(
t
)
→
Λ
h
/
μ
h, N
m
(
t
)
→
Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
) as t
→
∞.
Therefore the largest compact invariantD
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
L
˙
=
0 when R
0
<
1 is the singleton E
0 in D. By LaSalle’s invariant principle [28], E
0 is globally asymptotically stable for R
0
<
1 in D.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 (E
1) of the model is given by(15)E1=Sh∗∗,Eh∗∗,Ih∗∗,Rh∗∗,Sm∗∗,Em∗∗,Im∗∗.To derive E
1, we solve model (1) by equating it to zero.Substituting and solving forI
h
∗
∗ (as an expression of parameters only) through some algebraic manipulation give(16)Λh+ψb+τu2Ih∗∗ψ+μh+1-u1βϵϕNh1-u1μm+au1+pu3R0mIh∗∗1-u1λϵϕIh∗∗+Nhμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh+1-u4βϵϕNhw1-u4μm+au1+pu3R0mIh∗∗1-u4λϵϕIh∗∗+Nhwμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh-μhNhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh=0or(17)AIh∗∗2+BIh∗∗+C=0,where(18)
A
=
ψ
+
μ
h
Λ
h
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
+
ψ
b
+
τ
u
2
I
h
∗
∗
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
·
N
h
w
N
h
,
B
=
ψ
+
μ
h
1
-
u
1
β
ϵ
ϕ
1
-
u
1
·
μ
m
+
a
u
1
+
p
u
3
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
·
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
w
+
ψ
+
μ
h
1
-
u
4
β
ϵ
ϕ
1
-
u
4
μ
m
+
a
u
1
+
p
u
3
·
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
,
C
=
ψ
+
μ
h
μ
h
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
.We use the quadratic formula to find the roots of (17); that is,(19)Ih∗=-B±B2-4AC2A=-B+B2-4AC2Aor-B-B2-4AC2AIh∗=-B+B2-4AC2A=-B-B2-4AC2A=Φ.Hence(20)Sh∗∗=Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Eh∗∗=1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦNhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Ih∗∗=Φ=αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh≥0,Rh∗∗=b+τu2αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦμh+ψδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Sm∗∗=ΛmNh1-u1λϵϕΦ+μm+au1+pu3Nh,Em∗∗=1-u1αhλϵϕ1-u1μm+au1+pu3R0mΦΛmNhNhαm+μm+au1+pu31-u1λϵϕΦ+Nhμm+au1+pu31-u1λϵϕΦ+μm+au1+pu3Nh,Im∗∗=1-u1μm+au1+pu3R0mΦ1-u1λϵϕΦ+Nhμm+au1+pu3.From the quadratic equation (17) we analyze the possibility of multiple endemic equilibria. It is important to note that the coefficient A is always positive with B and C having different signs. All the negative terms in C are in the form of (
1
-
u
) where u is the control that is bounded by 1.It follows that(i)
there is a unique endemic equilibrium ifB
<
0 and C
=
0 or B
2
-
4
A
C
=
0,
(ii)
there is a unique endemic equilibrium ifC
<
0 (i.e., if R
0
>
1),
(iii)
there are two endemic equilibria ifC
>
0, B
<
0, and B
2
-
4
A
C
>
0,
(iv)
there are no endemic equilibria otherwise.Note that the hypothesis C
>
0 is equivalent to R
0
<
1.Thus the results of this section can be summarized in the following theorem.Theorem 3.
IfR
0
<
1, E
0 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 R
0
>
1, then E
0 is unstable and there exists a unique endemic equilibrium.Item (iii) indicates the possibility of backward bifurcation in model (1) when R
0
<
1. In the next section we will prove the occurrence of multiple equilibria for R
0
<
1 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 outbifurcation analysis of system (1) at R
0
=
1. We consider a transmission rate β as bifurcation parameter Ψ so that R
0
=
1.To apply the theory, we introduce dimensionless state variables into system (1).The system of (1) can be written as (21)dShdt=Λh+ψRh-μhSh-1-u1βϵϕShImNh-1-u4βϵϕShImNhw,dEhdt=1-u1βϵϕShImNh+1-u4βϵϕShImNhw-μhEh-αhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-μhRh-ψRh,dSmdt=Λm-1-u1λϵϕIhSmNh-μm+au1+pu3Sm,dEmdt=1-u1λϵϕIhSmNh-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im.We let x
1
=
S
h, x
2
=
E
h, x
3
=
I
h, x
4
=
R
h, x
5
=
S
m, x
6
=
E
m, and x
7
=
I
m
.Therefore system (1) in vector form can be written as(22)dXidt=HXi,where X
i
=
x
1
,
x
2
,
…
,
x
7
T and H
=
h
1
,
h
2
,
…
,
h
7
T are transposed matrices and N
h
=
Λ
h
/
μ
h with Ψ
∗
=
β, (23)dx1dt=Λh+ψx4-μhSh-1-u1Ψ∗ϵϕx7x1μhΛh-1-u4Ψ∗ϵϕx7x1μhΛh=h1,dx2dt=1-u1Ψ∗ϵϕx7x1μhΛh+1-u4Ψ∗ϵϕx7x1μhΛh-μh+αhx2=h2,dx3dt=αhx2-δh+μh+b+τu2x3=h3,dx4dt=b+τu2x3-μh+ψx4=h4,dx5dt=Λm-1-u1λϵϕx3x5μhΛh-μm+au1+pu3x5=h5,dx6dt=1-u1λϵϕx3x5μhΛh-αm+μm+au1+pu3x6=h6,dx7dt=αmx6-μm+au1+pu3x7=h7.Let Ψ
∗ be the bifurcation parameter; the system is linearized at disease-free equilibrium point when β
=
Ψ
∗ with R
0
=
1. That is, (24)R0=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,12=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,Ψ∗=Λhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2αhαmΛmμh1-u12ϕ2ϵλ+αhαmΛmμh1-u11-u4ϕ2ϵ2.The Jacobian matrix of (1) calculated at Ψ
∗ is given by (25)
-
μ
h
0
0
ψ
0
0
-
Ψ
∗
ϵ
ϕ
0
-
α
h
-
μ
h
0
0
0
0
Ψ
∗
ϵ
ϕ
0
α
h
-
δ
h
-
μ
h
-
b
-
τ
u
2
0
0
0
0
0
0
b
+
τ
u
2
-
μ
h
-
ψ
0
0
0
0
0
-
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
-
μ
m
+
a
u
1
+
p
u
3
0
0
0
0
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
0
0
0
0
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
.Centre manifold approach [29] is then applied.A right eigenvector associated with the eigenvalue zero isw
=
w
1
,
w
2
,
…
,
w
7. Solving the system we have the following right eigenvector:(26)w1=ψw4-Ψ∗ϵϕw7μh,w2=Ψ∗ϵϕw7αh+μh,w3=αhw2b+τu2+μh+δh,w4=b+τu2w3μh+ψ,w5=-1-u1λϵϕΛmμhw3Λhμm+au1+pu32,w6=1-u1λϵϕΛmμhw3Λhμm+au1+pu3αm+μm+au1+pu3,w7=αmw6μm+au1+pu3>0.The left eigenvectors satisfying v
·
w
=
1 are v
=
v
1
,
v
2
,
…
,
v
7. Solving the system, the left eigenvector will be given by(27)v1=0,v2=αhv3αh+μh,v3=v61-u1λϵϕΛmμhΛhμm+au1+pu3-δh-uh-b-τu2,v4=0,v5=0,v6=αmv7-μm-αm-au1-pu3,v7=Ψ∗ϵϕv2μm+au1+pu3.Computing for the sign of a and b as indicated in the theorem gives(28)
a
=
v
2
w
7
-
1
-
u
1
Ψ
∗
ϵ
ϕ
μ
h
-
1
-
u
4
Ψ
∗
ϵ
ϕ
μ
h
Λ
h
·
Ψ
∗
ϵ
w
7
α
h
+
μ
h
+
α
h
w
2
b
+
τ
u
2
+
μ
h
+
δ
h
+
v
6
w
3
-
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
h
·
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
m
w
3
Λ
h
μ
m
+
a
u
1
+
p
u
3
α
m
+
μ
m
+
a
u
1
+
p
u
3
+
α
m
w
3
μ
m
+
a
u
1
+
p
u
3
b
=
v
2
w
7
ϵ
ϕ
>
0so that b is always positive.Therefore the following result is established.Theorem 4.
Model (1) exhibits backward bifurcation at R
0
=
1 whenever a
>
0 and b
>
0 and R
0
<
1. Whenever a
<
0 and b
>
0, then model (1) exhibits a forward bifurcation at R
0
=
1.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 R
0
>
1 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 D. 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, E
1: (i) the uniqueness of E
1 in the interior of D (condition H
1); (ii) the existence of an absorbing compact set in the interior of D (condition H
2); and (iii) the fulfillment of a Bendixson criterion (i.e., inequality (A.6)).Theorem 5.
Ifa
<
0 and b
>
0 and R
0
>
1, then the endemic equilibrium of the malaria model (1) is globally asymptotically stable in the interior of D.Proof.
Following Li and Muldowney [30], for system (1), under the assumption of R
0
>
1, satisfies conditions H
1
-
(
H
2
), 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 c
>
0 such that any solutions (
S
h
t
,
I
h
t
,
I
m
t
) with (
S
h
0
,
I
h
0
,
I
m
0
) in the interior of D satisfy(29)min⁡limt→∞⁡inf⁡Sht,limt→∞⁡inf⁡Iht,limt→∞⁡inf⁡Imt.The uniform persistence together with boundedness of D is equivalent to the existence of a compact set in the interior of D which is absorbing for (4) as described by Hutson and Schmitt [32]. Thus, (
H
1
) is verified. Moreover, E
1 is the only equilibrium in the interior of D, so that (
H
2
) 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 S
h
,
I
h
,
I
m of reduced system (1), we get the second additive compound matrix J
2:(30)J2Sh,Ih,Im=-a11a121-u1ϵϕβSh+1-u4ϵϕβSh1-u1ϵϕλNm-Em-Im-a22-δh01-u1ϵϕβIm+1-u4ϵϕβIm-a33,where (31)a11=μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+δh+b+τu2,a12=1-u1ϵϕβSh+1-u4ϵϕβSh,a22=μh+δh+1-u1ϵϕβIm+1-u4ϵϕβIm+1-u1ϵϕλIh+μm+au1+pu3,a33=m+μm+au1+pu3+1-u1ϵϕλIh+1-u1ϵϕβIm+1-u4ϵϕβIm,where m
=
Λ
m
/
μ
m
+
a
u
1
+
p
u
3
/
Λ
h
/
μ
h
.
Choose now matrixP
(
x
)
=
P
S
h
,
I
h
,
I
m
=
diag
⁡
(
1
,
I
h
/
I
m
,
I
h
/
I
m
). Then P
f
P
-
1
=
diag
⁡
(
0
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
), and the matrix B
=
P
f
P
-
1
+
P
J
2
P
-
1 can be written in block form as(32)B=B11B12B21B22,where(33)B11=-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+b+τu2+δh,B12=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh,1-u1ϵϕβSh·ImIh+1-u4ϵϕβSh·ImIh,B21=1-u1ϵλϕImNm-Em+ImIhImT,B22=Ih˙Ih-Im˙Im-a2201-u1ϵϕβIm-1-u4ϵϕβIm-1-u1ϵϕλIhIh˙Ih-Im˙Im-a33.The vector norm · in R
+
3 is here chosen to be (34)x,y,z=max⁡x,y,z.Let σ
(
·
) denote the Lozinskii measure with respect to this norm. Using the method of estimating σ
(
·
) in [29], we have(35)σB≤sup⁡g1,g2=sup⁡σ1B11+B12,σ1B22+B21,where B
12 and B
21 are matrix norms with respect to the L vector norm and σ
1 denotes the Lozinskii measure with respect to L norm. Since B
11 is a scalar, its Lozinskii measure with respect to any norm in R
1 is equal to B
11.
Therefore(36)
σ
1
B
11
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
·
ϵ
ϕ
β
I
m
+
μ
h
+
b
+
τ
u
2
+
δ
h
,
σ
1
B
22
=
max
⁡
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
I
h
˙
I
h
-
I
m
˙
I
m
-
1
-
u
1
ϵ
ϕ
β
I
m
-
1
-
u
4
ϵ
ϕ
β
I
m
-
1
-
u
1
ϵ
ϕ
λ
I
h
-
μ
h
-
δ
h
-
b
-
τ
u
2
-
μ
m
-
a
u
1
-
p
u
3
,
σ
1
B
22
=
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
B
12
=
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
B
21
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
.Therefore(37)
g
1
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
μ
h
+
δ
h
+
b
+
τ
u
2
+
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
(38)
g
2
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
+
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
μ
m
+
a
u
1
+
p
u
3
+
1
-
u
1
ϵ
ϕ
λ
I
h
.We rewrite the last two equations of system (1) for I
h
˙ and I
m
˙ as (39)Ih˙Ih=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh-μh+δh+b+τu2,(40)Im˙Im=1-u1ϵϕλNm-Em-Im·IhIm-μm-au1-pu3.Substituting (39) into (37) and (40) into (38) we have(41)g1t=Ih˙Ih-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm.g2t=Ih˙Ih-μh+αh+1-u1ϵϕλIh.For the uniform persistence constant ε
>
0, there exists a time T
0
>
0 independent of x
(
0
)
∈
K, the compact absorbing set, such that(42)Iht>ε,Imt>ε;for t
>
T
0 we have(43)g1t≤Ih˙Ih-μh+αh+1-u1ϵϕβε+1-u4ϵϕβε,g2t≤Ih˙Ih-μh+αh+1-u1ϵϕλε.Relations (41) imply (44)σB≤Ih˙Ih-μfort>T0,where(45)
μ
=
min
⁡
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
ε
+
1
-
u
4
ϵ
ϕ
β
ε
,
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
λ
ε
.Along each solution S
h
t
,
I
h
t
,
I
m
(
t
) to (1) with S
h
0
,
I
h
0
,
I
m
(
0
)
∈
K where K is the compact absorbing set, we have, for t
>
T
0,(46)1t∫0tσBds≤1t∫0T0σBds+1tln⁡IhtIhT0-μt-T0t.Which implies q
-
2
<
μ
/
2
<
0. 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 Table4 is computed. The index measures the relative change in a variable with respect to relative changes in parameters.Table 4
Sensitivity indices (SI) ofR
0 to parameters for the malaria model.
Parameter
Sensitivity indices
Endemic
Seasonal
Epidemic
Low risk
μ
_
h
−0.0402531
−0.0402531
−0.0402531
−0.0402531
μ
_
m
−1.07211
−1.07211
−1.07211
−1.07211
α
_
h
0.00038817
0.00038817
0.00038817
0.00038817
α
_
m
0.22445
0.22445
0.22445
0.22445
λ
0.5
0.5
0.5
0.5
β
0.5
0.5
0.5
0.5
ϵ
1
1
1
1
Λ
_
h
−0.4987
−0.4987
−0.4987
−0.4987
Λ
_
m
0.5
0.5
0.5
0.5
b
−0.01818
−0.02048
−0.01639
−0.02563
τ
−0.322497
−0.322497
−0.322497
−0.322497
δ
_
h
−0.13695
−0.10336
−0.11321
−0.03508
ϕ
1
1
1
1Definition 6.
Following Chitnis et al. [33], the normalized forward sensitivity index of a variable, h, that depends on a parameter, l, is defined as ξ
l
=
(
l
/
h
)
×
(
∂
h
/
∂
l
).
The sensitivity index of the model parameters is given by(47)ξαhR0=μh2αh+μh,ξαmR0=pu3+μm+au12αm+pu3+au1+μm,ξμmR0=-μm2αm+3μm+au1+pu32αm+μm+au1+pu3μm+au1+pu3,ξδhR0=-δh2b+τu2+δh+μh,ξbR0=-b2b+τu2+δh+μh,ξμhR0=-μh2+αhδh+αhb+αhτu22μh+δh+b+τu2μh+αh.Sensitivity indices for the control parameters are given by(48)ξu1R0=-μ11-μ1,ξu2R0=-τμ22μh+ϵ+b+τμ2,ξu3R0=-p3pu3+3μm+au1+2αmu32pu3+au1+μmpu3+μm+αm,ξu4R0=-u41-u4.Through sensitivity analysis, it is observed that the most sensitive parameters to R
0 across all the settings were the mosquito’s natural death rate, μ
m, 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
R
0. Sensitivity indices are calculated using parameters in Table 5.Table 5
Parameter values for the full malaria model.
Parameter
Estimated value
Source
Endemic
Epidemic
Seasonal
Low risk
μ
h
0.00005447
0.00004644
0.00004281
0.00004566
KNBS (2009 Census estimates) [20]
μ
m
0.04
Estimated
α
h
0.07143
Estimated
α
m
0.0909
Chitnis et al. [21]
λ
0.42
Estimated
β
0.0655
Estimated
ϵ
0.2
Blayneh et al. [22]
ψ
0.01095
Estimated
Λ
h
0.00000575
0.00000575
0.00000548
0.00000438
KNBS (estimates based on 2009 Census) [20]
Λ
m
0.071
Niger and Gumel [23]
b
0.005
Chiyaka et al. [24]
τ
0.5
Assumed
δ
h
0.05
KNBS and ICF Macro [25]
p
0.25
Assumed
a
0.25
Assumed
ϕ
0.502
Blayneh et al. [22]
λ
h
0.00000149
0.00000123
0.00000445
0.00000226
Estimated
λ
h
w
0.00000247
0.00000203
0.00000693
0.00000328
Estimated
λ
m
0.00000048
0.00000394
0.00000143
0.00000073
Estimated
N
h
440169
535093
147832
290703
KNBS (2009 Census estimate) [20]
N
h
w
266343
324658
94857
200216
KNBS (2009 Census estimate) [20]
N
m
4401690
5350930
1478320
2907030
Estimated


## 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 t
≥
0.The feasible solutions set for model (1) given by (3)
D
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
S
h
,
S
m
>
0
,
E
h
,
I
h
,
R
h
,
E
m
,
I
m
≥
0
;
N
h
≤
Λ
h
μ
h
;
N
m
≤
Λ
m
μ
m
+
a
u
1
+
p
u
3is positively invariant and hence model (1) is biologically, epidemiologically meaningful and mathematically well posed in the domain D.System (1) has always a disease-free equilibrium given by(4)E0=Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗=Λhμh,0,0,0,Λmμm+au1+pu3,0,0.


## 2.1.2. Basic Reproduction Number

The matricesF and V for the new infection terms and the remaining transfer terms at disease-free equilibrium [26], respectively, are given by(5)F=0001-u1βϵϕ+1-u4βϵϕ000001-u1λϵϕμhΛmαhμm+au1+pu3000000,V=μh+αh000-αhδh+μh+b+τu20000αm+μm+au1+pu3000-αmμm+au1+pu3.It follows that the basic reproduction of model (1), denoted by R
0, is given by(6)R0=ρFV-1,R0=αh1-u1λϵϕΛmμh1-u1βϵϕαm+αh1-u1λϵϕΛmμh1-u4βϵϕαmμm+au1+pu3μh+αhδh+μh+b+τu2Λhαm+μm+au1+pu3μm+au1+pu3.


## 2.1.3. Stability Analysis of Disease-Free Equilibrium Point

(
1) Local Stability of Disease-Free Equilibrium PointTheorem 1.
The disease-free equilibrium point for system (1) is locally asymptotically stable if R
0
<
1.Proof.
The Jacobian matrix(
J
) of the malaria model (1) at the disease-free equilibrium point is given by(7)-αh+μh0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu20000b+τu2-ψ+μh0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu30000αm-μm+au1+pu3.The eigenvalues of the Jacobian matrix are the solutions of the characteristic equation(8)J-λI=0,-αh+μh+λ0001-u1βϵϕ+1-u4βϵϕαh-δh+μh+b+τu2+λ0000b+τu2-ψ+μh+λ0001-u1λϵϕΛmμhμm+au1+pu3Λh0-αm+μm+au1+pu3+λ0000αm-μm+au1+pu3+λ=0.Expanding the determinant into a characteristic equation we have(9)-δh+μhαh+μh+λδh+μh+b+τu2+λαm+μm+au1+pu3+λμm+au1+pu3+λ-1-u12λϵ2ϕ2Λmμhαmβαh+1-u41-u1λϵ2ϕ2Λmμhαmβαhμm+au1+pu3Λm=0.Hence we have(10)
A
1
=
μ
m
+
a
u
1
+
p
u
3
,
A
2
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
,
A
3
=
δ
h
+
μ
h
+
b
+
τ
u
2
,
A
4
=
α
h
+
μ
h
,
Q
=
1
-
u
1
2
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
+
1
-
u
4
1
-
u
1
λ
ϵ
2
ϕ
2
Λ
m
μ
h
α
m
β
α
h
·
μ
m
+
a
u
1
+
p
u
3
Λ
m
-
1
,
λ
+
A
1
λ
+
A
2
λ
+
A
3
λ
+
A
4
-
Q
=
0
,
λ
4
+
B
1
λ
3
+
B
2
λ
2
+
B
3
λ
+
B
4
=
0
,where(11)B1=A4+A3+A2+A1,B2=A4A3+A2+A1+A3A2+A1+A1A2,B3=A4A3A2+A4A3A1+A4A2A1+A3A2A1,B4=A4A3A2A1-Q.Thus, 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 R
0
<
1.Next, we study the global behavior of the disease-free equilibrium for system (1).(
2) Global Stability of Disease-Free Equilibrium PointTheorem 2.
The DFE,E
0, of system of (1) is globally asymptotically stable if R
0
<
1.Proof.
We consider the following Lyapunov function:(12)L=c1Eh+c2Ih+c3Em+c4Im,where(13)c1=αhμm+au1+pu3αh+μhb+τu2+δh+μh,c2=1μm+au1+pu3b+τu2+δh+μh,c3=11-u1ϵϕλΛm,c4=αm+μm+au1+pu31-u1ϵϕλΛmαm.Computing the derivative of L along the solution of the system of differential equation (1),(14)
L
˙
=
α
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
E
h
+
1
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
I
h
+
1
1
-
u
1
ϵ
ϕ
λ
Λ
m
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
I
m
,
L
˙
=
α
h
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
b
+
τ
u
2
+
δ
h
+
μ
h
E
h
+
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
μ
m
+
a
u
1
+
p
u
3
b
+
τ
u
2
+
δ
h
+
μ
h
I
h
+
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
E
m
+
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
I
m
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
1
β
ϵ
ϕ
α
m
+
α
h
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
1
-
u
4
β
ϵ
ϕ
α
m
μ
m
+
a
u
1
+
p
u
3
μ
h
+
α
h
δ
h
+
μ
h
+
b
+
τ
u
2
Λ
h
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
-
1
I
h
,
L
˙
=
α
m
+
μ
m
+
a
u
1
+
p
u
3
μ
m
+
a
u
1
+
p
u
3
+
δ
h
1
-
u
1
ϵ
ϕ
λ
Λ
m
α
m
R
0
-
1
I
h
≤
0
iff
R
0
≤
1
.Thus we have established that L
˙
≤
0 if R
0
<
1 and the equality L
˙
=
0 holds if and only if R
0
=
1 and E
h
=
I
h
=
E
m
=
I
m
=
0. If R
0
>
1 then L
˙
>
0 when S
h
(
t
) and S
m
(
t
) are sufficiently close to Λ
h
/
μ
h and Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
), respectively, except when E
h
=
I
h
=
E
m
=
I
m
=
0
.
On the boundary whenE
h
=
I
h
=
E
m
=
I
m
=
0, N
h
˙
t
=
Λ
h
-
μ
h
N
h and N
m
˙
t
=
Λ
m
-
(
μ
m
+
a
u
1
+
p
u
3
) and N
h
(
t
)
→
Λ
h
/
μ
h, N
m
(
t
)
→
Λ
m
/
(
μ
m
+
a
u
1
+
p
u
3
) as t
→
∞.
Therefore the largest compact invariantD
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m
∈
R
+
7
:
L
˙
=
0 when R
0
<
1 is the singleton E
0 in D. By LaSalle’s invariant principle [28], E
0 is globally asymptotically stable for R
0
<
1 in D.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 (E
1) of the model is given by(15)E1=Sh∗∗,Eh∗∗,Ih∗∗,Rh∗∗,Sm∗∗,Em∗∗,Im∗∗.To derive E
1, we solve model (1) by equating it to zero.Substituting and solving forI
h
∗
∗ (as an expression of parameters only) through some algebraic manipulation give(16)Λh+ψb+τu2Ih∗∗ψ+μh+1-u1βϵϕNh1-u1μm+au1+pu3R0mIh∗∗1-u1λϵϕIh∗∗+Nhμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh+1-u4βϵϕNhw1-u4μm+au1+pu3R0mIh∗∗1-u4λϵϕIh∗∗+Nhwμm+au1+pu3·Nhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh-μhNhwμh1-u1λϵϕIh∗∗+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh=0or(17)AIh∗∗2+BIh∗∗+C=0,where(18)
A
=
ψ
+
μ
h
Λ
h
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
+
ψ
b
+
τ
u
2
I
h
∗
∗
μ
h
N
h
w
R
0
2
1
-
u
1
2
·
μ
m
+
a
u
1
+
p
u
3
α
h
+
μ
h
δ
h
+
μ
h
+
b
+
τ
u
2
+
α
h
1
-
u
4
1
-
u
1
μ
m
+
a
u
1
+
p
u
3
β
ϵ
ϕ
R
0
m
μ
h
·
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
·
N
h
w
N
h
,
B
=
ψ
+
μ
h
1
-
u
1
β
ϵ
ϕ
1
-
u
1
·
μ
m
+
a
u
1
+
p
u
3
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
·
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
w
+
ψ
+
μ
h
1
-
u
4
β
ϵ
ϕ
1
-
u
4
μ
m
+
a
u
1
+
p
u
3
·
R
0
m
I
h
∗
∗
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
N
h
,
C
=
ψ
+
μ
h
μ
h
N
h
w
μ
h
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
Λ
h
N
h
w
μ
m
+
a
u
1
+
p
u
3
1
-
u
1
λ
ϵ
ϕ
I
h
∗
∗
+
N
h
μ
m
+
a
u
1
+
p
u
3
N
h
w
N
h
.We use the quadratic formula to find the roots of (17); that is,(19)Ih∗=-B±B2-4AC2A=-B+B2-4AC2Aor-B-B2-4AC2AIh∗=-B+B2-4AC2A=-B-B2-4AC2A=Φ.Hence(20)Sh∗∗=Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Eh∗∗=1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦNhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Ih∗∗=Φ=αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh≥0,Rh∗∗=b+τu2αh1-u1βϵϕ+1-u4βϵϕ1-u1μm+au1+pu3R0mΦμh+ψδh+μh+b+τu2NhwNhμh+αh1-u1λϵϕΦ+Nhμm+au1+pu3·Nhwμh1-u1λϵϕΦ+ΛhNhwμm+au1+pu3μhNhwR021-u12μm+au1+pu3αh+μhδh+μh+b+τu2+αh1-u41-u1μm+au1+pu3βϵϕR0mμh,Sm∗∗=ΛmNh1-u1λϵϕΦ+μm+au1+pu3Nh,Em∗∗=1-u1αhλϵϕ1-u1μm+au1+pu3R0mΦΛmNhNhαm+μm+au1+pu31-u1λϵϕΦ+Nhμm+au1+pu31-u1λϵϕΦ+μm+au1+pu3Nh,Im∗∗=1-u1μm+au1+pu3R0mΦ1-u1λϵϕΦ+Nhμm+au1+pu3.From the quadratic equation (17) we analyze the possibility of multiple endemic equilibria. It is important to note that the coefficient A is always positive with B and C having different signs. All the negative terms in C are in the form of (
1
-
u
) where u is the control that is bounded by 1.It follows that(i)
there is a unique endemic equilibrium ifB
<
0 and C
=
0 or B
2
-
4
A
C
=
0,
(ii)
there is a unique endemic equilibrium ifC
<
0 (i.e., if R
0
>
1),
(iii)
there are two endemic equilibria ifC
>
0, B
<
0, and B
2
-
4
A
C
>
0,
(iv)
there are no endemic equilibria otherwise.Note that the hypothesis C
>
0 is equivalent to R
0
<
1.Thus the results of this section can be summarized in the following theorem.Theorem 3.
IfR
0
<
1, E
0 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 R
0
>
1, then E
0 is unstable and there exists a unique endemic equilibrium.Item (iii) indicates the possibility of backward bifurcation in model (1) when R
0
<
1. In the next section we will prove the occurrence of multiple equilibria for R
0
<
1 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 outbifurcation analysis of system (1) at R
0
=
1. We consider a transmission rate β as bifurcation parameter Ψ so that R
0
=
1.To apply the theory, we introduce dimensionless state variables into system (1).The system of (1) can be written as (21)dShdt=Λh+ψRh-μhSh-1-u1βϵϕShImNh-1-u4βϵϕShImNhw,dEhdt=1-u1βϵϕShImNh+1-u4βϵϕShImNhw-μhEh-αhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-μhRh-ψRh,dSmdt=Λm-1-u1λϵϕIhSmNh-μm+au1+pu3Sm,dEmdt=1-u1λϵϕIhSmNh-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im.We let x
1
=
S
h, x
2
=
E
h, x
3
=
I
h, x
4
=
R
h, x
5
=
S
m, x
6
=
E
m, and x
7
=
I
m
.Therefore system (1) in vector form can be written as(22)dXidt=HXi,where X
i
=
x
1
,
x
2
,
…
,
x
7
T and H
=
h
1
,
h
2
,
…
,
h
7
T are transposed matrices and N
h
=
Λ
h
/
μ
h with Ψ
∗
=
β, (23)dx1dt=Λh+ψx4-μhSh-1-u1Ψ∗ϵϕx7x1μhΛh-1-u4Ψ∗ϵϕx7x1μhΛh=h1,dx2dt=1-u1Ψ∗ϵϕx7x1μhΛh+1-u4Ψ∗ϵϕx7x1μhΛh-μh+αhx2=h2,dx3dt=αhx2-δh+μh+b+τu2x3=h3,dx4dt=b+τu2x3-μh+ψx4=h4,dx5dt=Λm-1-u1λϵϕx3x5μhΛh-μm+au1+pu3x5=h5,dx6dt=1-u1λϵϕx3x5μhΛh-αm+μm+au1+pu3x6=h6,dx7dt=αmx6-μm+au1+pu3x7=h7.Let Ψ
∗ be the bifurcation parameter; the system is linearized at disease-free equilibrium point when β
=
Ψ
∗ with R
0
=
1. That is, (24)R0=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,12=αhαmΛmμh1-u12ϕ2ϵβλ+αhαmΛmμh1-u11-u4ϕ2ϵ2βΛhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2,Ψ∗=Λhpu3+au1+μm2μh+α1pu3+μm+au1+αmμh+δh+b+τu2αhαmΛmμh1-u12ϕ2ϵλ+αhαmΛmμh1-u11-u4ϕ2ϵ2.The Jacobian matrix of (1) calculated at Ψ
∗ is given by (25)
-
μ
h
0
0
ψ
0
0
-
Ψ
∗
ϵ
ϕ
0
-
α
h
-
μ
h
0
0
0
0
Ψ
∗
ϵ
ϕ
0
α
h
-
δ
h
-
μ
h
-
b
-
τ
u
2
0
0
0
0
0
0
b
+
τ
u
2
-
μ
h
-
ψ
0
0
0
0
0
-
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
-
μ
m
+
a
u
1
+
p
u
3
0
0
0
0
1
-
u
1
λ
ϵ
ϕ
Λ
m
μ
h
Λ
h
μ
m
+
a
u
1
+
p
u
3
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
0
0
0
0
0
0
-
α
m
-
μ
m
-
a
u
1
-
p
u
3
.Centre manifold approach [29] is then applied.A right eigenvector associated with the eigenvalue zero isw
=
w
1
,
w
2
,
…
,
w
7. Solving the system we have the following right eigenvector:(26)w1=ψw4-Ψ∗ϵϕw7μh,w2=Ψ∗ϵϕw7αh+μh,w3=αhw2b+τu2+μh+δh,w4=b+τu2w3μh+ψ,w5=-1-u1λϵϕΛmμhw3Λhμm+au1+pu32,w6=1-u1λϵϕΛmμhw3Λhμm+au1+pu3αm+μm+au1+pu3,w7=αmw6μm+au1+pu3>0.The left eigenvectors satisfying v
·
w
=
1 are v
=
v
1
,
v
2
,
…
,
v
7. Solving the system, the left eigenvector will be given by(27)v1=0,v2=αhv3αh+μh,v3=v61-u1λϵϕΛmμhΛhμm+au1+pu3-δh-uh-b-τu2,v4=0,v5=0,v6=αmv7-μm-αm-au1-pu3,v7=Ψ∗ϵϕv2μm+au1+pu3.Computing for the sign of a and b as indicated in the theorem gives(28)
a
=
v
2
w
7
-
1
-
u
1
Ψ
∗
ϵ
ϕ
μ
h
-
1
-
u
4
Ψ
∗
ϵ
ϕ
μ
h
Λ
h
·
Ψ
∗
ϵ
w
7
α
h
+
μ
h
+
α
h
w
2
b
+
τ
u
2
+
μ
h
+
δ
h
+
v
6
w
3
-
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
h
·
1
-
u
1
λ
ϵ
ϕ
μ
h
Λ
m
w
3
Λ
h
μ
m
+
a
u
1
+
p
u
3
α
m
+
μ
m
+
a
u
1
+
p
u
3
+
α
m
w
3
μ
m
+
a
u
1
+
p
u
3
b
=
v
2
w
7
ϵ
ϕ
>
0so that b is always positive.Therefore the following result is established.Theorem 4.
Model (1) exhibits backward bifurcation at R
0
=
1 whenever a
>
0 and b
>
0 and R
0
<
1. Whenever a
<
0 and b
>
0, then model (1) exhibits a forward bifurcation at R
0
=
1.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 R
0
>
1 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 D. 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, E
1: (i) the uniqueness of E
1 in the interior of D (condition H
1); (ii) the existence of an absorbing compact set in the interior of D (condition H
2); and (iii) the fulfillment of a Bendixson criterion (i.e., inequality (A.6)).Theorem 5.
Ifa
<
0 and b
>
0 and R
0
>
1, then the endemic equilibrium of the malaria model (1) is globally asymptotically stable in the interior of D.Proof.
Following Li and Muldowney [30], for system (1), under the assumption of R
0
>
1, satisfies conditions H
1
-
(
H
2
), 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 c
>
0 such that any solutions (
S
h
t
,
I
h
t
,
I
m
t
) with (
S
h
0
,
I
h
0
,
I
m
0
) in the interior of D satisfy(29)min⁡limt→∞⁡inf⁡Sht,limt→∞⁡inf⁡Iht,limt→∞⁡inf⁡Imt.The uniform persistence together with boundedness of D is equivalent to the existence of a compact set in the interior of D which is absorbing for (4) as described by Hutson and Schmitt [32]. Thus, (
H
1
) is verified. Moreover, E
1 is the only equilibrium in the interior of D, so that (
H
2
) 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 S
h
,
I
h
,
I
m of reduced system (1), we get the second additive compound matrix J
2:(30)J2Sh,Ih,Im=-a11a121-u1ϵϕβSh+1-u4ϵϕβSh1-u1ϵϕλNm-Em-Im-a22-δh01-u1ϵϕβIm+1-u4ϵϕβIm-a33,where (31)a11=μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+δh+b+τu2,a12=1-u1ϵϕβSh+1-u4ϵϕβSh,a22=μh+δh+1-u1ϵϕβIm+1-u4ϵϕβIm+1-u1ϵϕλIh+μm+au1+pu3,a33=m+μm+au1+pu3+1-u1ϵϕλIh+1-u1ϵϕβIm+1-u4ϵϕβIm,where m
=
Λ
m
/
μ
m
+
a
u
1
+
p
u
3
/
Λ
h
/
μ
h
.
Choose now matrixP
(
x
)
=
P
S
h
,
I
h
,
I
m
=
diag
⁡
(
1
,
I
h
/
I
m
,
I
h
/
I
m
). Then P
f
P
-
1
=
diag
⁡
(
0
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
,
I
h
˙
/
I
h
-
I
m
˙
/
I
m
), and the matrix B
=
P
f
P
-
1
+
P
J
2
P
-
1 can be written in block form as(32)B=B11B12B21B22,where(33)B11=-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm+μh+b+τu2+δh,B12=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh,1-u1ϵϕβSh·ImIh+1-u4ϵϕβSh·ImIh,B21=1-u1ϵλϕImNm-Em+ImIhImT,B22=Ih˙Ih-Im˙Im-a2201-u1ϵϕβIm-1-u4ϵϕβIm-1-u1ϵϕλIhIh˙Ih-Im˙Im-a33.The vector norm · in R
+
3 is here chosen to be (34)x,y,z=max⁡x,y,z.Let σ
(
·
) denote the Lozinskii measure with respect to this norm. Using the method of estimating σ
(
·
) in [29], we have(35)σB≤sup⁡g1,g2=sup⁡σ1B11+B12,σ1B22+B21,where B
12 and B
21 are matrix norms with respect to the L vector norm and σ
1 denotes the Lozinskii measure with respect to L norm. Since B
11 is a scalar, its Lozinskii measure with respect to any norm in R
1 is equal to B
11.
Therefore(36)
σ
1
B
11
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
·
ϵ
ϕ
β
I
m
+
μ
h
+
b
+
τ
u
2
+
δ
h
,
σ
1
B
22
=
max
⁡
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
I
h
˙
I
h
-
I
m
˙
I
m
-
1
-
u
1
ϵ
ϕ
β
I
m
-
1
-
u
4
ϵ
ϕ
β
I
m
-
1
-
u
1
ϵ
ϕ
λ
I
h
-
μ
h
-
δ
h
-
b
-
τ
u
2
-
μ
m
-
a
u
1
-
p
u
3
,
σ
1
B
22
=
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
1
-
u
1
ϵ
ϕ
λ
I
h
+
μ
m
+
a
u
1
+
p
u
3
,
B
12
=
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
B
21
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
.Therefore(37)
g
1
=
-
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
I
m
+
1
-
u
4
ϵ
ϕ
β
I
m
+
μ
h
+
δ
h
+
b
+
τ
u
2
+
1
-
u
1
ϵ
ϕ
β
S
h
+
1
-
u
4
ϵ
ϕ
β
S
h
·
I
m
I
h
,
(38)
g
2
=
1
-
u
1
ϵ
ϕ
λ
N
m
-
E
m
-
I
m
·
I
h
I
m
+
I
h
˙
I
h
-
I
m
˙
I
m
-
μ
h
+
α
h
+
μ
m
+
a
u
1
+
p
u
3
+
1
-
u
1
ϵ
ϕ
λ
I
h
.We rewrite the last two equations of system (1) for I
h
˙ and I
m
˙ as (39)Ih˙Ih=1-u1ϵϕβSh+1-u4ϵϕβSh·ImIh-μh+δh+b+τu2,(40)Im˙Im=1-u1ϵϕλNm-Em-Im·IhIm-μm-au1-pu3.Substituting (39) into (37) and (40) into (38) we have(41)g1t=Ih˙Ih-μh+αh+1-u1ϵϕβIm+1-u4ϵϕβIm.g2t=Ih˙Ih-μh+αh+1-u1ϵϕλIh.For the uniform persistence constant ε
>
0, there exists a time T
0
>
0 independent of x
(
0
)
∈
K, the compact absorbing set, such that(42)Iht>ε,Imt>ε;for t
>
T
0 we have(43)g1t≤Ih˙Ih-μh+αh+1-u1ϵϕβε+1-u4ϵϕβε,g2t≤Ih˙Ih-μh+αh+1-u1ϵϕλε.Relations (41) imply (44)σB≤Ih˙Ih-μfort>T0,where(45)
μ
=
min
⁡
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
β
ε
+
1
-
u
4
ϵ
ϕ
β
ε
,
μ
h
+
α
h
+
1
-
u
1
ϵ
ϕ
λ
ε
.Along each solution S
h
t
,
I
h
t
,
I
m
(
t
) to (1) with S
h
0
,
I
h
0
,
I
m
(
0
)
∈
K where K is the compact absorbing set, we have, for t
>
T
0,(46)1t∫0tσBds≤1t∫0T0σBds+1tln⁡IhtIhT0-μt-T0t.Which implies q
-
2
<
μ
/
2
<
0. 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 Table4 is computed. The index measures the relative change in a variable with respect to relative changes in parameters.Table 4
Sensitivity indices (SI) ofR
0 to parameters for the malaria model.
Parameter
Sensitivity indices
Endemic
Seasonal
Epidemic
Low risk
μ
_
h
−0.0402531
−0.0402531
−0.0402531
−0.0402531
μ
_
m
−1.07211
−1.07211
−1.07211
−1.07211
α
_
h
0.00038817
0.00038817
0.00038817
0.00038817
α
_
m
0.22445
0.22445
0.22445
0.22445
λ
0.5
0.5
0.5
0.5
β
0.5
0.5
0.5
0.5
ϵ
1
1
1
1
Λ
_
h
−0.4987
−0.4987
−0.4987
−0.4987
Λ
_
m
0.5
0.5
0.5
0.5
b
−0.01818
−0.02048
−0.01639
−0.02563
τ
−0.322497
−0.322497
−0.322497
−0.322497
δ
_
h
−0.13695
−0.10336
−0.11321
−0.03508
ϕ
1
1
1
1Definition 6.
Following Chitnis et al. [33], the normalized forward sensitivity index of a variable, h, that depends on a parameter, l, is defined as ξ
l
=
(
l
/
h
)
×
(
∂
h
/
∂
l
).
The sensitivity index of the model parameters is given by(47)ξαhR0=μh2αh+μh,ξαmR0=pu3+μm+au12αm+pu3+au1+μm,ξμmR0=-μm2αm+3μm+au1+pu32αm+μm+au1+pu3μm+au1+pu3,ξδhR0=-δh2b+τu2+δh+μh,ξbR0=-b2b+τu2+δh+μh,ξμhR0=-μh2+αhδh+αhb+αhτu22μh+δh+b+τu2μh+αh.Sensitivity indices for the control parameters are given by(48)ξu1R0=-μ11-μ1,ξu2R0=-τμ22μh+ϵ+b+τμ2,ξu3R0=-p3pu3+3μm+au1+2αmu32pu3+au1+μmpu3+μm+αm,ξu4R0=-u41-u4.Through sensitivity analysis, it is observed that the most sensitive parameters to R
0 across all the settings were the mosquito’s natural death rate, μ
m, 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
R
0. Sensitivity indices are calculated using parameters in Table 5.Table 5
Parameter values for the full malaria model.
Parameter
Estimated value
Source
Endemic
Epidemic
Seasonal
Low risk
μ
h
0.00005447
0.00004644
0.00004281
0.00004566
KNBS (2009 Census estimates) [20]
μ
m
0.04
Estimated
α
h
0.07143
Estimated
α
m
0.0909
Chitnis et al. [21]
λ
0.42
Estimated
β
0.0655
Estimated
ϵ
0.2
Blayneh et al. [22]
ψ
0.01095
Estimated
Λ
h
0.00000575
0.00000575
0.00000548
0.00000438
KNBS (estimates based on 2009 Census) [20]
Λ
m
0.071
Niger and Gumel [23]
b
0.005
Chiyaka et al. [24]
τ
0.5
Assumed
δ
h
0.05
KNBS and ICF Macro [25]
p
0.25
Assumed
a
0.25
Assumed
ϕ
0.502
Blayneh et al. [22]
λ
h
0.00000149
0.00000123
0.00000445
0.00000226
Estimated
λ
h
w
0.00000247
0.00000203
0.00000693
0.00000328
Estimated
λ
m
0.00000048
0.00000394
0.00000143
0.00000073
Estimated
N
h
440169
535093
147832
290703
KNBS (2009 Census estimate) [20]
N
h
w
266343
324658
94857
200216
KNBS (2009 Census estimate) [20]
N
m
4401690
5350930
1478320
2907030
Estimated


## 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 variablesu
t
=
u
1
,
u
2
,
u
3
,
u
4
∈
U relative to the state variables S
h, E
h, I
h, R
h, S
m, E
m, I
m where control variables are bounded and measured with(49)
U
=
u
1
,
u
2
,
u
3
,
u
4
∈
U
is  Lebesgue  measurable  on
0,1
,
0
≤
u
i
t
≤
1
,
t
∈
0
,
T
,
i
=
1,2
,
3,4
.We define the objective function as(50)
J
u
1
,
u
2
,
u
3
,
u
4
=
∫
0
T
A
1
N
m
+
A
2
I
h
+
A
3
E
h
+
1
2
B
1
u
1
2
+
B
2
u
2
2
+
B
3
u
3
2
+
B
4
u
4
2
d
tsubject to(51)dShdt=Λh+ψRh-1-u1λhSh-1-u4λhwSh-μhSh,dEhdt=1-u1λhSh+1-u4λhwSh-αh+μhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-ψ+μhRh,dSmdt=Λm-1-u1λmSm-μm+au1+pu3Sm,dEmdt=1-u1λmSm-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im,Sh0≥0,Eh0≥0,Ih0≥0,Rh0≥0,Sm0≥0,Em0≥0,Im0≥0.In the objective function T is the final time and quantities A
1, A
2, and A
3 are weights constants of the total mosquito population, infected individuals, and exposed individuals, respectively, while B
1, B
2, B
3, and B
4 are weight constants for use with ITNs, treatment effort, IRS, and IPTp efforts, respectively. The total mosquito population N
m
=
S
m
+
E
m
+
I
m is part of the objective function because it is affected by the use of IRS and ITNs. In addition, E
h and I
h 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 B
1
u
1
2, treatment of infected individuals is B
2
u
2
2, spraying of houses with IRS is B
3
u
3
2, and preventive method of IPTp is B
4
u
4
2. A linear function has been chosen for the cost incurred by exposed individuals A
3
E
h, infected individuals, A
2
I
h, and the mosquito population, A
1
N
m. A quadratic form is used for the cost on the controls B
1
u
1
2, B
2
u
2
2, B
3
u
3
2, and B
4
u
4
2, such that the terms (
1
/
2
)
B
1
u
1
2, (
1
/
2
)
B
2
u
2
2, (
1
/
2
)
B
3
u
3
2, and (
1
/
2
)
B
4
u
4
2 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 populationN
m, the number of exposed humans E
h
(
t
), and infected humans I
h
(
t
) while minimizing the cost of control u
1
t, u
2
t, u
3
t, and u
4
(
t
). 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 [
0
,
1
].We seek an optimal controlu
1
∗, u
2
∗, u
3
∗, and u
4
∗ such that(52)Ju1∗,u2∗,u3∗,u4∗=minu1,u2,u3,u4∈U⁡Ju1,u2,u3,u4.Pontryagin’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, L, and Hamiltonian, H, with respect to u
1, u
2, u
3, and u
4.


### 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:(H
1):
the set of controls and corresponding state variables is nonempty.
(H
2):
the control setU is convex and closed.
(H
3):
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 u with coefficients depending on time and state.
(H
4):
there exist constantsc
1
,
c
2
>
0 and β
>
1 such that the integrand L
(
y
,
u
,
t
) of the objective functional J is convex and satisfies(53)Ly,u,t≥c1u12+u22+u32+u42β/2-c2.In 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, (
H
1
) is satisfied. The solutions are bounded and hence the control set satisfies (
H
2
) as well. The system is bilinear in u
1
,
u
2
,
u
3
,
u
4 and hence the right hand side of state system (1) satisfies condition (
H
3
) (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 U is closed and convex. The integrand of the objective cost function J expressed by (51) is a convex function of u
1
,
u
2
,
u
3
,
u
4 on the control set U. 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 ξ
1, ξ
2 and a constant ε
>
1 such that (54)A1Nm+A2Ih+A3Eh+12B1u12+B2u22+B3u32+B4u42≥ξ1u12+u22+u32+u42ε/2-ξ2,where ξ
1, ξ
2
>
0, A
1, A
2, B
1, B
2, B
3, B
4
>
0, and ε
>
1
.This concludes existence of an optimal control since the state variables are bounded. Hence we have the following theorem.Theorem 7.
Given the objective functionalJ
u
1
,
u
2
,
u
3
,
u
4
=
∫
0
T
A
1
N
m
+
A
2
I
h
+
A
3
E
h
+
(
1
/
2
)
(
B
1
u
1
2
+
B
2
u
2
2
+
B
3
u
3
2
+
B
4
u
4
2
)
d
t, where U
=
u
1
,
u
2
,
u
3
,
u
4
∣
0
≤
u
1
,
u
2
,
u
3
,
u
4
≤
1
,
0
≤
u
i
t
≤
1
,
t
∈
0
,
T
,
i
=
1,2
,
3,4 subject to (1) with initial conditions, then there exists an optimal control u
∗
=
u
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗ such that J
u
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗
=
min
U
⁡
J
u
1
,
u
2
,
u
3
,
u
4.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 by(55)LIh,Eh,Nm,u1,u2,u3,u4=A1Nm+A2Ih+A3Eh+12B1u12+B2u22+B3u32+B4u42.We seek to find the minimal value of the Lagrangian. To do this, we define the Hamiltonian H for the control problem which consists of the integrand of the objective functional (Lagrangian, L) and the inner product of the right hand sides of the state equations and the costate variables or adjoint variables λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7 as(56)H=L+λ1dShdt+λ2dEhdt+λ3dIhdt+λ4dRhdt+λ5dSmdt+λ6dEmdt+λ7dImdt.Taking X
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m, U
=
u
1
,
u
2
,
u
3
,
u
4, and λ
=
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7 we obtain the Hamiltonian given by(57)
H
X
,
U
,
λ
=
L
I
h
,
E
h
,
N
m
,
u
1
,
u
2
,
u
3
,
u
4
+
λ
1
Λ
h
+
ψ
R
h
-
1
-
u
1
λ
h
S
h
-
1
-
u
4
λ
h
w
S
h
-
μ
h
S
h
+
λ
2
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
E
h
+
λ
3
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
I
h
+
λ
4
b
+
τ
u
2
I
h
-
ψ
+
μ
h
R
h
+
λ
5
Λ
h
-
1
-
u
1
λ
m
S
m
-
μ
m
+
a
u
1
+
p
u
3
S
m
+
λ
6
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
E
m
+
λ
7
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
I
m
.


### 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.Ifu
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗ are optimal solution of an optimal control problem, then there exists a nontrivial vector function λ
t
=
λ
1
t
,
λ
2
t
,
…
,
λ
n
(
t
) satisfying the following conditions.The state equation is(58)dxdt=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂λ.The optimality condition is(59)0=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂uand the adjoint equation is(60)dλdt=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂x.Now, we apply the necessary conditions to the Hamiltonian.Theorem 8.
Given the optimal controlsu
1
∗, u
2
∗, u
3
∗, u
4
∗ and solutions S
h
∗, E
h
∗, I
h
∗, R
h
∗, S
m
∗, E
m
∗, I
m
∗ of the corresponding state system (1), there exist adjoint variables λ
1, λ
2, λ
3, λ
4, λ
5, λ
6, λ
7 satisfying(61)-dλ1dt=∂H∂Sh=1-u1λhλ2-λ1+1-u4λhwλ2-λ1-μhλ1,-dλ2dt=∂H∂Eh=αhλ3-λ2-μhλ2+A3,-dλ3dt=∂H∂Ih=b+τu2λ4-b+τu2+μh+δhλ3+A2-1-u1λϵϕSmNhλ5+1-u1λϵϕSmNhλ6,-dλ4dt=∂H∂Rh=ψλ1-μh+ψλ4,-dλ5dt=∂H∂Sm=1-u1λ6-λ5λm-μm+au1+pu3λ5+A1,-dλ6dt=∂H∂Em=αmλ7-λ6-μm+au1+pu3λ6+A1,-dλ7dt=∂H∂Im=-μm+au1+pu3λ7+A1+-1-u1βϵϕShNh-1-u4βϵϕShNhwλ1+1-u1βϵϕShNh+1-u4βϵϕShNhwλ2with transversality conditions(62)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=0.Furthermore u
1
∗, u
2
∗, u
3
∗, u
4
∗ are represented by(63)u1∗=max⁡0,min⁡1,λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=max⁡0,min⁡1,τλ3-λ4Ih∗B2,u3∗=max⁡0,min⁡1,pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=max⁡0,min⁡1,λ2-λ1λhwSh∗B4.Proof.
To determine the adjoint equations and the transversality conditions we use the HamiltonianH. The Hamiltonian function, H, is differentiated with respect to S
h, E
h, I
h, R
h, S
m, E
m, and I
m. The adjoint/costate equation is given by(64)-dλ1dt=∂H∂Sh=1-u1λhλ2-λ1+1-u4λhwλ2-λ1-μhλ1,-dλ2dt=∂H∂Eh=αhλ3-λ2-μhλ2+A3,-dλ3dt=∂H∂Ih=b+τu2λ4-b+τu2+μh+δhλ3+A2-1-u1λϵϕSmNhλ5+1-u1λϵϕSmNhλ6,-dλ4dt=∂H∂Rh=ψλ1-μh+ψλ4,-dλ5dt=∂H∂Sm=1-u1λ6-λ5λm-μm+au1+pu3λ5+A1,-dλ6dt=∂H∂Em=αmλ7-λ6-μm+au1+pu3λ6+A1,-dλ7dt=∂H∂Im=-μm+au1+pu3λ7+A1+-1-u1βϵϕShNh-1-u4βϵϕShNhwλ1+1-u1βϵϕShNh+1-u4βϵϕShNhwλ2with transversality conditions(65)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=0.In order to minimize Hamiltonian,H, with respect to the controls at the optimal controls, H is differentiated with respect to u
1, u
2, u
3, and u
4 on the set U, and the solution for the optimal control point is obtained after equating to zero. This is the optimality condition.Solving∂
H
/
∂
u
1
=
0, ∂
H
/
∂
u
2
=
0, ∂
H
/
∂
u
3
=
0, and ∂
H
/
∂
u
4
=
0, evaluating at the optimal control on the interior of the control set, where 0
<
u
i
<
1, for i
=
1,2
,
3,4, and letting S
h
=
S
h
∗, E
h
=
E
h
∗, I
h
=
I
h
∗, R
h
=
R
h
∗, S
m
=
S
m
∗, E
m
=
E
m
∗, and I
m
=
I
m
∗ yield (66)∂H∂u1=B1u1∗+λ1λhSh∗-λ2λhSh∗+λ5λmSm∗-aSm∗λ5-aEm∗λ6-aIm∗λ7=0,∂H∂u2=B2u2∗-τλ3Ih∗+τλ4Ih∗=0,∂H∂u3=B3u3∗-pλ5Sm∗-pλ6Em∗-pλ7Im∗=0,∂H∂u4=B4u4∗+λhwλ2Sh∗-λ1λhwSh∗=0for which(67)u1∗=λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=τλ3-λ4Ih∗B2,u3∗=pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=λ2-λ1λhwSh∗B4.By applying the boundary condition of each control, the solution of (67) becomes (68)u1∗=max⁡0,min⁡1,λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=max⁡0,min⁡1,τλ3-λ4Ih∗B2,u3∗=max⁡0,min⁡1,pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=max⁡0,min⁡1,λ2-λ1λhwSh∗B4.The optimality system is comprised of state system (1), adjoint system (61), initial conditions at t
=
0, 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 u
1, u
2, u
3, and u
4, respectively, are positive indicates that the optimal problem is a minimum at controls u
1
∗, u
2
∗, u
3
∗, and u
4
∗.The optimality system is solved using the forward-backward fourth-order Runge-Kutta scheme inR 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 u
1, u
2, u
3, and u
4 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].


## 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 variablesu
t
=
u
1
,
u
2
,
u
3
,
u
4
∈
U relative to the state variables S
h, E
h, I
h, R
h, S
m, E
m, I
m where control variables are bounded and measured with(49)
U
=
u
1
,
u
2
,
u
3
,
u
4
∈
U
is  Lebesgue  measurable  on
0,1
,
0
≤
u
i
t
≤
1
,
t
∈
0
,
T
,
i
=
1,2
,
3,4
.We define the objective function as(50)
J
u
1
,
u
2
,
u
3
,
u
4
=
∫
0
T
A
1
N
m
+
A
2
I
h
+
A
3
E
h
+
1
2
B
1
u
1
2
+
B
2
u
2
2
+
B
3
u
3
2
+
B
4
u
4
2
d
tsubject to(51)dShdt=Λh+ψRh-1-u1λhSh-1-u4λhwSh-μhSh,dEhdt=1-u1λhSh+1-u4λhwSh-αh+μhEh,dIhdt=αhEh-δh+μhIh-b+τu2Ih,dRhdt=b+τu2Ih-ψ+μhRh,dSmdt=Λm-1-u1λmSm-μm+au1+pu3Sm,dEmdt=1-u1λmSm-αmEm-μm+au1+pu3Em,dImdt=αmEm-μm+au1+pu3Im,Sh0≥0,Eh0≥0,Ih0≥0,Rh0≥0,Sm0≥0,Em0≥0,Im0≥0.In the objective function T is the final time and quantities A
1, A
2, and A
3 are weights constants of the total mosquito population, infected individuals, and exposed individuals, respectively, while B
1, B
2, B
3, and B
4 are weight constants for use with ITNs, treatment effort, IRS, and IPTp efforts, respectively. The total mosquito population N
m
=
S
m
+
E
m
+
I
m is part of the objective function because it is affected by the use of IRS and ITNs. In addition, E
h and I
h 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 B
1
u
1
2, treatment of infected individuals is B
2
u
2
2, spraying of houses with IRS is B
3
u
3
2, and preventive method of IPTp is B
4
u
4
2. A linear function has been chosen for the cost incurred by exposed individuals A
3
E
h, infected individuals, A
2
I
h, and the mosquito population, A
1
N
m. A quadratic form is used for the cost on the controls B
1
u
1
2, B
2
u
2
2, B
3
u
3
2, and B
4
u
4
2, such that the terms (
1
/
2
)
B
1
u
1
2, (
1
/
2
)
B
2
u
2
2, (
1
/
2
)
B
3
u
3
2, and (
1
/
2
)
B
4
u
4
2 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 populationN
m, the number of exposed humans E
h
(
t
), and infected humans I
h
(
t
) while minimizing the cost of control u
1
t, u
2
t, u
3
t, and u
4
(
t
). 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 [
0
,
1
].We seek an optimal controlu
1
∗, u
2
∗, u
3
∗, and u
4
∗ such that(52)Ju1∗,u2∗,u3∗,u4∗=minu1,u2,u3,u4∈U⁡Ju1,u2,u3,u4.Pontryagin’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, L, and Hamiltonian, H, with respect to u
1, u
2, u
3, and u
4.


## 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:(H
1):
the set of controls and corresponding state variables is nonempty.
(H
2):
the control setU is convex and closed.
(H
3):
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 u with coefficients depending on time and state.
(H
4):
there exist constantsc
1
,
c
2
>
0 and β
>
1 such that the integrand L
(
y
,
u
,
t
) of the objective functional J is convex and satisfies(53)Ly,u,t≥c1u12+u22+u32+u42β/2-c2.In 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, (
H
1
) is satisfied. The solutions are bounded and hence the control set satisfies (
H
2
) as well. The system is bilinear in u
1
,
u
2
,
u
3
,
u
4 and hence the right hand side of state system (1) satisfies condition (
H
3
) (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 U is closed and convex. The integrand of the objective cost function J expressed by (51) is a convex function of u
1
,
u
2
,
u
3
,
u
4 on the control set U. 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 ξ
1, ξ
2 and a constant ε
>
1 such that (54)A1Nm+A2Ih+A3Eh+12B1u12+B2u22+B3u32+B4u42≥ξ1u12+u22+u32+u42ε/2-ξ2,where ξ
1, ξ
2
>
0, A
1, A
2, B
1, B
2, B
3, B
4
>
0, and ε
>
1
.This concludes existence of an optimal control since the state variables are bounded. Hence we have the following theorem.Theorem 7.
Given the objective functionalJ
u
1
,
u
2
,
u
3
,
u
4
=
∫
0
T
A
1
N
m
+
A
2
I
h
+
A
3
E
h
+
(
1
/
2
)
(
B
1
u
1
2
+
B
2
u
2
2
+
B
3
u
3
2
+
B
4
u
4
2
)
d
t, where U
=
u
1
,
u
2
,
u
3
,
u
4
∣
0
≤
u
1
,
u
2
,
u
3
,
u
4
≤
1
,
0
≤
u
i
t
≤
1
,
t
∈
0
,
T
,
i
=
1,2
,
3,4 subject to (1) with initial conditions, then there exists an optimal control u
∗
=
u
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗ such that J
u
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗
=
min
U
⁡
J
u
1
,
u
2
,
u
3
,
u
4.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 by(55)LIh,Eh,Nm,u1,u2,u3,u4=A1Nm+A2Ih+A3Eh+12B1u12+B2u22+B3u32+B4u42.We seek to find the minimal value of the Lagrangian. To do this, we define the Hamiltonian H for the control problem which consists of the integrand of the objective functional (Lagrangian, L) and the inner product of the right hand sides of the state equations and the costate variables or adjoint variables λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7 as(56)H=L+λ1dShdt+λ2dEhdt+λ3dIhdt+λ4dRhdt+λ5dSmdt+λ6dEmdt+λ7dImdt.Taking X
=
S
h
,
E
h
,
I
h
,
R
h
,
S
m
,
E
m
,
I
m, U
=
u
1
,
u
2
,
u
3
,
u
4, and λ
=
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7 we obtain the Hamiltonian given by(57)
H
X
,
U
,
λ
=
L
I
h
,
E
h
,
N
m
,
u
1
,
u
2
,
u
3
,
u
4
+
λ
1
Λ
h
+
ψ
R
h
-
1
-
u
1
λ
h
S
h
-
1
-
u
4
λ
h
w
S
h
-
μ
h
S
h
+
λ
2
1
-
u
1
λ
h
S
h
+
1
-
u
4
λ
h
w
S
h
-
α
h
+
μ
h
E
h
+
λ
3
α
h
E
h
-
δ
h
+
μ
h
I
h
-
b
+
τ
u
2
I
h
+
λ
4
b
+
τ
u
2
I
h
-
ψ
+
μ
h
R
h
+
λ
5
Λ
h
-
1
-
u
1
λ
m
S
m
-
μ
m
+
a
u
1
+
p
u
3
S
m
+
λ
6
1
-
u
1
λ
m
S
m
-
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
E
m
+
λ
7
α
m
E
m
-
μ
m
+
a
u
1
+
p
u
3
I
m
.


## 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.Ifu
1
∗
,
u
2
∗
,
u
3
∗
,
u
4
∗ are optimal solution of an optimal control problem, then there exists a nontrivial vector function λ
t
=
λ
1
t
,
λ
2
t
,
…
,
λ
n
(
t
) satisfying the following conditions.The state equation is(58)dxdt=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂λ.The optimality condition is(59)0=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂uand the adjoint equation is(60)dλdt=∂Ht,u1∗,u2∗,u3∗,u4∗,λt∂x.Now, we apply the necessary conditions to the Hamiltonian.Theorem 8.
Given the optimal controlsu
1
∗, u
2
∗, u
3
∗, u
4
∗ and solutions S
h
∗, E
h
∗, I
h
∗, R
h
∗, S
m
∗, E
m
∗, I
m
∗ of the corresponding state system (1), there exist adjoint variables λ
1, λ
2, λ
3, λ
4, λ
5, λ
6, λ
7 satisfying(61)-dλ1dt=∂H∂Sh=1-u1λhλ2-λ1+1-u4λhwλ2-λ1-μhλ1,-dλ2dt=∂H∂Eh=αhλ3-λ2-μhλ2+A3,-dλ3dt=∂H∂Ih=b+τu2λ4-b+τu2+μh+δhλ3+A2-1-u1λϵϕSmNhλ5+1-u1λϵϕSmNhλ6,-dλ4dt=∂H∂Rh=ψλ1-μh+ψλ4,-dλ5dt=∂H∂Sm=1-u1λ6-λ5λm-μm+au1+pu3λ5+A1,-dλ6dt=∂H∂Em=αmλ7-λ6-μm+au1+pu3λ6+A1,-dλ7dt=∂H∂Im=-μm+au1+pu3λ7+A1+-1-u1βϵϕShNh-1-u4βϵϕShNhwλ1+1-u1βϵϕShNh+1-u4βϵϕShNhwλ2with transversality conditions(62)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=0.Furthermore u
1
∗, u
2
∗, u
3
∗, u
4
∗ are represented by(63)u1∗=max⁡0,min⁡1,λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=max⁡0,min⁡1,τλ3-λ4Ih∗B2,u3∗=max⁡0,min⁡1,pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=max⁡0,min⁡1,λ2-λ1λhwSh∗B4.Proof.
To determine the adjoint equations and the transversality conditions we use the HamiltonianH. The Hamiltonian function, H, is differentiated with respect to S
h, E
h, I
h, R
h, S
m, E
m, and I
m. The adjoint/costate equation is given by(64)-dλ1dt=∂H∂Sh=1-u1λhλ2-λ1+1-u4λhwλ2-λ1-μhλ1,-dλ2dt=∂H∂Eh=αhλ3-λ2-μhλ2+A3,-dλ3dt=∂H∂Ih=b+τu2λ4-b+τu2+μh+δhλ3+A2-1-u1λϵϕSmNhλ5+1-u1λϵϕSmNhλ6,-dλ4dt=∂H∂Rh=ψλ1-μh+ψλ4,-dλ5dt=∂H∂Sm=1-u1λ6-λ5λm-μm+au1+pu3λ5+A1,-dλ6dt=∂H∂Em=αmλ7-λ6-μm+au1+pu3λ6+A1,-dλ7dt=∂H∂Im=-μm+au1+pu3λ7+A1+-1-u1βϵϕShNh-1-u4βϵϕShNhwλ1+1-u1βϵϕShNh+1-u4βϵϕShNhwλ2with transversality conditions(65)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=0.In order to minimize Hamiltonian,H, with respect to the controls at the optimal controls, H is differentiated with respect to u
1, u
2, u
3, and u
4 on the set U, and the solution for the optimal control point is obtained after equating to zero. This is the optimality condition.Solving∂
H
/
∂
u
1
=
0, ∂
H
/
∂
u
2
=
0, ∂
H
/
∂
u
3
=
0, and ∂
H
/
∂
u
4
=
0, evaluating at the optimal control on the interior of the control set, where 0
<
u
i
<
1, for i
=
1,2
,
3,4, and letting S
h
=
S
h
∗, E
h
=
E
h
∗, I
h
=
I
h
∗, R
h
=
R
h
∗, S
m
=
S
m
∗, E
m
=
E
m
∗, and I
m
=
I
m
∗ yield (66)∂H∂u1=B1u1∗+λ1λhSh∗-λ2λhSh∗+λ5λmSm∗-aSm∗λ5-aEm∗λ6-aIm∗λ7=0,∂H∂u2=B2u2∗-τλ3Ih∗+τλ4Ih∗=0,∂H∂u3=B3u3∗-pλ5Sm∗-pλ6Em∗-pλ7Im∗=0,∂H∂u4=B4u4∗+λhwλ2Sh∗-λ1λhwSh∗=0for which(67)u1∗=λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=τλ3-λ4Ih∗B2,u3∗=pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=λ2-λ1λhwSh∗B4.By applying the boundary condition of each control, the solution of (67) becomes (68)u1∗=max⁡0,min⁡1,λ2-λ1λhSh∗+λ6-λ5λmSm∗+aSm∗λ5+aEm∗λ6+aIm∗λ7B1,u2∗=max⁡0,min⁡1,τλ3-λ4Ih∗B2,u3∗=max⁡0,min⁡1,pλ5Sm∗+λ6Em∗+λ7Im∗B3,u4∗=max⁡0,min⁡1,λ2-λ1λhwSh∗B4.The optimality system is comprised of state system (1), adjoint system (61), initial conditions at t
=
0, 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 u
1, u
2, u
3, and u
4, respectively, are positive indicates that the optimal problem is a minimum at controls u
1
∗, u
2
∗, u
3
∗, and u
4
∗.The optimality system is solved using the forward-backward fourth-order Runge-Kutta scheme inR 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 u
1, u
2, u
3, and u
4 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 asu
1
=
0.0904, u
2
=
0.165, u
3
=
0.076, and u
4
=
0.035 and the cost of intervention is for u
1
=
$
2.5
-
5, u
2
=
$
2.5, u
3
=
$
1.5, and u
4
=
$
2.5 (White et al. [38] and Hansen et al. [39]). The initial state variables are constant across all the epidemiological zones and are chosen as S
h
0
=
700, E
h
0
=
250, I
h
0
=
30, R
h
0
=
30, S
m
0
=
5000, E
m
0
=
500, and I
m
0
=
100. The following weight factors were also fixed for the different epidemiological scenarios as B
1
=
20, B
2
=
65, B
3
=
10, B
4
=
10, A
3
=
100, A
2
=
92, and A
1
=
20.


### 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 (I
h) as shown in Figure 2.Figure 2
Simulations showing the dynamics of human population with intervention strategies for the endemic setting.
(a)
(b)
(c)
(d)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 Figures3–6 show a significant difference in I
h and I
m with the control strategy compared to I
h and I
m without the control strategy. It is observed that the control strategy leads to decrease in the number of infected humans (I
h). The uncontrollable case leads to a decrease in the number of infected mosquitoes (I
m), 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.Figure 3
Simulations of the model showing the effect of treatment on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 4
Simulations of the model showing the effect of ITNs and treatment on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 5
Simulations of the model showing the effect of ITNs, treatment, and IRS on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 6
Simulations of the model showing the effect of ITNs, treatment, IRS, and IPTp on the spread of malaria for the endemic setting.
(a)
(b)
(c)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 (u
1) on ITNs is used to optimize the objective function J, while the control on treatment (u
2), the control on IRS (u
3), and control on IPTp (u
4) are set to zero.(b) Optimal Treatment. Only the control (u
2) on treatment is used to optimize the objective function J, while the control on ITNs (u
1), the control on IRS (u
3), and control on IPTp (u
4) are set to zero.(c) Optimal IRS. Only the control (u
3) on IRS is used to optimize the objective function J, while the control on treatment (u
2), the control on ITNs (u
1), and control on IPTp (u
4) are set to zero.(d) Optimal IPTp. Only the control (u
4) on IPTp is used to optimize the objective function J, while the control on treatment (u
2), the control on IRS (u
3), and control on ITNs (u
1) 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 (u
1) and the treatment (u
2) are used to optimize the objective function J while setting the control on IRS (u
3) and control on IPTp (u
4) to zero.(b) Optimal ITN and IRS. With this strategy, the control on ITNs (u
1) and the IRS (u
3) are used to optimize the objective function J while setting the control on treatment (u
2) and control on IPTp (u
4) to zero.(c) Optimal ITN and IPTp. With this strategy, the control on ITNs (u
1) and IPTp (u
4) are used to optimize the objective function J while setting the control on treatment (u
2) and control on IRS (u
3) to zero.(d) Optimal Treatment and IRS. With this strategy, the control on treatment (u
2) and the IRS (u
3) are used to optimize the objective function J while setting the control on ITNs (u
1) and control on IPTp (u
4) to zero.(e) Optimal Treatment and IPTp. With this strategy, the control on treatment (u
2) and the IPTp (u
4) are used to optimize the objective function J while setting the control on IRS (u
3) and control on ITNs (u
1) to zero.(f) Optimal IRS and IPTp. With this strategy, the control on IRS (u
3) and the IPTp (u
4) are used to optimize the objective function J while setting the control on treatment (u
2) and control on ITNs (u
1) 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 (u
1), treatment control (u
2), and IRS control (u
3) are used to optimize the objective function J, while IPTp control (u
4) is set to zero.(b) Optimal ITN, Treatment, and IPTp. In this case ITNs control (u
1), treatment control (u
2), and IPTp control (u
4) are used to optimize the objective function J, while IRS control (u
3) is set to zero.(c) Optimal ITN, IRS, and IPTp. In this case ITNs control (u
1), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J, while treatment control (u
2) is set to zero.(d) Optimal Treatment, IRS, and IPTp. In this case treatment control (u
2), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J, while ITNs control (u
1) 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 (u
1), treatment control (u
2), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J.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.


## 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 (I
h) as shown in Figure 2.Figure 2
Simulations showing the dynamics of human population with intervention strategies for the endemic setting.
(a)
(b)
(c)
(d)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 Figures3–6 show a significant difference in I
h and I
m with the control strategy compared to I
h and I
m without the control strategy. It is observed that the control strategy leads to decrease in the number of infected humans (I
h). The uncontrollable case leads to a decrease in the number of infected mosquitoes (I
m), 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.Figure 3
Simulations of the model showing the effect of treatment on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 4
Simulations of the model showing the effect of ITNs and treatment on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 5
Simulations of the model showing the effect of ITNs, treatment, and IRS on the spread of malaria for the endemic setting.
(a)
(b)
(c)Figure 6
Simulations of the model showing the effect of ITNs, treatment, IRS, and IPTp on the spread of malaria for the endemic setting.
(a)
(b)
(c)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 (u
1) on ITNs is used to optimize the objective function J, while the control on treatment (u
2), the control on IRS (u
3), and control on IPTp (u
4) are set to zero.(b) Optimal Treatment. Only the control (u
2) on treatment is used to optimize the objective function J, while the control on ITNs (u
1), the control on IRS (u
3), and control on IPTp (u
4) are set to zero.(c) Optimal IRS. Only the control (u
3) on IRS is used to optimize the objective function J, while the control on treatment (u
2), the control on ITNs (u
1), and control on IPTp (u
4) are set to zero.(d) Optimal IPTp. Only the control (u
4) on IPTp is used to optimize the objective function J, while the control on treatment (u
2), the control on IRS (u
3), and control on ITNs (u
1) 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 (u
1) and the treatment (u
2) are used to optimize the objective function J while setting the control on IRS (u
3) and control on IPTp (u
4) to zero.(b) Optimal ITN and IRS. With this strategy, the control on ITNs (u
1) and the IRS (u
3) are used to optimize the objective function J while setting the control on treatment (u
2) and control on IPTp (u
4) to zero.(c) Optimal ITN and IPTp. With this strategy, the control on ITNs (u
1) and IPTp (u
4) are used to optimize the objective function J while setting the control on treatment (u
2) and control on IRS (u
3) to zero.(d) Optimal Treatment and IRS. With this strategy, the control on treatment (u
2) and the IRS (u
3) are used to optimize the objective function J while setting the control on ITNs (u
1) and control on IPTp (u
4) to zero.(e) Optimal Treatment and IPTp. With this strategy, the control on treatment (u
2) and the IPTp (u
4) are used to optimize the objective function J while setting the control on IRS (u
3) and control on ITNs (u
1) to zero.(f) Optimal IRS and IPTp. With this strategy, the control on IRS (u
3) and the IPTp (u
4) are used to optimize the objective function J while setting the control on treatment (u
2) and control on ITNs (u
1) 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 (u
1), treatment control (u
2), and IRS control (u
3) are used to optimize the objective function J, while IPTp control (u
4) is set to zero.(b) Optimal ITN, Treatment, and IPTp. In this case ITNs control (u
1), treatment control (u
2), and IPTp control (u
4) are used to optimize the objective function J, while IRS control (u
3) is set to zero.(c) Optimal ITN, IRS, and IPTp. In this case ITNs control (u
1), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J, while treatment control (u
2) is set to zero.(d) Optimal Treatment, IRS, and IPTp. In this case treatment control (u
2), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J, while ITNs control (u
1) 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 (u
1), treatment control (u
2), IRS control (u
3), and IPTp control (u
4) are used to optimize the objective function J.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 ifR
0
≤
1 in D. If R
0
>
1, the unique endemic equilibrium exists and is globally asymptotically stable in D. The model exhibited backward bifurcation backward bifurcation at R
0
=
1 implying the existence of multiple endemic equilibria for R
0
<
1. The existence of multiple endemic equilibria emphasizes the fact that R
0
<
1 is not sufficient to eradicate disease from the population and the need to lower R
0 much below one to make the disease-free equilibrium stable globally. This behavior has important public health implications because it means that bringing R
0 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 att
=
0, 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 theR 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 freeR 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 ifR
0
≤
1, the disease-free equilibrium is globally asymptotically stable in D. If R
0
>
1, the unique endemic equilibrium exists and is globally asymptotically stable. The model also exhibits backward bifurcation at R
0
=
1.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.



---

*Source: 8013574-2016-09-01.xml*